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