1! 2!--------------------------------------------------------------- 3subroutine green(y,lam,e,dvy,chi,vpot,ze2) 4!--------------------------------------------------------------- 5 ! 6 use kinds, only: DP 7 use ld1inc, only: grid 8 use radial_grids, only: ndmx, series 9 implicit none 10 ! 11 ! I/O variables 12 ! 13 integer :: lam ! l angular momentum 14 real(DP) :: y(ndmx), chi(ndmx), dvy(ndmx), vpot(ndmx) 15 real(DP) :: e, ze2 16 ! 17 ! local variables 18 ! 19 integer :: i, l1, ncross, imatch 20 real(DP) :: f(ndmx), g(ndmx) 21 real(DP) :: work(ndmx), int_0_inf_dr 22 real(DP) :: a(0:3), b(0:3), c0, c1, c2, c3, c4, b0e 23 real(DP) :: rr1, rr2, r1, r2 24 real(DP) :: ddx12, sqlhf, xl1, x4l6, x6l12, x8l20 25 real(DP) :: gi, gim1 26 real(DP) :: fac 27! set up constants and initialize 28 ddx12=grid%dx*grid%dx/12.0 29 l1=lam+1 30 sqlhf=(DBLE(lam)+.5)**2 31 xl1=l1 32 x4l6=4*lam+6 33 x6l12=6*lam+12 34 x8l20=8*lam+20 35! series developement of the potential and r.h.s. term near the origin 36 do i=1,4 37 y(i)=vpot(i)-ze2/grid%r(i) 38 end do 39 call series(y,grid%r,grid%r2,b) 40 do i=1,4 41 y(i)= dvy(i)/grid%r(i)**l1 42 end do 43 call series(y,grid%r,grid%r2,a) 44! 45! determine the position of the last change of sign of the f-function 46! f < 0 (approximatively) means classically allowed region 47! f > 0 " " " forbidden " 48! set up the f-, g- and y-functions 49! 50 f(1)=ddx12*(grid%r2(1)*(vpot(1)-e)+sqlhf) 51 do i=2,grid%mesh 52 f(i)=ddx12*(grid%r2(i)*(vpot(i)-e)+sqlhf) 53 if( f(i) .ne. sign(f(i),f(i-1)) ) imatch=i 54 end do 55 56 do i=1,grid%mesh 57 f(i)=1-f(i) 58 g(i)=ddx12*dvy(i)*grid%r(i)*grid%sqr(i) 59 y(i)=0.d0 60 end do 61! wave-function in the first two points by series developement 62 b0e=b(0)-e 63 c0 = 1.0d0 64 c1=0.5*ze2*c0/xl1 65 c2=(c1*ze2+b0e*c0+a(0))/x4l6 66 c3=(c2*ze2+c1*b0e+b(1)*c0+a(1))/x6l12 67 c4=(c3*ze2+c2*b0e+c1*b(1)+b(2)*c0+a(2))/x8l20 68 r1=grid%r(1) 69 r2=grid%r(2) 70 rr1=(c0+r1*(c1+r1*(c2+r1*(c3+r1*c4))))*r1**l1 71 rr2=(c0+r2*(c1+r2*(c2+r2*(c3+r2*c4))))*r2**l1 72 73 y(1)=rr1/grid%sqr(1) 74 y(2)=rr2/grid%sqr(2) 75! 76! new set up of the g- function 77 gim1 = g(1) 78 do i=2,grid%mesh-1 79 gi=g(i) 80 g(i)=gim1+10.0d0*gi+g(i+1) 81 gim1=gi 82 end do 83! outward integration: numerov's algorithm 84 call outward(y,f,g,grid%mesh,imatch,ncross) 85 86! inward integration: froese's algorithm 87 call inward(y,f,g,grid%mesh,imatch) 88 89!-orthogonalize to the 0^th order eigenvector 90 do i=1,grid%mesh 91 y(i) = y(i)*grid%sqr(i) 92 work(i)=y(i)*chi(i) 93 end do 94 95 fac = int_0_inf_dr(work,grid,grid%mesh,2*lam+2) 96 do i=1,grid%mesh 97 y(i) = y(i) - fac*chi(i) 98 end do 99! 100 return 101 102end subroutine green 103