TABLE OF CONTENTS
ABINIT/setvtr [ Functions ]
NAME
setvtr
FUNCTION
Set up the trial potential and some energy terms
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (XG, GMR, FJ, MT, EB, SPr) 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
[add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy atindx1(dtset%natom)=index table for atoms, inverse of atindx dtset <type(dataset_type)>=all input variables in this dataset | if =0,1 no xc kernel, =2 spin-averaged (LDA) kernel | densfor_pred=govern the choice of preconditioner for the SCF cycle | iscf=determines the way the SCF cycle is handled | natom=number of atoms in cell. | nspden=number of spin-density components | qprtrb(3)= integer wavevector of possible perturbing potential | in basis of reciprocal lattice translations | typat(natom)=type integer for each atom in cell | vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero, | perturbing potential is added of the form | V(G)=(vprtrb(1)+I*vprtrb(2))/2 at the values G=qprtrb and | (vprtrb(1)-I*vprtrb(2))/2 at G=-qprtrb (integers) | for each type of atom, from psp (used in norm-conserving only) gmet(3,3)=metric tensor for G vecs (in bohr**-2) gprimd(3,3)=dimensional reciprocal space primitive translations gsqcut=cutoff on (k+G)^2 (bohr^-2) (sphere for density and potential) istep=step number in the main loop of scfcv mgfft=maximum size of 1D FFTs moved_rhor=1 if the density was moved just before mpi_enreg=information about MPI parallelization nattyp(ntypat)=number of atoms of each type in cell. nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nhat(nfft,nspden*usepaw)= -PAW only- compensation density nhatgr(nfft,nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise nkxc=second dimension of the array kxc ntypat=number of types of atoms in unit cell. n1xccc=dimension of xccc1d ; 0 if no XC core correction is used n3xccc=dimension of the xccc3d array (0 or nfft). optene=>0 if some additional energies have to be computed pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data ph1d(2,3*(2*mgfft+1)*natom)=phase (structure factor) information. psps <type(pseudopotential_type)>=variables related to pseudopotentials rhog(2,nfft)=Fourier transform of electron density rhor(nfft,nspden)=electron density in electrons/bohr**3. | definition for spin components: | case of nspden = 2 | rhor(:,1) => rho_up + rho_dwn | rhor(:,2) => rho_up | case of nspden = 4 | rhor(:,1) => rho_upup + rho_dwndwn | rhor(:,2:4) => {m_x,m_y,m_z} rmet(3,3)=real space metric (bohr**2) rprimd(3,3)=dimensional primitive translations in real space (bohr) [taug(2,nfftf*dtset%usekden)]=array for Fourier transform of kinetic energy density [taur(nfftf,nspden*dtset%usekden)]=array for kinetic energy density ucvol = unit cell volume (bohr^3) usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
energies <type(energies_type)>=all part of total energy. | e_xc=exchange-correlation energy (hartree) | In case of hybrid compensation algorithm: | e_hybcomp_E0=energy compensation term for hybrid exchange-correlation energy (hartree) at fixed density | e_hybcomp_v0=potential compensation term for hybrid exchange-correlation energy (hartree) at fixed density | e_hybcomp_v=potential compensation term for hybrid exchange-correlation energy (hartree) at self-consistent density ==== if optene==2 or 4 | e_localpsp=local psp energy (hartree) ==== if dtset%icoulomb == 0 | e_ewald=Ewald energy (hartree) ==== if optene>=1 | e_hartree=Hartree part of total energy (hartree) ==== if optene==3 or 4 | e_xcdc=exchange-correlation double-counting energy (hartree) ==== if dtset%vdw_xc == 5 or 6 or 7 | e_vdw_dftd=Dispersion energy from DFT-D Van der Waals correction (hartree) grchempottn(3,natom)=grads of spatially-varying chemical energy (hartree) grewtn(3,natom)=grads of Ewald energy (hartree) grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D2 dispersion (hartree) kxc(nfft,nkxc)=exchange-correlation kernel, will be computed if nkxc/=0 . see routine rhotoxc for a more complete description strsxc(6)=xc contribution to stress tensor (hartree/bohr^3) vxcavg=mean of the vxc potential
SIDE EFFECTS
[electronpositron <type(electronpositron_type)>]=quantities for the electron-positron annihilation (optional argument) moved_atm_inside=1 if the atomic positions were moved inside the SCF loop. vhartr(nfft)=Hartree potential (Hartree) vpsp(nfft)=local psp (Hartree) vtrial(nfft,nspden)= trial potential (Hartree) vxc(nfft,nspden)= xc potential (Hartree) [vxc_hybcomp(nfft,nspden)= compensation xc potential (Hartree) in case of hybrids] Optional output i.e. difference between the hybrid Vxc at fixed density and the auxiliary Vxc at fixed density [vxctau(nfftf,dtset%nspden,4*dtset%usekden)]=derivative of XC energy density with respect to kinetic energy density (metaGGA cases) (optional output) xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
NOTES
In case of PAW calculations:
All computations are done on the fine FFT grid.
All variables (nfft,ngfft,mgfft) refer to this fine FFT grid.
All arrays (densities/potentials...) are computed on this fine FFT grid.
! Developpers have to be careful when introducing others arrays:
they have to be stored on the fine FFT grid.
In case of norm-conserving calculations the FFT grid is the usual FFT grid.
PARENTS
bethe_salpeter,scfcv,screening,sigma
CHILDREN
atm2fft,denspot_set_history,dotprod_vn,ewald,hartre,ionion_realspace ionion_surface,jellium,mag_constr,mkcore,mkcore_alt,mkcore_wvl,mklocl psolver_rhohxc,rhohxcpositron,rhotoxc,spatialchempot,timab,vdw_dftd2 vdw_dftd3,wvl_psitohpsi,wvl_vtrial_abi2big,xcdata_init,xchybrid_ncpp_cc xred2xcart
SOURCE
132 #if defined HAVE_CONFIG_H 133 #include "config.h" 134 #endif 135 136 #include "abi_common.h" 137 138 subroutine setvtr(atindx1,dtset,energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcut,& 139 & istep,kxc,mgfft,moved_atm_inside,moved_rhor,mpi_enreg,& 140 & nattyp,nfft,ngfft,ngrvdw,nhat,nhatgr,nhatgrdim,nkxc,ntypat,n1xccc,n3xccc,& 141 & optene,pawrad,pawtab,ph1d,psps,rhog,rhor,rmet,rprimd,strsxc,& 142 & ucvol,usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,wvl,xccc3d,xred,& 143 & electronpositron,taug,taur,vxc_hybcomp,vxctau,add_tfw) ! optionals arguments 144 145 use defs_basis 146 use defs_datatypes 147 use defs_abitypes 148 use defs_wvltypes 149 use m_profiling_abi 150 use m_errors 151 use m_abi2big 152 use m_xmpi 153 use m_xcdata 154 155 use m_geometry, only : xred2xcart 156 use m_cgtools, only : dotprod_vn 157 use m_ewald, only : ewald 158 use m_energies, only : energies_type 159 use m_electronpositron, only : electronpositron_type,electronpositron_calctype 160 use libxc_functionals, only : libxc_functionals_is_hybrid 161 use m_pawrad, only : pawrad_type 162 use m_pawtab, only : pawtab_type 163 164 #if defined HAVE_BIGDFT 165 use BigDFT_API, only: denspot_set_history 166 #endif 167 168 !This section has been created automatically by the script Abilint (TD). 169 !Do not modify the following lines by hand. 170 #undef ABI_FUNC 171 #define ABI_FUNC 'setvtr' 172 use interfaces_18_timing 173 use interfaces_56_xc 174 use interfaces_62_poisson 175 use interfaces_62_wvl_wfs 176 use interfaces_64_psp 177 use interfaces_67_common, except_this_one => setvtr 178 !End of the abilint section 179 180 implicit none 181 182 !Arguments ------------------------------------ 183 !scalars 184 integer,intent(in) :: istep,mgfft,n1xccc,n3xccc,nfft,ngrvdw,nhatgrdim,nkxc,ntypat 185 integer,intent(in) :: optene,usexcnhat 186 integer,intent(inout) :: moved_atm_inside,moved_rhor 187 logical,intent(in),optional :: add_tfw 188 real(dp),intent(in) :: gsqcut,ucvol 189 real(dp),intent(out) :: vxcavg 190 type(MPI_type),intent(in) :: mpi_enreg 191 type(dataset_type),intent(inout) :: dtset 192 type(electronpositron_type),pointer,optional :: electronpositron 193 type(energies_type),intent(inout) :: energies 194 type(pseudopotential_type),intent(in) :: psps 195 type(wvl_data), intent(inout) :: wvl 196 !arrays 197 integer, intent(in) :: atindx1(dtset%natom),nattyp(ntypat),ngfft(18) 198 real(dp),intent(in) :: gmet(3,3),gprimd(3,3) 199 real(dp),intent(in) :: nhat(nfft,dtset%nspden*psps%usepaw) 200 real(dp),intent(in) :: nhatgr(:,:,:) !(nfft,dtset%nspden,3*nhatgrdim) 201 real(dp),intent(in) :: rhog(2,nfft) 202 real(dp),intent(in),optional :: taug(2,nfft*dtset%usekden) 203 real(dp),intent(in) :: rmet(3,3),rprimd(3,3) 204 real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom) 205 real(dp),intent(inout) :: rhor(nfft,dtset%nspden),vhartr(nfft),vpsp(nfft) 206 real(dp),intent(inout),optional :: taur(nfft,dtset%nspden*dtset%usekden) 207 real(dp),intent(inout) :: vtrial(nfft,dtset%nspden),vxc(nfft,dtset%nspden) 208 real(dp),intent(out),optional :: vxctau(nfft,dtset%nspden,4*dtset%usekden) 209 real(dp),intent(out),optional :: vxc_hybcomp(:,:) ! (nfft,nspden) 210 real(dp),intent(inout) :: xccc3d(n3xccc) 211 real(dp),intent(in) :: xred(3,dtset%natom) 212 real(dp),intent(out) :: grchempottn(3,dtset%natom) 213 real(dp),intent(out) :: grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(nfft,nkxc),strsxc(6) 214 type(pawtab_type),intent(in) :: pawtab(ntypat*dtset%usepaw) 215 type(pawrad_type),intent(in) :: pawrad(ntypat*dtset%usepaw) 216 217 !Local variables------------------------------- 218 !scalars 219 integer :: coredens_method,mpi_comm_sphgrid,nk3xc 220 integer :: iatom,ifft,ipositron,ispden,nfftot 221 integer :: optatm,optdyfr,opteltfr,optgr,option,option_eff,optn,optn2,optstr,optv,vloc_method 222 real(dp) :: doti,e_xcdc_vxctau,ebb,ebn,evxc,ucvol_local,rpnrm 223 logical :: add_tfw_,is_hybrid_ncpp,non_magnetic_xc,with_vxctau,wvlbigdft 224 real(dp), allocatable :: xcart(:,:) 225 character(len=500) :: message 226 type(xcdata_type) :: xcdata,xcdatahyb 227 !arrays 228 real(dp),parameter :: identity(1:4)=(/1._dp,1._dp,0._dp,0._dp/) 229 real(dp) :: dummy6(6),tsec(2) 230 real(dp) :: grewtn_fake(3,1) 231 real(dp) :: dummy_in(0) 232 real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0) 233 real(dp) :: strn_dummy6(6), strv_dummy6(6) 234 real(dp) :: vzeeman(4) 235 real(dp),allocatable :: grtn(:,:),dyfr_dum(:,:,:),gr_dum(:,:) 236 real(dp),allocatable :: rhojellg(:,:),rhojellr(:),rhowk(:,:),vjell(:) 237 real(dp),allocatable :: Vmagconstr(:,:),rhog_dum(:,:) 238 239 ! ********************************************************************* 240 241 call timab(91,1,tsec) 242 243 ! Initialise non_magnetic_xc for rhohxc 244 non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14) 245 246 !Check that usekden is not 0 if want to use vxctau 247 with_vxctau = (present(vxctau).and.present(taur).and.(dtset%usekden/=0)) 248 249 !Check if we're in hybrid norm conserving pseudopotential with a core correction 250 is_hybrid_ncpp=(dtset%usepaw==0 .and. n3xccc/=0 .and. & 251 & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid())) 252 253 !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed 254 wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1) 255 256 !Get size of FFT grid 257 nfftot=PRODUCT(ngfft(1:3)) 258 259 !mpi 260 mpi_comm_sphgrid=mpi_enreg%comm_fft 261 if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl 262 263 !Test electron-positron case 264 ipositron=0;if (present(electronpositron)) ipositron=electronpositron_calctype(electronpositron) 265 266 !Test addition of Weiszacker gradient correction to Thomas-Fermi kin energy 267 add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw 268 269 !Get Ewald energy and Ewald forces, as well as vdW-DFTD energy and forces, and chemical potential energy and forces. 270 !------------------------------------------------------------------------------------------------------------------- 271 call timab(5,1,tsec) 272 if (ipositron/=1) then 273 if (dtset%icoulomb == 0 .or. (dtset%usewvl == 0 .and. dtset%icoulomb == 2)) then 274 ! Periodic system, need to compute energy and forces due to replica and 275 ! to correct the shift in potential calculation. 276 call ewald(energies%e_ewald,gmet,grewtn,dtset%natom,ntypat,rmet,dtset%typat,ucvol,xred,psps%ziontypat) 277 ! For a periodic system bearing a finite charge, the monopole correction to the 278 ! energy is relevant. 279 ! See Leslie and Gillan, JOURNAL OF PHYSICS C-SOLID STATE PHYSICS 18, 973 (1985) 280 if(abs(dtset%charge)>tol8) then 281 call ewald(energies%e_monopole,gmet,grewtn_fake,1,1,rmet,(/1/),ucvol,(/0.0_dp,0.0_dp,0.0_dp/),(/dtset%charge/)) 282 energies%e_monopole=-energies%e_monopole 283 end if 284 else if (dtset%icoulomb == 1) then 285 ! In a non periodic system (real space computation), the G=0 divergence 286 ! doesn't occur and ewald is not needed. Only the ion/ion interaction 287 ! energy is relevant and used as Ewald energy and gradient. 288 call ionion_realSpace(dtset, energies%e_ewald, grewtn, rprimd, xred, psps%ziontypat) 289 else if (dtset%icoulomb == 2) then 290 call ionion_surface(dtset, energies%e_ewald, grewtn, mpi_enreg%me_wvl, mpi_enreg%nproc_wvl, rprimd, & 291 & wvl%descr, wvl%den, xred) 292 end if 293 if (dtset%nzchempot>0) then 294 call spatialchempot(energies%e_chempot,dtset%chempot,grchempottn,dtset%natom,ntypat,dtset%nzchempot,dtset%typat,xred) 295 end if 296 if (dtset%vdw_xc==5.and.ngrvdw==dtset%natom) then 297 call vdw_dftd2(energies%e_vdw_dftd,dtset%ixc,dtset%natom,ntypat,1,dtset%typat,rprimd,& 298 & dtset%vdw_tol,xred,psps%znucltypat,fred_vdw_dftd2=grvdw) 299 end if 300 if ((dtset%vdw_xc==6.or.dtset%vdw_xc==7).and.ngrvdw==dtset%natom) then 301 call vdw_dftd3(energies%e_vdw_dftd,dtset%ixc,dtset%natom,& 302 & ntypat,1,dtset%typat,rprimd,dtset%vdw_xc,dtset%vdw_tol,dtset%vdw_tol_3bt,& 303 & xred,psps%znucltypat,fred_vdw_dftd3=grvdw) 304 end if 305 else 306 energies%e_ewald=zero 307 energies%e_chempot=zero 308 grchempottn=zero 309 grewtn=zero 310 energies%e_vdw_dftd=zero 311 if (ngrvdw>0) grvdw=zero 312 end if 313 call timab(5,2,tsec) 314 315 !Compute parts of total energy depending on potentials 316 !-------------------------------------------------------------- 317 if (dtset%usewvl == 0) then 318 ucvol_local = ucvol 319 #if defined HAVE_BIGDFT 320 else 321 ! We need to tune the volume when wavelets are used because, not all FFT points are used. 322 ! ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3) 323 ucvol_local = product(wvl%den%denspot%dpbox%hgrids) * real(product(wvl%den%denspot%dpbox%ndims), dp) 324 #endif 325 end if 326 327 !Determine by which method the local ionic potential and/or the pseudo core charge density 328 ! have to be computed 329 !Local ionic potential: 330 ! Method 1: PAW 331 ! Method 2: Norm-conserving PP, icoulomb>0, wavelets 332 vloc_method=1;if (psps%usepaw==0) vloc_method=2 333 if (dtset%icoulomb>0) vloc_method=2 334 if (psps%usewvl==1) vloc_method=2 335 !Pseudo core charge density: 336 ! Method 1: PAW, nc_xccc_gspace 337 ! Method 2: Norm-conserving PP, wavelets 338 coredens_method=1;if (psps%usepaw==0) coredens_method=2 339 if (psps%nc_xccc_gspace==1) coredens_method=1 340 if (psps%nc_xccc_gspace==0) coredens_method=2 341 if (psps%usewvl==1) coredens_method=2 342 343 !Local ionic potential and/or pseudo core charge by method 1 344 if (vloc_method==1.or.coredens_method==1) then 345 call timab(552,1,tsec) 346 optv=0;if (vloc_method==1) optv=1 347 optn=0;if (coredens_method==1) optn=n3xccc/nfft 348 optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=1 349 call atm2fft(atindx1,xccc3d,vpsp,dummy_out1,dummy_out2,dummy_out3,dummy_in,& 350 & gmet,gprimd,dummy_out4,dummy_out5,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,& 351 & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,& 352 & dummy_in,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,& 353 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 354 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 355 call timab(552,2,tsec) 356 end if 357 358 !Local ionic potential by method 2 359 if (vloc_method==2) then 360 option=1 361 ABI_ALLOCATE(gr_dum,(3,dtset%natom)) 362 ABI_ALLOCATE(dyfr_dum,(3,3,dtset%natom)) 363 ABI_ALLOCATE(rhog_dum,(2,nfft)) 364 call mklocl(dtset,dyfr_dum,energies%e_localpsp,gmet,gprimd,& 365 & gr_dum,gsqcut,dummy6,mgfft,mpi_enreg,dtset%natom,nattyp,& 366 & nfft,ngfft,dtset%nspden,ntypat,option,pawtab,ph1d,psps,& 367 & dtset%qprtrb,rhog_dum,rhor,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred) 368 ABI_DEALLOCATE(gr_dum) 369 ABI_DEALLOCATE(dyfr_dum) 370 ABI_DEALLOCATE(rhog_dum) 371 end if 372 373 !3D pseudo core electron density xccc3d by method 2 374 if (coredens_method==2.and.n1xccc/=0) then 375 call timab(91,2,tsec) 376 call timab(92,1,tsec) 377 option=1 378 ABI_ALLOCATE(gr_dum,(3,dtset%natom)) 379 ABI_ALLOCATE(dyfr_dum,(3,3,dtset%natom)) 380 if (psps%usewvl==0.and.psps%usepaw==0.and.dtset%icoulomb==0) then 381 call mkcore(dummy6,dyfr_dum,gr_dum,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,& 382 & ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,dtset%typat,ucvol,& 383 & vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred) 384 else if (psps%usewvl==0.and.(psps%usepaw==1.or.dtset%icoulomb==1)) then 385 call mkcore_alt(atindx1,dummy6,dyfr_dum,gr_dum,dtset%icoulomb,mpi_enreg,dtset%natom,& 386 & nfft,dtset%nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,& 387 & ucvol,vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred,pawrad,pawtab,psps%usepaw) 388 else if (psps%usewvl==1.and.psps%usepaw==1) then 389 #if defined HAVE_BIGDFT 390 ! call mkcore_wvl_old(atindx1,dummy6,dyfr_dum,wvl%descr%atoms%astruct%geocode,gr_dum,wvl%descr%h,& 391 ! & dtset%natom,nattyp,nfft,wvl%den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl,:),& 392 ! & dtset%nspden,ntypat,wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,wvl%descr%Glr%d%n2,& 393 ! & wvl%descr%Glr%d%n2i,wvl%descr%Glr%d%n3,wvl%den%denspot%dpbox%n3pi,n3xccc,option,& 394 ! & pawrad,pawtab,psps%gth_params%psppar,rprimd,ucvol_local,vxc,xccc3d,xred,& 395 ! & mpi_comm_wvl=mpi_enreg%comm_wvl) 396 call mkcore_wvl(atindx1,dummy6,gr_dum,dtset%natom,nattyp,nfft,dtset%nspden,ntypat,& 397 & n1xccc,n3xccc,option,pawrad,pawtab,rprimd,vxc,psps%xccc1d,xccc3d,& 398 & psps%xcccrc,xred,wvl%den,wvl%descr,mpi_comm_wvl=mpi_enreg%comm_wvl) 399 #endif 400 end if 401 ABI_DEALLOCATE(gr_dum) 402 ABI_DEALLOCATE(dyfr_dum) 403 call timab(92,2,tsec) 404 call timab(91,1,tsec) 405 end if 406 407 !Adds the jellium potential to the local part of ionic potential 408 if (dtset%jellslab/=0) then 409 ABI_ALLOCATE(vjell,(nfft)) 410 ABI_ALLOCATE(rhojellg,(2,nfft)) 411 ABI_ALLOCATE(rhojellr,(nfft)) 412 option=1 413 call jellium(gmet,gsqcut,mpi_enreg,nfft,ngfft,dtset%nspden,option,mpi_enreg%paral_kgb,& 414 & dtset%slabwsrad,rhojellg,rhojellr,rprimd,vjell,dtset%slabzbeg,dtset%slabzend) 415 ! Compute background-background energy 416 call dotprod_vn(1,rhojellr,ebb,doti,nfft,nfftot,1,1,vjell,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid) 417 ebb=half*ebb 418 ! Compute electrostatic energy between background and nuclei before adding vjell to vpsp 419 call dotprod_vn(1,rhojellr,ebn,doti,nfft,nfftot,1,1,vpsp,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid) 420 ! Update e_ewald with ebb and ebn 421 energies%e_ewald=energies%e_ewald+ebb+ebn 422 ! Compute gradient of ebn wrt tn 423 ! This is not yet coded for usewvl or icoulomb=1 424 if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then 425 write(message,'(3a)')& 426 & 'The computation of forces due to jellium background',ch10,& 427 & 'has to be verified in the PAW formalism.' 428 MSG_WARNING(message) 429 430 ABI_ALLOCATE(grtn,(3,dtset%natom)) 431 optatm=0;optdyfr=0;opteltfr=0;optgr=1;optstr=0;optv=1;optn=0;optn2=1 432 call atm2fft(atindx1,dummy_out1,vpsp,dummy_out2,dummy_out3,dummy_out4,dummy_in,& 433 & gmet,gprimd,dummy_out5,grtn,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,& 434 & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,& 435 & rhojellg,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,& 436 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 437 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 438 439 ! Update grewtn with gradient of ebn wrt tn 440 do iatom=1,dtset%natom 441 grewtn(1:3,iatom)=grewtn(1:3,iatom)+grtn(1:3,iatom) 442 end do 443 ABI_DEALLOCATE(grtn) 444 else ! of usepaw==1 445 option=2 446 ABI_ALLOCATE(dyfr_dum,(3,3,dtset%natom)) 447 ABI_ALLOCATE(grtn,(3,dtset%natom)) 448 call mklocl(dtset,dyfr_dum,energies%e_localpsp,gmet,gprimd,& 449 & grtn,gsqcut,dummy6,mgfft,mpi_enreg,dtset%natom,nattyp,& 450 & nfft,ngfft,1,ntypat,option,pawtab,ph1d,psps,dtset%qprtrb,rhojellg,& 451 & rhojellr,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred) 452 ! Update grewtn with gradient of ebn wrt tn (reestablish order of atoms) 453 do iatom=1,dtset%natom 454 grewtn(1:3,atindx1(iatom))=grewtn(1:3,atindx1(iatom))+grtn(1:3,iatom) 455 end do 456 ABI_DEALLOCATE(dyfr_dum) 457 ABI_DEALLOCATE(grtn) 458 end if ! of usepaw==1 459 vpsp(:)=vpsp(:)+vjell(:) 460 ABI_DEALLOCATE(vjell) 461 ABI_DEALLOCATE(rhojellg) 462 ABI_DEALLOCATE(rhojellr) 463 end if 464 465 !Additional stuff for electron-positron calculation 466 !Compute the electronic/positronic local (Hartree) potential 467 if (ipositron==1) vpsp=-vpsp 468 469 !If we are at the initialisation, or 470 !if the atom positions has changed and the non-linear core correction 471 !is included, or the rhor has changed, one needs to compute the xc stuff. 472 !One needs also to compute the Hartree stuff if the density changed, 473 !or at initialisation. 474 !-------------------------------------------------------------- 475 476 if(istep==1 .or. n1xccc/=0 .or. moved_rhor==1 .or. dtset%positron<0 .or. mod(dtset%fockoptmix,100)==11) then 477 478 option=0 479 if(istep==1 .or. moved_rhor==1 .or. dtset%positron<0 .or. mod(dtset%fockoptmix,100)==11) option=1 480 if (nkxc>0) option=2 481 if (dtset%iscf==-1) option=-2 482 if (ipositron/=1) then 483 if (dtset%icoulomb == 0 .and. dtset%usewvl == 0) then 484 if(option/=0 .and. option/=10)then 485 call hartre(1,gsqcut,psps%usepaw,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog,rprimd,vhartr) 486 end if 487 call xcdata_init(xcdata,dtset=dtset) 488 if(mod(dtset%fockoptmix,100)==11)then 489 xcdatahyb=xcdata 490 ! Setup the auxiliary xc functional information 491 call xcdata_init(xcdata,dtset=dtset,auxc_ixc=0,ixc=dtset%auxc_ixc) 492 end if 493 ! Use the periodic solver to compute Hxc 494 nk3xc=1 495 if (ipositron==0) then 496 497 ! Compute energies%e_xc and associated quantities 498 if(.not.is_hybrid_ncpp .or. mod(dtset%fockoptmix,100)==11)then 499 ! Not yet able to deal fully with the full XC kernel in case of GGA + spin 500 option_eff=option 501 if(xcdata%xclevel==2.and.(nkxc==3-2*mod(xcdata%nspden,2))) option_eff=12 502 call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,& 503 & nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,& 504 & option_eff,dtset%paral_kgb,rhor,rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,& 505 & taug=taug,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_) 506 else 507 ! Only when is_hybrid_ncpp, and moreover, the xc functional is not the auxiliary xc functional, then call xchybrid_ncpp_cc 508 call xchybrid_ncpp_cc(dtset,energies%e_xc,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,& 509 & strsxc,vxcavg,xccc3d,vxc=vxc) 510 end if 511 512 ! Possibly compute energies%e_hybcomp_E0 513 if(mod(dtset%fockoptmix,100)==11)then 514 ! This call to rhotoxc uses the hybrid xc functional 515 if(.not.is_hybrid_ncpp)then 516 ! Not yet able to deal fully with the full XC kernel in case of GGA + spin 517 option_eff=option 518 if(xcdata%xclevel==2.and.(nkxc==3-2*mod(xcdata%nspden,2))) option_eff=12 519 call rhotoxc(energies%e_hybcomp_E0,kxc,mpi_enreg,nfft,ngfft,& 520 & nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,& 521 & option,dtset%paral_kgb,rhor,rprimd,strsxc,usexcnhat,vxc_hybcomp,vxcavg,xccc3d,xcdatahyb,& 522 & taug=taug,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_) 523 else 524 call xchybrid_ncpp_cc(dtset,energies%e_hybcomp_E0,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,& 525 & strsxc,vxcavg,xccc3d,vxc=vxc_hybcomp) 526 end if 527 528 ! Combine hybrid and auxiliary quantities 529 energies%e_xc=energies%e_xc*dtset%auxc_scal 530 energies%e_hybcomp_E0=energies%e_hybcomp_E0-energies%e_xc 531 vxc(:,:)=vxc(:,:)*dtset%auxc_scal 532 vxc_hybcomp(:,:)=vxc_hybcomp(:,:)-vxc(:,:) 533 534 end if 535 536 else if (ipositron==2) then 537 ! Not yet able to deal fully with the full XC kernel in case of GGA + spin 538 option_eff=option 539 if(xcdata%xclevel==2.and.(nkxc==3-2*mod(xcdata%nspden,2))) option_eff=12 540 call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,& 541 & nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,& 542 & option,dtset%paral_kgb,rhor,rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,& 543 & taug=taug,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_,& 544 & electronpositron=electronpositron) 545 end if 546 547 elseif(.not. wvlbigdft) then 548 ! Use the free boundary solver 549 call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, & 550 & dtset%icoulomb, dtset%ixc, & 551 & mpi_enreg, nfft, ngfft,& 552 & nhat,psps%usepaw,& 553 & dtset%nscforder,dtset%nspden,n3xccc,rhor,rprimd, & 554 & usexcnhat,psps%usepaw,dtset%usewvl,vhartr, vxc, & 555 & vxcavg,wvl%descr,wvl%den,wvl%e,& 556 & xccc3d,dtset%xclevel,dtset%xc_denpos) 557 end if 558 else 559 energies%e_xc=zero 560 call rhohxcpositron(electronpositron,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat,nkxc,dtset%nspden,n3xccc,& 561 & dtset%paral_kgb,rhor,strsxc,ucvol,usexcnhat,psps%usepaw,vhartr,vxc,vxcavg,xccc3d,dtset%xc_denpos) 562 end if 563 if (ipositron/=0) then 564 if (optene>=1) then 565 call dotprod_vn(1,rhor,electronpositron%e_hartree,doti,& 566 & nfft,nfftot,1,1,electronpositron%vha_ep,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid) 567 end if 568 vhartr=vhartr+electronpositron%vha_ep 569 end if 570 end if 571 572 !Compute the trial potential 573 !------------------------------------------------------------- 574 if (.not. wvlbigdft) then 575 ! Now, compute trial Hxc potential. Local psp potential will be added later. 576 if(moved_atm_inside==0 .or.dtset%iscf>=10) then 577 578 ! Compute starting Hxc potential. 579 ! Multiply by identity, should not change anything if nspden /= 4 580 do ispden=1,dtset%nspden 581 vtrial(:,ispden)=vhartr(:)*identity(ispden)+vxc(:,ispden) 582 end do 583 584 else 585 586 ! One should be here only when moved_atm_inside==1 587 ! The (H)xc now added corrects the previous one. 588 if(dtset%densfor_pred==1)then 589 ! xc was substracted off. This should be rationalized later 590 do ispden=1,dtset%nspden 591 vtrial(:,ispden)=vtrial(:,ispden)+vxc(:,ispden) 592 end do 593 else if(abs(dtset%densfor_pred)==2.or.abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)then 594 ! Hxc was substracted off. This should be rationalized later 595 do ispden=1,dtset%nspden 596 vtrial(:,ispden)=vtrial(:,ispden)+vhartr(:)*identity(ispden)+vxc(:,ispden) 597 end do 598 end if 599 end if 600 601 ! Adds the local part of the potential 602 if ((moved_atm_inside==0).or.(dtset%densfor_pred/=3)) then 603 do ispden=1,min(2,dtset%nspden) 604 do ifft=1,nfft 605 vtrial(ifft,ispden)=vtrial(ifft,ispden)+vpsp(ifft) 606 end do 607 end do 608 end if 609 610 ! Adds the compensating vxc for hybrids 611 if(mod(dtset%fockoptmix,100)==11)then 612 vtrial(:,:)=vtrial(:,:)+vxc_hybcomp(:,:) 613 end if 614 615 if(dtset%usewvl==1) then 616 call wvl_vtrial_abi2big(1,vtrial,wvl%den) 617 end if 618 619 else 620 621 ! Compute with covering comms the different part of the potential. 622 #if defined HAVE_BIGDFT 623 ! Copy e_ewald. 624 wvl%e%energs%eion = energies%e_ewald 625 ! Setup the mixing, if necessary 626 call denspot_set_history(wvl%den%denspot,dtset%iscf,dtset%nsppol, & 627 & wvl%den%denspot%dpbox%ndims(1),wvl%den%denspot%dpbox%ndims(2)) 628 #endif 629 ABI_ALLOCATE(xcart,(3, dtset%natom)) 630 call xred2xcart(dtset%natom, rprimd, xcart, xred) 631 call wvl_psitohpsi(dtset%diemix,energies%e_exactX, energies%e_xc, energies%e_hartree, & 632 & energies%e_kinetic, energies%e_localpsp, energies%e_nonlocalpsp, energies%e_sicdc, & 633 & istep, 1, dtset%iscf, mpi_enreg%me_wvl, dtset%natom, dtset%nfft, mpi_enreg%nproc_wvl, dtset%nspden, & 634 & rpnrm, .true.,evxc, wvl,.true., xcart, strsxc,& 635 & vtrial, vxc) 636 if (optene==3.or.optene==4) energies%e_xcdc=evxc 637 ABI_DEALLOCATE(xcart) 638 639 end if 640 641 !Add the zeeman field to vtrial 642 if (any(abs(dtset%zeemanfield(:))>tol8)) then 643 vzeeman(:) = zero ! vzeeman_ij = -1/2*sigma_ij^alpha*B_alpha 644 if(dtset%nspden==2)then 645 vzeeman(1) = -half*dtset%zeemanfield(3) ! v_dwndwn = -1/2*B_z 646 vzeeman(2) = half*dtset%zeemanfield(3) ! v_upup = 1/2*B_z 647 do ifft=1,nfft 648 vtrial(ifft,1) = vtrial(ifft,1) + vzeeman(1) !SPr: added 1st component 649 vtrial(ifft,2) = vtrial(ifft,2) + vzeeman(2) 650 end do !ifft 651 end if 652 if(dtset%nspden==4)then 653 vzeeman(1)=-half*dtset%zeemanfield(3) ! v_dwndwn => v_11 654 vzeeman(2)= half*dtset%zeemanfield(3) ! v_upup => v_22 655 vzeeman(3)=-half*dtset%zeemanfield(1) ! Re(v_dwnup) = Re(v_updwn) => Re(v_12) 656 vzeeman(4)= half*dtset%zeemanfield(2) ! Im(v_dwnup) =-Im(v_dwnup) => Im(v_12) 657 do ispden=1,dtset%nspden 658 do ifft=1,nfft 659 vtrial(ifft,ispden) = vtrial(ifft,ispden) + vzeeman(ispden) 660 end do 661 end do 662 end if 663 end if 664 665 !Compute the constrained potential for the magnetic moments 666 if (dtset%magconon==1.or.dtset%magconon==2) then 667 ABI_ALLOCATE(Vmagconstr, (nfft,dtset%nspden)) 668 Vmagconstr = zero 669 call mag_constr(dtset%natom,dtset%spinat,dtset%nspden,dtset%magconon,dtset%magcon_lambda,rprimd, & 670 & mpi_enreg,nfft,ngfft,dtset%ntypat,dtset%ratsph,rhor,dtset%typat,Vmagconstr,xred) 671 if(dtset%nspden==4)then 672 do ispden=1,dtset%nspden ! (SPr: both components should be used? EB: Yes it should be the case, corrected now) 673 do ifft=1,nfft 674 vtrial(ifft,ispden) = vtrial(ifft,ispden) + Vmagconstr(ifft,ispden) 675 end do !ifft 676 end do !ispden 677 else if(dtset%nspden==2)then 678 do ifft=1,nfft 679 ! TODO : MJV: check that magnetic constraint works also for nspden 2 or add input variable condition 680 ! EB: ispden=2 is rho_up only: to be tested 681 ! SPr: for ispden=2, both components should be used (e.g. see definition for vzeeman)? 682 vtrial(ifft,1) = vtrial(ifft,1) + Vmagconstr(ifft,1) !SPr: added the first component here 683 vtrial(ifft,2) = vtrial(ifft,2) + Vmagconstr(ifft,2) 684 end do !ifft 685 end if 686 ABI_DEALLOCATE(Vmagconstr) 687 end if 688 689 !Compute parts of total energy depending on potentials 690 !-------------------------------------------------------------- 691 692 !For icoulomb==0 or usewvl Ehartree is calculated in psolver_rhohxc(). 693 !For PAW we recalculate this since nhat was not taken into account 694 !in psolver_rhohxc: E_H= int v_H (n+nhat) dr 695 696 if (optene>=1 .and. .not. wvlbigdft .and. (dtset%icoulomb==0 .or. dtset%usepaw==1 ) ) then 697 ! Compute Hartree energy ehart 698 ! Already available in the Psolver case through psolver_rhohxc(). 699 if (ipositron/=1) then 700 call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol_local,& 701 & mpi_comm_sphgrid=mpi_comm_sphgrid) 702 if (ipositron==0) energies%e_hartree = half * energies%e_hartree 703 if (ipositron==2) energies%e_hartree = half * (energies%e_hartree-electronpositron%e_hartree) 704 else 705 energies%e_hartree=zero 706 end if 707 end if 708 709 if(mod(dtset%fockoptmix,100)==11)then 710 if (.not. wvlbigdft) then 711 call dotprod_vn(1,rhor,energies%e_hybcomp_v0,doti,nfft,nfftot,1,1,vxc_hybcomp,ucvol_local,& 712 & mpi_comm_sphgrid=mpi_comm_sphgrid) 713 energies%e_hybcomp_v=energies%e_hybcomp_v0 714 end if 715 end if 716 717 if (optene==2.or.optene==4 .and. .not. wvlbigdft) then 718 ! Compute local psp energy eei 719 call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfft,nfftot,1,1,vpsp,ucvol_local,& 720 & mpi_comm_sphgrid=mpi_comm_sphgrid) 721 end if 722 723 if (optene==3.or.optene==4 .and. .not. wvlbigdft) then 724 ! Compute double-counting XC energy enxcdc 725 if (ipositron/=1) then 726 if (dtset%usepaw==0.or.usexcnhat/=0) then 727 call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,& 728 & mpi_comm_sphgrid=mpi_comm_sphgrid) 729 if (with_vxctau) then 730 call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfft,nfftot,dtset%nspden,1,vxctau(:,:,1),& 731 & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid) 732 energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau 733 end if 734 else 735 ABI_ALLOCATE(rhowk,(nfft,dtset%nspden)) 736 rhowk=rhor-nhat 737 call dotprod_vn(1,rhowk,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,& 738 & mpi_comm_sphgrid=mpi_comm_sphgrid) 739 ABI_DEALLOCATE(rhowk) 740 end if 741 if (ipositron==2) energies%e_xcdc=energies%e_xcdc-electronpositron%e_xcdc 742 else 743 energies%e_xcdc=zero 744 end if 745 end if 746 747 !-------------------------------------------------------------- 748 749 !The initialisation for the new atomic positions has been done 750 moved_atm_inside=0 751 752 call timab(91,2,tsec) 753 754 end subroutine setvtr