1!
2! Copyright (C) 2004 PWSCF group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!
9!---------------------------------------------------------------
10         subroutine intref(lam,e,mesh,grid,vpot,ze2,chi)
11!---------------------------------------------------------------
12!
13!  numerical integration of the radial schroedinger equation
14!  computing logarithmic derivatives
15!  thresh dermines the absolute accuracy for the eigenvalue
16!
17!
18!
19      use kinds, only : DP
20      use radial_grids, only: radial_grid_type, series
21      implicit none
22      type(radial_grid_type), intent(in):: grid
23      integer :: &
24              mesh,    &     ! the mesh size
25              lam           ! the angular momentum
26      real(DP) :: &
27              vpot(mesh), &   ! the local potential
28              chi(mesh),  &   ! the solution
29              ze2,       &   ! the nuclear charge in Ry units
30              e             ! the eigenvalue
31
32      integer :: &
33              ierr,  &      ! used to control allocation
34              n             ! generic counter
35
36      real(DP) :: &
37              lamsq,   &     ! combined angular momentum
38              b(0:3),c(4)      ! used for starting guess of the solution
39
40      real(DP),allocatable :: &
41              al(:)      ! the known part of the differential equation
42
43      if (mesh.gt.grid%mesh) &
44          call errore('intref','mesh dimension is too large',1)
45
46      allocate(al(mesh),stat=ierr)
47
48      do n=1,4
49         al(n)=vpot(n)-ze2/grid%r(n)
50      enddo
51      call series(al,grid%r,grid%r2,b)
52
53
54      lamsq=(lam+0.5_DP)**2
55!
56!     b) find the value of solution s in the first two points
57!
58      call start_scheq( lam, e, b, grid, ze2, chi )
59
60      do n=1,mesh
61         al(n)=( (vpot(n)-e)*grid%r2(n) + lamsq )*grid%dx**2/12.0_DP
62         al(n)=1.0_DP-al(n)
63      enddo
64!
65!     Integrate forward the equation:
66!     c) integrate the equation from 0 to matching radius
67!
68      do n=2,mesh-1
69         chi(n+1)=((12.0_DP-10.0_DP*al(n))*chi(n) -al(n-1)*chi(n-1))/al(n+1)
70      enddo
71!
72!     set the correct r dependence
73!
74      do n=1,mesh
75         chi(n)= chi(n)*grid%sqr(n)
76      enddo
77
78
79      deallocate(al)
80      return
81      end subroutine intref
82