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!--------------------------------------------------------------------- 10subroutine integrate_outward (lam,jam,e,mesh,ndm,grid,f, & 11 b,y,beta,ddd,qq,nbeta,nwfx,lls,jjs,ikk,ik) 12 !--------------------------------------------------------------------- 13 ! 14 ! Integrate the wavefunction from 0 to r(ik) 15 ! generalized separable or US pseudopotentials are allowed 16 ! This routine assumes that y countains already the 17 ! correct values in the first two points 18 ! 19 use kinds, only : DP 20 use radial_grids, only : radial_grid_type 21 implicit none 22 type(radial_grid_type):: grid 23 integer :: & 24 lam, & ! l angular momentum 25 mesh, & ! size of radial mesh 26 ndm, & ! maximum radial mesh 27 nbeta, & ! number of beta function 28 nwfx, & ! maximum number of beta functions 29 lls(nbeta),&! for each beta the angular momentum 30 ikk(nbeta),&! for each beta the integration point 31 ik ! the last integration point 32 33 real(DP) :: & 34 e, & ! output eigenvalue 35 jam, & ! j angular momentum 36 f(mesh), & ! the f function 37 b(0:3), & ! the taylor expansion of the potential 38 y(mesh), & ! the output solution 39 jjs(nwfx), & ! the j angular momentum 40 beta(ndm,nwfx),& ! the beta functions 41 ddd(nwfx,nwfx),qq(nwfx,nwfx) ! parameters for computing B_ij 42 43 integer :: & 44 nst, & ! the exponential around the origin 45 n, & ! counter on mesh points 46 iib,jjb, & ! counter on beta with correct lam 47 ierr, & ! used to control allocation 48 ib,jb, & ! counter on beta 49 info ! info on exit of LAPACK subroutines 50 51 integer, allocatable :: iwork(:) ! auxiliary space 52 53 real(DP) :: & 54 b0e, & ! the expansion of the known part 55 ddx12, & ! the deltax entering the equations 56 x4l6, & ! auxiliary for small r expansion 57 j1(4),d(4),& ! auxiliary for starting values of chi 58 delta,xc(4),& ! auxiliary for starting values of eta 59 int_0_inf_dr ! the integral function 60 61 real(DP), allocatable :: & 62 el(:), & ! auxiliary for integration 63 cm(:,:), &! the linear system 64 bm(:), & ! the known part of the linear system 65 c(:), & ! the chi functions 66 coef(:), & ! the solution of the linear system 67 eta(:,:) ! the partial solution of the nonomogeneous 68 69 70 allocate(c(ik), stat=ierr) 71 allocate(el(ik), stat=ierr) 72 allocate(cm(nbeta,nbeta), stat=ierr) 73 allocate(bm(nbeta), stat=ierr) 74 allocate(coef(nbeta), stat=ierr) 75 allocate(iwork(nbeta), stat=ierr) 76 allocate(eta(ik,nbeta), stat=ierr) 77 78 if (mesh.ne.grid%mesh) call errore('integrate_outward','mesh dimension is not as expected',1) 79 ddx12=grid%dx*grid%dx/12.0_DP 80 b0e=b(0)-e 81 x4l6=4*lam+6 82 nst=(lam+1)*2 83 ! 84 ! first solve the homogeneous equation 85 ! 86 ! f is the original function 87 ! of the form 1 + h^2/12 * d2Rdr2 88 do n=2,ik-1 89 y(n+1)=( 12.0_DP*y(n) - 10.0_DP* f(n)*y(n) - f(n-1)*y(n-1) )/f(n+1) 90 enddo 91 ! 92 ! for each beta function with correct angular momentum 93 ! solve the inhomogeneous equation 94 ! 95 iib=0 96 jjb=0 97 do ib=1,nbeta 98 ! 99 ! set up the known part 100 ! 101 if (lls(ib).eq.lam.and.jjs(ib).eq.jam) then 102 iib=iib+1 103 c=0.0_DP 104 do jb=1,nbeta 105 if (lls(jb).eq.lam.and.jjs(jb).eq.jam) then 106 do n=1,ikk(jb) 107 c(n)= c(n)+(ddd(jb,ib) -e*qq(jb,ib))*beta(n,jb) 108 enddo 109 endif 110 enddo 111 ! 112 ! compute the starting values of the solutions 113 ! 114 do n=1,4 115 j1(n)=c(n)/grid%r(n)**(lam+1) 116 enddo 117 call seriesbes(j1,grid%r,grid%r2,4,d) 118 delta=b0e**2+x4l6*b(2) 119 xc(1)=(-d(1)*b0e-x4l6*d(3))/delta 120 xc(3)=(-b0e*d(3)+d(1)*b(2))/delta 121 xc(2)=0.0_DP 122 xc(4)=0.0_DP 123 do n=1,3 124 eta(n,iib)=grid%r(n)**(lam+1)*(xc(1)+grid%r2(n)*xc(3))/grid%sqr(n) 125 enddo 126 127 do n=1,ik 128 c(n)=c(n)*grid%r2(n)/grid%sqr(n) 129 enddo 130 ! 131 ! solve the inhomogeneous equation 132 ! 133 do n=3,ik-1 134 eta(n+1,iib)=((12.0_DP-10.0_DP*f(n))*eta(n,iib) & 135 -f(n-1)*eta(n-1,iib) & 136 + ddx12*(10.0_DP*c(n)+c(n-1)+c(n+1)) )/f(n+1) 137 enddo 138 ! 139 ! compute the coefficents of the linear system 140 ! 141 jjb=0 142 do jb=1,nbeta 143 if (lls(jb).eq.lam.and.jjs(jb).eq.jam) then 144 jjb=jjb+1 145 do n=1,min(ik,ikk(jb)) 146 el(n)=beta(n,jb)*eta(n,iib)*grid%sqr(n) 147 enddo 148 cm(jjb,iib)=-int_0_inf_dr(el,grid,min(ik,ikk(jb)),nst) 149 endif 150 enddo 151 152 do n=1,min(ik,ikk(ib)) 153 el(n)=beta(n,ib)*y(n)*grid%sqr(n) 154 enddo 155 bm(iib)=int_0_inf_dr(el,grid,min(ik,ikk(ib)),nst) 156 cm(iib,iib)=1.0_DP+cm(iib,iib) 157 endif 158 enddo 159 if (iib.ne.jjb) call errore('integrate_outward','jjb.ne.iib',1) 160 161 if (iib.gt.0) then 162 call dcopy(iib,bm,1,coef,1) 163 164 call DGESV(iib,1,cm,nbeta,iwork,coef,nbeta,info) 165 166 if (info /= 0) call errore('integrate_outward', & 167 & 'problems solving the linear system',info) 168 169 do ib=1,iib 170 do n=1,ik 171 y(n)= y(n)+coef(ib)*eta(n,ib) 172 enddo 173 enddo 174 endif 175 176 deallocate(eta) 177 deallocate(iwork) 178 deallocate(coef) 179 deallocate(bm) 180 deallocate(cm) 181 deallocate(el) 182 deallocate(c) 183 184 return 185end subroutine integrate_outward 186