TABLE OF CONTENTS
ABINIT/psp2lo [ Functions ]
NAME
psp2lo
FUNCTION
Treat local part of Goedecker-Teter-Hutter pseudopotentials (pspcod=2), as well as Hartwigsen-Goedecker-Hutter pseudopotentials (pspcod=3)
COPYRIGHT
Copyright (C) 1998-2018 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
cc1,2,3,4=parameters from analytic pseudopotential form
mqgrid=number of grid points in q from 0 to qmax.
qgrid(mqgrid)=values of q (or G) on grid from 0 to qmax (bohr^-1)
if vlspl_recipSpace is .true. else values of r on grid from
0 to 2pi / qmax * mqgrid_ff (bohr).
rloc=local pseudopotential core radius (bohr)
vlspl_recipSpace= .true. if computation of vlspl is done in reciprocal space
zion=valence charge of atom
parameters for local potential: rloc,c1,c2,c3,c4
OUTPUT
dvloc(mqgrid)=dVloc(r)/dr (only allocated if vlspl_recipSpace is false).
epsatm=$4\pi\int[r^2 (v(r)+\frac{Zv}{r} dr]$
q2vq(mqgrid)&=&q^2 v(q) \nonumber \\
&=&-Zv/\pi
+q^2 4\pi\int[(\frac{\sin(2\pi qr)}{2\pi qr})(r^2 v(r)+r Zv)dr]\nonumber\\
&=&\exp(-K^2*rloc^2/2) \nonumber \\
&& *(-\frac{zion}{\pi}+(\frac{K^2*rloc^3}{\sqrt{2*\pi}}*
(c1+c2*(3-(rloc*K)^2) \nonumber \\
&& +c3*(15-10(rloc*K)^2+(rloc*K)^4) \nonumber \\
&& +c4*(105-105*(rloc*K)^2+21*(rloc*K)^4-(rloc*K)^6)) \nonumber
for GTH vloc with $K=(2\pi q)$.
yp1,ypn=derivative of q^2 v(q) wrt q at q=0 and q=qmax
(needed for spline fitter).
PARENTS
psp10in,psp2in,psp3in
CHILDREN
wrtout
SOURCE
53 #if defined HAVE_CONFIG_H 54 #include "config.h" 55 #endif 56 57 #include "abi_common.h" 58 59 60 subroutine psp2lo(cc1,cc2,cc3,cc4,dvloc,epsatm,mqgrid,qgrid,q2vq,& 61 & rloc,vlspl_recipSpace,yp1,ypn,zion) 62 63 use defs_basis 64 use m_profiling_abi 65 66 use m_special_funcs, only : abi_derfc 67 68 !This section has been created automatically by the script Abilint (TD). 69 !Do not modify the following lines by hand. 70 #undef ABI_FUNC 71 #define ABI_FUNC 'psp2lo' 72 use interfaces_14_hidewrite 73 !End of the abilint section 74 75 implicit none 76 77 !Arguments ------------------------------------ 78 !scalars 79 integer,intent(in) :: mqgrid 80 real(dp),intent(in) :: cc1,cc2,cc3,cc4,rloc,zion 81 real(dp),intent(out) :: epsatm,yp1,ypn 82 logical,intent(in) :: vlspl_recipSpace 83 !arrays 84 real(dp),intent(in) :: qgrid(mqgrid) 85 real(dp),intent(out) :: dvloc(mqgrid),q2vq(mqgrid) 86 87 !Local variables------------------------------- 88 !scalars 89 integer :: iqgrid 90 real(dp) :: erfValue,gaussValue,polyValue,qmax,rq,rq2 91 character(len=500) :: message 92 93 ! ************************************************************************* 94 95 !Compute epsatm = lim(q->0) [Vloc(q) + zion/(Pi*q^2)] 96 epsatm=2.d0*pi*rloc**2*zion+(2.d0*pi)**(1.5d0)*rloc**3*& 97 & (cc1+3.d0*cc2+15.d0*cc3+105.d0*cc4) 98 99 !If vlspl_recipSpace is .true., we compute V(q)*q^2 in reciprocal space, 100 !else we compute V(r) in real space. 101 if (vlspl_recipSpace) then 102 write(message, '(a)' ) '- Local part computed in reciprocal space.' 103 call wrtout(ab_out,message,'COLL') 104 call wrtout(std_out, message,'COLL') 105 106 ! d(q^2*V(q))/d(q) at q=0 and q=qmax 107 qmax=qgrid(mqgrid) 108 rq2=(2.d0*pi*qmax*rloc)**2 109 yp1=0.d0 110 ypn= (2.d0*pi*qmax*rloc**2)*exp(-0.5d0*rq2)* & 111 & (2.d0*zion + sqrt(2.d0*pi)*rloc*& 112 & (cc1*(2.d0-rq2) + cc2*(6.d0-7.d0*rq2+rq2**2) +& 113 & cc3*(30.d0-55.d0*rq2+16.d0*rq2**2-rq2**3) +& 114 & cc4*(210.d0-525.d0*rq2+231.d0*rq2**2-29.d0*rq2**3+rq2**4))) 115 ! ypn has been tested against Maple-derived expression. 116 117 ! Compute q^2*vloc(q) on uniform grid 118 do iqgrid=1,mqgrid 119 rq2=(2.d0*pi*qgrid(iqgrid)*rloc)**2 120 q2vq(iqgrid)=exp(-0.5d0*rq2)*(-zion/pi+rq2*(rloc/sqrt(2.d0*pi)) *& 121 & ( cc1 + cc2*(3.d0-rq2) + cc3*(15.d0-10.d0*rq2+rq2**2) +& 122 & cc4*(105.d0-rq2*(105.d0-rq2*(21.d0-rq2))) )) 123 end do 124 else 125 write(message, '(a)' ) '- Local part computed in real space.' 126 call wrtout(ab_out,message,'COLL') 127 call wrtout(std_out, message,'COLL') 128 129 ! Compute derivatives for splines computations 130 yp1 = 0.d0 131 rq2 = (qgrid(mqgrid) / rloc) ** 2 132 erfValue = abi_derfc(sqrt(0.5d0 * rq2)) 133 ypn = - 2.0d0 * zion / sqrt(2.d0 * pi) / qgrid(mqgrid) / rloc 134 ypn = ypn - rq2 * (cc1 + cc2 * rq2 + cc3 * rq2 ** 2 + cc4 * rq2 ** 3) / qgrid(mqgrid) 135 ypn = ypn + (2.d0 * cc2 * rq2 + 4.d0 * cc3 * rq2 ** 2 + 6.d0 * cc4 * rq2 ** 3) / qgrid(mqgrid) 136 ypn = ypn * exp(-0.5d0 * rq2) 137 ypn = ypn + zion / qgrid(mqgrid) ** 2 * erfValue 138 ! Note that ypn has been calculated on a full-proof a4 paper sheet. 139 140 ! Compute local potential and its first derivatives. 141 do iqgrid = 1, mqgrid, 1 142 rq2 = (qgrid(iqgrid) / rloc) ** 2 143 ! Compute erf() part 144 ! Case r = 0 145 gaussValue = exp(-0.5d0 * rq2) 146 if (qgrid(iqgrid) == 0.d0) then 147 q2vq(iqgrid) = -zion / rloc * sqrt(2.d0 / pi) 148 dvloc(iqgrid) = 0.d0 149 else 150 erfValue = abi_derfc(sqrt(0.5d0 * rq2)) 151 q2vq(iqgrid) = -zion / qgrid(iqgrid) * (1.0d0 - erfValue) 152 dvloc(iqgrid) = - sqrt(2.d0 / pi) * zion * gaussValue / (qgrid(iqgrid) * rloc) - & 153 & q2vq(iqgrid) / qgrid(iqgrid) 154 end if 155 ! Add the gaussian part 156 polyValue = cc1 + cc2 * rq2 + cc3 * rq2 ** 2 + cc4 * rq2 ** 3 157 q2vq(iqgrid) = q2vq(iqgrid) + gaussValue * polyValue 158 rq = qgrid(iqgrid) / rloc 159 dvloc(iqgrid) = dvloc(iqgrid) - qgrid(iqgrid) / rloc ** 2 * gaussValue * polyValue + & 160 & gaussValue * (2.0d0 * cc2 * rq / rloc + 3.0d0 * cc3 * rq ** 3 / rloc + & 161 & 6.0d0 * cc4 * rq ** 5 / rloc) 162 end do 163 164 write(message, '(a,f12.7,a,a,f12.7,a,a,a,f12.7)' ) & 165 & ' | dr spline step is : ', qgrid(2), ch10, & 166 & ' | r > ', qgrid(mqgrid) ,' is set to 0.', ch10, & 167 & ' | last non-nul potential value is : ', q2vq(mqgrid) 168 call wrtout(ab_out,message,'COLL') 169 call wrtout(std_out, message,'COLL') 170 end if 171 172 end subroutine psp2lo