TABLE OF CONTENTS
ABINIT/psp5lo [ Functions ]
NAME
psp5lo
FUNCTION
Compute sine transform to transform from V(r) to q^2 V(q). Computes integrals on logarithmic grid using related uniform grid in exponent and corrected trapezoidal integration.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, FrD) 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
al=spacing in exponent for radial atomic grid. mmax=number of radial r grid points (logarithmic atomic grid). mqgrid=number of grid points in q from 0 to qmax. qgrid(mqgrid)=q grid values (bohr**-1). rad(mmax)=r grid values (bohr). vloc(mmax)=V(r) on radial grid. zion=nominal valence charge of atom.
OUTPUT
epsatm=$ 4\pi\int[r^2 (V(r)+\frac{Zv}{r}dr]$. q2vq(mqgrid) =q^2 V(q) = -\frac{Zv}{\pi} + q^2 4\pi\int[(\frac{\sin(2\pi q r)}{2\pi q r})(r^2 V(r)+r Zv)dr]. yp1,ypn=derivative of q^2 V(q) wrt q at q=0 and q=qmax (needed for spline fitter).
PARENTS
psp5in,psp6in
CHILDREN
ctrap
SOURCE
46 #if defined HAVE_CONFIG_H 47 #include "config.h" 48 #endif 49 50 #include "abi_common.h" 51 52 subroutine psp5lo(al,epsatm,mmax,mqgrid,qgrid,q2vq,rad,& 53 & vloc,yp1,ypn,zion) 54 55 use defs_basis 56 use m_profiling_abi 57 58 use m_numeric_tools, only : ctrap 59 60 !This section has been created automatically by the script Abilint (TD). 61 !Do not modify the following lines by hand. 62 #undef ABI_FUNC 63 #define ABI_FUNC 'psp5lo' 64 !End of the abilint section 65 66 implicit none 67 68 !Arguments---------------------------------------------------------- 69 !scalars 70 integer,intent(in) :: mmax,mqgrid 71 real(dp),intent(in) :: al,zion 72 real(dp),intent(out) :: epsatm,yp1,ypn 73 !arrays 74 real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax) 75 real(dp),intent(out) :: q2vq(mqgrid) 76 77 !Local variables------------------------------- 78 !scalars 79 integer :: iq,ir 80 real(dp),parameter :: scale=10.0d0 81 real(dp) :: arg,result,rmtoin,test,ztor1 82 !arrays 83 real(dp),allocatable :: work(:) 84 85 ! ************************************************************************* 86 87 ABI_ALLOCATE(work,(mmax)) 88 89 !Do q=0 separately (compute epsatm) 90 !Do integral from 0 to r1 91 ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2 92 93 !Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$ 94 !with extra factor of r to convert to uniform grid in exponent 95 do ir=1,mmax 96 ! First handle tail region 97 test=vloc(ir)+zion/rad(ir) 98 ! DEBUG 99 ! write(std_out,*)ir,rad(ir),test 100 ! ENDDEBUG 101 ! Ignore small contributions, or impose a cut-off in the case 102 ! the pseudopotential data are in single precision. 103 ! (it is indeed expected that vloc is very close to zero beyond 20, 104 ! so a value larger than 2.0d-8 is considered anomalous) 105 if (abs(test)<1.0d-20 .or. (rad(ir)>20.0d0 .and. abs(test)>2.0d-8) ) then 106 work(ir)=zero 107 else 108 work(ir)=(rad(ir)*rad(ir))*(rad(ir)*vloc(ir)+zion) 109 end if 110 end do 111 !DEBUG 112 !write(std_out,*)' psp5lo : stop ' 113 !stop 114 !ENDDEBUG 115 116 !Do integral from r(1) to r(max) 117 call ctrap(mmax,work,al,result) 118 !Do integral from r(mmax) to infinity 119 !compute decay length lambda at r(mmax) 120 !$\lambda=-\log((rad(im1)*vloc(im1)+zion)$/ & 121 !$(rad(imat)*vloc(imat)+zion))/(rad(im1)-rad(imat))$ 122 !rmtoin=$(rad(mmax)*vloc(mmax)+zion)*(rad(mmax)+1.d0/\lambda)/\lambda$ 123 !Due to inability to fit exponential decay to r*V(r)+Zv 124 !in tail, NO TAIL CORRECTION IS APPLIED 125 !(numerical trouble might be removed if atomic code is 126 !cleaned up in tail region) 127 rmtoin=0.0d0 128 129 epsatm=4.d0*pi*(result+ztor1+rmtoin) 130 131 q2vq(1)=-zion/pi 132 133 !Loop over q values 134 do iq=2,mqgrid 135 arg=2.d0*pi*qgrid(iq) 136 ! ztor1=$ -Zv/\pi+2q \int_0^{r1}[\sin(2\pi q r)(rV(r)+Zv) dr]$ 137 ztor1=(vloc(1)*sin(arg*rad(1))/arg-(rad(1)*vloc(1)+zion)* & 138 & cos(arg*rad(1)) )/pi 139 140 ! set up integrand 141 do ir=1,mmax 142 test=vloc(ir)+zion/rad(ir) 143 ! Ignore contributions within decade of machine precision 144 if ((scale+abs(test)).eq.scale) then 145 work(ir)=zero 146 else 147 work(ir)=rad(ir)*sin(arg*rad(ir))*(rad(ir)*vloc(ir)+zion) 148 end if 149 end do 150 ! do integral from r(1) to r(mmax) 151 call ctrap(mmax,work,al,result) 152 153 ! do integral from r(mmax) to infinity 154 ! rmtoin=(r(mmax)*vr(mmax)+zion)*(lambda*sin(arg*r(mmax))+ 155 ! arg*cos(arg*r(mmax)))/(arg**2+lambda**2) 156 ! See comment above; no tail correction 157 rmtoin=0.0d0 158 159 ! store q^2 v(q) 160 q2vq(iq)=ztor1+2.d0*qgrid(iq)*(result+rmtoin) 161 162 end do 163 164 !Compute derivatives of q^2 v(q) at ends of interval 165 yp1=0.0d0 166 !ypn=$ 2\int_0^\infty[(\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r)(r V(r)+Z) dr]$ 167 !integral from 0 to r1 168 arg=2.0d0*pi*qgrid(mqgrid) 169 ztor1=zion*rad(1)*sin(arg*rad(1)) 170 ztor1=ztor1+ 3.d0*rad(1)*vloc(1)*cos(arg*rad(1))/arg + & 171 & (rad(1)**2-1.0d0/arg**2)*vloc(1)*sin(arg*rad(1)) 172 !integral from r(mmax) to infinity is overkill; ignore 173 !set up integrand 174 do ir=1,mmax 175 test=vloc(ir)+zion/rad(ir) 176 ! Ignore contributions within decade of machine precision 177 if ((scale+abs(test)).eq.scale) then 178 work(ir)=0.0d0 179 else 180 work(ir)=rad(ir)*(sin(arg*rad(ir))+arg*rad(ir)*cos(arg*rad(ir))) * & 181 & (rad(ir)*vloc(ir)+zion) 182 end if 183 end do 184 call ctrap(mmax,work,al,result) 185 ypn=2.0d0 * (ztor1 + result) 186 187 ABI_DEALLOCATE(work) 188 189 end subroutine psp5lo