TABLE OF CONTENTS
ABINIT/setup_positron [ Functions ]
NAME
setup_positron
FUNCTION
Do various initializations for the positron lifetime calculation
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (GJ,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
atindx(natom)=index table for atoms (see gstate.f) atindx1(natom)=index table for atoms, inverse of atindx dtefield <type(efield_type)> = variables related to Berry phase dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables for this dataset ecore=core psp energy (part of total energy) (hartree) etotal=current value of total energy fock <type(fock_type)>= quantities to calculate Fock exact exchange forces_needed=if >0 forces are needed fred(3,natom)=forces in reduced coordinates gprimd(3,3)=dimensional primitive translations for reciprocal space gmet(3,3)=reciprocal space metric grchempottn(3,natom)=d(E_chemical_potential)/d(xred) (hartree) grewtn(3,natom)=d(Ewald)/d(xred) (hartree) grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D dispersion (hartree) gsqcut=cutoff value on G**2 for sphere inside fft box hdr <type(hdr_type)>=the header of wf, den and pot files ifirst_gs= 0 if we are in a single ground-state calculation or in the first ground-state calculation of a structural minimization/dynamics indsym(4,nsym,natom)=index showing transformation of atom labels under symmetry operations (computed in symatm) istep=index of the number of steps in the routine scfcv istep_mix=index of the number of steps for the SCF mixing (can be <istep) kg(3,mpw*mkmem)=reduced (integer) coordinates of G vecs in basis sphere kxc(nfft,nkxc)=exchange-correlation kernel, needed only if nkxc>0 maxfor=maximum absolute value of fcart (forces) mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mgfft=maximum size of 1D FFTs mpi_enreg=informations about MPI parallelization my_natom=number of atoms treated by current processor n3xccc=dimension of the xccc3d array (0 or nfftf). nattyp(ntypat)= # atoms of each type. nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc nhat(nfftf,nspden*usepaw)= -PAW only- compensation density nkxc=second dimension of the array kxc, see rhotoxc.f for a description npwarr(nkpt)=number of planewaves in basis and on boundary for each k nvresid(nfftf,nspden)=array for the residual of the density/potential optres=0 if the potential residual has to be used for forces corrections =1 if the density residual has to be used for forces corrections paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawang <type(pawang_type)>=paw angular mesh and related data pawfgr(my_natom*usepaw) <type(pawfgr_type)>=fine grid parameters and related data pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfftf+1)*natom)=one-dimensional structure factor information (fine FFT grid) ph1dc(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases (coarse FFT grid) psps <type(pseudopotential_type)>=variables related to pseudopotentials rprimd(3,3)=dimensional primitive translations in real space (bohr) stress_needed=if >0 stresses are needed strsxc(6)=xc correction to stress symrec(3,3,nsym)=symmetries in reciprocal space, reduced coordinates ucvol=unit cell volume in bohr**3. usecprj= 1 if cprj array is stored in memory vhartr(nfftf)=array for holding Hartree potential vpsp(nfftf)=array for holding local psp vxc(nfftf,nspden)=exchange-correlation potential (hartree) in real space xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3 xred(3,natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
SIDE EFFECTS
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation energies <type(energies_type)>=all part of total energy. cg(2,mcg)=wavefunctions cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors: cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector. eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) occ(mband*nkpt*nsppol)=occupation number for each band at each k point pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data rhog(2,nfft)=Fourier transform of total electron/positron density rhor(nfft,nspden)=total electron/positron density (el/bohr**3)
PARENTS
scfcv
CHILDREN
energies_copy,energies_init,forstr,fourdp,getcut,hartre,initrhoij initro,pawcprj_alloc,pawcprj_copy,pawcprj_free,pawmknhat,pawrhoij_alloc pawrhoij_copy,pawrhoij_free,read_rhor,wrtout
SOURCE
104 #if defined HAVE_CONFIG_H 105 #include "config.h" 106 #endif 107 108 #include "abi_common.h" 109 110 subroutine setup_positron(atindx,atindx1,cg,cprj,dtefield,dtfil,dtset,ecore,eigen,etotal,electronpositron,& 111 & energies,fock,forces_needed,fred,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcut,hdr,ifirst_gs,indsym,istep,istep_mix,kg,& 112 & kxc,maxfor,mcg,mcprj,mgfft,mpi_enreg,my_natom,n3xccc,nattyp,nfft,ngfft,ngrvdw,nhat,nkxc,npwarr,nvresid,occ,optres,& 113 & paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,ph1d,ph1dc,psps,rhog,rhor,& 114 & rprimd,stress_needed,strsxc,symrec,ucvol,usecprj,vhartr,vpsp,vxc,& 115 & xccc3d,xred,ylm,ylmgr) 116 117 use defs_basis 118 use defs_datatypes 119 use defs_abitypes 120 use m_efield 121 use m_errors 122 use m_profiling_abi 123 use m_energies 124 use m_wffile 125 use m_electronpositron 126 use m_hdr 127 128 use m_ioarr, only : ioarr, read_rhor 129 use m_pawang, only : pawang_type 130 use m_pawrad, only : pawrad_type 131 use m_pawtab, only : pawtab_type 132 use m_paw_ij, only : paw_ij_type 133 use m_pawfgrtab,only : pawfgrtab_type 134 use m_pawrhoij, only : pawrhoij_type, pawrhoij_copy, pawrhoij_free, pawrhoij_alloc 135 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy 136 use m_pawfgr, only : pawfgr_type 137 use m_fock, only : fock_type 138 use m_kg, only : getcut 139 use defs_wvltypes, only : wvl_data 140 141 !This section has been created automatically by the script Abilint (TD). 142 !Do not modify the following lines by hand. 143 #undef ABI_FUNC 144 #define ABI_FUNC 'setup_positron' 145 use interfaces_14_hidewrite 146 use interfaces_53_ffts 147 use interfaces_56_xc 148 use interfaces_65_paw 149 use interfaces_67_common, except_this_one => setup_positron 150 !End of the abilint section 151 152 implicit none 153 154 !Arguments ------------------------------------ 155 !scalars 156 integer,intent(in) :: forces_needed,ifirst_gs,istep,mcg,mcprj,mgfft,my_natom,n3xccc,nfft 157 integer,intent(in) :: ngrvdw,nkxc,optres,stress_needed,usecprj 158 integer,intent(inout) :: istep_mix 159 real(dp),intent(in) :: ecore,etotal,gsqcut,maxfor,ucvol 160 type(efield_type),intent(in) :: dtefield 161 type(datafiles_type),intent(in) :: dtfil 162 type(dataset_type),intent(in) :: dtset 163 type(electronpositron_type),pointer :: electronpositron 164 type(energies_type),intent(inout) :: energies 165 type(hdr_type),intent(inout) :: hdr 166 type(MPI_type),intent(inout) :: mpi_enreg 167 type(pawang_type),intent(in) :: pawang 168 type(pawfgr_type),intent(in) :: pawfgr 169 type(pseudopotential_type), intent(in) :: psps 170 type(fock_type),pointer, intent(inout) :: fock 171 !arrays 172 integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),indsym(4,dtset%nsym,dtset%natom) 173 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),nattyp(dtset%natom),ngfft(18) 174 integer,intent(in) :: npwarr(dtset%nkpt),symrec(3,3,dtset%nsym) 175 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),grchempottn(3,dtset%natom) 176 real(dp),intent(in) :: grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(nfft,nkxc) 177 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),ph1dc(2,(3*(2*dtset%mgfft+1)*dtset%natom)*dtset%usepaw) 178 real(dp),intent(in) :: rprimd(3,3),strsxc(6),vhartr(nfft),vpsp(nfft),vxc(nfft,dtset%nspden) 179 real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 180 real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm) 181 real(dp),intent(inout) :: cg(2,mcg) 182 real(dp),intent(inout) :: nhat(nfft,dtset%nspden*dtset%usepaw) 183 real(dp),intent(inout) :: nvresid(nfft,dtset%nspden) 184 real(dp),intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol),fred(3,dtset%natom) 185 real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 186 real(dp),intent(inout) :: rhog(2,nfft),rhor(nfft,dtset%nspden) 187 real(dp),intent(inout) :: xccc3d(n3xccc),xred(3,dtset%natom) 188 type(pawcprj_type),intent(inout) :: cprj(dtset%natom,mcprj*usecprj) 189 type(paw_ij_type),intent(in) :: paw_ij(my_natom*dtset%usepaw) 190 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*dtset%usepaw) 191 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw) 192 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw) 193 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*dtset%usepaw) 194 195 !Local variables------------------------------- 196 !scalars 197 integer,parameter :: cplex1=1 198 integer :: history_level,iatom,iband,icalctype,icalctype0,icg,ifft,ikpt 199 integer :: iocc,ireadwf,ispden,isppol,n3xccc0,nocc,optfor,optstr,rdwrpaw,comm_cell 200 logical,parameter :: always_restart=.false. ! Set to true to restart by a pure electronic step at each new atomic structure 201 logical :: need_scocc,new_calctype 202 real(dp) :: boxcut_dum,diffor_dum,ecut_eff,eigtmp,etotal_read,gsqcut_eff,maxfor_dum 203 real(dp) :: maxocc,nelect,occlast,occtmp,rhotmp 204 character(len=69) :: TypeCalcStrg 205 character(len=500) :: message 206 character(len=fnlen) :: fname 207 type(energies_type) :: energies_tmp 208 type(wvl_data) :: wvl 209 type(hdr_type) :: hdr_den 210 !arrays 211 integer,allocatable :: nlmn(:) 212 real(dp) :: cgtmp(2) 213 real(dp),parameter :: qphon(3)=(/zero,zero,zero/) 214 real(dp),allocatable :: favg_dum(:),fcart_dum(:,:),forold_dum(:,:),fred_tmp(:,:) 215 real(dp),allocatable :: gresid_dum(:,:),grhf_dum(:,:),grxc_dum(:,:) 216 real(dp),allocatable :: rhog_ep(:,:),scocc(:),str_tmp(:),synlgr_dum(:,:) 217 real(dp) :: nhatgr(0,0,0) 218 type(pawcprj_type),allocatable :: cprj_tmp(:,:) 219 type(pawrhoij_type),allocatable :: pawrhoij_tmp(:) 220 221 ! ************************************************************************* 222 223 !Compatibility tests 224 if (dtset%positron==0) then 225 MSG_BUG('Not valid for dtset%positron=0!') 226 end if 227 228 if (istep>1.and.nfft/=electronpositron%nfft) then 229 MSG_BUG('Invalid value for nfft!') 230 end if 231 232 if (dtset%usewvl==1) then 233 MSG_BUG('Not valid for wavelets!') 234 end if 235 236 if (dtset%positron==1) then 237 do isppol=1,dtset%nsppol 238 do ikpt=1,dtset%nkpt 239 if (dtset%nband(ikpt+dtset%nkpt*(isppol-1))/=dtset%nband(1)) then 240 message = "dtset%positron needs nband to be the same at each k-point !" 241 MSG_ERROR(message) 242 end if 243 end do 244 end do 245 end if 246 247 comm_cell = mpi_enreg%comm_cell 248 249 !----------------------------------------------------------------------- 250 !Compute new value for calctype (type of electron-positron calculation) 251 !----------------------------------------------------------------------- 252 icalctype0=electronpositron%calctype 253 254 new_calctype=.false. 255 if (dtset%positron==1.or.dtset%positron==2) then 256 if (istep==1) new_calctype=.true. 257 else if (dtset%positron<0) then 258 if (ifirst_gs/=0.and.istep==1.and.(.not.always_restart)) new_calctype=.true. 259 if (electronpositron%scf_converged) new_calctype=.true. 260 end if 261 262 !Comment: 263 !history_level=-1: not used 264 !history_level= 0: rhor from scratch, rhor_ep from scratch or read 265 !history_level= 1: rhor in memory, rhor_ep from scratch or read 266 !history_level= 2: rhor_ep <-rhor, rhor from scratch 267 !history_level= 3: rhor_ep <-> rhor 268 !history_level= 4: rhor in memory, rhor_ep in memory 269 history_level=-1 270 if (dtset%positron==1.or.dtset%positron==2) then 271 if (ifirst_gs==0.and.istep==1) history_level=0 272 if (ifirst_gs/=0.and.istep==1) history_level=4 273 else if (dtset%positron<0) then 274 if (.not.electronpositron%scf_converged) then 275 if (ifirst_gs/=0.and.istep==1.and.(.not.always_restart)) history_level=4 276 else if (electronpositron%scf_converged) then 277 if (icalctype0==0) history_level=2 278 if (icalctype0> 0) history_level=3 279 end if 280 end if 281 282 electronpositron%calctype=icalctype0 283 if (dtset%positron==1.or.dtset%positron==2) then 284 electronpositron%calctype=dtset%positron 285 else if (dtset%positron<0) then 286 if (electronpositron%scf_converged) then 287 if (icalctype0==0) electronpositron%calctype=1 288 if (icalctype0>0 ) electronpositron%calctype=3-electronpositron%calctype 289 else if (ifirst_gs/=0.and.istep==1) then 290 if (always_restart) then 291 electronpositron%calctype=0 292 else 293 electronpositron%calctype=2 294 ! if (electronpositron%particle==EP_POSITRON) electronpositron%calctype=1 295 end if 296 end if 297 end if 298 299 electronpositron%scf_converged=.false. 300 if (new_calctype) electronpositron%istep=electronpositron%istep+1 301 if (istep==1) electronpositron%istep=1 302 ireadwf=dtfil%ireadwf;if (electronpositron%istep>1) ireadwf=0 303 304 !============================================ 305 !The following lines occur only when the type 306 !of electron-positron calculation changes 307 !============================================ 308 if (new_calctype) then 309 310 ! Reset some indexes 311 if (electronpositron%calctype==0) then 312 electronpositron%particle=EP_NOTHING 313 else if (electronpositron%calctype==1) then 314 electronpositron%particle=EP_ELECTRON 315 else if (electronpositron%calctype==2) then 316 electronpositron%particle=EP_POSITRON 317 end if 318 electronpositron%has_pos_ham=mod(electronpositron%calctype,2) 319 electronpositron%istep_scf=1 320 istep_mix=1 321 322 ! ----------------------------------------------------------------------------------------- 323 ! Update forces and stresses 324 ! If electronpositron%calctype==1: fred_ep/stress_ep are the electronic fred/stress 325 ! If electronpositron%calctype==2: fred_ep/stress_ep are the positronic fred/stress 326 ! ----------------------------------------------------------------------------------------- 327 if (history_level==2.or.history_level==3) then 328 optstr=0;optfor=0 329 if (allocated(electronpositron%stress_ep)) optstr=stress_needed 330 if (allocated(electronpositron%fred_ep).and.forces_needed==2) optfor=1 331 if (optfor>0.or.optstr>0) then 332 ABI_ALLOCATE(favg_dum,(3)) 333 ABI_ALLOCATE(fcart_dum,(3,dtset%natom)) 334 ABI_ALLOCATE(forold_dum,(3,dtset%natom)) 335 ABI_ALLOCATE(gresid_dum,(3,dtset%natom)) 336 ABI_ALLOCATE(grhf_dum,(3,dtset%natom)) 337 ABI_ALLOCATE(grxc_dum,(3,dtset%natom)) 338 ABI_ALLOCATE(synlgr_dum,(3,dtset%natom)) 339 ABI_ALLOCATE(fred_tmp,(3,dtset%natom)) 340 ABI_ALLOCATE(str_tmp,(6)) 341 forold_dum=zero;n3xccc0=n3xccc 342 icalctype=electronpositron%calctype;electronpositron%calctype=-icalctype0 343 if (electronpositron%calctype==0) electronpositron%calctype=-100 344 if (electronpositron%calctype==-1) n3xccc0=0 ! Note: if calctype=-1, previous calculation was positron 345 call forstr(atindx1,cg,cprj,diffor_dum,dtefield,dtset,eigen,electronpositron,energies,& 346 & favg_dum,fcart_dum,fock,forold_dum,fred_tmp,grchempottn,gresid_dum,grewtn,grhf_dum,grvdw,grxc_dum,gsqcut,& 347 & indsym,kg,kxc,maxfor_dum,mcg,mcprj,mgfft,mpi_enreg,my_natom,n3xccc0,nattyp,nfft,ngfft,& 348 & ngrvdw,nhat,nkxc,npwarr,dtset%ntypat,nvresid,occ,optfor,optres,paw_ij,pawang,pawfgr,& 349 & pawfgrtab,pawrad,pawrhoij,pawtab,ph1dc,ph1d,psps,rhog,rhor,rprimd,optstr,strsxc,str_tmp,symrec,& 350 & synlgr_dum,ucvol,usecprj,vhartr,vpsp,vxc,wvl,xccc3d,xred,ylm,ylmgr,0.0_dp) 351 electronpositron%calctype=icalctype 352 if (optfor>0) electronpositron%fred_ep(:,:)=fred_tmp(:,:) 353 if (optstr>0) electronpositron%stress_ep(:)=str_tmp(:) 354 ABI_DEALLOCATE(favg_dum) 355 ABI_DEALLOCATE(fcart_dum) 356 ABI_DEALLOCATE(forold_dum) 357 ABI_DEALLOCATE(gresid_dum) 358 ABI_DEALLOCATE(grhf_dum) 359 ABI_DEALLOCATE(grxc_dum) 360 ABI_DEALLOCATE(synlgr_dum) 361 ABI_DEALLOCATE(fred_tmp) 362 ABI_DEALLOCATE(str_tmp) 363 end if 364 if (optfor==0.and.forces_needed>0.and.allocated(electronpositron%fred_ep)) then 365 electronpositron%fred_ep(:,:)=fred(:,:)-electronpositron%fred_ep(:,:) 366 end if 367 end if 368 369 ! ---------------------------------------------------------------------------------------------------- 370 ! Initialize/Update densities 371 ! If electronpositron%calctype==1: rhor is the positronic density, rhor_ep is the electronic density 372 ! If electronpositron%calctype==2: rhor is the electronic density, rhor_ep is the positronic density 373 ! --------------------------------------------------------------------------------------------------- 374 ABI_ALLOCATE(rhog_ep,(2,nfft)) 375 376 ! ===== PREVIOUS DENSITY RHOR_EP: 377 if (history_level==0.or.history_level==1) then 378 ! ----- Read from disk 379 if (dtset%positron>0) then 380 rdwrpaw=dtset%usepaw 381 fname=trim(dtfil%fildensin);if (dtset%positron==2) fname=trim(dtfil%fildensin)//'_POSITRON' 382 call read_rhor(trim(fname), cplex1, dtset%nspden, nfft, ngfft, rdwrpaw, mpi_enreg, electronpositron%rhor_ep, & 383 hdr_den, electronpositron%pawrhoij_ep, comm_cell, check_hdr=hdr) 384 etotal_read = hdr_den%etot; call hdr_free(hdr_den) 385 call fourdp(1,rhog_ep,electronpositron%rhor_ep,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 386 if (dtset%usepaw==1.and.allocated(electronpositron%nhat_ep)) then 387 call pawmknhat(occtmp,1,0,0,0,0,gprimd,my_natom,dtset%natom,nfft,ngfft,0,& 388 & dtset%nspden,dtset%ntypat,pawang,pawfgrtab,nhatgr,electronpositron%nhat_ep,& 389 & electronpositron%pawrhoij_ep,electronpositron%pawrhoij_ep,pawtab,& 390 & qphon,rprimd,ucvol,dtset%usewvl,xred,distribfft=mpi_enreg%distribfft,& 391 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 392 & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0) 393 end if 394 end if 395 ! ----- Electronic from scratch 396 if (dtset%positron<0.and.electronpositron%calctype==1) then 397 ecut_eff=dtset%pawecutdg*(dtset%dilatmx)**2 398 call getcut(boxcut_dum,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,std_out,qphon,ngfft) 399 call initro(atindx,dtset%densty,gmet,gsqcut_eff,dtset%usepaw,mgfft,mpi_enreg,& 400 & psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,dtset%ntypat,& 401 & dtset%paral_kgb,psps,pawtab,ph1d,psps%qgrid_vl,rhog_ep,electronpositron%rhor_ep,& 402 & dtset%spinat,ucvol,dtset%usepaw,dtset%ziontypat,dtset%znucl) 403 if (dtset%usepaw==1) then 404 if (size(electronpositron%pawrhoij_ep)>0) then 405 ABI_DATATYPE_ALLOCATE(pawrhoij_tmp,(my_natom)) 406 call initrhoij(electronpositron%pawrhoij_ep(1)%cplex,dtset%lexexch,& 407 & dtset%lpawu,my_natom,dtset%natom,dtset%nspden,& 408 & electronpositron%pawrhoij_ep(1)%nspinor,dtset%nsppol,& 409 & dtset%ntypat,pawrhoij_tmp,dtset%pawspnorb,pawtab,dtset%spinat,dtset%typat,& 410 & ngrhoij=electronpositron%pawrhoij_ep(1)%ngrhoij,& 411 & nlmnmix=electronpositron%pawrhoij_ep(1)%lmnmix_sz,& 412 & use_rhoij_=electronpositron%pawrhoij_ep(1)%use_rhoij_,& 413 & use_rhoijres=electronpositron%pawrhoij_ep(1)%use_rhoijres,& 414 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 415 if (electronpositron%pawrhoij_ep(1)%lmnmix_sz>0) then 416 do iatom=1,my_natom 417 pawrhoij_tmp(iatom)%kpawmix(:)=electronpositron%pawrhoij_ep(iatom)%kpawmix(:) 418 end do 419 end if 420 call pawrhoij_copy(pawrhoij_tmp,electronpositron%pawrhoij_ep) 421 call pawrhoij_free(pawrhoij_tmp) 422 ABI_DATATYPE_DEALLOCATE(pawrhoij_tmp) 423 end if 424 if (allocated(electronpositron%nhat_ep)) then 425 call pawmknhat(occtmp,1,0,0,0,0,gprimd,my_natom,dtset%natom,nfft,ngfft,0,& 426 & dtset%nspden,dtset%ntypat,pawang,pawfgrtab,nhatgr,electronpositron%nhat_ep,& 427 & electronpositron%pawrhoij_ep,electronpositron%pawrhoij_ep,pawtab,qphon,rprimd,& 428 & ucvol,dtset%usewvl,xred,distribfft=mpi_enreg%distribfft,& 429 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 430 & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0) 431 end if 432 end if 433 end if 434 ! ----- Positronic from scratch 435 if (dtset%positron<0.and.electronpositron%calctype==2) then 436 electronpositron%rhor_ep(:,1)=one/ucvol 437 if (dtset%nspden>=2) electronpositron%rhor_ep(:,2)=half/ucvol 438 if (dtset%nspden==4) electronpositron%rhor_ep(:,3:4)=zero 439 rhog_ep=zero;rhog_ep(1,1)=one/ucvol 440 if (dtset%usepaw==1) then 441 do iatom=1,dtset%natom 442 electronpositron%pawrhoij_ep(iatom)%rhoijp(:,:)=zero 443 electronpositron%pawrhoij_ep(iatom)%nrhoijsel=0 444 end do 445 if (allocated(electronpositron%nhat_ep)) electronpositron%nhat_ep=zero 446 end if 447 end if 448 end if 449 ! ----- Deduced from rhor in memory 450 if (history_level==2) then 451 electronpositron%rhor_ep(:,:)=rhor(:,:) 452 rhog_ep(:,:)=rhog(:,:) 453 if (dtset%usepaw==1) then 454 call pawrhoij_copy(pawrhoij,electronpositron%pawrhoij_ep) 455 if (allocated(electronpositron%nhat_ep)) electronpositron%nhat_ep(:,:)=nhat(:,:) 456 end if 457 end if 458 459 ! ===== CURRENT DENSITY RHOR: 460 if (history_level==0.or.history_level==2) then 461 if (ireadwf==0) then 462 ! ----- Positronic from scratch 463 if (electronpositron%calctype==1) then 464 rhor(:,1)=one/ucvol 465 if (dtset%nspden>=2) rhor(:,2)=half/ucvol 466 if (dtset%nspden==4) rhor(:,3:4)=zero 467 rhog=zero;rhog(1,1)=one/ucvol 468 if (dtset%usepaw==1) then 469 do iatom=1,my_natom 470 pawrhoij(iatom)%rhoijp(:,:)=zero 471 pawrhoij(iatom)%nrhoijsel=0 472 end do 473 nhat(:,:)=zero 474 end if 475 end if 476 ! ----- Electronic from scratch 477 if (electronpositron%calctype==2) then 478 ecut_eff=dtset%pawecutdg*(dtset%dilatmx)**2 479 call getcut(boxcut_dum,ecut_eff,gmet,gsqcut_eff,dtset%iboxcut,std_out,qphon,ngfft) 480 call initro(atindx,dtset%densty,gmet,gsqcut_eff,dtset%usepaw,mgfft,mpi_enreg,& 481 & psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%nspden,dtset%ntypat,& 482 & dtset%paral_kgb,psps,pawtab,ph1d,psps%qgrid_vl,rhog,rhor,dtset%spinat,ucvol,& 483 & dtset%usepaw,dtset%ziontypat,dtset%znucl) 484 485 if (dtset%usepaw==1) then 486 if (size(pawrhoij)>0) then 487 ABI_DATATYPE_ALLOCATE(pawrhoij_tmp,(my_natom)) 488 call initrhoij(pawrhoij(1)%cplex,dtset%lexexch,dtset%lpawu,& 489 & my_natom,dtset%natom,dtset%nspden,pawrhoij(1)%nspinor,dtset%nsppol,& 490 & dtset%ntypat,pawrhoij_tmp,dtset%pawspnorb,pawtab,dtset%spinat,& 491 & dtset%typat,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,& 492 & use_rhoij_=pawrhoij(1)%use_rhoij_,use_rhoijres=pawrhoij(1)%use_rhoijres,& 493 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 494 do iatom=1,my_natom 495 pawrhoij_tmp(iatom)%kpawmix(:)=pawrhoij(iatom)%kpawmix(:) 496 end do 497 call pawrhoij_copy(pawrhoij_tmp,pawrhoij) 498 call pawrhoij_free(pawrhoij_tmp) 499 ABI_DATATYPE_DEALLOCATE(pawrhoij_tmp) 500 end if 501 call pawmknhat(occtmp,1,0,0,0,0,gprimd,my_natom,dtset%natom,nfft,ngfft,0,& 502 & dtset%nspden,dtset%ntypat,pawang,pawfgrtab,nhatgr,nhat,& 503 & pawrhoij,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred, & 504 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,& 505 & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,& 506 & me_g0=mpi_enreg%me_g0,distribfft=mpi_enreg%distribfft) 507 end if 508 509 end if 510 end if 511 end if 512 513 ! ===== EXCHANGE POSITRONIC AND ELECTRONIC DENSITY (CURRENT AND PREVIOUS) 514 if (history_level==3) then 515 do ispden=1,dtset%nspden 516 do ifft=1,nfft 517 rhotmp=rhor(ifft,ispden) 518 rhor(ifft,ispden)=electronpositron%rhor_ep(ifft,ispden) 519 electronpositron%rhor_ep(ifft,ispden)=rhotmp 520 end do 521 end do 522 rhog_ep(:,:)=rhog 523 call fourdp(1,rhog,rhor,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 524 ! If PAW, exchange "positronic" and "electronic" rhoij 525 if (dtset%usepaw==1) then 526 if (size(pawrhoij)>0.and.size(electronpositron%pawrhoij_ep)>0) then 527 ABI_DATATYPE_ALLOCATE(pawrhoij_tmp,(my_natom)) 528 call pawrhoij_alloc(pawrhoij_tmp,pawrhoij(1)%cplex,pawrhoij(1)%nspden,& 529 & pawrhoij(1)%nspinor,pawrhoij(1)%nsppol,dtset%typat,& 530 & pawtab=pawtab,ngrhoij=pawrhoij(1)%ngrhoij,nlmnmix=pawrhoij(1)%lmnmix_sz,& 531 & use_rhoij_=pawrhoij(1)%use_rhoij_,use_rhoijres=pawrhoij(1)%use_rhoijres, & 532 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 533 call pawrhoij_copy(pawrhoij,pawrhoij_tmp) 534 call pawrhoij_copy(electronpositron%pawrhoij_ep,pawrhoij) 535 call pawrhoij_copy(pawrhoij_tmp,electronpositron%pawrhoij_ep) 536 call pawrhoij_free(pawrhoij_tmp) 537 ABI_DATATYPE_DEALLOCATE(pawrhoij_tmp) 538 end if 539 if (allocated(electronpositron%nhat_ep)) then 540 do ispden=1,dtset%nspden 541 do ifft=1,nfft 542 rhotmp=nhat(ifft,ispden) 543 nhat(ifft,ispden)=electronpositron%nhat_ep(ifft,ispden) 544 electronpositron%nhat_ep(ifft,ispden)=rhotmp 545 end do 546 end do 547 end if 548 end if 549 end if 550 551 ! ===== COMPUTE HARTREE POTENTIAL ASSOCIATED TO RHOR_EP 552 if (history_level==4) then 553 call fourdp(1,rhog_ep,electronpositron%rhor_ep,-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,0) 554 end if 555 if (history_level/=-1) then 556 call hartre(1,gsqcut,dtset%usepaw,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog_ep,rprimd,& 557 & electronpositron%vha_ep) 558 electronpositron%vha_ep=-electronpositron%vha_ep 559 else 560 electronpositron%vha_ep=zero 561 end if 562 ABI_DEALLOCATE(rhog_ep) 563 564 ! ---------------------------------------------------------------------- 565 ! Initialize/Update energies 566 ! ---------------------------------------------------------------------- 567 electronpositron%etotal_prev=etotal 568 electronpositron%maxfor_prev=maxfor 569 570 ! Inits/exchange news energies 571 ! Retrieve energy of non-evolving particle(s) 572 if (history_level== 0) then 573 call energies_init(energies) 574 call energies_init(electronpositron%energies_ep) 575 if (dtset%positron>0) energies%e0_electronpositron=etotal_read 576 if (dtset%positron<0) energies%e0_electronpositron=zero 577 else if (history_level== 1) then 578 call energies_init(electronpositron%energies_ep) 579 if (dtset%positron>0) energies%e0_electronpositron=etotal_read 580 else if (history_level== 2) then 581 call energies_copy(energies,electronpositron%energies_ep) 582 call energies_init(energies) 583 energies%e0_electronpositron=electronpositron%e0 584 else if (history_level== 3) then 585 call energies_copy(electronpositron%energies_ep,energies_tmp) 586 call energies_copy(energies,electronpositron%energies_ep) 587 call energies_copy(energies_tmp,energies) 588 energies%e0_electronpositron=electronpositron%e0 589 ! else if (history_level== 4) then 590 end if 591 592 ! Adjust core psps energy 593 if (electronpositron%calctype/=1) energies%e_corepsp=ecore/ucvol 594 595 ! ----------------------------------------------------------------------------------------- 596 ! Update wavefunctions 597 ! If electronpositron%calctype==1: cg are the positronic WFs, cg_ep are the electronic WFs 598 ! If electronpositron%calctype==2: cg are the electronic WFs, cg_ep are the positronic WFs 599 ! ----------------------------------------------------------------------------------------- 600 if (electronpositron%dimcg>0.or.electronpositron%dimcprj>0) then 601 602 if (history_level==0.or.history_level==1) then 603 electronpositron%cg_ep=zero 604 end if 605 606 if (history_level==2) then 607 if (electronpositron%dimcg>0) then 608 do icg=1,electronpositron%dimcg 609 electronpositron%cg_ep(1:2,icg)=cg(1:2,icg) 610 end do 611 end if 612 if (dtset%usepaw==1.and.electronpositron%dimcprj>0) then 613 call pawcprj_copy(cprj,electronpositron%cprj_ep) 614 end if 615 end if 616 617 if (history_level==3) then 618 if (electronpositron%dimcg>0) then 619 do icg=1,electronpositron%dimcg 620 cgtmp(1:2)=electronpositron%cg_ep(1:2,icg) 621 electronpositron%cg_ep(1:2,icg)=cg(1:2,icg) 622 cg(1:2,icg)=cgtmp(1:2) 623 end do 624 end if 625 if (dtset%usepaw==1.and.electronpositron%dimcprj>0) then 626 ABI_ALLOCATE(nlmn,(dtset%natom)) 627 ABI_DATATYPE_ALLOCATE(cprj_tmp,(dtset%natom,electronpositron%dimcprj)) 628 do iatom=1,dtset%natom 629 nlmn(iatom)=cprj(iatom,1)%nlmn 630 end do 631 call pawcprj_alloc(cprj_tmp,cprj(1,1)%ncpgr,nlmn) 632 ABI_DEALLOCATE(nlmn) 633 call pawcprj_copy(electronpositron%cprj_ep,cprj_tmp) 634 call pawcprj_copy(cprj,electronpositron%cprj_ep) 635 call pawcprj_copy(cprj_tmp,cprj) 636 call pawcprj_free(cprj_tmp) 637 ABI_DATATYPE_DEALLOCATE(cprj_tmp) 638 end if 639 end if 640 641 end if ! dimcg>0 or dimcprj>0 642 643 ! ----------------------------------------------------------------------------------------------------------- 644 ! Initialize/Update occupations 645 ! If electronpositron%calctype==1: occ are the positronic occupations, occ_ep are the electronic occupations 646 ! If electronpositron%calctype==2: occ are the electronic occupations, occ_ep are the positronic occupations 647 ! ----------------------------------------------------------------------------------------------------------- 648 ! When needed, precompute electronic occupations with semiconductor occupancies 649 need_scocc=.false. 650 if (electronpositron%dimocc>0.and.electronpositron%calctype==1.and. & 651 & (history_level==0.or.history_level==1)) need_scocc=.true. 652 if (electronpositron%calctype==2.and.ireadwf==0.and. & 653 & (history_level==0.or.history_level==2.or. & 654 & (history_level==3.and.electronpositron%dimocc==0))) need_scocc=.true. 655 if (need_scocc) then 656 nelect=-dtset%charge 657 do iatom=1,dtset%natom 658 nelect=nelect+dtset%ziontypat(dtset%typat(iatom)) 659 end do 660 maxocc=two/real(dtset%nsppol*dtset%nspinor,dp) 661 nocc=(nelect-tol8)/maxocc + 1 662 nocc=min(nocc,dtset%nband(1)*dtset%nsppol) 663 occlast=nelect-maxocc*(nocc-1) 664 ABI_ALLOCATE(scocc,(dtset%nband(1)*dtset%nsppol)) 665 scocc=zero 666 if (1<nocc) scocc(1:nocc-1)=maxocc 667 if (1<=nocc) scocc(nocc)=occlast 668 end if 669 670 ! ===== PREVIOUS OCCUPATIONS OCC_EP: 671 if (electronpositron%dimocc>0) then 672 if (history_level==0.or.history_level==1) then 673 ! ----- Electronic from scratch 674 if (electronpositron%calctype==1) then 675 ! Initialize electronic occupations with semiconductor occupancies 676 do ikpt=1,dtset%nkpt 677 do iband=1,dtset%nband(1) 678 do isppol=1,dtset%nsppol 679 electronpositron%occ_ep(iband+dtset%nband(1)*(ikpt-1+dtset%nkpt*(isppol-1)))=& 680 & scocc(isppol+dtset%nsppol*(iband-1)) 681 end do 682 end do 683 end do 684 end if 685 ! ----- Positronic from scratch 686 if (electronpositron%calctype==1) then 687 ! Initialize positronic occupations with only one positron (or less) 688 electronpositron%occ_ep(:)=zero 689 isppol=1;iocc=1 690 do ikpt=1,dtset%nkpt 691 electronpositron%occ_ep(iocc)=electronpositron%posocc 692 iocc=iocc+dtset%nband(ikpt+dtset%nkpt*(isppol-1)) 693 end do 694 end if 695 end if 696 ! ----- Deduced from occ in memory 697 if (history_level==2) then 698 electronpositron%occ_ep(:)=occ(:) 699 end if 700 end if ! dimocc>0 701 702 ! ===== CURRENT OCCUPATIONS OCC: 703 if (history_level==0.or.history_level==2.or.(history_level==3.and.electronpositron%dimocc==0)) then 704 if (ireadwf==0) then 705 ! ----- Positronic from scratch 706 if (electronpositron%calctype==1) then 707 ! Initialize positronic occupations with only one positron (or less) 708 occ(:)=zero 709 isppol=1;iocc=1 710 do ikpt=1,dtset%nkpt 711 occ(iocc)=electronpositron%posocc 712 iocc=iocc+dtset%nband(ikpt+dtset%nkpt*(isppol-1)) 713 end do 714 end if 715 ! ----- Electronic from scratch 716 if (electronpositron%calctype==2) then 717 ! Initialize electronic occupations with semiconductor occupancies 718 do ikpt=1,dtset%nkpt 719 do iband=1,dtset%nband(1) 720 do isppol=1,dtset%nsppol 721 occ(iband+dtset%nband(1)*(ikpt-1+dtset%nkpt*(isppol-1)))=& 722 & scocc(isppol+dtset%nsppol*(iband-1)) 723 end do 724 end do 725 end do 726 end if 727 end if 728 end if 729 730 ! ===== EXCHANGE POSITRONIC AND ELECTRONIC OCCUPATIONS (CURRENT AND PREVIOUS) 731 if (history_level==3.and.electronpositron%dimocc>0) then 732 do iocc=1,electronpositron%dimocc 733 occtmp=occ(iocc) 734 occ(iocc)=electronpositron%occ_ep(iocc) 735 electronpositron%occ_ep(iocc)=occtmp 736 end do 737 end if 738 739 if (need_scocc) then 740 ABI_DEALLOCATE(scocc) 741 end if 742 743 ! ----------------------------------------------------------------------------------------------------------- 744 ! Initialize/Update eigen energies 745 ! If electronpositron%calctype==1: eigen are the positronic eigen E, eigen_ep are the electronic eigen E 746 ! If electronpositron%calctype==2: eigen are the electronic eigen E, eigen_ep are the positronic eigen E 747 ! ----------------------------------------------------------------------------------------------------------- 748 749 ! ===== PREVIOUS EIGEN ENERGIES EIGEN_EP: 750 if (electronpositron%dimeigen>0) then 751 if (history_level==0.or.history_level==1) then 752 ! ----- Electronic or positronic from scratch 753 electronpositron%eigen_ep(:)=zero 754 end if 755 ! ----- Deduced from eigen in memory 756 if (history_level==2) then 757 electronpositron%eigen_ep(:)=eigen(:) 758 end if 759 end if ! dimeigen>0 760 761 ! ===== CURRENT EIGEN ENERGIES EIGEN: 762 if (history_level==0.or.history_level==2.or.(history_level==3.and.electronpositron%dimeigen==0)) then 763 if (ireadwf==0) then 764 ! ----- Electronic or positronic from scratch 765 eigen(:)=zero 766 end if 767 end if 768 769 ! ===== EXCHANGE POSITRONIC AND ELECTRONIC EIGEN ENERGIES (CURRENT AND PREVIOUS) 770 if (history_level==3.and.electronpositron%dimeigen>0) then 771 do iocc=1,electronpositron%dimeigen 772 eigtmp=eigen(iocc) 773 eigen(iocc)=electronpositron%eigen_ep(iocc) 774 electronpositron%eigen_ep(iocc)=eigtmp 775 end do 776 end if 777 778 ! ============================================= 779 end if ! the type of e-p calculation changes 780 !============================================= 781 782 !------------------------------------------------------------------ 783 !Write messages 784 !------------------------------------------------------------------ 785 if (istep_mix==1.and.dtset%positron/=0) then 786 ! Log message 787 if (electronpositron%calctype==0) then 788 message = 'Were are now performing an electronic ground-state calculation...' 789 else if (electronpositron%calctype==1) then 790 message = 'Were are now performing a positronic ground-state calculation...' 791 else if (electronpositron%calctype==2) then 792 message = 'Were are now performing an electronic ground-state calculation in presence of a positron...' 793 end if 794 MSG_COMMENT(message) 795 ! Output message 796 if (dtset%positron<0) then 797 if (electronpositron%calctype==0) then 798 TypeCalcStrg='ELECTRONIC GROUND-STATE CALCULATION' 799 else if (electronpositron%calctype==1) then 800 TypeCalcStrg='POSITRONIC GROUND-STATE CALCULATION IN PRESENCE OF ELECTRONS AND IONS' 801 else if (electronpositron%calctype==2) then 802 TypeCalcStrg='ELECTRONIC GROUND-STATE CALCULATION IN PRESENCE OF A POSITRON' 803 end if 804 if (istep>1) then 805 write(message,'(2a,i3,2a)') ch10,'TC-DFT STEP ',electronpositron%istep,' - ',trim(TypeCalcStrg) 806 else 807 write(message,'(a,i3,2a)') 'TC-DFT STEP ',electronpositron%istep,' - ',trim(TypeCalcStrg) 808 end if 809 call wrtout(ab_out,message,'COLL') 810 end if 811 end if 812 813 end subroutine setup_positron