TABLE OF CONTENTS
ABINIT/psp1nl [ Functions ]
NAME
psp1nl
FUNCTION
Make Kleinman-Bylander form factors f_l(q) for each l from 0 to lmax; Vloc is assumed local potential.
COPYRIGHT
Copyright (C) 1998-2017 ABINIT group (DCA, XG, GMR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
dr(mmax)=inverse of grid spacing for radial grid
lloc=angular momentum of local channel (avoid doing integrals for this l)
lmax=maximum ang momentum for which nonlocal form factor is desired.
mmax=number of radial grid points for atomic grid
mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
mqgrid=number of grid points for q grid
qgrid(mqgrid)=values at which form factors are returned
rad(mmax)=radial grid values
vloc(mmax)=local pseudopotential on radial grid
vpspll(mmax,lmax+1)=nonlocal pseudopotentials for each l on radial grid
wfll(mmax,lmax+1)=reference state wavefunctions on radial grid
wksincos(mmax,2,2)=contains sine and cosine of 2*pi*r(:)*dq and 2*pi*r(:)*q
at input : wksincos(:,1,1)=sine of 2*pi*r(:)*dq
wksincos(:,2,1)=cosine of 2*pi*r(:)*dq
wksincos(:,:,2) is not initialized, will be used inside the routine
OUTPUT
ekb(mpsang)=Kleinman-Bylander energy, {{\ \begin{equation} \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]} {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]} \end{equation} }} for each l ffspl(mqgrid,2,mpsang)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum
NOTES
u_l(r) is reference state wavefunction (input as wfll);
j_l(q) is a spherical Bessel function;
dV_l(r) = vpsp_l(r)-vloc(r) for angular momentum l;
f_l(q) =$ \int_0^{rmax}[j_l(2\pi q r) u_l(r) dV_l(r) r dr]/\sqrt{dvms}$
where dvms=$\displaystyle \int_0^{rmax}[(u_l(r) dV_l(r))^2 dr]$ is the mean
square value of the nonlocal correction for angular momentum l.
E_KB = $\displaystyle \frac{dvms}{\int_0^{rmax}[(u_l(r))^2 dV_l(r) dr]}$.
This is the eigenvalue of the Kleinman-Bylander operator and sets
the energy scale of the nonlocal psp corrections.
Bessel functions replaced by besj, which accomodates args near 0.
PARENTS
psp1in
CHILDREN
besjm,der_int,sincos,spline
SOURCE
64 #if defined HAVE_CONFIG_H 65 #include "config.h" 66 #endif 67 68 #include "abi_common.h" 69 70 subroutine psp1nl(dr,ekb,ffspl,lloc,lmax,mmax,mpsang,mqgrid,& 71 & qgrid,rad,vloc,vpspll,wfll,wksincos) 72 73 use defs_basis 74 use m_errors 75 use m_profiling_abi 76 use m_splines 77 78 use m_special_funcs, only : besjm 79 80 !This section has been created automatically by the script Abilint (TD). 81 !Do not modify the following lines by hand. 82 #undef ABI_FUNC 83 #define ABI_FUNC 'psp1nl' 84 use interfaces_64_psp, except_this_one => psp1nl 85 !End of the abilint section 86 87 implicit none 88 89 !Arguments ------------------------------------ 90 !scalars 91 integer,intent(in) :: lloc,lmax,mmax,mpsang,mqgrid 92 !arrays 93 real(dp),intent(in) :: dr(mmax),qgrid(mqgrid),rad(mmax),vloc(mmax) 94 real(dp),intent(in) :: vpspll(mmax,mpsang),wfll(mmax,mpsang) 95 real(dp),intent(inout) :: wksincos(mmax,2,2) 96 real(dp),intent(out) :: ekb(mpsang),ffspl(mqgrid,2,mpsang) 97 98 !Local variables------------------------------- 99 !scalars 100 integer,parameter :: dpsang=5 101 integer :: iq,ir,irmax,lp1 102 real(dp) :: dvwf,result,test,tpiq,yp1,ypn 103 character(len=500) :: message 104 !arrays 105 real(dp) :: ckb(dpsang),dvms(dpsang),eta(dpsang),renorm(dpsang) 106 real(dp),allocatable :: besjx(:),work1(:),work2(:),work3(:),work4(:),work5(:) 107 real(dp),allocatable :: work_spl(:) 108 109 ! ************************************************************************* 110 111 !DEBUG 112 !write(std_out,*)' psp1nl : enter' 113 !stop 114 !ENDDEBUG 115 116 !Zero out Kleinman-Bylander energies ekb 117 ekb(:)=0.0d0 118 !Zero out eta and other parameters too (so 0 s show up in output later) 119 eta(:)=0.0d0 120 dvms(:)=0.0d0 121 ckb(:)=0.0d0 122 123 !Allow for no nonlocal correction (lmax=-1) 124 if (lmax/=-1) then 125 126 ! Check that lmax is within allowed range 127 if (lmax<0.or.lmax>3) then 128 write(message, '(a,i12,a,a,a,a,a,a,a)' )& 129 & 'lmax=',lmax,' is not an allowed value.',ch10,& 130 & 'Allowed values are -1 for no nonlocal correction or else',ch10,& 131 & '0, 1, 2, or 3 for maximum l nonlocal correction.',ch10,& 132 & 'Action: check the input atomic psp data file for lmax.' 133 MSG_ERROR(message) 134 end if 135 136 ! Compute normalizing integrals eta=<dV> and mean square 137 ! nonlocal psp correction dvms=<dV^2> 138 ! "dvwf" consistently refers to dV(r)*wf(r) where dV=nonlocal correction 139 140 ABI_ALLOCATE(work1,(mmax+1)) 141 ABI_ALLOCATE(work2,(mmax+1)) 142 ABI_ALLOCATE(work_spl,(mqgrid)) 143 ABI_ALLOCATE(work5,(mmax)) 144 ABI_ALLOCATE(besjx,(mmax)) 145 146 do lp1=1,lmax+1 147 148 ! Only do the work if nonlocal correction is nonzero 149 if (lp1 /= lloc+1) then 150 151 ! integrand for 0 to r(mmax) 152 do ir=1,mmax 153 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1) 154 work1(ir)=wfll(ir,lp1)*dvwf 155 end do 156 157 ! do integral 158 ! first need derivative of function; note use of 159 ! shifted indices to accomodate Mike Teter s choice of 0:mmax-1 160 call der_int(work1,work2,rad,dr,mmax-1,result) 161 eta(lp1)=result 162 163 ! DEBUG 164 ! write(std_out,*)' psp1nl : write eta(lp1)' 165 ! write(std_out,*)result 166 ! do ir=1,mmax,61 167 ! write(std_out,*)vpspll(ir,lp1),vloc(ir),wfll(ir,lp1) 168 ! end do 169 ! write(std_out,*) 170 ! do ir=1,mmax,61 171 ! write(std_out,*)work1(ir),rad(ir),dr(ir) 172 ! end do 173 ! ENDDEBUG 174 175 do ir=1,mmax 176 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1) 177 work1(ir)=dvwf**2 178 end do 179 call der_int(work1,work2,rad,dr,mmax-1,result) 180 181 dvms(lp1)=result 182 183 ! If dvms is not 0 for any given angular momentum l, 184 ! compute Xavier Gonze s definition of the Kleinman-Bylander 185 ! energy E_KB = dvms/eta. In this case also renormalize 186 ! the projection operator to u_KB(r)=$u_l(r) dV(r)/\sqrt{dvms}$. 187 ! This means dvwf gets multiplied by the normalization factor 188 ! "renorm"=$1/\sqrt{dvms}$ as seen below. 189 ! With dvwf=dV(r)*wf(r) for wf(r)=``radial'' wf, the integrand 190 ! for each angular momentum l is 191 ! Bessel_l(2 $\pi$ q r) * wf(r) * dV(r) * r; 192 ! NOTE presence of extra r in integrand. 193 194 if (dvms(lp1)/=0.0d0) then 195 ekb(lp1)=dvms(lp1)/eta(lp1) 196 renorm(lp1)=1.0d0/sqrt(dvms(lp1)) 197 ! ckb is Kleinman-Bylander "cosine" (Xavier Gonze) 198 ckb(lp1)=eta(lp1)/sqrt(dvms(lp1)) 199 else 200 ekb(lp1)=0.0d0 201 end if 202 end if 203 end do 204 205 ! Loop on angular momenta 206 do lp1=1,lmax+1 207 208 ! Compute form factor if ekb(lp1) not 0 209 if (ekb(lp1)/=0.0d0) then 210 211 ! do q=0 separately, non-zero if l=0 212 if(lp1==1)then 213 do ir=1,mmax 214 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 215 work1(ir)=rad(ir)*dvwf 216 end do 217 call der_int(work1,work2,rad,dr,mmax-1,result) 218 ffspl(1,1,lp1)=result 219 else 220 ! For l non-zero, f(q=0) vanishes ! 221 ffspl(1,1,lp1)=0.0d0 222 end if 223 224 ! Prepare loop over q values 225 irmax=mmax+1 226 do ir=mmax,2,-1 227 test=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)*rad(ir) 228 work5(ir)=test 229 work1(ir)=0.0d0 230 ! Will ignore tail within decade of machine precision 231 if ((10.0d0+abs(test))==10.0d0 .and. irmax==ir+1) then 232 irmax=ir 233 end if 234 end do 235 ! Increase irmax a bit 236 irmax=irmax+4 237 ! Ask irmax to be lower than mmax 238 if(irmax>mmax-1)irmax=mmax-1 239 240 ABI_ALLOCATE(work3,(irmax-1)) 241 ABI_ALLOCATE(work4,(irmax-1)) 242 243 ! Loop over q values 244 do iq=2,mqgrid 245 tpiq=two_pi*qgrid(iq) 246 call sincos(iq,irmax,mmax,wksincos,rad,tpiq) 247 work3(:)=wksincos(2:irmax,2,2) !Temporary array (Intel compiler compatibility) 248 work4(:)=wksincos(2:irmax,1,2) !Temporary array (Intel compiler compatibility) 249 250 ! Handle r=0 separately 251 work1(1)=0.0d0 252 call besjm(tpiq,besjx(2:irmax),work3,(lp1-1),irmax-1,work4,rad(2:irmax)) 253 do ir=2,irmax 254 work1(ir)=besjx(ir)*work5(ir) 255 end do 256 ! do integral 257 call der_int(work1,work2,rad,dr,irmax,result) 258 ffspl(iq,1,lp1)=result 259 end do 260 261 ! Compute yp1=derivative of f(q) at q=0 262 if(lp1/=2)then 263 ! For l/=1, yp1=0 264 yp1=0.0d0 265 else 266 ! For l=1, yp1=Int [2 Pi r^2 wf(r) dV(r)]/3 267 do ir=1,irmax 268 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 269 work1(ir)=(two_pi*rad(ir)**2)*dvwf/3.0d0 270 end do 271 call der_int(work1,work2,rad,dr,irmax,result) 272 yp1=result 273 end if 274 275 ! Compute ypn=derivative of f(q) at q=qgrid(mqgrid) 276 tpiq=two_pi*qgrid(mqgrid) 277 ! Treat ir=1, r=0, separately 278 work1(1)=0.0d0 279 ! Here, must distinguish l==0 from others 280 if(lp1==1)then 281 ! l==0 : ypn=$\int [2\pi r (-bes1(2\pi r q)) wf(r) dV(r) r dr]$ 282 ! The sine and cosine of the last point were computed in the previous loop 283 ! So, there is no need to call sincos. Note that the rank of besj is 1. 284 call besjm(tpiq,besjx(2:irmax),work3,1,irmax-1,work4,rad(2:irmax)) 285 do ir=2,irmax 286 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 287 work1(ir)=-besjx(ir)*two_pi*rad(ir)*rad(ir)*dvwf 288 end do 289 else 290 ! l==1 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_0(x)-(2/x)j_1(x)) dr]$ 291 ! l==2 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_1(x)-(3/x)j_2(x)) dr]$ 292 ! l==3 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_2(x)-(4/x)j_3(x)) dr]$ 293 ! The sine and cosine of the last point were computed in the previous loop 294 ! Store first previously computed value with besj of order l, then use 295 ! besj of order l-1 (=lp1-2) 296 work1(2:irmax)=besjx(2:irmax) 297 call besjm(tpiq,besjx(2:irmax),work3,(lp1-2),irmax-1,work4,rad(2:irmax)) 298 do ir=2,irmax 299 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 300 work1(ir)=(two_pi*rad(ir)**2)*dvwf*& 301 & ( besjx(ir) - ( dble(lp1)*work1(ir)/(tpiq*rad(ir)) ) ) 302 end do 303 end if 304 ! work1 is ready for integration 305 call der_int(work1,work2,rad,dr,irmax,result) 306 ypn=result 307 308 ! Fit spline to get second derivatives by spline fit 309 call spline(qgrid,ffspl(:,1,lp1),mqgrid,yp1,ypn,& 310 & ffspl(:,2,lp1)) 311 312 ABI_DEALLOCATE(work3) 313 ABI_DEALLOCATE(work4) 314 315 else 316 317 ! KB energy is zero, put nonlocal correction at l=0 to 0 318 ffspl(:,:,lp1)=0.0d0 319 320 end if 321 322 ! End loop on angular momenta 323 end do 324 325 ABI_DEALLOCATE(work1) 326 ABI_DEALLOCATE(work2) 327 ABI_DEALLOCATE(work_spl) 328 ABI_DEALLOCATE(work5) 329 ABI_DEALLOCATE(besjx) 330 331 ! End of lmax/=-1 condition 332 end if 333 334 end subroutine psp1nl