TABLE OF CONTENTS
ABINIT/pspatm [ Functions ]
NAME
pspatm
FUNCTION
Open atomic pseudopotential data file for a given atom, read the three first lines, make some checks, then call appropriate subroutine for the reading of V(r) and wf R(r) data for each angular momentum, and subsequent Fourier and Bessel function transforms for local and nonlocal potentials. Close psp file at end. Handles pseudopotential files produced by (pspcod=1 or 4) Teter code, or from the Goedecker-Teter-Hutter paper (pspcod=2), or from the Hartwigsen-Goedecker-Hutter paper (pspcod=3 or 10) or "Phoney pseudopotentials" (Hamman grid in real space) (pspcod=5) or "Troullier-Martins pseudopotentials" from the FHI (pspcod=6) or "XML format" (pspcod=9) or "UPF PWSCF format" (pspcod=11)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, FrD, AF, DRH) 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
dq= spacing of the q-grid dtset <type(dataset_type)>=all input variables in this dataset | ixc=exchange-correlation choice from main routine data file | pawxcdev=choice of XC development in PAW formalism | usexcnhat_orig=choice for use of nhat in Vxc in PAW formalism | xclevel= XC functional level ipsp=id in the array of the currently read pseudo.
OUTPUT
ekb(dimekb)=
->NORM-CONSERVING PSPS ONLY (pspcod/=7):
(Real) Kleinman-Bylander energies (hartree)
{{\ \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 number of basis functions (l,n) (dimekb=lnmax)
If any, spin-orbit components begin at l=mpsang+1
epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
or i=lmn (if useylm=1)
pawrad <type(pawrad_type)>=paw radial mesh and related data
pawtab <type(pawtab_type)>=paw tabulated starting data
vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
ffspl(mqgrid_ff,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
second derivative from spline fit for each angular momentum and
each projector; if any, spin-orbit components begin at l=mpsang+1
xcccrc=XC core correction cutoff radius (bohr) from psp file
xccc1d(n1xccc*(1-usepaw),6)=1D core charge function and five derivatives, from psp file (used in NC only)
nctab=<nctab_t>
has_tvale=True if the pseudo provides the valence density (used in NC only)
tvalespl(mqgrid_vl(1-usepaw),2)=the pseudo valence density and 2nd derivative in reciprocal space on a regular grid
(used in NC only)
SIDE EFFECTS
Input/Output : psps <type(pseudopotential_type)>=at output, values depending on the read pseudo are set. | dimekb(IN)=dimension of ekb (see module defs_psp.f) | filpsp(IN)=name of formatted external file containing atomic psp data. | 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 | mpssoang(IN)= 1+maximum (spin*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 for Vloc (array vlspl) | n1xccc(IN)=dimension of xccc1d ; 0 if no XC core correction is used | optnlxccc(IN)=option for nl XC core correction | positron(IN)=0 if electron GS calculation | 1 if positron GS calculation | 2 if electron GS calculation in presence of the positron | pspso(INOUT)=spin-orbit characteristics, govern the content of ffspl and ekb | if =0 : this input requires NO spin-orbit characteristics of the psp | if =2 : this input requires HGH characteristics of the psp | if =3 : this input requires HFN characteristics of the psp | if =1 : this input will be changed at output to 1, 2, 3, according | to the intrinsic characteristics of the psp file | 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 | usepaw(IN)= 0 for non paw calculation; =1 for paw calculation | 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. | znuclpsp(IN)=atomic number of atom as specified in input file to main routine
NOTES
Format expected for the three first lines of pseudopotentials (1) title (character) line (2) znucl,zion,pspdat (3) pspcod,pspxc,lmax,lloc,mmax,r2well Dimensions of form factors and Vloc q grids must be the same in Norm-Conserving case
PARENTS
pspini
CHILDREN
nctab_eval_tcorespl,pawpsp_17in,pawpsp_7in,pawpsp_bcast pawpsp_read_header_xml,pawpsp_read_pawheader,pawpsp_wvl,psp10in,psp1in psp2in,psp3in,psp5in,psp6in,psp8in,psp9in,psp_dump_outputs psxml2abheader,test_xml_xmlpaw_upf,timab,upf2abinit,wrtout wvl_descr_psp_fill
SOURCE
117 #if defined HAVE_CONFIG_H 118 #include "config.h" 119 #endif 120 121 #include "abi_common.h" 122 123 subroutine pspatm(dq,dtset,dtfil,ekb,epsatm,ffspl,indlmn,ipsp,pawrad,pawtab,& 124 & psps,vlspl,dvlspl,xcccrc,xccc1d,nctab,comm_mpi) 125 126 use defs_basis 127 use defs_datatypes 128 use defs_abitypes 129 use m_profiling_abi 130 use m_errors 131 use m_xmpi 132 use m_psps 133 134 use m_io_tools, only : open_file 135 use m_pawrad, only : pawrad_type 136 use m_pawtab, only : pawtab_type 137 use m_pawpsp, only : pawpsp_bcast, pawpsp_read_pawheader, pawpsp_read_header_xml,& 138 & pawpsp_header_type, pawpsp_wvl, pawpsp_7in, pawpsp_17in 139 use m_pawxmlps, only : paw_setup, ipsp2xml 140 use m_psxml2ab 141 #if defined HAVE_BIGDFT 142 use BigDFT_API, only : dictionary, atomic_info, dict_init, dict_free, UNINITIALIZED 143 #endif 144 145 !This section has been created automatically by the script Abilint (TD). 146 !Do not modify the following lines by hand. 147 #undef ABI_FUNC 148 #define ABI_FUNC 'pspatm' 149 use interfaces_14_hidewrite 150 use interfaces_18_timing 151 use interfaces_43_wvl_wrappers 152 use interfaces_64_psp, except_this_one => pspatm 153 !End of the abilint section 154 155 implicit none 156 157 !Arguments --------------------------------------------- 158 !scalars 159 integer,intent(in) :: ipsp 160 integer, optional,intent(in) :: comm_mpi 161 real(dp),intent(in) :: dq 162 real(dp),intent(out) :: epsatm,xcccrc 163 type(dataset_type),intent(in) :: dtset 164 type(datafiles_type),intent(in) :: dtfil 165 type(pawrad_type),intent(inout) :: pawrad 166 type(pawtab_type),intent(inout) :: pawtab 167 type(nctab_t),intent(inout) :: nctab 168 type(pseudopotential_type),intent(inout) :: psps 169 !arrays 170 integer,intent(out) :: indlmn(6,psps%lmnmax) 171 real(dp),intent(out) :: dvlspl(psps%mqgrid_vl,2) 172 real(dp),intent(inout) :: ekb(psps%dimekb*(1-psps%usepaw)) 173 real(dp),intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) 174 real(dp),intent(out) :: vlspl(psps%mqgrid_vl,2) 175 real(dp),intent(inout) :: xccc1d(psps%n1xccc*(1-psps%usepaw),6) 176 177 !Local variables --------------------------------------- 178 !scalars 179 integer :: ii,il,ilmn,iln,iln0,lloc,lmax,me,mmax 180 integer :: paral_mode,pspcod,pspdat,pspxc,useupf,usexml,xmlpaw,unt 181 real(dp) :: maxrad,qchrg,r2well,zion,znucl 182 character(len=500) :: message,errmsg 183 character(len=fnlen) :: title 184 character(len=fnlen) :: filnam 185 type(pawpsp_header_type):: pawpsp_header 186 !arrays 187 integer,allocatable :: nproj(:) 188 real(dp) :: tsec(2) 189 real(dp),allocatable :: e990(:),e999(:),ekb1(:),ekb2(:),epspsp(:),rcpsp(:) 190 real(dp),allocatable :: rms(:) 191 #if defined HAVE_PSML 192 !! usexml= 0 for non xml ps format ; =1 for xml ps format 193 character(len=3) :: atmsymb 194 character(len=30) :: creator 195 type(pspheader_type) :: psphead 196 #endif 197 198 ! ****************************************************************************** 199 200 !paral_mode defines how we access to the psp file 201 ! paral_mode=0: all processes access to the file (sequentially) 202 ! paral_mode=1: only proc. 0 access to the file and then broadcast 203 paral_mode=0 204 if (present(comm_mpi)) then 205 if (psps%usepaw==1.and.xmpi_comm_size(comm_mpi)>1) paral_mode=1 206 end if 207 me=0;if (paral_mode==1) me=xmpi_comm_rank(comm_mpi) 208 209 if (paral_mode == 1) then 210 ABI_CHECK(psps%usepaw==1, "paral_mode==1 is only compatible with PAW, see call to pawpsp_bcast below") 211 end if 212 213 nctab%has_tvale = .False.; nctab%has_tcore = .False. 214 215 if (me==0) then 216 ! Dimensions of form factors and Vloc q grids must be the same in Norm-Conserving case 217 if (psps%usepaw==0 .and. psps%mqgrid_ff/=psps%mqgrid_vl) then 218 write(message, '(a,a,a,a,a)' )& 219 & 'Dimension of q-grid for nl form factors (mqgrid_ff)',ch10,& 220 & 'is different from dimension of q-grid for Vloc (mqgrid_vl) !',ch10,& 221 & 'This is not allowed for norm-conserving psp.' 222 MSG_ERROR(message) 223 end if 224 225 write(message, '(a,t38,a)' )'- pspatm: opening atomic psp file',trim(psps%filpsp(ipsp)) 226 call wrtout(ab_out, message,'COLL') 227 call wrtout(std_out, message,'COLL') 228 !write(message, "(2a)")"- md5: ",trim(psps%md5_pseudos(ipsps)) 229 230 ! Check if the file pseudopotential file is written in (XML| XML-PAW | UPF) 231 call test_xml_xmlpaw_upf(psps%filpsp(ipsp), usexml, xmlpaw, useupf) 232 233 ! ---------------------------------------------------------------------------- 234 ! allocate nproj here: can be read in now for UPF 235 ABI_ALLOCATE(nproj,(psps%mpssoang)) 236 nproj(:)=0 237 238 if (usexml /= 1 .and. useupf /= 1) then 239 240 ! Open the atomic data file, and read the three first lines 241 ! These three first lines have a similar format in all allowed psp files 242 243 ! Open atomic data file (note: formatted input file) 244 if (open_file(psps%filpsp(ipsp), message, unit=tmp_unit, form='formatted', status='old') /= 0) then 245 MSG_ERROR(message) 246 end if 247 rewind (unit=tmp_unit,err=10,iomsg=errmsg) 248 249 ! Read and write some description of file from first line (character data) 250 read (tmp_unit,'(a)',err=10,iomsg=errmsg) title 251 write(message, '(a,a)' ) '- ',trim(title) 252 call wrtout(ab_out,message,'COLL') 253 call wrtout(std_out, message,'COLL') 254 255 ! Read and write more data describing psp parameters 256 read (tmp_unit,*,err=10,iomsg=errmsg) znucl,zion,pspdat 257 write(message,'(a,f9.5,f10.5,2x,i8,t47,a)')'-',znucl,zion,pspdat,'znucl, zion, pspdat' 258 call wrtout(ab_out,message,'COLL') 259 call wrtout(std_out,message,'COLL') 260 261 read (tmp_unit,*,err=10,iomsg=errmsg) pspcod,pspxc,lmax,lloc,mmax,r2well 262 if(pspxc<0) then 263 write(message, '(i5,i8,2i5,i10,f10.5,t47,a)' ) & 264 & pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well' 265 else 266 write(message, '(4i5,i10,f10.5,t47,a)' ) & 267 & pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well' 268 end if 269 call wrtout(ab_out,message,'COLL') 270 call wrtout(std_out, message,'COLL') 271 272 else if (usexml == 1 .and. xmlpaw == 0) then 273 274 ! the following is probably useless - already read in everything in inpspheads 275 #if defined HAVE_PSML 276 ! write(message,'(a,a)') & 277 !& '- pspatm: Reading pseudopotential header in XML form from ', trim(psps%filpsp(ipsp)) 278 ! call wrtout(ab_out,message,'COLL') 279 ! call wrtout(std_out, message,'COLL') 280 281 call psxml2abheader( psps%filpsp(ipsp), psphead, atmsymb, creator, 0 ) 282 znucl = psphead%znuclpsp 283 zion = psphead%zionpsp 284 pspdat = psphead%pspdat 285 pspcod = psphead%pspcod 286 pspxc = psphead%pspxc 287 lmax = psphead%lmax 288 !lloc = 0 ! does this mean s? in psml case the local potential can be different from any l channel 289 lloc = -1 290 mmax = -1 291 r2well = 0 292 293 write(message,'(a,1x,a3,3x,a)') "-",atmsymb,trim(creator) 294 call wrtout(ab_out,message,'COLL') 295 call wrtout(std_out,message,'COLL') 296 write(message,'(a,f9.5,f10.5,2x,i8,t47,a)')'-',znucl,zion,pspdat,'znucl, zion, pspdat' 297 call wrtout(ab_out,message,'COLL') 298 call wrtout(std_out,message,'COLL') 299 if(pspxc<0) then 300 write(message, '(i5,i8,2i5,i10,f10.5,t47,a)' ) & 301 & pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well' 302 else 303 write(message, '(4i5,i10,f10.5,t47,a)' ) & 304 & pspcod,pspxc,lmax,lloc,mmax,r2well,'pspcod,pspxc,lmax,lloc,mmax,r2well' 305 end if 306 call wrtout(ab_out,message,'COLL') 307 call wrtout(std_out,message,'COLL') 308 #else 309 write(message,'(a,a)') & 310 & 'ABINIT is not compiled with XML support for reading this type of pseudopotential ', & 311 & trim(psps%filpsp(ipsp)) 312 MSG_BUG(message) 313 #endif 314 ! END useless 315 else if (usexml == 1 .and. xmlpaw == 1) then 316 write(message,'(a,a)') & 317 & '- pspatm : Reading pseudopotential header in XML form from ', trim(psps%filpsp(ipsp)) 318 call wrtout(ab_out,message,'COLL') 319 call wrtout(std_out, message,'COLL') 320 321 !PENDING: These routines will be added to libpaw: 322 ! Return header informations 323 call pawpsp_read_header_xml(lloc,lmax,pspcod,& 324 & pspxc,paw_setup(ipsp2xml(ipsp)),r2well,zion,znucl) 325 ! Fill in pawpsp_header object: 326 call pawpsp_read_pawheader(pawpsp_header%basis_size,& 327 & lmax,pawpsp_header%lmn_size,& 328 & pawpsp_header%l_size,pawpsp_header%mesh_size,& 329 & pawpsp_header%pawver,paw_setup(ipsp2xml(ipsp)),& 330 & pawpsp_header%rpaw,pawpsp_header%rshp,pawpsp_header%shape_type) 331 332 else if (useupf == 1) then 333 if (psps%usepaw /= 0) then 334 MSG_ERROR("UPF format not allowed with PAW (USPP part not read yet)") 335 end if 336 337 pspcod = 11 338 r2well = 0 339 340 ! should initialize znucl,zion,pspxc,lmax,lloc,mmax 341 call upf2abinit (psps%filpsp(ipsp), znucl, zion, pspxc, lmax, lloc, mmax, & 342 & psps, epsatm, xcccrc, indlmn, ekb, ffspl, nproj, vlspl, xccc1d) 343 344 else 345 MSG_ERROR("You should not be here! erroneous type or pseudopotential input") 346 end if 347 348 ! ------------------------------------------------------------------------------ 349 ! Check data for consistency against main routine input 350 351 ! Does required spin-orbit characteristics agree with format 352 ! TODO: in case of pspcod 5 (phoney) and 8 (oncvpsp) this is not specific enough. 353 ! they can be non-SOC as well. 354 ! HGH is ok - can always turn SOC on or off. 355 ! write(std_out,*) pspso 356 if((pspcod/=3).and.(pspcod/=5).and.(pspcod/=8).and.(pspcod/=10))then 357 ! If pspso requires internal characteristics, set it to 1 for non-HGH psps 358 if(psps%pspso(ipsp)==1) psps%pspso(ipsp)=0 359 if(psps%pspso(ipsp)/=0)then 360 write(message, '(3a,i0,3a)' )& 361 & 'Pseudopotential file cannot give spin-orbit characteristics,',ch10,& 362 & 'while pspso(itypat)= ',psps%pspso(ipsp),'.',ch10,& 363 & 'Action: check your pseudopotential and input files for consistency.' 364 MSG_ERROR(message) 365 end if 366 end if 367 368 ! Does nuclear charge znuclpsp agree with psp input znucl 369 !write(std_out,*)znucl,ipsp,psps%znuclpsp(ipsp) 370 !MGNAG: v5[66] gives NAG in %znuclpsp if -nan 371 if (abs(psps%znuclpsp(ipsp)-znucl)>tol8) then 372 write(message, '(a,f10.5,2a,f10.5,5a)' )& 373 & 'Pseudopotential file znucl=',znucl,ch10,& 374 & 'does not equal input znuclpsp=',psps%znuclpsp(ipsp),' better than 1e-08 .',ch10,& 375 & 'znucl is read from the psp file in pspatm, while',ch10,& 376 & 'znuclpsp is read in iofn2.' 377 MSG_BUG(message) 378 end if 379 380 ! Is the highest angular momentum within limits? 381 ! Recall mpsang is 1+highest l for nonlocal correction. 382 ! Nonlocal corrections for s, p, d, and f are supported. 383 if (lmax+1>psps%mpsang) then 384 write(message, '(a,i0,a,i0,a,a)' )& 385 & 'input lmax+1= ',lmax+1,' exceeds mpsang= ',psps%mpsang,ch10,& 386 & 'indicates input lmax too large for dimensions.' 387 MSG_BUG(message) 388 end if 389 390 ! Check several choices for ixc against pspxc 391 ! ixc is from ABINIT code; pspxc is from atomic psp file 392 if (dtset%ixc==0) then 393 MSG_WARNING('Note that input ixc=0 => no xc is being used.') 394 else if(dtset%ixc/=pspxc) then 395 write(message, '(a,i0,3a,i0,9a)' )& 396 & 'Pseudopotential file pspxc= ',pspxc,',',ch10,& 397 & 'not equal to input ixc= ',dtset%ixc,'.',ch10,& 398 & 'These parameters must agree to get the same xc ',ch10,& 399 & 'in ABINIT code as in psp construction.',ch10,& 400 & 'Action: check psp design or input file.',ch10,& 401 & 'Assume experienced user. Execution will continue.' 402 MSG_WARNING(message) 403 end if 404 405 if (lloc>lmax .and. pspcod/=4 .and. pspcod/=8 .and. pspcod/=10) then 406 write(message, '(a,2i12,a,a,a,a)' )& 407 & 'lloc,lmax=',lloc,lmax,ch10,& 408 & 'chosen l of local psp exceeds range from input data.',ch10,& 409 & 'Action: check pseudopotential input file.' 410 MSG_ERROR(message) 411 end if 412 413 ! Does the pspcod agree with type of calculation (paw or not)? 414 if (((pspcod/=7.and.pspcod/=17).and.psps%usepaw==1).or.((pspcod==7.or.pspcod==17).and.psps%usepaw==0)) then 415 write(message, '(a,i2,a,a,i0,a)' )& 416 & 'In reading atomic psp file, finds pspcod= ',pspcod,ch10,& 417 & 'This is not an allowed value with usepaw= ',psps%usepaw,'.' 418 MSG_BUG(message) 419 end if 420 421 if (.not.psps%vlspl_recipSpace .and. & 422 & (pspcod /= 2 .and. pspcod /= 3 .and. pspcod /= 10 .and. pspcod /= 7)) then 423 ! The following "if" statement can substitute the one just before once libBigDFT 424 ! has been upgraded to include pspcod 10 425 ! if (.not.psps%vlspl_recipSpace .and. (pspcod /= 2 .and. pspcod /= 3 .and. pspcod /= 10)) then 426 write(message, '(a,i2,a,a)' )& 427 & 'In reading atomic psp file, finds pspcod=',pspcod,ch10,& 428 & 'This is not an allowed value with real space computation.' 429 MSG_BUG(message) 430 end if 431 432 ! MJV 16/6/2009 added pspcod 11 for upf format 433 if( pspcod<1 .or. (pspcod>11.and.pspcod/=17) ) then 434 write(message, '(a,i0,4a)' )& 435 & 'In reading atomic psp file, finds pspcod= ',pspcod,ch10,& 436 & 'This is not an allowed value. Allowed values are 1 to 11 .',ch10,& 437 & 'Action: check pseudopotential input file.' 438 MSG_ERROR(message) 439 end if 440 441 ! ----------------------------------------------------------------------- 442 ! Set various terms to 0 in case not defined below 443 ABI_ALLOCATE(e990,(psps%mpssoang)) 444 ABI_ALLOCATE(e999,(psps%mpssoang)) 445 ABI_ALLOCATE(rcpsp,(psps%mpssoang)) 446 ABI_ALLOCATE(rms,(psps%mpssoang)) 447 ABI_ALLOCATE(epspsp,(psps%mpssoang)) 448 ABI_ALLOCATE(ekb1,(psps%mpssoang)) 449 ABI_ALLOCATE(ekb2,(psps%mpssoang)) 450 e990(:)=zero ;e999(:)=zero 451 rcpsp(:)=zero;rms(:)=zero 452 ekb1(:)=zero ;ekb2(:)=zero 453 epspsp(:)=zero 454 qchrg=0 455 456 ! ---------------------------------------------------------------------- 457 if(pspcod==1 .or. pspcod==4)then 458 459 ! Teter pseudopotential (pspcod=1 or 4) 460 call psp1in(dq,ekb,ekb1,ekb2,epsatm,epspsp,& 461 & e990,e999,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,& 462 & mmax,psps%mpsang,psps%mqgrid_ff,nproj,psps%n1xccc,pspcod,qchrg,psps%qgrid_ff,& 463 & rcpsp,rms,psps%useylm,vlspl,xcccrc,xccc1d,zion,psps%znuclpsp(ipsp)) 464 465 else if (pspcod==2)then 466 467 ! GTH pseudopotential 468 call psp2in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps,vlspl,dvlspl,zion) 469 xccc1d(:,:)=0.0d0 ; qchrg=0.0d0 ; xcccrc=0.0d0 470 471 else if (pspcod==3)then 472 473 ! HGH pseudopotential 474 call psp3in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps, psps%pspso(ipsp), & 475 & vlspl,zion) 476 xccc1d(:,:)=0.0d0 ; qchrg=0.0d0 ; xcccrc=0.0d0 477 478 else if (pspcod==5)then 479 480 ! Old phoney pseudopotentials 481 call psp5in(ekb,ekb1,ekb2,epsatm,epspsp,& 482 & e990,e999,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,& 483 & mmax,psps%mpsang,psps%mpssoang,psps%mqgrid_ff,nproj,psps%n1xccc,psps%pspso(ipsp),qchrg,psps%qgrid_ff,& 484 & rcpsp,rms,psps%useylm,vlspl,xcccrc,xccc1d,zion,psps%znuclpsp(ipsp)) 485 486 else if (pspcod==6)then 487 488 ! FHI pseudopotentials 489 call psp6in(ekb,epsatm,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,mmax,& 490 & psps%mpsang,psps%mqgrid_ff,nproj,psps%n1xccc,psps%optnlxccc,psps%positron,qchrg,psps%qgrid_ff,psps%useylm,vlspl,& 491 & xcccrc,xccc1d,zion,psps%znuclpsp(ipsp)) 492 493 else if (pspcod==7)then 494 ! PAW "pseudopotentials" 495 call pawpsp_7in(epsatm,ffspl,dtset%icoulomb,dtset%ixc,& 496 & lmax,psps%lnmax,mmax,psps%mqgrid_ff,psps%mqgrid_vl,& 497 & pawrad,pawtab,dtset%pawxcdev,psps%qgrid_ff,psps%qgrid_vl,& 498 & dtset%usewvl,dtset%usexcnhat_orig,vlspl,xcccrc,dtset%xclevel,& 499 & dtset%xc_denpos,zion,psps%znuclpsp(ipsp)) 500 501 else if (pspcod==8)then 502 503 ! DRH pseudopotentials 504 call psp8in(ekb,epsatm,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,mmax,& 505 & psps%mpsang,psps%mpssoang,psps%mqgrid_ff,psps%mqgrid_vl,nproj,psps%n1xccc,psps%pspso(ipsp),& 506 & qchrg,psps%qgrid_ff,psps%qgrid_vl,psps%useylm,vlspl,xcccrc,xccc1d,zion,psps%znuclpsp(ipsp),nctab,maxrad) 507 508 #if defined DEV_YP_DEBUG_PSP 509 call psp_dump_outputs("DBG",pspcod,psps%lmnmax,psps%lnmax,psps%mpssoang, & 510 & psps%mqgrid_ff,psps%n1xccc,mmax,maxrad,epsatm,qchrg,xcccrc,nctab, & 511 & indlmn,nproj,ekb,ffspl,vlspl,xccc1d) 512 #endif 513 514 else if (pspcod==9)then 515 516 #if defined HAVE_PSML 517 call psp9in(psps%filpsp(ipsp),ekb,epsatm,ffspl,indlmn,lloc,lmax,psps%lmnmax,psps%lnmax,mmax,& 518 & psps%mpsang,psps%mpssoang,psps%mqgrid_ff,psps%mqgrid_vl,nproj,psps%n1xccc, & 519 & psps%pspso(ipsp),qchrg,psps%qgrid_ff,psps%qgrid_vl,psps%useylm,vlspl,& 520 & xcccrc,xccc1d,zion,psps%znuclpsp(ipsp),nctab,maxrad) 521 522 #if defined DEV_YP_DEBUG_PSP 523 call psp_dump_outputs("DBG",pspcod,psps%lmnmax,psps%lnmax,psps%mpssoang, & 524 & psps%mqgrid_ff,psps%n1xccc,mmax,maxrad,epsatm,qchrg,xcccrc,nctab, & 525 & indlmn,nproj,ekb,ffspl,vlspl,xccc1d) 526 #endif 527 #else 528 write(message,'(2a)') & 529 & 'ABINIT is not compiled with XML support for reading this type of pseudopotential ', & 530 & trim(psps%filpsp(ipsp)) 531 MSG_BUG(message) 532 #endif 533 534 else if (pspcod==10)then 535 536 ! HGH pseudopotential, full h/k matrix read 537 call psp10in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps, psps%pspso(ipsp), & 538 & vlspl,zion) 539 xccc1d(:,:)=0.0d0 ; qchrg=0.0d0 ; xcccrc=0.0d0 540 541 ! NB for pspcod 11 the reading has already been done above. 542 else if (pspcod==17)then 543 ! PAW XML "pseudopotentials" 544 call pawpsp_17in(epsatm,ffspl,dtset%icoulomb,ipsp,dtset%ixc,lmax,& 545 & psps%lnmax,mmax,psps%mqgrid_ff,psps%mqgrid_vl,pawpsp_header,pawrad,pawtab,& 546 & dtset%pawxcdev,psps%qgrid_ff,psps%qgrid_vl,dtset%usewvl,& 547 & dtset%usexcnhat_orig,vlspl,xcccrc,& 548 & dtset%xclevel,dtset%xc_denpos,zion,psps%znuclpsp(ipsp)) 549 end if 550 551 close (unit=tmp_unit) 552 553 ! ---------------------------------------------------------------------- 554 if (pspcod==2 .or. pspcod==3 .or. pspcod==10)then 555 write(message, '(a,a,a,a,a,a,a,a,a,a)' )ch10,& 556 & ' pspatm : COMMENT -',ch10,& 557 & ' the projectors are not normalized,',ch10,& 558 & ' so that the KB energies are not consistent with ',ch10,& 559 & ' definition in PRB44, 8503 (1991). ',ch10,& 560 & ' However, this does not influence the results obtained hereafter.' 561 call wrtout(ab_out,message,'COLL') 562 call wrtout(std_out,message,'COLL') 563 ! The following lines are added to keep backward compatibilty 564 maxrad=zero 565 #if defined HAVE_BIGDFT 566 do ii=1,size(psps%gth_params%psppar,1)-1 ! psppar first dim begins at 0 567 if (psps%gth_params%psppar(ii,0,ipsp)/=zero) maxrad=max(maxrad,psps%gth_params%psppar(ii,0,ipsp)) 568 end do 569 if (abs(maxrad)<=tol12) then 570 psps%gth_params%radii_cf(ipsp,3)=zero 571 else 572 psps%gth_params%radii_cf(ipsp,3)=max( & 573 & min(dtset%wvl_crmult*psps%gth_params%radii_cf(ipsp,1),15._dp*maxrad)/dtset%wvl_frmult, & 574 & psps%gth_params%radii_cf(ipsp,2)) 575 end if 576 #endif 577 end if 578 579 if (pspcod/=7.and.pspcod/=17) then 580 write(message, '(a,f14.8,a,a)' ) ' pspatm : epsatm=',epsatm,ch10,& 581 & ' --- l ekb(1:nproj) -->' 582 call wrtout(ab_out,message,'COLL') 583 call wrtout(std_out, message,'COLL') 584 iln0=0 585 do ilmn=1,psps%lmnmax 586 iln=indlmn(5,ilmn) 587 if (iln>iln0) then 588 il=indlmn(1,ilmn) 589 if (indlmn(6,ilmn)==1) then 590 iln0=iln0+nproj(il+1) 591 !if (dtset%optdriver == RUNL_SIGMA) then 592 ! do ii=0,nproj(il+1)-1 593 ! ekb(iln+ii) = zero 594 ! end do 595 !end if 596 write(message, '(13x,i1,4f12.6)' ) il,(ekb(iln+ii),ii=0,nproj(il+1)-1) 597 else 598 iln0=iln0+nproj(il+psps%mpsang) 599 write(message, '(2x,a,i1,4f12.6)' ) 'spin-orbit ',il,(ekb(iln+ii),ii=0,nproj(il+psps%mpsang)-1) 600 end if 601 call wrtout(ab_out,message,'COLL') 602 call wrtout(std_out,message,'COLL') 603 end if 604 end do 605 end if 606 607 ! NC: Evalute spline-fit of the model core charge in reciprocal space. 608 ! TODO: Be careful, because we will be using the PAW part in which tcore is always avaiable! 609 ! Should add a test with 2 NC pseudos: one with NLCC and the other without! 610 if (psps%usepaw == 0) then 611 call nctab_eval_tcorespl(nctab, psps%n1xccc, xcccrc, xccc1d, psps%mqgrid_vl, psps%qgrid_vl) 612 end if 613 614 write(message,'(3a)') ' pspatm: atomic psp has been read ',' and splines computed',ch10 615 call wrtout(ab_out,message,'COLL') 616 call wrtout(std_out,message,'COLL') 617 618 ABI_DEALLOCATE(e990) 619 ABI_DEALLOCATE(e999) 620 ABI_DEALLOCATE(rcpsp) 621 ABI_DEALLOCATE(rms) 622 ABI_DEALLOCATE(ekb1) 623 ABI_DEALLOCATE(ekb2) 624 ABI_DEALLOCATE(epspsp) 625 ABI_DEALLOCATE(nproj) 626 627 if (dtset%prtvol > 9 .and. psps%usepaw==0 .and. psps%lmnmax>3) then 628 629 write (filnam, '(a,i0,a)') trim(dtfil%fnameabo_pspdata), ipsp, ".dat" 630 if (open_file(filnam, message, newunit=unt) /= 0) then 631 MSG_ERROR(message) 632 end if 633 write (unt,*) '# Pseudopotential data in reciprocal space as used by ABINIT' 634 write (unt,'(a)', ADVANCE='NO') '# index vlocal ' 635 if (psps%lnmax > 0) & 636 & write (unt,'(a,I3)', ADVANCE='NO') ' 1st proj(l=', indlmn(1,1) 637 if (psps%lnmax > 1) & 638 & write (unt,'(a,I3)', ADVANCE='NO') ') 2nd(l=', indlmn(1,2) 639 if (psps%lnmax > 2) & 640 & write (unt,'(a,I3,a)', ADVANCE='NO') ') 3rd(l=', indlmn(1,3), ')' 641 write (unt,*) 642 643 do ii = 1, psps%mqgrid_vl 644 write(unt, '(I5,E24.16)', ADVANCE='NO') ii, vlspl(ii,1) 645 if (psps%lnmax > 0) write(unt, '(E24.16)', ADVANCE='NO') ffspl(ii,1,1) 646 if (psps%lnmax > 1) write(unt, '(E24.16)', ADVANCE='NO') ffspl(ii,1,2) 647 if (psps%lnmax > 2) write(unt, '(E24.16)', ADVANCE='NO') ffspl(ii,1,3) 648 write(unt, *) 649 end do 650 close(unt) 651 652 write (filnam, '(a,i0,a)') trim(dtfil%fnameabo_nlcc_derivs), ipsp, ".dat" 653 if (open_file(filnam, message, newunit=unt) /= 0) then 654 MSG_ERROR(message) 655 end if 656 write (unt,*) '# Non-linear core corrections' 657 write (unt,*) '# r, pseudocharge, 1st, 2nd, 3rd, 4th, 5th derivatives' 658 do ii = 1, psps%n1xccc 659 write (unt,*) xcccrc*(ii-1)/(psps%n1xccc-1), xccc1d(ii,1), xccc1d(ii,2), & 660 & xccc1d(ii,3), xccc1d(ii,4), & 661 & xccc1d(ii,5), xccc1d(ii,6) 662 end do 663 close(unt) 664 end if 665 666 end if !me=0 667 668 if (paral_mode==1) then 669 call timab(48,1,tsec) 670 call pawpsp_bcast(comm_mpi,epsatm,ffspl,pawrad,pawtab,vlspl,xcccrc) 671 call timab(48,2,tsec) 672 end if 673 674 if (psps%usepaw==1) then 675 indlmn(:,:)=0 676 indlmn(1:6,1:pawtab%lmn_size)=pawtab%indlmn(1:6,1:pawtab%lmn_size) 677 end if 678 679 !-------------------------------------------------------------------- 680 !WVL+PAW: 681 if (dtset%usepaw==1 .and. (dtset%icoulomb /= 0 .or. dtset%usewvl==1)) then 682 #if defined HAVE_BIGDFT 683 psps%gth_params%psppar(:,:,ipsp) = UNINITIALIZED(1._dp) 684 psps%gth_params%radii_cf(ipsp,:) = UNINITIALIZED(1._dp) 685 call wvl_descr_psp_fill(psps%gth_params, ipsp, psps%pspxc(1), int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), 0) 686 #endif 687 688 ! The following lines are added to keep backward compatibilty 689 maxrad=zero 690 #if defined HAVE_BIGDFT 691 do ii=1,size(psps%gth_params%psppar,1)-1 ! psppar first dim begins at 0 692 if (psps%gth_params%psppar(ii,0,ipsp)/=zero) maxrad=max(maxrad,psps%gth_params%psppar(ii,0,ipsp)) 693 end do 694 if (abs(maxrad)<=tol12) then 695 !== MT COMMENT 696 ! Damien wants to activate this (in order to directly compare to bigDFT): 697 psps%gth_params%radii_cf(ipsp,3)= psps%gth_params%radii_cf(ipsp,2) 698 ! But, this changes strongly file references. 699 ! So, I keep this, waiting for Tonatiuh s validation 700 psps%gth_params%radii_cf(ipsp,3)= & 701 & (psps%gth_params%radii_cf(ipsp,1)+psps%gth_params%radii_cf(ipsp,2))*half 702 !== MT COMMENT 703 else 704 psps%gth_params%radii_cf(ipsp,3)=max( & 705 & min(dtset%wvl_crmult*psps%gth_params%radii_cf(ipsp,1),15._dp*maxrad)/dtset%wvl_frmult, & 706 & psps%gth_params%radii_cf(ipsp,2)) 707 end if 708 if(present(comm_mpi)) then 709 call pawpsp_wvl(psps%filpsp(ipsp),pawrad,pawtab,dtset%usewvl,dtset%wvl_ngauss,comm_mpi) 710 else 711 call pawpsp_wvl(psps%filpsp(ipsp),pawrad,pawtab,dtset%usewvl,dtset%wvl_ngauss) 712 end if 713 #endif 714 end if 715 716 !end of WVL+PAW section 717 !---------------------------------------------------- 718 719 return 720 721 ! Handle IO error 722 10 continue 723 MSG_ERROR(errmsg) 724 725 end subroutine pspatm