atomic.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module atomic_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use logrid_oct_m
  use messages_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use sort_oct_m
  use xc_f03_lib_m
#ifdef HAVE_LIBXC_FUNCS
  use xc_f03_funcs_m
#endif
  use iso_c_binding
  implicit none
 
  private
  public ::        &
    atomxc,        &
    atomhxc,       &
    vhrtre, egofv
 
  public ::            &
    valconf_t,         &
    valconf_copy,      &
    valconf_sort_by_l, &
    write_valconf,     &
    read_valconf,      &
    valconf_unpolarized_to_polarized, &
    valconf_string_length
 
  character(len=1), parameter :: &
    spec_notation(0:3) = (/ 's', 'p', 'd', 'f' /)
 
  integer, parameter :: VALCONF_STRING_LENGTH = 80
 
  type valconf_t
    ! Components are public by default
    integer           :: z = 0
    character(len=3)  :: symbol = ""
    integer           :: type = 0
    integer           :: p = 0
    integer           :: n(15) = 0
    integer           :: l(15) = 0
    real(real64)      :: occ(15,2) = m_zero 
    real(real64)      :: j(15) = m_zero     
  end type valconf_t
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine valconf_copy(cout, cin)
    type(valconf_t), intent(out) :: cout
    type(valconf_t), intent(in)  :: cin
 
    push_sub(valconf_copy)
 
    cout%z      = cin%z
    cout%symbol = cin%symbol
    cout%type   = cin%type
    cout%p      = cin%p
    cout%n      = cin%n
    cout%l      = cin%l
    cout%occ    = cin%occ
    cout%j      = cin%j
 
    pop_sub(valconf_copy)
  end subroutine valconf_copy
 
  subroutine valconf_sort_by_l(c, lmax)
    type(valconf_t), intent(inout)  :: c
    integer,         intent(in)     :: lmax
 
    integer :: ll
    integer, allocatable :: order(:)
    type(valconf_t) :: tmp_c
 
    push_sub(valconf_sort_by_l)
 
    call valconf_copy(tmp_c, c)
 
    ! In order to work in the following, we need to sort the occupations by angular momentum
    safe_allocate(order(1:lmax+1))
    call sort(c%l(1:lmax+1), order)
    do ll = 0, lmax
      c%occ(ll+1, 1:2) = tmp_c%occ(order(ll+1), 1:2)
      c%n(ll+1) = tmp_c%n(order(ll+1))
      c%j(ll+1) = tmp_c%j(order(ll+1))
    end do
    safe_deallocate_a(order)
 
    pop_sub(valconf_sort_by_l)
  end subroutine valconf_sort_by_l
 
  ! ---------------------------------------------------------
  subroutine write_valconf(c, s)
    type(valconf_t), intent(in) :: c
    character(len=VALCONF_STRING_LENGTH), intent(out) :: s
 
    integer :: j
 
    push_sub(write_valconf)
 
    write(s,'(i2,1x,a2,i1,1x,i1,a1,6(i1,a1,f6.3,a1))') c%z, c%symbol, c%type, c%p, ':',&
      (c%n(j),spec_notation(c%l(j)),c%occ(j,1),',',j=1,c%p)
 
    pop_sub(write_valconf)
  end subroutine write_valconf
 
 
  ! ---------------------------------------------------------
  subroutine read_valconf(namespace, s, c)
    type(namespace_t),        intent(in)    :: namespace
    character(len=VALCONF_STRING_LENGTH), intent(in) :: s
    type(valconf_t), intent(out) :: c
 
    integer :: j
    character(len=1) :: lvalues(1:6)
 
    push_sub(read_valconf)
 
    read (s,'(i2,1x,a2,i1,1x,i1,1x,6(i1,a1,f6.3,1x))') c%z, c%symbol, c%type, c%p,&
      (c%n(j),lvalues(j),c%occ(j,1),j=1,c%p)
    do j = 1, c%p
      select case (lvalues(j))
      case ('s')
        c%l(j) = 0
      case ('p')
        c%l(j) = 1
      case ('d')
        c%l(j) = 2
      case ('f')
        c%l(j) = 3
      case default
        message(1) = 'read_valconf.'
        call messages_fatal(1, namespace=namespace)
      end select
    end do
 
    c%j = m_zero
 
    pop_sub(read_valconf)
  end subroutine read_valconf
 
  !In case of spin-polarized calculations, we properly distribute the electrons
  subroutine valconf_unpolarized_to_polarized(conf)
    type(valconf_t), intent(inout) :: conf
 
    integer :: ii, ll
    real(real64)   :: xx
 
    push_sub(valconf_unpolarized_to_polarized)
 
    do ii = 1, conf%p
      ll = conf%l(ii)
      xx = conf%occ(ii, 1)
      conf%occ(ii, 1) = min(xx, real(2*ll+1, real64) )
      conf%occ(ii, 2) = xx - conf%occ(ii,1)
    end do
 
    pop_sub(valconf_unpolarized_to_polarized)
  end subroutine valconf_unpolarized_to_polarized
 
 
  ! ---------------------------------------------------------
  subroutine atomhxc(namespace, functl, g, nspin, dens, v, extra)
    type(namespace_t), intent(in)  :: namespace
    character(len=*),  intent(in)  :: functl
    type(logrid_t),    intent(in)  :: g
    integer,           intent(in)  :: nspin
    real(real64),      intent(in)  :: dens(g%nrval, nspin)
    real(real64),      intent(out) :: v(g%nrval, nspin)
    real(real64),      intent(in), optional :: extra(g%nrval)
 
    character(len=5) :: xcfunc, xcauth
    integer :: is, ir
    real(real64), allocatable :: xc(:, :), ve(:, :), rho(:, :)
    real(real64) :: r2, ex, ec, dx, dc
 
    push_sub(atomhxc)
 
    safe_allocate( ve(1:g%nrval, 1:nspin))
    safe_allocate( xc(1:g%nrval, 1:nspin))
    safe_allocate(rho(1:g%nrval, 1:nspin))
    ve = m_zero
    xc = m_zero
    rho = m_zero
 
    ! To calculate the Hartree term, we put all the density in one variable.
    do is = 1, nspin
      rho(:, 1) = rho(:, 1) + dens(:, is)
    end do
    ve = m_zero
    do is = 1, nspin
      call vhrtre(rho(:, 1), ve(:, is), g%rofi, g%drdi, g%s, g%nrval, g%a)
    end do
    ve(1, 1:nspin) = ve(2, 1:nspin)
 
    select case (functl)
    case ('LDA')
      xcfunc = 'LDA'
      xcauth = 'PZ'
    case ('GGA')
      xcfunc = 'GGA'
      xcauth = 'PBE'
    case default
      message(1) = 'Internal Error in atomhxc: unknown functl'
      call messages_fatal(1, namespace=namespace)
    end select
 
    rho(:, :) = dens(:, :)
    do is = 1, nspin
      if (present(extra)) rho(:, is) = rho(:, is) + extra(:)/nspin
      rho(2:g%nrval, is) = rho(2:g%nrval, is)/(m_four*m_pi*g%r2ofi(2:g%nrval))
    end do
    r2 = g%rofi(2)/(g%rofi(3)-g%rofi(2))
    rho(1, 1:nspin) = rho(2, 1:nspin) - (rho(3, 1:nspin)-rho(2, 1:nspin))*r2
 
    do is = 1, nspin
      do ir = 1, g%nrval
        if (rho(ir, is) < m_epsilon) rho(ir, is) = m_zero
      end do
    end do
 
    call atomxc(namespace, xcfunc, xcauth, g%nrval, g%nrval, g%rofi, nspin, rho, ex, ec, dx, dc, xc)
    v = ve + xc
 
    safe_deallocate_a(ve)
    safe_deallocate_a(xc)
    safe_deallocate_a(rho)
 
    pop_sub(atomhxc)
  end subroutine atomhxc
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine atomxc(namespace, FUNCTL, AUTHOR, NR, MAXR, RMESH, NSPIN, DENS, EX, EC, DX, DC, V_XC)
    type(namespace_t), intent(in) :: namespace
    character(len=*), intent(in) :: functl, author
    integer, intent(in) :: nr, maxr
    real(real64), intent(in) :: rmesh(maxr)
    integer, intent(in) :: nspin
    real(real64), intent(in) :: dens(maxr,nspin)
    real(real64), intent(out) :: dc, dx, ec, ex, v_xc(maxr,nspin)
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Internal parameters                                                         !
! NN     : order of the numerical derivatives: the number of radial           !
!          points used is 2*NN+1                                              !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    integer, parameter :: nn = 5
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Fix energy unit:  EUNIT=1.0 => Hartrees,                                    !
!                   EUNIT=0.5 => Rydbergs,                                    !
!                   EUNIT=0.03674903 => eV                                    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    real(real64), parameter :: eunit = m_half
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! DVMIN is added to differential of volume to avoid division by zero          !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    real(real64), parameter :: dvmin = 1.0e-12_real64
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Local variables and arrays                                                  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    logical :: gga
    integer :: in, in1, in2, ir, is, jn
    real(real64) :: aux(maxr), d(nspin),                              &
      dgdm(-nn:nn), dgidfj(-nn:nn), drdm, dvol,                      &
      f1, f2, gd(3,nspin), sigma(3), vxsigma(3), vcsigma(3),   &
      epsx(1), epsc(1),                                              &
      dexdd(nspin), dexdgd(3,nspin), decdd(nspin), decdgd(3,nspin)
    type(xc_f03_func_t) :: x_conf, c_conf
 
    push_sub(atomxc)
 
    ! sanity check
    assert(nspin == 1 .or. nspin == 2)
 
    ! Set GGA switch
    if (functl .eq. 'LDA' .or. functl .eq. 'lda' .or.                              &
      functl .eq. 'LSD' .or. functl .eq. 'lsd') then
      gga = .false.
    else if (functl .eq. 'GGA' .or. functl .eq. 'gga') then
      gga = .true.
    else
      gga = .false.
      write(message(1),'(a,a)') 'atomxc: Unknown functional ', functl
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! initialize xc functional
    if (gga) then
      call xc_f03_func_init(x_conf, xc_gga_x_pbe, nspin)
      call xc_f03_func_init(c_conf, xc_gga_c_pbe, nspin)
    else
      call xc_f03_func_init(x_conf, xc_lda_x, nspin)
      if (author .eq. 'CA' .or. author .eq. 'ca' .or.                                &
        author .eq. 'PZ' .or. author .eq. 'pz') then
        call xc_f03_func_init(c_conf, xc_lda_c_pz, nspin)
      else if (author .eq. 'PW92' .or. author .eq. 'pw92') then
        call xc_f03_func_init(c_conf, xc_lda_c_pw, nspin)
      else
        write(message(1),'(a,a)') 'LDAXC: Unknown author ', author
        call messages_fatal(1, namespace=namespace)
      end if
    end if
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Initialize output                                                           !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    ex = m_zero
    ec = m_zero
    dx = m_zero
    dc = m_zero
    do is = 1,nspin
      do ir = 1,nr
        v_xc(ir,is) = m_zero
      end do
    end do
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Loop on mesh points                                                         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    do ir = 1,nr
 
      ! Find interval of neighbour points to calculate derivatives
      in1 = max(  1, ir-nn ) - ir
      in2 = min( nr, ir+nn ) - ir
 
      ! Find weights of numerical derivative from Lagrange
      ! interpolation formula
      do in = in1,in2
        if (in == 0) then
          dgdm(in) = 0
          do jn = in1,in2
            if (jn /= 0) dgdm(in) = dgdm(in) + m_one / (0 - jn)
          end do
        else
          f1 = 1
          f2 = 1
          do jn = in1,in2
            if (jn /= in .and. jn /= 0) f1 = f1 * (0  - jn)
            if (jn /= in)               f2 = f2 * (in - jn)
          end do
          dgdm(in) = f1 / f2
        end if
      end do
 
      ! Find dr/dmesh
      drdm = 0
      do in = in1,in2
        drdm = drdm + rmesh(ir+in) * dgdm(in)
      end do
 
      ! Find differential of volume. Use trapezoidal integration rule
      dvol = 4 * m_pi * rmesh(ir)**2 * drdm
 
      ! DVMIN is a small number added to avoid a division by zero
      dvol = safe_tol(dvol, dvmin)
      if (ir == 1 .or. ir == nr) dvol = dvol / 2
      if (gga) aux(ir) = dvol
 
      ! Find the weights for the derivative d(gradF(i))/d(F(j)), of
      ! the gradient at point i with respect to the value at point j
      if (gga) then
        do in = in1,in2
          dgidfj(in) = dgdm(in) / drdm
        end do
      end if
 
      ! Find density and gradient of density at this point
      do is = 1,nspin
        d(is) = dens(ir,is)
      end do
      if (gga) then
        do is = 1,nspin
          gd(:,is) = m_zero
          do in = in1,in2
            gd(3,is) = gd(3,is) + dgidfj(in) * dens(ir+in,is)
          end do
        end do
      else
        gd = m_zero
      end if
 
      ! The xc_f03_gga routines need as input these combinations of
      ! the gradient of the densities.
      sigma(1) = gd(3, 1)*gd(3, 1)
      if (nspin > 1) then
        sigma(2) = gd(3, 1) * gd(3, 2)
        sigma(3) = gd(3, 2) * gd(3, 2)
      end if
 
      ! Find exchange and correlation energy densities and their
      ! derivatives with respect to density and density gradient
      if (gga) then
        call xc_f03_gga_exc_vxc(x_conf, int(1, c_size_t), d, sigma, epsx, dexdd, vxsigma)
        call xc_f03_gga_exc_vxc(c_conf, int(1, c_size_t), d, sigma, epsc, decdd, vcsigma)
      else
        call xc_f03_lda_exc_vxc(x_conf, int(1, c_size_t), d, epsx, dexdd)
        call xc_f03_lda_exc_vxc(c_conf, int(1, c_size_t), d, epsc, decdd)
      end if
 
      ! The derivatives of the exchange and correlation energies per particle with
      ! respect to the density gradient have to be recovered from the derivatives
      ! with respect to the sigma values defined above.
      if (gga) then
        dexdgd(:, 1) = m_two * vxsigma(1) * gd(:, 1)
        decdgd(:, 1) = m_two * vcsigma(1) * gd(:, 1)
        if (nspin == 2) then
          dexdgd(:, 1) = dexdgd(:, 1) + vxsigma(2)*gd(:, 2)
          dexdgd(:, 2) = m_two * vxsigma(3) * gd(:, 2) + vxsigma(2) * gd(:, 1)
          decdgd(:, 1) = decdgd(:, 1) + vcsigma(2)*gd(:, 2)
          decdgd(:, 2) = m_two * vcsigma(3) * gd(:, 2) + vcsigma(2) * gd(:, 1)
        end if
      end if
 
      ! Add contributions to exchange-correlation energy and its
      ! derivatives with respect to density at all points
      do is = 1, nspin
        ex = ex + dvol * d(is) * epsx(1)
        ec = ec + dvol * d(is) * epsc(1)
        dx = dx + dvol * d(is) * (epsx(1) - dexdd(is))
        dc = dc + dvol * d(is) * (epsc(1) - decdd(is))
        if (gga) then
          v_xc(ir,is) = v_xc(ir,is) + dvol * (dexdd(is) + decdd(is))
          do in = in1,in2
            dx = dx - dvol * dens(ir+in,is) * dexdgd(3,is) * dgidfj(in)
            dc = dc - dvol * dens(ir+in,is) * decdgd(3,is) * dgidfj(in)
            v_xc(ir+in,is) = v_xc(ir+in,is) + dvol *                         &
              (dexdgd(3,is) + decdgd(3,is)) * dgidfj(in)
          end do
        else
          v_xc(ir,is) = dexdd(is) + decdd(is)
        end if
      end do
 
    end do
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Divide by volume element to obtain the potential (per electron)             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    if (gga) then
      do is = 1,nspin
        do ir = 1,nr
          dvol = aux(ir)
          v_xc(ir,is) = v_xc(ir,is) / dvol
        end do
      end do
    end if
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Divide by energy unit                                                       !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    ex = ex / eunit
    ec = ec / eunit
    dx = dx / eunit
    dc = dc / eunit
    do is = 1,nspin
      do ir = 1,nr
        v_xc(ir,is) = v_xc(ir,is) / eunit
      end do
    end do
 
    call xc_f03_func_end(x_conf)
    call xc_f03_func_end(c_conf)
 
    pop_sub(atomxc)
  end subroutine atomxc
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine vhrtre(rho, v, r, drdi, srdrdi, nr, a)
    real(real64),   intent(in)  :: rho(:)
    real(real64),   intent(out) :: v(:)
    real(real64),   intent(in)  :: r(:), drdi(:), srdrdi(:)
    integer, intent(in)  :: nr
    real(real64),   intent(in)  :: a
 
    real(real64) :: x,y,q,a2by4,ybyq,qbyy,qpartc,v0,qt,dz,t,beta,dv
    integer :: nrm1,nrm2,ir
 
    push_sub(vhrtre)
 
    nrm1 = nr - 1
    nrm2 = nr - 2
    a2by4 = a*a/m_four
    ybyq = m_one - a*a/48.0_real64
    qbyy = m_one/ybyq
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  SIMPSON`S RULE IS USED TO PERFORM TWO INTEGRALS OVER THE ELECTRON          !
!  DENSITY.  THE TOTAL CHARGE QT IS USED TO FIX THE POTENTIAL AT R=R(NR)      !
!  AND V0 (THE INTEGRAL OF THE ELECTRON DENSITY DIVIDED BY R) FIXES           !
!  THE ELECTROSTATIC POTENTIAL AT THE ORIGIN                                  !
!  Modified to be consistent with pseudopotential generation (use the         !
!  trapezoidal rule for integration). A. Rubio. Jan. 2000                     !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    v0 = m_zero
    qt = m_zero
    do ir = 2, nrm1, 2
      dz = drdi(ir)*rho(ir)
      qt = qt + dz
      v0 = v0 + dz/r(ir)
    end do
 
    do ir=3, nrm2, 2
      dz = drdi(ir)*rho(ir)
      qt = qt + dz
      v0 = v0 + dz/r(ir)
    end do
    dz = drdi(nr)*rho(nr)
 
    qt = (qt + qt + dz)*m_half
    v0 = (v0 + v0 + dz/r(nr))
    v(1) = m_two*v0
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  THE ELECTROSTATIC POTENTIAL AT R=0 IS SET EQUAL TO                         !
!                       THE AVERAGE VALUE OF RHO(R)/R                         !
!  BEGIN CONSTRUCTION OF THE POTENTIAL AT FINITE                              !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    ir   = 2
    t    = srdrdi(ir)/r(ir)
    beta = drdi(ir)*t*rho(ir)
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  THE NEXT 4 STATEMENTS INDICATE THAT WE FIRST FIND THE PARTICULAR           !
!  SOLUTION TO THE INHOMOGENEOUS EQN. FOR WHICH Q(2)=0, WE THEN               !
!  ADD TO THIS PARTICULAR SOLUTION A SOLUTION OF THE HOMOGENEOUS EQN.         !
!  (A CONSTANT IN V OR A Q PROPORTIONAL TO R)                                 !
!  WHICH WHEN DIVIDED BY R IN GOING FROM Q TO V GIVES                         !
!  THE POTENTIAL THE DESIRED COULOMB TAIL OUTSIDE THE ELECTRON DENSITY.       !
!  THE SIGNIFICANCE OF THE SOLUTION VANISHING AT THE SECOND RADIAL            !
!  MESH POINT IS THAT, SINCE ALL REGULAR SOLUTIONS OF THE EQUATION            !
!  FOR Q=R*V VANISH AT THE ORIGIN, THE KNOWLEDGE OF THE SOLUTION              !
!  VALUE AT THE SECOND MESH POINT PROVIDES THE TWO SOLUTION VALUES            !
!  REQUIRED TO START THE NUMEROV PROCEDURE.                                   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    x = m_zero
    y = m_zero
    q = (y - beta / 12.0_real64) * qbyy
    v(ir) = m_two * t * q
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  BEGINNING OF THE NUMEROV ALGORITHM                                         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    do ir = 2, nr
      x  = x + a2by4 * q - beta
      y  = y + x
      t  = srdrdi(ir) / r(ir)
      beta = t * drdi(ir) * rho(ir)
      q = (y - beta / 12.0_real64) * qbyy
      v(ir) = m_two * t * q
    end do
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  END OF THE NUMEROV ALGORITHM                                               !
!                                                                             !
!  WE HAVE NOW FOUND A PARTICULAR SOLUTION TO THE INHOMOGENEOUS EQN.          !
!  FOR WHICH Q(R) AT THE SECOND RADIAL MESH POINT EQUALS ZERO.                !
!  NOTE THAT ALL REGULAR SOLUTIONS TO THE EQUATION FOR Q=R*V                  !
!  VANISH AT THE ORIGIN.                                                      !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    qpartc = r(nr)*v(nr)/m_two
    dz = qt - qpartc
    dv = m_two*dz/r(nr)
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  THE LOOP FOLLOWING ADDS THE CONSTANT SOLUTION OF THE HOMOGENEOUS           !
!  EQN TO THE PARTICULAR SOLUTION OF THE INHOMOGENEOUS EQN.                   !
!  NOTE THAT V(1) IS CONSTRUCTED INDEPENDENTLY                                !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    do ir = 2, nr
      v(ir) = v(ir) + dv
    end do
 
    pop_sub(vhrtre)
  end subroutine vhrtre
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! FROM NOW, ON, ALL SUBROUTINES ARE IN DOUBLE PRECISION, NO MATTER IF THE CODE
! IS COMPILED WITH SINGLE PRECISION OR NOT. APPARENTLY DOUBLE PRECISION IS
! NEEDED FOR THESE PROCEDURES.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine egofv(namespace, h, s, n, e, g, y, l, z, a, b, rmax, nprin, nnode, dr, ierr)
    type(namespace_t), intent(in)    :: namespace
    real(real64),          intent(in)    :: h(:), s(:)
    integer,           intent(in)    :: n
    real(real64),          intent(inout) :: e
    real(real64),          intent(out)   :: g(:), y(:)
    integer,           intent(in)    :: l
    real(real64),          intent(in)    :: z, a, b, rmax
    integer,           intent(in)    :: nprin, nnode
    real(real64),          intent(in)    :: dr
    integer,           intent(out)   :: ierr
 
    integer :: i,ncor,n1,n2,niter,nt
    real(real64)               :: e1, e2, del, de, et, t
    real(real64), parameter    :: tol = 1.0e-5_real64
 
    push_sub(egofv)
 
    ncor = nprin-l-1
    n1 = nnode
    n2 = nnode-1
    e1 = e
    e2 = e
    ierr = 1
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  the labels 1 and 2 refer to the bisection process, defining the            !
!  range in which the desired solution is located.  the initial               !
!  settings of n1, n2, e1 and e2 are not consistent with the bisection        !
!  algorithm; they are set to consistent values when the desired              !
!  energy interval has been located.                                          !
!                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    nt = 0
    del = m_half
    de  = m_zero
 
    do niter = 0, 40
 
      et = e + de
      ! the following line is the fundamental "bisection"
      e = m_half*(e1+e2)
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  the following concatenation of logical ors ensures that node               !
!  counting is used unless the previous integration of the radial             !
!  eq produced both the correct number of nodes and a sensible                !
!  prediction for the energy.                                                 !
!                                                                             !
!     sensible means that et must be greater than e1 and less than e2         !
!     correct number of nodes means that nt = nnode or nnode-1.               !
!                                                                             !
!     leaving e set to its bisection value, and transfering to                !
!     the call to yofe means that we are performing bisection,                !
!     whereas setting e to et is use of the variational estimate.             !
!                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 
      if (.not. (et <= e1 .or. et >= e2 .or. nt < nnode-1 .or. nt > nnode)) then
        e = et
        if (abs(de) < tol) then
          ierr = 0
          exit
        end if
      end if
      call yofe(e,de,dr,rmax,h,s,y,n,l,ncor,nt,z,a,b)
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  yofe integrates the schro eq.; now the bisection logic                     !
!                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
      if (nt < nnode) then
        !  too few nodes; set e1 and n1
        e1 = e
        n1 = nt
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  at this point, we have just set the bottom of the bisection range;         !
!  if the top is also set, we proceed.  if the top of the range has not       !
!  been set, it means that we have yet to find an e greater than the          !
!  desired energy.  the upper end of the range is extended.                   !
!                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
        if (n2 >= nnode) cycle
        del=del*m_two
        e2=e1+del
        cycle
      end if
 
      !  too many nodes; set e2 and n2
      e2=e
      n2=nt
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  at this point, we have just set the top of the bisection range;            !
!  if the top is also set, we proceed.  if the top of the range has           !
!  not been set, it means that we have yet to find an e less than the         !
!  desired energy.  the lower end of the range is extended.                   !
!                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
      if (n1  <  nnode) cycle
      del=del*m_two
      e1=e2-del
 
    end do
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  the Numerov method uses a transformation of the radial wavefunction.       !
!  that we call "y".  having located the eigenenergy, we transform            !
!  y to "g", from which the density is easily constructed.                    !
!  finally, the call to "nrmlzg" normalizes g to one electron.                !
!                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    if (ierr /= 1) then
      g(1) = m_zero
      do i = 2, n
        t = h(i) - e * s(i)
        g(i) = y(i) / (m_one - t / 12._real64)
      end do
      call nrmlzg(namespace, g, s, n)
    end if
 
    pop_sub(egofv)
  end subroutine egofv
 
 
  subroutine yofe(e,de,dr,rmax,h,s,y,nmax,l,ncor,nnode,z,a,b)
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    real(real64), intent(inout) :: e
    real(real64), intent(out)   :: de
    real(real64), intent(in)    :: dr, rmax
    real(real64), intent(in)    :: h(:), s(:)
    real(real64), intent(out)   :: y(:)
    integer,  intent(in)    :: nmax, l
    integer,  intent(in)    :: ncor
    integer,  intent(out)   :: nnode
    real(real64), intent(in)    :: z, a, b
 
    real(real64) :: y2, g, gsg, x, gin, gsgin, xin
    integer :: n,knk,nndin,i
    real(real64) :: zdr, yn, ratio, t
 
    ! No PUSH SUB, called too often.
 
    zdr = z*a*b
 
    do n = nmax, 1, -1
      if (h(n)-e*s(n) < m_one) then
        knk = n
        exit
      end if
      y(n) = m_zero
    end do
 
    call bcorgn(e,h,s,l,zdr,y2)
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  bcorgn computes y2, which embodies the boundary condition                  !
!  satisfied by the radial wavefunction at the origin                         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    call numout(e,h,s,y,ncor,knk,nnode,y2,g,gsg,x)
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  the outward integration is now complete                                    !
!                                                                             !
!     we first decide if the kinetic energy is sufficiently non               !
!     negative to permit use of the numerov eq at rmax.  if                   !
!     it is not, then zero-value boundary condition is used                   !
!                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 
    yn = m_zero
    if (.not. (n < nmax .or. abs(dr) > 1.e3_real64)) then
      call bcrmax(e,dr,rmax,h,s,n,yn,a,b)
    end if
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  bcrmax computes yn, which embodies the boundary condition                  !
!  satisfied by the radial wavefunction at rmax                               !
!                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    ! n should be larger than 1
    n = max(n, 2)
    knk = max(knk, 1)
    call numin(e,h,s,y,n,nndin,yn,gin,gsgin,xin,knk)
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  numin performs the inward integration by the Numerov method                !
!                                                                             !
!  the energy increment is now evaluated from the kink in psi                 !
!                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 
    ratio = g / gin
    xin = xin * ratio
    gsg = gsg + gsgin * ratio * ratio
    t = h(knk) - e * s(knk)
    de = g * (x + xin + t * g) / gsg
    nnode = nnode + nndin
    if (de < m_zero) nnode = nnode + 1
 
    do i = knk, n
      y(i) = y(i) * ratio
    end do
 
  end subroutine yofe
 
 
  subroutine nrmlzg(namespace, g, s, n)
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    type(namespace_t), intent(in)    :: namespace
    real(real64),          intent(inout) :: g(:)
    real(real64),          intent(in)    :: s(:)
    integer,           intent(in)    :: n
 
    integer :: nm1,nm2,i
    real(real64) :: norm, srnrm
 
    push_sub(nrmlzg)
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  determine the norm of g(i) using Simpson`s rule                            !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 
    if (mod(n,2) /= 1) then
      write(message(1),'(a,i6)') ' nrmlzg: n should be odd. n =', n
      call messages_warning(1, namespace=namespace)
    end if
 
    norm = m_zero
    nm1 = n - 1
    do i = 2,nm1,2
      norm=norm + g(i)*s(i)*g(i)
    end do
    norm = norm * m_two
    nm2  = n - 2
    do i = 3,nm2,2
      norm=norm + g(i)*s(i)*g(i)
    end do
    norm = norm * m_two
    nm2  = n - 2
    do i = 1, n, nm1
      norm = norm + g(i) * s(i) * g(i)
    end do
    norm = norm / m_three
    srnrm = sqrt(norm)
    do i = 1, n
      g(i) = g(i) / srnrm
    end do
 
    pop_sub(nrmlzg)
  end subroutine nrmlzg
 
 
  subroutine bcorgn(e,h,s,l,zdr,y2)
    real(real64), intent(in) :: e
    real(real64), intent(in) :: h(:), s(:)
    integer,     intent(in) :: l
    real(real64), intent(in) :: zdr
    real(real64), intent(out) :: y2
 
    real(real64) :: t2, t3, d2, c0, c1, c2
 
    ! no PUSH SUB, called too often
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   the quantity called d(i) in the program is actually the inverse           !
!   of the diagonal of the tri-diagonal Numerov matrix                        !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    t2 = h(2) - e * s(2)
    d2 = -((24._real64 + 10.0_real64 * t2) / (12._real64 - t2))
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  The following section deals with the fact that the independent             !
!  variable "y" in the Numerov equation is not zero at the origin             !
!  for l less than 2.                                                         !
!  The l=0 solution g vanishes, but the first and second                      !
!  derivatives are finite, making the Numerov variable y finite.              !
!  The l=1 solution g vanishes, and gprime also vanishes, but                 !
!  The second derivative gpp is finite making y finite. For l > 1,            !
!  g and its first two derivatives vanish, making y zero.                     !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    if (l  <  2) then
      if (l  <=  0) then
        c0=zdr/6.0_real64
        c0=c0/(m_one-0.75_real64*zdr)
      else
        c0=m_one/12._real64
        c0=(-c0)*8.0_real64/m_three
      end if
 
      c1=c0*(12._real64+13._real64*t2)/(12._real64-t2)
      t3=h(3)-e*s(3)
      c2=(-m_half)*c0*(24._real64-t3)/(12._real64-t3)
      d2=(d2-c1)/(m_one-c2)
    end if
    y2=(-m_one)/d2
 
    if (l  <  2) then
      if (l  <=  0) then
        c0=zdr/6.0_real64
        c0=c0/(m_one-0.75_real64*zdr)
      else
        c0=m_one/12._real64
        c0=(-c0)*8.0_real64/m_three
      end if
      c1=c0*(12._real64+13._real64*t2)/(12._real64-t2)
      t3=h(3)-e*s(3)
      c2=(-m_half)*c0*(24._real64-t3)/(12._real64-t3)
      d2=(d2-c1)/(m_one-c2)
    end if
    y2=(-m_one)/d2
 
  end subroutine bcorgn
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 22.7.85                                                                     !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine bcrmax(e, dr, rmax, h, s, n, yn, a, b)
    real(real64), intent(in)  :: e, dr, rmax
    real(real64), intent(in)  :: h(:), s(:)
    integer,  intent(in)  :: n
    real(real64), intent(out) :: yn
    real(real64), intent(in)  :: a, b
 
    real(real64) :: tnm1, tn, tnp1, beta, dg, c1, c2, c3, dn
 
    push_sub(bcrmax)
 
    tnm1=h(n-1)-e*s(n-1)
    tn  =h(n  )-e*s(n  )
    tnp1=h(n+1)-e*s(n+1)
    beta=m_one+b/rmax
    dg=a*beta*(dr+m_one-m_half/beta)
 
    c2=24._real64*dg/(12._real64-tn)
    dn=-((24._real64+10._real64*tn)/(12._real64-tn))
 
    c1= (m_one-tnm1/6.0_real64)/(m_one-tnm1/12._real64)
    c3=-((m_one-tnp1/6.0_real64)/(m_one-tnp1/12._real64))
    yn=-((m_one-c1/c3)/(dn-c2/c3))
 
    pop_sub(bcrmax)
  end subroutine bcrmax
 
 
  subroutine numin(e,h,s,y,n,nnode,yn,g,gsg,x,knk)
    real(real64), intent(in)    :: e, h(:), s(:)
    real(real64), intent(inout) :: y(:)
    integer,  intent(in)    :: n
    integer,  intent(out)   :: nnode
    real(real64), intent(in)    :: yn
    real(real64), intent(out)   :: g, gsg, x
    integer,  intent(in)    :: knk
 
    integer :: i
    real(real64) :: t
 
    ! no PUSH SUB, called too often
 
    y(n)=yn
    t=h(n)-e*s(n)
    g=y(n)/(m_one-t/12._real64)
    gsg=g*s(n)*g
    i=n-1
    y(i)=m_one
    t=h(i)-e*s(i)
    g=y(i)/(m_one-t/12._real64)
    gsg=gsg+g*s(i)*g
    x=y(i)-y(n)
    nnode=0
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  begin the inward integration by the Numerov method                         !
!                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    assert(knk > 0)
 
    do i = n - 2, knk, -1
      x = x + t * g
      y(i) = y(i+1) + x
      ! y(i) * y(i+1) can lead to overflows. This is now prevented
      ! by comparing only the signs
      if (sign(m_one,y(i)) * sign(m_one,y(i+1)) < m_zero) nnode = nnode + 1
      t = h(i) - e * s(i)
      g = y(i) / (m_one - t / 12._real64)
      gsg = gsg + g * s(i) * g
    end do
 
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  The kink radius is defined as the point where                              !
!  psi first turns downward.  this usually means at the outermost             !
!  maximum                                                                    !
!                                                                             !
!  the inward integration is now complete                                     !
!                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
  end subroutine numin
 
 
  subroutine numout(e, h, s, y, ncor, knk, nnode, y2, g, gsg, x)
    real(real64), intent(in)    :: e, h(:), s(:)
    real(real64), intent(inout) :: y(:)
    integer,  intent(in)    :: ncor
    integer,  intent(inout) :: knk
    integer,  intent(out)   :: nnode
    real(real64), intent(in)    :: y2
    real(real64), intent(out)   :: g, gsg, x
 
    integer :: i, nm4
    real(real64) :: t, xl
 
    ! no PUSH SUB, called too often
 
    y(1)=m_zero
    y(2)=y2
    t=h(2)-e*s(2)
    g=y(2)/(m_one-t/12._real64)
    gsg=g*s(2)*g
    y(3)=m_one
    t=h(3)-e*s(3)
    g=y(3)/(m_one-t/12._real64)
    gsg=gsg+g*s(3)*g
    x=y(3)-y(2)
    nnode=0
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  begin the outward integration by the Numerov method                        !
!                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
    nm4 = knk - 4
    knk = nm4
    do i = 4, nm4
      xl = x
      x = x + t * g
      y(i) = y(i-1) + x
      ! y(i) * y(i-1) can lead to overflows. This is now prevented
      ! by comparing only the signs
      if (sign(m_one,y(i)) * sign(m_one,y(i-1)) < m_zero) nnode = nnode + 1
      t = h(i) - e * s(i)
      g = y(i) / (m_one - t / 12._real64)
      gsg = gsg + g * s(i) * g
      if (.not. (nnode < ncor .or. xl*x > m_zero)) then
        knk = i
        exit
      end if
    end do
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                             !
!  the outward integration is now complete                                    !
!                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  end subroutine numout
 
end module atomic_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: