TABLE OF CONTENTS
ABINIT/psp2in [ Functions ]
NAME
psp2in
FUNCTION
Initialize pspcod=2 pseudopotentials (GTH format): continue to read the file, then compute the corresponding local and non-local potentials.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, MT) 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
dtset <type(dataset_type)>=all input variables in this dataset ipsp=id in the array of the pseudo-potential. lmax=value of lmax mentioned at the second line of the psp file zion=nominal valence of atom as specified in psp file
OUTPUT
ekb(lnmax)=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,n) epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+Zv/r) dr]$ (hartree) ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum and each projector indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0) or i=lmn (if useylm=1) nproj(mpsang)=number of projection functions for each angular momentum vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit dvlspl(mqgrid_vl,2)=dVloc(r)/dr and second derivatives from spline fit (only allocated if vlspl_recipSpace is false.
SIDE EFFECTS
Input/output
lmax : at input =value of lmax mentioned at the second line of the psp file
at output= 1
psps <type(pseudopotential_type)>=at output, values depending on the read
pseudo are set.
| lmnmax(IN)=if useylm=1, max number of (l,m,n) comp. over all type of psps
| =if useylm=0, max number of (l,n) comp. over all type of psps
| lnmax(IN)=max. number of (l,n) components over all type of psps
| angular momentum of nonlocal pseudopotential
| mpsang(IN)= 1+maximum angular momentum for nonlocal pseudopotentials
| mqgrid_ff(IN)=dimension of q (or G) grid for nl form factors (array ffspl)
| mqgrid_vl(IN)=dimension of q (or G) grid or r grid (if vlspl_recipSpace = .false.)
| qgrid_ff(mqgrid_ff)(IN)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors
| qgrid_vl(mqgrid_vl)(IN)=values of q on grid from 0 to qmax (bohr^-1) for Vloc
| if vlspl_recipSpace is .true. else values of r on grid from
| 0 to 2pi / qmax * mqgrid_ff (bohr).
| useylm(IN)=governs the way the nonlocal operator is to be applied:
| 1=using Ylm, 0=using Legendre polynomials
| vlspl_recipSpace(IN)=.true. if pseudo are expressed in reciprocal space.
| gth_params(OUT)=store GTH coefficients and parameters.
PARENTS
pspatm
CHILDREN
psp2lo,psp2nl,spline,wrtout,wvl_descr_psp_fill
SOURCE
72 #if defined HAVE_CONFIG_H 73 #include "config.h" 74 #endif 75 76 #include "abi_common.h" 77 78 subroutine psp2in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps,vlspl,dvlspl,zion) 79 80 use defs_basis 81 use defs_datatypes 82 use defs_abitypes 83 use m_splines 84 use m_profiling_abi 85 use m_errors 86 #if defined HAVE_BIGDFT 87 use BigDFT_API, only: atomic_info 88 #endif 89 90 !This section has been created automatically by the script Abilint (TD). 91 !Do not modify the following lines by hand. 92 #undef ABI_FUNC 93 #define ABI_FUNC 'psp2in' 94 use interfaces_14_hidewrite 95 use interfaces_43_wvl_wrappers 96 use interfaces_64_psp, except_this_one => psp2in 97 !End of the abilint section 98 99 implicit none 100 101 !Arguments ------------------------------------ 102 !scalars 103 integer,intent(in) :: ipsp,lmax 104 real(dp),intent(in) :: zion 105 real(dp),intent(out) :: epsatm 106 type(dataset_type),intent(in) :: dtset 107 type(pseudopotential_type),intent(inout) :: psps 108 !arrays 109 integer,intent(out) :: indlmn(6,psps%lmnmax),nproj(psps%mpsang) 110 real(dp),intent(out) :: dvlspl(psps%mqgrid_vl,2),ekb(psps%lnmax) 111 real(dp),intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) !vz_i 112 real(dp),intent(out) :: vlspl(psps%mqgrid_vl,2) 113 114 !Local variables------------------------------- 115 !scalars 116 integer :: iln,index,ipsang,kk,ll,mm 117 real(dp) :: cc1,cc2,cc3,cc4,h1p,h1s,h2s,rloc,rrp,rrs 118 real(dp) :: yp1,ypn 119 character(len=500) :: message,errmsg 120 !arrays 121 real(dp),allocatable :: work_space(:),work_spl(:) 122 real(dp),allocatable :: dvloc(:) 123 124 ! *************************************************************************** 125 126 !Set various terms to 0 in case not defined below 127 !GTH values 128 rloc=0.d0 129 cc1=0.d0 130 cc2=0.d0 131 cc3=0.d0 132 cc4=0.d0 133 rrs=0.d0 134 h1s=0.d0 135 h2s=0.d0 136 rrp=0.d0 137 h1p=0.d0 138 nproj(1:psps%mpsang)=0 139 140 !Read and write different lines of the pseudopotential file 141 read (tmp_unit,*, err=10, iomsg=errmsg) rloc,cc1,cc2,cc3,cc4 142 write(message, '(a,f12.7)' ) ' rloc=',rloc 143 call wrtout(ab_out,message,'COLL') 144 call wrtout(std_out, message,'COLL') 145 write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )& 146 & ' cc1=',cc1,'; cc2=',cc2,'; cc3=',cc3,'; cc4=',cc4 147 call wrtout(ab_out,message,'COLL') 148 call wrtout(std_out, message,'COLL') 149 150 read (tmp_unit,*, err=10, iomsg=errmsg) rrs,h1s,h2s 151 write(message, '(a,f12.7,a,f12.7,a,f12.7)' )& 152 & ' rrs=',rrs,'; h1s=',h1s,'; h2s=',h2s 153 call wrtout(ab_out,message,'COLL') 154 call wrtout(std_out, message,'COLL') 155 156 read (tmp_unit,*, err=10, iomsg=errmsg) rrp,h1p 157 write(message, '(a,f12.7,a,f12.7)' )& 158 & ' rrp=',rrp,'; h1p=',h1p 159 call wrtout(ab_out,message,'COLL') 160 call wrtout(std_out, message,'COLL') 161 162 !Store the coefficients. 163 psps%gth_params%set(ipsp) = .true. 164 psps%gth_params%psppar(0, :, ipsp) = (/ rloc, cc1, cc2, cc3, cc4, 0.d0, 0.d0 /) 165 psps%gth_params%psppar(1, :, ipsp) = (/ rrs, h1s, h2s, 0.d0, 0.d0, 0.d0, 0.d0 /) 166 psps%gth_params%psppar(2, :, ipsp) = (/ rrp, h1p, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0 /) 167 if (dtset%usewvl == 1) then 168 call wvl_descr_psp_fill(psps%gth_params, ipsp, 0, int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), tmp_unit) 169 end if 170 171 if (abs(h1s)>1.d-08) nproj(1)=1 172 if (abs(h2s)>1.d-08) nproj(1)=2 173 174 if (abs(h1p)>1.d-08) then 175 if(psps%mpsang<2)then 176 write(message, '(a,es12.4,a,a,a,i2,a)' )& 177 & 'With non-zero h1p (=',h1p,'), mpsang should be at least 2,',ch10,& 178 & 'while mpsang=',psps%mpsang,'.' 179 MSG_ERROR(message) 180 end if 181 nproj(2)=1 182 if (lmax<1) then 183 write(message, '(a,i5,a,e12.4,a,a,a,a)' )& 184 & 'Input lmax=',lmax,' disagree with input h1p=',h1p,'.',& 185 & 'Your pseudopotential is incoherent.',ch10,& 186 & 'Action : correct your pseudopotential file.' 187 MSG_ERROR(message) 188 end if 189 end if 190 191 !Initialize array indlmn array giving l,m,n,lm,ln,s for i=lmn 192 index=0;iln=0;indlmn(:,:)=0 193 do ipsang=1,lmax+1 194 if(nproj(ipsang)>0)then 195 ll=ipsang-1 196 do kk=1,nproj(ipsang) 197 iln=iln+1 198 do mm=1,2*ll*psps%useylm+1 199 index=index+1 200 indlmn(1,index)=ll 201 indlmn(2,index)=mm-ll*psps%useylm-1 202 indlmn(3,index)=kk 203 indlmn(4,index)=ll*ll+(1-psps%useylm)*ll+mm 204 indlmn(5,index)=iln 205 indlmn(6,index)=1 206 end do 207 end do 208 end if 209 end do 210 211 !First, the local potential -- 212 !compute q^2V(q) or V(r) 213 !MJV NOTE: psp2lo should never be called with dvspl unallocated, which 214 !is possible unless .not.psps%vlspl_recipSpace 215 ABI_ALLOCATE(dvloc,(psps%mqgrid_vl)) 216 call psp2lo(cc1,cc2,cc3,cc4,dvloc,epsatm,psps%mqgrid_vl,psps%qgrid_vl,& 217 & vlspl(:,1),rloc,psps%vlspl_recipSpace,yp1,ypn,zion) 218 219 !DEBUG 220 !write(std_out,*)' psp2in : after psp2lo ' 221 !stop 222 !ENDDEBUG 223 224 !Fit spline to (q^2)V(q) or V(r) 225 ABI_ALLOCATE(work_space,(psps%mqgrid_vl)) 226 ABI_ALLOCATE(work_spl,(psps%mqgrid_vl)) 227 call spline (psps%qgrid_vl,vlspl(:,1),psps%mqgrid_vl,yp1,ypn,work_spl) 228 vlspl(:,2)=work_spl(:) 229 if (.not.psps%vlspl_recipSpace) then 230 dvlspl(:,1) = dvloc 231 call spline (psps%qgrid_vl,dvlspl(:,1),psps%mqgrid_vl,yp1,ypn,work_spl) 232 dvlspl(:,2)=work_spl(:) 233 end if 234 ABI_DEALLOCATE(work_space) 235 ABI_DEALLOCATE(work_spl) 236 ABI_DEALLOCATE(dvloc) 237 238 239 !Second, compute KB energies and form factors and fit splines 240 ekb(:)=0.0d0 241 !First check if any nonlocal projectors are being used 242 if (maxval(nproj(1:lmax+1))>0) then 243 call psp2nl(ekb,ffspl,h1p,h1s,h2s,psps%lnmax,psps%mqgrid_ff,psps%qgrid_ff,rrp,rrs) 244 end if 245 246 return 247 248 ! Handle IO error 249 10 continue 250 MSG_ERROR(errmsg) 251 252 end subroutine psp2in
ABINIT/psp2nl [ Functions ]
NAME
psp2nl
FUNCTION
Goedecker-Teter-Hutter nonlocal pseudopotential (from preprint of 1996). Uses Gaussians for fully nonlocal form, analytic expressions.
INPUTS
h1p=factor defining strength of 1st projector for l=1 channel h1s=factor defining strength of 1st projector for l=0 channel h2s=factor defining strength of 2nd projector for l=0 channel lnmax=max. number of (l,n) components over all type of psps mqgrid=number of grid points for qgrid qgrid(mqgrid)=array of |G| values rrp=core radius for p channel (bohr) rrs=core radius for s channel (bohr)
OUTPUT
ekb(lnmax)=Kleinman-Bylander energy ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum and each projector
PARENTS
psp2in
CHILDREN
spline
SOURCE
287 subroutine psp2nl(ekb,ffspl,h1p,h1s,h2s,lnmax,mqgrid,qgrid,rrp,rrs) 288 289 use defs_basis 290 use m_splines 291 use m_profiling_abi 292 293 !This section has been created automatically by the script Abilint (TD). 294 !Do not modify the following lines by hand. 295 #undef ABI_FUNC 296 #define ABI_FUNC 'psp2nl' 297 !End of the abilint section 298 299 implicit none 300 301 !Arguments ------------------------------------ 302 !scalars 303 integer,intent(in) :: lnmax,mqgrid 304 real(dp),intent(in) :: h1p,h1s,h2s,rrp,rrs 305 !arrays 306 real(dp),intent(in) :: qgrid(mqgrid) 307 real(dp),intent(inout) :: ekb(lnmax),ffspl(mqgrid,2,lnmax) !vz_i 308 309 !Local variables------------------------------- 310 !scalars 311 integer :: iln,iqgrid 312 real(dp) :: qmax,yp1,ypn 313 !arrays 314 real(dp),allocatable :: work(:) 315 316 ! ************************************************************************* 317 318 ABI_ALLOCATE(work,(mqgrid)) 319 320 !Kleinman-Bylander energies ekb were set to zero in calling program 321 322 !Compute KB energies 323 iln=0 324 if (abs(h1s)>1.d-12) then 325 iln=iln+1 326 ekb(iln)=h1s*32.d0*rrs**3*(pi**(2.5d0)/(4.d0*pi)**2) 327 end if 328 if (abs(h2s)>1.d-12) then 329 iln=iln+1 330 ekb(iln) =h2s*(128.d0/15.d0)*rrs**3*(pi**(2.5d0)/(4.d0*pi)**2) 331 end if 332 if (abs(h1p)>1.d-12) then 333 iln=iln+1 334 ekb(iln)=h1p*(64.d0/3.d0)*rrp**5*(pi**(2.5d0)/(4.d0*pi)**2) 335 end if 336 337 !Compute KB form factor 338 iln=0 339 340 !l=0 first projector 341 if (abs(h1s)>1.d-12) then 342 iln=iln+1 343 do iqgrid=1,mqgrid 344 ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) 345 end do 346 ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid) 347 yp1=0.d0 348 qmax=qgrid(mqgrid) 349 ypn=-4.d0*pi**2*qmax*rrs**2*exp(-0.5d0*(two_pi*qmax*rrs)**2) 350 ! Fit spline to get second derivatives by spline fit 351 call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln)) 352 ! else 353 ! or else put first projector nonlocal correction at l=0 to 0 354 ! ffspl(:,:,iln)=0.0d0 355 end if 356 357 !l=0 second projector 358 if (abs(h2s)>1.d-12) then 359 iln=iln+1 360 do iqgrid=1,mqgrid 361 ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) * & 362 & (3.d0-(two_pi*qgrid(iqgrid)*rrs)**2) 363 end do 364 ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid) 365 yp1=0.d0 366 qmax=qgrid(mqgrid) 367 ypn=4.d0*pi**2*qmax*rrs**2*exp(-0.5d0*(two_pi*qmax*rrs)**2) * & 368 & (-5.d0+(two_pi*qmax*rrs)**2) 369 ! Fit spline to get second derivatives by spline fit 370 call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln)) 371 ! else if(mproj>=2)then 372 ! or else put second projector nonlocal correction at l=0 to 0 373 ! ffspl(:,:,iln)=0.0d0 374 end if 375 376 !l=1 first projector 377 if (abs(h1p)>1.d-12) then 378 iln=iln+1 379 do iqgrid=1,mqgrid 380 ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * & 381 & (two_pi*qgrid(iqgrid)) 382 end do 383 ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid) 384 yp1=two_pi 385 qmax=qgrid(mqgrid) 386 ypn=-two_pi*((two_pi*qmax*rrp)**2-1.d0) * exp(-0.5d0*(two_pi*qmax*rrp)**2) 387 ! Fit spline to get second derivatives by spline fit 388 call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln)) 389 ! else if(mpsang>=2)then 390 ! or else put first projector l=1 nonlocal correction to 0 391 ! ffspl(:,:,iln)=0.0d0 392 end if 393 394 ABI_DEALLOCATE(work) 395 396 end subroutine psp2nl