TABLE OF CONTENTS
ABINIT/pspcor [ Functions ]
NAME
pspcor
FUNCTION
Compute ecore pseudoion-pseudoion correction energy from epsatm for different types of atoms in unit cell.
INPUTS
natom=number of atoms in cell ntypat=number of types of atoms typat(natom)=integer label of 'typat' for each atom in cell epsatm(ntypat)=pseudoatom energy for each type of atom zion(ntypat)=valence charge on each type of atom in cell
OUTPUT
ecore=resulting psion-psion energy in Hartrees
PARENTS
pspini
CHILDREN
SOURCE
629 subroutine pspcor(ecore,epsatm,natom,ntypat,typat,zion) 630 631 use defs_basis 632 use m_profiling_abi 633 634 !This section has been created automatically by the script Abilint (TD). 635 !Do not modify the following lines by hand. 636 #undef ABI_FUNC 637 #define ABI_FUNC 'pspcor' 638 !End of the abilint section 639 640 implicit none 641 642 !Arguments ------------------------------------ 643 !scalars 644 integer,intent(in) :: natom,ntypat 645 real(dp),intent(out) :: ecore 646 !arrays 647 integer,intent(in) :: typat(natom) 648 real(dp),intent(in) :: epsatm(ntypat),zion(ntypat) 649 650 !Local variables------------------------------- 651 !scalars 652 integer :: ia 653 real(dp) :: charge,esum 654 655 ! ************************************************************************* 656 657 charge = 0.d0 658 esum = 0.d0 659 do ia=1,natom 660 ! compute pseudocharge: 661 charge=charge+zion(typat(ia)) 662 ! add pseudocore energies together: 663 esum = esum + epsatm(typat(ia)) 664 end do 665 666 ecore=charge*esum 667 668 end subroutine pspcor
ABINIT/pspini [ Functions ]
NAME
pspini
FUNCTION
Looping over atom types 1 ... ntypat, read pseudopotential data filename, then call pspatm for each psp. Might combine the psps to generate pseudoatoms, thanks to alchemy. Also compute ecore=[Sum(i) zion(i)] * [Sum(i) epsatm(i)] by calling pspcor.
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 | iscf=parameter controlling scf or non-scf calculations | ixc=exchange-correlation choice as input to main routine | natom=number of atoms in unit cell | pawxcdev=choice of XC development in PAW formalism | prtvol= control output volume | typat(natom)=type (integer) for each atom | main routine, for each type of atom gsqcut=cutoff for G^2 based on ecut for basis sphere (bohr^-2) gsqcutdg=PAW only - cutoff for G^2 based on ecutdg (fine grid) for basis sphere (bohr^-2) rprimd(3,3)=dimensional primitive translations in real space (bohr) used to estimate real space mesh (if necessary)
OUTPUT
ecore=total psp core correction energy*ucvol (hartree*bohr^3) pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data gencond=general condition for new computation of pseudopotentials (if gencond=1, new psps have been re-computed)
SIDE EFFECTS
psps <type(pseudopotential_type)>=at output, psps is completely initialized At the input, it is already partially or completely initialized.
NOTES
The interplay with the multi-dataset mode is interesting : the pseudopotentials are independent of the dataset, but the largest q vector, the spin-orbit characteristics, the use of Ylm as well as ixc play a role in the set up of pseudopotentials (ixc plays a very minor role, however). So, the pseudopotential data ought not be recomputed when gsqcut, gsqcutdg, mqgrid_ff, mqgrid_vl, npspso, ixc, dimekb and useylm do not change. In many cases, this routine is also called just to write the psp line of the header, without reading again the psp. This psp line is constant throughout run.
PARENTS
bethe_salpeter,eph,gstate,nonlinear,respfn,screening,sigma,wfk_analyze
CHILDREN
nctab_free,nctab_init,nctab_mixalch,pawtab_set_flags,pspatm,pspcor psps_ncwrite,psps_print,timab,wrtout,xmpi_sum
SOURCE
65 #if defined HAVE_CONFIG_H 66 #include "config.h" 67 #endif 68 69 #include "abi_common.h" 70 71 subroutine pspini(dtset,dtfil,ecore,gencond,gsqcut,gsqcutdg,pawrad,pawtab,psps,rprimd,comm_mpi) 72 73 use defs_basis 74 use defs_datatypes 75 use defs_abitypes 76 use m_errors 77 use m_profiling_abi 78 use m_xmpi 79 80 use m_psps, only : psps_print, psps_ncwrite, nctab_init, nctab_free, nctab_mixalch 81 use m_pawrad, only : pawrad_type 82 use m_pawtab, only : pawtab_type, pawtab_set_flags 83 84 !This section has been created automatically by the script Abilint (TD). 85 !Do not modify the following lines by hand. 86 #undef ABI_FUNC 87 #define ABI_FUNC 'pspini' 88 use interfaces_14_hidewrite 89 use interfaces_18_timing 90 use interfaces_64_psp, except_this_one => pspini 91 !End of the abilint section 92 93 implicit none 94 95 !Arguments ------------------------------------ 96 !scalars 97 integer, optional,intent(in) :: comm_mpi 98 integer,intent(out) :: gencond 99 real(dp),intent(in) :: gsqcut,gsqcutdg 100 real(dp),intent(out) :: ecore 101 type(dataset_type),intent(in) :: dtset 102 type(datafiles_type),intent(in) :: dtfil 103 !arrays 104 real(dp),intent(in) :: rprimd(3,3) 105 !no_abirules 106 type(pseudopotential_type), target,intent(inout) :: psps 107 type(pawrad_type), intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 108 type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 109 110 !Local variables------------------------------- 111 !scalars 112 integer,parameter :: npspmax=50 113 integer,save :: dimekb_old=0,ifirst=1,ixc_old=-1,lmnmax_old=0,lnmax_old=0 114 integer,save :: mpssoang_old=0,mqgridff_old=0,mqgridvl_old=0,optnlxccc_old=-1 115 integer,save :: paw_size_old=-1,pawxcdev_old=-1,positron_old=-2,usepaw_old=-1 116 integer,save :: usexcnhat_old=-1,usewvl_old=-1,useylm_old=-1 117 integer :: comm_mpi_,ierr,ii,ilang,ilmn,ilmn0,iproj,ipsp,ipspalch 118 integer :: ispin,itypalch,itypat,mtypalch,npsp,npspalch,ntypalch 119 integer :: ntypat,ntyppure,paw_size 120 logical :: has_kij,has_tproj,has_tvale,has_nabla,has_shapefncg,has_vminushalf,has_wvl 121 real(dp),save :: ecore_old=zero,gsqcut_old=zero,gsqcutdg_old=zero 122 real(dp) :: dq,epsatm_psp,qmax,rmax,xcccrc 123 character(len=500) :: message 124 type(pawrad_type) :: pawrad_dum 125 type(pawtab_type) :: pawtab_dum 126 type(nctab_t) :: nctab_dum 127 type(nctab_t),pointer :: nctab_ptr 128 !arrays 129 integer :: paw_options(9) 130 integer,save :: paw_options_old(9)=(/-1,-1,-1,-1,-1,-1,-1,-1,-1/) 131 integer,save :: pspso_old(npspmax),pspso_zero(npspmax) 132 integer,allocatable :: indlmn_alch(:,:,:),new_pspso(:) 133 integer,pointer :: indlmn(:,:) 134 real(dp),save :: epsatm(npspmax) 135 real(dp) :: tsec(2) 136 real(dp),allocatable :: dvlspl(:,:),dvlspl_alch(:,:,:),ekb(:),ekb_alch(:,:) 137 real(dp),allocatable :: epsatm_alch(:),ffspl(:,:,:),ffspl_alch(:,:,:,:) 138 real(dp),allocatable :: vlspl(:,:),vlspl_alch(:,:,:),xccc1d(:,:) 139 real(dp),allocatable :: xccc1d_alch(:,:,:),xcccrc_alch(:) 140 type(nctab_t),target,allocatable :: nctab_alch(:) 141 142 ! ************************************************************************* 143 144 DBG_ENTER("COLL") 145 146 !Keep track of time spent in this subroutine 147 call timab(15,1,tsec) 148 149 !------------------------------------------------------------- 150 ! Some initializations 151 !------------------------------------------------------------- 152 153 !Useful sizes 154 ntypat=psps%ntypat 155 mtypalch=psps%mtypalch 156 npsp=psps%npsp 157 if (npsp>npspmax) then 158 MSG_BUG("npsp>npspmax in pspini !") 159 end if 160 161 !Size of grids for atomic data represented in reciprocal space 162 163 !Set up q grids, make qmax 20% larger than largest expected: 164 qmax=1.2d0 * sqrt(gsqcut) 165 !ffnl is always computed in reciprocal space 166 dq=qmax/(one*(psps%mqgrid_ff-1)) 167 do ii=1,psps%mqgrid_ff 168 psps%qgrid_ff(ii)=(ii-1)*dq 169 end do 170 if (psps%usepaw==1) qmax=1.2d0 * sqrt(gsqcutdg) 171 !If vlspl is computed in real space, qgrid contains a real space mesh 172 !the max is taken as the biggest distance in the box. 173 if (psps%vlspl_recipSpace) then 174 dq=qmax/(one*(psps%mqgrid_vl-1)) 175 else 176 rmax = (rprimd(1, 1) + rprimd(1, 2) + rprimd(1, 3)) ** 2 177 rmax = rmax + (rprimd(2, 1) + rprimd(2, 2) + rprimd(2, 3)) ** 2 178 rmax = rmax + (rprimd(3, 1) + rprimd(3, 2) + rprimd(3, 3)) ** 2 179 rmax = sqrt(rmax) 180 dq = rmax /(one*(psps%mqgrid_vl-1)) 181 end if 182 do ii=1,psps%mqgrid_vl 183 psps%qgrid_vl(ii)=(ii-1)*dq 184 end do 185 186 !Determine whether new optional data requests have changed 187 paw_options=0;paw_size=0 188 if (psps%usepaw==1) then 189 paw_size=size(pawtab) 190 has_kij=(dtset%positron/=0) 191 has_tvale=.true. ! Will be modified later (depending on PAW dataset format) 192 has_nabla=.false. 193 has_shapefncg=(dtset%optdriver==RUNL_GSTATE.and.((dtset%iprcel>=20.and.dtset%iprcel<70).or.dtset%iprcel>=80)) 194 has_wvl=(dtset%usewvl==1.or.dtset%icoulomb/=0) 195 has_tproj=(dtset%usewvl==1) ! projectors will be free at the end of the psp reading 196 has_vminushalf=(maxval(dtset%ldaminushalf)==1) 197 if (has_kij) paw_options(1)=1 198 if (has_tvale) paw_options(2)=1 199 if (has_nabla) paw_options(5)=1 200 if (has_shapefncg) paw_options(6)=1 201 if (has_wvl) paw_options(7)=1 202 if (has_tproj) paw_options(8)=1 203 if (has_vminushalf)paw_options(9)=1 204 !if (dtset%prtvclmb /= 0) then 205 paw_options(3) = 1 206 paw_options(4) = 1 207 !end if 208 end if 209 210 !Determine whether the spin-orbit characteristic has changed 211 !Do not forget that the SO is not consistent with alchemy presently 212 ABI_ALLOCATE(new_pspso,(npsp)) 213 if (ifirst==1) pspso_old(:)=-1 214 if (ifirst==1) pspso_zero(:)=-1 215 do ipsp=1,npsp 216 new_pspso(ipsp)=1 217 ! No new characteristics if it is equal to the old one, 218 ! or, if it is one, the old one is equal to the intrinsic characteristic one. 219 if (psps%pspso(ipsp)==pspso_old(ipsp).or. & 220 & (psps%pspso(ipsp)==1.and.pspso_old(ipsp)==pspso_zero(ipsp))) then 221 new_pspso(ipsp)=0 222 end if 223 ! Prepare the saving of the intrinsic pseudopotential characteristics 224 if(psps%pspso(ipsp)==1) pspso_zero(ipsp)=0 225 end do 226 227 !Compute the general condition for new computation of pseudopotentials 228 gencond=0 229 if( ixc_old /= dtset%ixc & 230 & .or. mqgridff_old /= psps%mqgrid_ff & 231 & .or. mqgridvl_old /= psps%mqgrid_vl & 232 & .or. mpssoang_old /= psps%mpssoang & 233 & .or. abs(gsqcut_old-gsqcut)>1.0d-10 & 234 & .or. (psps%usepaw==1.and.abs(gsqcutdg_old-gsqcutdg)>1.0d-10) & 235 & .or. dimekb_old /= psps%dimekb & 236 & .or. lmnmax_old /= psps%lmnmax & 237 & .or. lnmax_old /= psps%lnmax & 238 & .or. optnlxccc_old /= psps%optnlxccc & 239 & .or. usepaw_old /= psps%usepaw & 240 & .or. useylm_old /= psps%useylm & 241 & .or. pawxcdev_old /= dtset%pawxcdev & 242 & .or. positron_old /= dtset%positron & 243 & .or. usewvl_old /= dtset%usewvl & 244 & .or. paw_size_old /= paw_size & 245 & .or. usexcnhat_old/=dtset%usexcnhat_orig & 246 & .or. any(paw_options_old(:)/=paw_options(:)) & 247 & .or. sum(new_pspso(:))/=0 & 248 & .or. mtypalch>0 & 249 & .or. (dtset%usewvl==1.and.psps%usepaw==1)& 250 & ) gencond=1 251 252 if (present(comm_mpi).and.psps%usepaw==1) then 253 if(xmpi_comm_size(comm_mpi)>1)then 254 call xmpi_sum(gencond,comm_mpi,ierr) 255 end if 256 if (gencond/=0) gencond=1 257 end if 258 ABI_DEALLOCATE(new_pspso) 259 260 !------------------------------------------------------------- 261 ! Following section is only reached when new computation 262 ! of pseudopotentials is needed 263 !------------------------------------------------------------- 264 265 if (gencond==1) then 266 267 write(message, '(a,a)' ) ch10,& 268 & '--- Pseudopotential description ------------------------------------------------' 269 call wrtout(ab_out,message,'COLL') 270 271 ABI_ALLOCATE(ekb,(psps%dimekb*(1-psps%usepaw))) 272 ABI_ALLOCATE(xccc1d,(psps%n1xccc*(1-psps%usepaw),6)) 273 ABI_ALLOCATE(ffspl,(psps%mqgrid_ff,2,psps%lnmax)) 274 ABI_ALLOCATE(vlspl,(psps%mqgrid_vl,2)) 275 if (.not.psps%vlspl_recipSpace) then 276 ABI_ALLOCATE(dvlspl,(psps%mqgrid_vl,2)) 277 else 278 ABI_ALLOCATE(dvlspl,(0,0)) 279 end if 280 281 ! PAW: reset flags for optional data 282 if (psps%usepaw==1) then 283 call pawtab_set_flags(pawtab,has_kij=paw_options(1),has_tvale=paw_options(2),& 284 & has_vhnzc=paw_options(3),has_vhtnzc=paw_options(4),& 285 & has_nabla=paw_options(5),has_shapefncg=paw_options(6),& 286 & has_wvl=paw_options(7),has_tproj=paw_options(8)) 287 ! the following have to be included in pawtab_set_flags 288 do ipsp=1,psps%ntypat 289 pawtab(ipsp)%has_vminushalf=dtset%ldaminushalf(ipsp) 290 end do 291 end if 292 293 ! Read atomic pseudopotential data and get transforms 294 ! for each atom type : two cases, alchemy or not. 295 296 ! No alchemical pseudoatom, in all datasets, npsp=ntypat 297 if(mtypalch==0)then 298 299 do ipsp=1,npsp 300 301 xcccrc=zero 302 ekb(:)=zero;ffspl(:,:,:)=zero;vlspl(:,:)=zero 303 if (.not.psps%vlspl_recipSpace) dvlspl(:, :)=zero 304 if (psps%usepaw==0) xccc1d(:,:)=zero 305 indlmn=>psps%indlmn(:,:,ipsp) 306 indlmn(:,:)=0 307 308 write(message, '(a,i4,a,t38,a)' ) & 309 & '- pspini: atom type',ipsp,' psp file is',trim(psps%filpsp(ipsp)) 310 call wrtout(ab_out,message,'COLL') 311 call wrtout(std_out,message,'COLL') 312 313 ! Read atomic psp V(r) and wf(r) to get local and nonlocal psp: 314 ! Cannot use the same call in case of bound checking, because of pawrad/pawtab 315 if(psps%usepaw==0)then 316 call pspatm(dq,dtset,dtfil,ekb,epsatm(ipsp),ffspl,indlmn,ipsp,& 317 & pawrad_dum,pawtab_dum,psps,vlspl,dvlspl,xcccrc,xccc1d,psps%nctab(ipsp)) 318 psps%ekb(:,ipsp)=ekb(:) 319 psps%xccc1d(:,:,ipsp)=xccc1d(:,:) 320 else 321 comm_mpi_=xmpi_comm_self;if (present(comm_mpi)) comm_mpi_=comm_mpi 322 call pspatm(dq,dtset,dtfil,ekb,epsatm(ipsp),ffspl,indlmn,ipsp,& 323 & pawrad(ipsp),pawtab(ipsp),psps,vlspl,dvlspl,xcccrc,xccc1d,nctab_dum,& 324 & comm_mpi=comm_mpi_) 325 end if 326 327 ! Copy data to psps datastructure. 328 psps%xcccrc(ipsp)=xcccrc 329 psps%znucltypat(ipsp)=psps%znuclpsp(ipsp) 330 psps%ffspl(:,:,:,ipsp)=ffspl(:,:,:) 331 psps%vlspl(:,:,ipsp)=vlspl(:,:) 332 if (.not.psps%vlspl_recipSpace) psps%dvlspl(:, :, ipsp) = dvlspl(:, :) 333 end do ! ipsp 334 335 else ! if mtypalch/=0 336 337 npspalch=psps%npspalch 338 ntyppure=npsp-npspalch 339 ntypalch=psps%ntypalch 340 ABI_ALLOCATE(epsatm_alch,(npspalch)) 341 ABI_ALLOCATE(ekb_alch,(psps%dimekb,npspalch*(1-psps%usepaw))) 342 ABI_ALLOCATE(ffspl_alch,(psps%mqgrid_ff,2,psps%lnmax,npspalch)) 343 ABI_ALLOCATE(xccc1d_alch,(psps%n1xccc*(1-psps%usepaw),6,npspalch)) 344 ABI_ALLOCATE(xcccrc_alch,(npspalch)) 345 ABI_ALLOCATE(vlspl_alch,(psps%mqgrid_vl,2,npspalch)) 346 if (.not.psps%vlspl_recipSpace) then 347 ABI_ALLOCATE(dvlspl_alch,(psps%mqgrid_vl,2,npspalch)) 348 end if 349 ABI_ALLOCATE(indlmn,(6,psps%lmnmax)) 350 ABI_ALLOCATE(indlmn_alch,(6,psps%lmnmax,npspalch)) 351 352 ! Allocate NC tables used for mixing. 353 if (psps%usepaw == 0) then 354 ABI_DT_MALLOC(nctab_alch, (npspalch)) 355 do ipspalch=1,npspalch 356 call nctab_init(nctab_alch(ipspalch), psps%mqgrid_vl, .False., .False.) 357 end do 358 end if 359 360 do ipsp=1,npsp 361 write(message, '(a,i4,a,t38,a)' ) & 362 & '- pspini: atom type',ipsp,' psp file is',trim(psps%filpsp(ipsp)) 363 call wrtout(ab_out,message,'COLL') 364 365 xcccrc=zero 366 ekb(:)=zero;ffspl(:,:,:)=zero;vlspl(:,:)=zero 367 if (.not.psps%vlspl_recipSpace) dvlspl(:, :)=zero 368 if (psps%usepaw==0) xccc1d(:,:)=zero 369 indlmn(:,:)=0 370 371 ! Read atomic psp V(r) and wf(r) to get local and nonlocal psp: 372 if (psps%usepaw==0) then 373 if (ipsp <= ntyppure) then 374 ! Store data in nctab if pure atom. 375 nctab_ptr => psps%nctab(ipsp) 376 else 377 ! Store data in nctab_alch (to be mixed afterwards). 378 nctab_ptr => nctab_alch(ipsp-ntyppure) 379 end if 380 381 call pspatm(dq,dtset,dtfil,ekb,epsatm_psp,ffspl,indlmn,ipsp,& 382 & pawrad_dum,pawtab_dum,psps,vlspl,dvlspl,xcccrc,xccc1d,nctab_ptr) 383 384 else if (psps%usepaw==1) then 385 comm_mpi_=xmpi_comm_self;if (present(comm_mpi)) comm_mpi_=comm_mpi 386 call pspatm(dq,dtset,dtfil,ekb,epsatm_psp,ffspl,indlmn,ipsp,& 387 & pawrad(ipsp),pawtab(ipsp),psps,vlspl,dvlspl,xcccrc,xccc1d,nctab_dum,& 388 & comm_mpi=comm_mpi_) 389 end if 390 391 if (ipsp<=ntyppure) then 392 ! Pure pseudopotentials, leading to pure pseudoatoms 393 epsatm(ipsp)=epsatm_psp 394 psps%znucltypat(ipsp)=psps%znuclpsp(ipsp) 395 if (psps%usepaw==0) psps%ekb(:,ipsp)=ekb(:) 396 psps%ffspl(:,:,:,ipsp)=ffspl(:,:,:) 397 psps%vlspl(:,:,ipsp)=vlspl(:,:) 398 if (.not.psps%vlspl_recipSpace) psps%dvlspl(:, :, ipsp)=dvlspl(:, :) 399 if (psps%usepaw==0) psps%xccc1d(:,:,ipsp)=xccc1d(:,:) 400 psps%xcccrc(ipsp)=xcccrc 401 psps%indlmn(:,:,ipsp)=indlmn(:,:) 402 403 else 404 ! Pseudopotentials for alchemical generation 405 ipspalch=ipsp-ntyppure 406 epsatm_alch(ipspalch)=epsatm_psp 407 ffspl_alch(:,:,:,ipspalch)=ffspl(:,:,:) 408 vlspl_alch(:,:,ipspalch)=vlspl(:,:) 409 if (.not.psps%vlspl_recipSpace) dvlspl_alch(:,:,ipspalch)=dvlspl(:,:) 410 if (psps%usepaw==0) then 411 ekb_alch(:,ipspalch)=ekb(:) 412 xccc1d_alch(:,:,ipspalch)=xccc1d(:,:) 413 end if 414 xcccrc_alch(ipspalch)=xcccrc 415 indlmn_alch(:,:,ipspalch)=indlmn(:,:) 416 ! write(std_out,'(a,6i4)' )' pspini : indlmn_alch(:,1,ipspalch)=',indlmn_alch(:,1,ipspalch) 417 ! write(std_out,'(a,6i4)' )' pspini : indlmn_alch(:,2,ipspalch)=',indlmn_alch(:,2,ipspalch) 418 end if 419 420 end do ! ipsp 421 422 ! Generate data for alchemical pseudos. 423 do itypalch=1,ntypalch 424 itypat=itypalch+ntyppure 425 psps%znucltypat(itypat)=200.0+itypalch ! Convention for alchemical pseudoatoms 426 vlspl(:,:)=zero 427 if (.not.psps%vlspl_recipSpace) dvlspl(:, :) = zero 428 epsatm(itypat)=zero 429 xcccrc=zero 430 if (psps%usepaw==0) xccc1d(:,:)=zero 431 432 ! Here, linear combination of the quantities 433 ! MG: FIXME I think that the mixing of xcccrc is wrong when the xxccrc are different! 434 ! but this is minor bug since alchemical pseudos should not have XCCC (?) 435 do ipspalch=1,npspalch 436 epsatm(itypat) = epsatm(itypat) + epsatm_alch(ipspalch) * psps%mixalch(ipspalch,itypalch) 437 vlspl(:,:) = vlspl(:,:) + vlspl_alch(:,:,ipspalch) * psps%mixalch(ipspalch,itypalch) 438 if (.not.psps%vlspl_recipSpace) then 439 dvlspl(:,:) = dvlspl(:,:) + dvlspl_alch(:,:,ipspalch) * psps%mixalch(ipspalch,itypalch) 440 end if 441 xcccrc = xcccrc + xcccrc_alch(ipspalch) * psps%mixalch(ipspalch,itypalch) 442 if (psps%usepaw==0) then 443 xccc1d(:,:) = xccc1d(:,:) + xccc1d_alch(:,:,ipspalch) * psps%mixalch(ipspalch,itypalch) 444 end if 445 end do ! ipspalch 446 447 psps%vlspl(:,:,itypat)=vlspl(:,:) 448 if (.not.psps%vlspl_recipSpace) psps%dvlspl(:, :, itypat) = dvlspl(:, :) 449 if (psps%usepaw==0) then 450 psps%xccc1d(:,:,itypat)=xccc1d(:,:) 451 end if 452 psps%xcccrc(itypat)=xcccrc 453 454 if (abs(xcccrc) > tol6) then 455 write(std_out, *)"xcccrc", xcccrc 456 MSG_WARNING("Alchemical pseudopotential with nlcc!") 457 end if 458 459 ! Combine the different non-local projectors : for the scalar part then 460 ! the spin-orbit part, treat the different angular momenta 461 ! WARNING : this coding does not work for PAW 462 ilmn=0 463 psps%indlmn(:,:,itypat)=0 464 do ispin=1,2 465 do ilang=0,3 466 if(ispin==2 .and. ilang==0)cycle 467 iproj=0 468 do ipspalch=1,npspalch 469 if(abs(psps%mixalch(ipspalch,itypalch))>tol10)then 470 do ilmn0=1,psps%lmnmax 471 if(indlmn_alch(5,ilmn0,ipspalch)/=0)then 472 if(indlmn_alch(6,ilmn0,ipspalch)==ispin)then 473 if(indlmn_alch(1,ilmn0,ipspalch)==ilang)then 474 ilmn=ilmn+1 ! increment the counter 475 iproj=iproj+1 ! increment the counter, this does not work for PAW 476 if(ilmn>psps%lmnmax)then 477 MSG_BUG('Problem with the alchemical pseudopotentials : ilmn>lmnmax.') 478 end if 479 psps%indlmn(1,ilmn,itypat)=ilang 480 psps%indlmn(2,ilmn,itypat)=indlmn_alch(2,ilmn0,ipspalch) 481 psps%indlmn(3,ilmn,itypat)=iproj ! This does not work for PAW 482 psps%indlmn(4,ilmn,itypat)=ilmn ! This does not work for PAW 483 psps%indlmn(5,ilmn,itypat)=ilmn 484 psps%indlmn(6,ilmn,itypat)=ispin 485 ! The two lines below do not work for PAW 486 if (psps%usepaw==0) then 487 psps%ekb(ilmn,itypat)=psps%mixalch(ipspalch,itypalch) *ekb_alch(ilmn0,ipspalch) 488 end if 489 psps%ffspl(:,:,ilmn,itypat)=ffspl_alch(:,:,ilmn0,ipspalch) 490 end if ! ilang is OK 491 end if ! ispin is OK 492 end if ! ilmn0 exist 493 end do ! ilmn0 494 end if ! mixalch>tol10 495 end do ! ipspalch 496 end do ! ilang 497 end do ! ispin 498 499 end do ! itypalch 500 501 ABI_DEALLOCATE(epsatm_alch) 502 ABI_DEALLOCATE(ekb_alch) 503 ABI_DEALLOCATE(ffspl_alch) 504 ABI_DEALLOCATE(xccc1d_alch) 505 ABI_DEALLOCATE(xcccrc_alch) 506 ABI_DEALLOCATE(vlspl_alch) 507 if (.not.psps%vlspl_recipSpace) then 508 ABI_DEALLOCATE(dvlspl_alch) 509 end if 510 ABI_DEALLOCATE(indlmn_alch) 511 ABI_DEALLOCATE(indlmn) 512 513 ! Mix NC tables. 514 if (psps%usepaw == 0) then 515 call nctab_mixalch(nctab_alch, npspalch, ntypalch, psps%algalch, psps%mixalch, psps%nctab(ntyppure+1:)) 516 do ipspalch=1,npspalch 517 call nctab_free(nctab_alch(ipspalch)) 518 end do 519 ABI_DT_FREE(nctab_alch) 520 end if 521 end if ! mtypalch 522 523 ABI_DEALLOCATE(ekb) 524 ABI_DEALLOCATE(ffspl) 525 ABI_DEALLOCATE(vlspl) 526 ABI_DEALLOCATE(xccc1d) 527 528 if (.not.psps%vlspl_recipSpace) then 529 ABI_DEALLOCATE(dvlspl) 530 end if 531 532 end if ! End condition of new computation needed 533 534 !------------------------------------------------------------- 535 ! Following section is always executed 536 !------------------------------------------------------------- 537 !One should move this section of code outside of pspini, 538 !but epsatm is needed, so should be in the psp datastructure. 539 !Compute pseudo correction energy. Will differ from an already 540 !computed one if the number of atom differ ... 541 call pspcor(ecore,epsatm,dtset%natom,ntypat,dtset%typat,psps%ziontypat) 542 if(abs(ecore_old-ecore)>tol8*abs(ecore_old+ecore))then 543 write(message, '(2x,es15.8,t50,a)' ) ecore,'ecore*ucvol(ha*bohr**3)' 544 ! ecore is useless if iscf<=0, but at least it has been initialized 545 if(dtset%iscf>=0)then 546 call wrtout(ab_out,message,'COLL') 547 end if 548 call wrtout(std_out,message,'COLL') 549 end if 550 551 !End of pseudopotential output section 552 write(message, '(2a)' )& 553 & '--------------------------------------------------------------------------------',ch10 554 call wrtout(ab_out,message,'COLL') 555 556 !------------------------------------------------------------- 557 ! Keep track of this call to the routine 558 !------------------------------------------------------------- 559 560 if (ifirst==1) ifirst=0 561 562 mqgridff_old=psps%mqgrid_ff 563 mqgridvl_old=psps%mqgrid_vl 564 mpssoang_old=psps%mpssoang 565 ixc_old=dtset%ixc 566 gsqcut_old=gsqcut;if (psps%usepaw==1) gsqcutdg_old=gsqcutdg 567 lmnmax_old=psps%lmnmax 568 lnmax_old=psps%lnmax 569 optnlxccc_old=psps%optnlxccc 570 usepaw_old=psps%usepaw 571 dimekb_old=psps%dimekb 572 useylm_old=psps%useylm 573 pawxcdev_old=dtset%pawxcdev 574 positron_old=dtset%positron 575 usewvl_old = dtset%usewvl 576 usexcnhat_old=dtset%usexcnhat_orig 577 paw_size_old=paw_size 578 ecore_old=ecore 579 paw_options_old(:)=paw_options(:) 580 581 do ipsp=1,npsp 582 pspso_old(ipsp)=psps%pspso(ipsp) 583 if(pspso_zero(ipsp)==0)pspso_zero(ipsp)=psps%pspso(ipsp) 584 end do 585 psps%mproj = maxval(psps%indlmn(3,:,:)) 586 587 if (gencond == 1) call psps_print(psps,std_out,dtset%prtvol) 588 589 ! Write the PSPS.nc file and exit here if requested by the user. 590 if (abs(dtset%prtpsps) == 1) then 591 if (xmpi_comm_rank(xmpi_world) == 0) call psps_ncwrite(psps, trim(dtfil%filnam_ds(4))//"_PSPS.nc") 592 if (dtset%prtpsps == -1) then 593 MSG_ERROR_NODUMP("prtpsps == -1 ==> aborting now") 594 end if 595 end if 596 597 call timab(15,2,tsec) 598 599 DBG_EXIT("COLL") 600 601 end subroutine pspini