TABLE OF CONTENTS
- ABINIT/ionion_realspace
- ABINIT/ionion_surface
- ABINIT/m_setvtr
- m_setvtr/setvtr
- m_setvtr/spatialchempot
ABINIT/ionion_realspace [ Functions ]
NAME
ionion_realspace
FUNCTION
Compute the ion/ion interaction energies and forces in real space case. Use ewald() instead if computations are done in reciprocal space since it also includes the correction for the shift done in potentials calculations and includes replica interactions.
INPUTS
dtset <type(dataset_type)>=all input variables in this dataset rmet(3,3)=metric tensor in real space (bohr^2) xred(3,natom)=relative coords of atoms in unit cell (dimensionless) zion(ntypat)=charge on each type of atom (real number)
OUTPUT
eew=final ion/ion energy in hartrees grewtn(3,natom)=grads of ion/ion wrt xred(3,natom), hartrees.
SOURCE
970 subroutine ionion_realSpace(dtset, eew, grewtn, rprimd, xred, zion) 971 972 !Arguments ------------------------------------ 973 !scalars 974 real(dp),intent(out) :: eew 975 type(dataset_type),intent(in) :: dtset 976 !arrays 977 real(dp),intent(in) :: rprimd(3,3),zion(dtset%ntypat) 978 real(dp),intent(in) :: xred(3,dtset%natom) 979 real(dp),intent(out) :: grewtn(3,dtset%natom) 980 981 !Local variables------------------------------- 982 !scalars 983 integer :: ia1,ia2,iatom,igeo 984 real(dp) :: r 985 !arrays 986 real(dp),allocatable :: grew_cart(:,:),xcart(:,:) 987 988 ! ************************************************************************* 989 990 !Store xcart for each atom 991 ABI_MALLOC(xcart,(3, dtset%natom)) 992 call xred2xcart(dtset%natom, rprimd, xcart, xred) 993 994 !Summing the interaction between ions. 995 eew = 0._dp 996 do ia1 = 1, dtset%natom, 1 997 do ia2 = ia1 + 1, dtset%natom, 1 998 r = sqrt((xcart(1, ia1) - xcart(1, ia2)) ** 2 + & 999 & (xcart(2, ia1) - xcart(2, ia2)) ** 2 + & 1000 & (xcart(3, ia1) - xcart(3, ia2)) ** 2) 1001 eew = eew + zion(dtset%typat(ia1)) * zion(dtset%typat(ia2)) / r 1002 end do 1003 end do 1004 1005 !Allocate temporary array to store cartesian gradients. 1006 ABI_MALLOC(grew_cart,(3, dtset%natom)) 1007 1008 !Summing the forces for each atom 1009 do ia1 = 1, dtset%natom, 1 1010 grew_cart(:, ia1) = 0._dp 1011 do ia2 = 1, dtset%natom, 1 1012 if (ia1 /= ia2) then 1013 r = (xcart(1, ia1) - xcart(1, ia2)) ** 2 + & 1014 & (xcart(2, ia1) - xcart(2, ia2)) ** 2 + & 1015 & (xcart(3, ia1) - xcart(3, ia2)) ** 2 1016 do igeo = 1, 3, 1 1017 grew_cart(igeo, ia1) = grew_cart(igeo, ia1) - (xcart(igeo, ia1) - xcart(igeo, ia2)) * & 1018 & zion(dtset%typat(ia1)) * zion(dtset%typat(ia2)) / (r ** 1.5_dp) 1019 end do 1020 end if 1021 end do 1022 end do 1023 1024 ABI_FREE(xcart) 1025 1026 !Transform cartesian gradients to reduced gradients. 1027 do iatom = 1, dtset%natom, 1 1028 do igeo = 1, 3, 1 1029 grewtn(igeo, iatom) = rprimd(1, igeo) * grew_cart(1, iatom) + & 1030 & rprimd(2, igeo) * grew_cart(2, iatom) + & 1031 & rprimd(3, igeo) * grew_cart(3, iatom) 1032 end do 1033 end do 1034 ABI_FREE(grew_cart) 1035 1036 end subroutine ionion_realSpace
ABINIT/ionion_surface [ Functions ]
NAME
ionion_surface
FUNCTION
Compute the ion/ion interaction energies and forces in real space case. Use ewald() instead if computations are done in reciprocal space since it also includes the correction for the shift done in potentials calculations and includes replica interactions.
INPUTS
dtset <type(dataset_type)>=all input variables in this dataset rmet(3,3)=metric tensor in real space (bohr^2) xred(3,natom)=relative coords of atoms in unit cell (dimensionless) zion(ntypat)=charge on each type of atom (real number)
OUTPUT
eew=final ion/ion energy in hartrees grewtn(3,natom)=grads of ion/ion wrt xred(3,natom), hartrees.
SOURCE
1062 subroutine ionion_surface(dtset, eew, grewtn, me, nproc, rprimd, wvl, wvl_den, xred) 1063 1064 #if defined HAVE_BIGDFT 1065 use BigDFT_API, only: IonicEnergyandForces 1066 #endif 1067 1068 !Arguments ------------------------------------ 1069 !scalars 1070 integer, intent(in) :: me, nproc 1071 real(dp),intent(out) :: eew 1072 type(dataset_type),intent(in) :: dtset 1073 type(wvl_internal_type), intent(in) :: wvl 1074 type(wvl_denspot_type), intent(inout) :: wvl_den 1075 !arrays 1076 real(dp),intent(in) :: rprimd(3,3) 1077 real(dp),intent(in) :: xred(3,dtset%natom) 1078 real(dp),intent(out) :: grewtn(3,dtset%natom) 1079 1080 !Local variables------------------------------- 1081 !scalars 1082 integer :: dispersion, iatom, igeo 1083 real(dp) :: psoffset 1084 !arrays 1085 real(dp),allocatable :: xcart(:,:) 1086 real(dp),pointer :: grew_cart(:,:),fdisp(:,:) 1087 #if defined HAVE_BIGDFT 1088 real(dp) :: edisp 1089 real(dp) :: ewaldstr(6) 1090 #endif 1091 1092 ! ************************************************************************* 1093 1094 !Store xcart for each atom 1095 ABI_MALLOC(xcart,(3, dtset%natom)) 1096 call xred2xcart(dtset%natom, rprimd, xcart, xred) 1097 1098 nullify(fdisp) 1099 nullify(grew_cart) 1100 dispersion = 0 1101 psoffset = 0._dp 1102 #if defined HAVE_BIGDFT 1103 call IonicEnergyandForces(me, nproc, wvl_den%denspot%dpbox,& 1104 & wvl%atoms, dtset%efield, xcart, & 1105 & eew, grew_cart, dispersion, edisp, fdisp,& 1106 & ewaldstr,wvl%Glr%d%n1,wvl%Glr%d%n2,wvl%Glr%d%n3,& 1107 & wvl_den%denspot%V_ext, wvl_den%denspot%pkernel,psoffset) 1108 1109 if (associated(fdisp)) then 1110 ABI_FREE(fdisp) 1111 end if 1112 #endif 1113 1114 ABI_FREE(xcart) 1115 1116 !Transform cartesian gradients to reduced gradients. 1117 do iatom = 1, dtset%natom, 1 1118 do igeo = 1, 3, 1 1119 grewtn(igeo, iatom) = -rprimd(1, igeo) * grew_cart(1, iatom) - & 1120 & rprimd(2, igeo) * grew_cart(2, iatom) - & 1121 & rprimd(3, igeo) * grew_cart(3, iatom) 1122 end do 1123 end do 1124 if (associated(grew_cart)) then 1125 ABI_FREE(grew_cart) 1126 end if 1127 1128 #if !defined HAVE_BIGDFT 1129 if (.false.) write(std_out,*) me,nproc,wvl%h(1),wvl_den%symObj 1130 #endif 1131 1132 end subroutine ionion_surface
ABINIT/m_setvtr [ Modules ]
NAME
m_setvtr
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 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 .
SOURCE
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_setvtr 22 23 use defs_basis 24 use defs_wvltypes 25 use m_abicore 26 use m_errors 27 use m_abi2big 28 use m_xmpi 29 use m_xcdata 30 use m_dtset 31 32 use defs_datatypes, only : pseudopotential_type 33 use defs_abitypes, only : MPI_type 34 use m_time, only : timab 35 use m_geometry, only : xred2xcart 36 use m_cgtools, only : dotprod_vn 37 use m_ewald, only : ewald 38 use m_energies, only : energies_type 39 use m_electronpositron, only : electronpositron_type, electronpositron_calctype, rhohxcpositron 40 use libxc_functionals, only : libxc_functionals_is_hybrid 41 use m_pawang, only : pawang_type 42 use m_pawrad, only : pawrad_type 43 use m_pawrhoij, only : pawrhoij_type 44 use m_pawtab, only : pawtab_type 45 use m_jellium, only : jellium 46 use m_spacepar, only : hartre 47 use m_dens, only : constrained_dft_t,constrained_dft_ini,constrained_dft_free,mag_penalty 48 use m_vdw_dftd2, only : vdw_dftd2 49 use m_vdw_dftd3, only : vdw_dftd3 50 use m_atm2fft, only : atm2fft 51 use m_rhotoxc, only : rhotoxc 52 use m_mklocl, only : mklocl 53 use m_xchybrid, only : xchybrid_ncpp_cc 54 use m_mkcore, only : mkcore, mkcore_alt 55 use m_psolver, only : psolver_rhohxc 56 use m_wvl_psi, only : wvl_psitohpsi 57 use m_mkcore_wvl, only : mkcore_wvl 58 use m_xc_tb09, only : xc_tb09_update_c 59 60 #if defined HAVE_BIGDFT 61 use BigDFT_API, only: denspot_set_history 62 #endif 63 64 implicit none 65 66 private
m_setvtr/setvtr [ Functions ]
[ Top ] [ m_setvtr ] [ Functions ]
NAME
setvtr
FUNCTION
Set up the trial potential and some energy terms
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 pawang <type(pawang_type)> =paw angular mesh and related data pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom) 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) 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 [taur(nfftf,nspden*dtset%usekden)]=array for kinetic energy density
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
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) [electronpositron <type(electronpositron_type)>]=quantities for the electron-positron annihilation (optional argument) [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 [xcctau3d(n3xccc*usekden)]=3D core electron kinetic energy 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.
Developers 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.
SOURCE
189 subroutine setvtr(atindx1,dtset,energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcut,& 190 & istep,kxc,mgfft,moved_atm_inside,moved_rhor,mpi_enreg,& 191 & nattyp,nfft,ngfft,ngrvdw,nhat,nhatgr,nhatgrdim,nkxc,ntypat,n1xccc,n3xccc,& 192 & optene,pawang,pawrad,pawrhoij,pawtab,ph1d,psps,rhog,rhor,rmet,rprimd,strsxc,& 193 & ucvol,usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,wvl,xccc3d,xred,& 194 & electronpositron,taur,vxc_hybcomp,vxctau,add_tfw,xcctau3d,calc_ewald) ! optionals arguments 195 196 !Arguments ------------------------------------ 197 !scalars 198 integer,intent(in) :: istep,mgfft,n1xccc,n3xccc,nfft,ngrvdw,nhatgrdim,nkxc,ntypat 199 integer,intent(in) :: optene,usexcnhat 200 integer,intent(inout) :: moved_atm_inside,moved_rhor 201 logical,intent(in),optional :: add_tfw 202 logical,intent(in),optional :: calc_ewald 203 real(dp),intent(in) :: gsqcut,ucvol 204 real(dp),intent(out) :: vxcavg 205 type(MPI_type),intent(in) :: mpi_enreg 206 type(dataset_type),intent(inout) :: dtset 207 type(electronpositron_type),pointer,optional :: electronpositron 208 type(energies_type),intent(inout) :: energies 209 type(pawang_type),intent(in) :: pawang 210 type(pseudopotential_type),intent(in) :: psps 211 type(wvl_data), intent(inout) :: wvl 212 !arrays 213 integer, intent(in) :: atindx1(dtset%natom),nattyp(ntypat),ngfft(18) 214 real(dp),intent(in) :: gmet(3,3),gprimd(3,3) 215 real(dp),intent(in) :: nhat(nfft,dtset%nspden*psps%usepaw) 216 real(dp),intent(in) :: nhatgr(:,:,:) !(nfft,dtset%nspden,3*nhatgrdim) 217 real(dp),intent(in) :: rhog(2,nfft) 218 real(dp),intent(inout) :: rmet(3,3),rprimd(3,3) 219 real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom) 220 real(dp),intent(inout) :: rhor(nfft,dtset%nspden),vhartr(nfft),vpsp(nfft) 221 real(dp),intent(inout),optional :: taur(nfft,dtset%nspden*dtset%usekden) 222 real(dp),intent(inout) :: vtrial(nfft,dtset%nspden),vxc(nfft,dtset%nspden) 223 real(dp),intent(out),optional :: vxctau(nfft,dtset%nspden,4*dtset%usekden) 224 real(dp),intent(out),optional :: vxc_hybcomp(:,:) ! (nfft,nspden) 225 real(dp),intent(inout) :: xccc3d(n3xccc) 226 real(dp),intent(inout),optional ::xcctau3d(n3xccc*dtset%usekden) 227 real(dp),intent(in) :: xred(3,dtset%natom) 228 real(dp),intent(out) :: grchempottn(3,dtset%natom) 229 real(dp),intent(out) :: grewtn(3,dtset%natom),grvdw(3,ngrvdw),kxc(nfft,nkxc),strsxc(6) 230 type(pawrhoij_type),intent(in) :: pawrhoij(:) 231 type(pawrad_type),intent(in) :: pawrad(ntypat*dtset%usepaw) 232 type(pawtab_type),intent(in) :: pawtab(ntypat*dtset%usepaw) 233 234 !Local variables------------------------------- 235 !scalars 236 integer :: coredens_method,coretau_method,mpi_comm_sphgrid,nk3xc 237 integer :: iatom,ifft,ipositron,ispden,nfftot 238 integer :: optatm,optdyfr,opteltfr,optgr,option,option_eff,optn,optn2,optstr,optv,vloc_method 239 real(dp) :: doti,e_xcdc_vxctau,ebb,ebn,evxc,ucvol_local,rpnrm 240 logical :: add_tfw_,is_hybrid_ncpp,non_magnetic_xc,with_vxctau,wvlbigdft,lewald 241 real(dp), allocatable :: xcart(:,:) 242 character(len=500) :: message 243 type(constrained_dft_t) :: constrained_dft 244 type(xcdata_type) :: xcdata,xcdatahyb 245 !arrays 246 real(dp),parameter :: identity(1:4)=(/1._dp,1._dp,0._dp,0._dp/) 247 real(dp) :: dummy6(6),tsec(2) 248 real(dp) :: grewtn_fake(3,1) 249 real(dp) :: dummy_in(0) 250 real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0) 251 real(dp) :: strn_dummy6(6), strv_dummy6(6) 252 real(dp) :: vzeeman(4) 253 real(dp),allocatable :: grtn(:,:),dyfr_dum(:,:,:),gr_dum(:,:) 254 real(dp),allocatable :: rhojellg(:,:),rhojellr(:),rhowk(:,:),vjell(:) 255 real(dp),allocatable :: v_constr_dft_r(:,:),rhog_dum(:,:) 256 257 ! ********************************************************************* 258 259 call timab(91,1,tsec) 260 261 !Check that usekden is not 0 if want to use vxctau 262 with_vxctau = (present(vxctau).and.present(taur).and.(dtset%usekden/=0)) 263 264 !Check if we're in hybrid norm conserving pseudopotential with a core correction 265 is_hybrid_ncpp=(dtset%usepaw==0 .and. n3xccc/=0 .and. & 266 & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid())) 267 268 !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed 269 wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1) 270 271 !Get size of FFT grid 272 nfftot=PRODUCT(ngfft(1:3)) 273 274 !mpi 275 mpi_comm_sphgrid=mpi_enreg%comm_fft 276 if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl 277 278 !Test electron-positron case 279 ipositron=0;if (present(electronpositron)) ipositron=electronpositron_calctype(electronpositron) 280 281 !Test addition of Weiszacker gradient correction to Thomas-Fermi kin energy 282 add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw 283 284 !Get Ewald energy and Ewald forces, as well as vdW-DFTD energy and forces, and chemical potential energy and forces. 285 !------------------------------------------------------------------------------------------------------------------- 286 call timab(5,1,tsec) 287 if (ipositron/=1) then 288 lewald = .true. 289 if (present(calc_ewald)) lewald=calc_ewald 290 if (lewald) then 291 if (dtset%icoulomb == 0 .or. (dtset%usewvl == 0 .and. dtset%icoulomb == 2)) then 292 ! Periodic system, need to compute energy and forces due to replica and 293 ! to correct the shift in potential calculation. 294 call ewald(energies%e_ewald,gmet,grewtn,gsqcut,dtset%icutcoul,dtset%natom,ngfft,dtset%nkpt,ntypat,& 295 &dtset%rcut,rmet,rprimd,dtset%typat,ucvol,dtset%vcutgeo,xred,psps%ziontypat) 296 ! For a periodic system bearing a finite charge, the monopole correction to the 297 ! energy is relevant. 298 ! See Leslie and Gillan, JOURNAL OF PHYSICS C-SOLID STATE PHYSICS 18, 973 (1985) 299 if(abs(dtset%cellcharge(1))>tol8) then 300 call ewald(energies%e_monopole,gmet,grewtn_fake,gsqcut,dtset%icutcoul,1,ngfft,dtset%nkpt,1,& 301 &dtset%rcut,rmet,rprimd,(/1/),ucvol,dtset%vcutgeo,(/0.0_dp,0.0_dp,0.0_dp/),(/dtset%cellcharge(1)/)) 302 energies%e_monopole=-energies%e_monopole 303 end if 304 else if (dtset%icoulomb == 1) then 305 ! In a non periodic system (real space computation), the G=0 divergence 306 ! doesn't occur and ewald is not needed. Only the ion/ion interaction 307 ! energy is relevant and used as Ewald energy and gradient. 308 call ionion_realSpace(dtset, energies%e_ewald, grewtn, rprimd, xred, psps%ziontypat) 309 else if (dtset%icoulomb == 2) then 310 call ionion_surface(dtset, energies%e_ewald, grewtn, mpi_enreg%me_wvl, mpi_enreg%nproc_wvl, rprimd, & 311 & wvl%descr, wvl%den, xred) 312 end if 313 end if 314 if (dtset%nzchempot>0) then 315 call spatialchempot(energies%e_chempot,dtset%chempot,grchempottn,dtset%natom,ntypat,dtset%nzchempot,dtset%typat,xred) 316 end if 317 if (dtset%vdw_xc==5.and.ngrvdw==dtset%natom) then 318 call vdw_dftd2(energies%e_vdw_dftd,dtset%ixc,dtset%natom,ntypat,1,dtset%typat,rprimd,& 319 & dtset%vdw_tol,xred,psps%znucltypat,gred_vdw_dftd2=grvdw) 320 end if 321 if ((dtset%vdw_xc==6.or.dtset%vdw_xc==7).and.ngrvdw==dtset%natom) then 322 call vdw_dftd3(energies%e_vdw_dftd,dtset%ixc,dtset%natom,& 323 & ntypat,1,dtset%typat,rprimd,dtset%vdw_xc,dtset%vdw_tol,dtset%vdw_tol_3bt,& 324 & xred,psps%znucltypat,gred_vdw_dftd3=grvdw) 325 end if 326 else 327 energies%e_ewald=zero 328 energies%e_chempot=zero 329 grchempottn=zero 330 grewtn=zero 331 energies%e_vdw_dftd=zero 332 if (ngrvdw>0) grvdw=zero 333 end if 334 call timab(5,2,tsec) 335 336 !Compute parts of total energy depending on potentials 337 !-------------------------------------------------------------- 338 if (dtset%usewvl == 0) then 339 ucvol_local = ucvol 340 #if defined HAVE_BIGDFT 341 else 342 ! We need to tune the volume when wavelets are used because, not all FFT points are used. 343 ! ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3) 344 ucvol_local = product(wvl%den%denspot%dpbox%hgrids) * real(product(wvl%den%denspot%dpbox%ndims), dp) 345 #endif 346 end if 347 348 !Determine by which method the local ionic potential and/or the pseudo core charge density 349 ! have to be computed 350 !Local ionic potential: 351 ! Method 1: PAW 352 ! Method 2: Norm-conserving PP, icoulomb>0, wavelets 353 vloc_method=1;if (psps%usepaw==0) vloc_method=2 354 if (dtset%icoulomb>0) vloc_method=2 355 if (psps%usewvl==1) vloc_method=2 356 !Pseudo core charge density: 357 ! Method 1: PAW, nc_xccc_gspace 358 ! Method 2: Norm-conserving PP, wavelets 359 coredens_method=1;if (psps%usepaw==0) coredens_method=2 360 if (psps%nc_xccc_gspace==1) coredens_method=1 361 if (psps%nc_xccc_gspace==0) coredens_method=2 362 if (psps%usewvl==1) coredens_method=2 363 coretau_method=0 364 if (dtset%usekden==1.and.psps%usepaw==1) then 365 coretau_method=1;if (psps%nc_xccc_gspace==0) coretau_method=2 366 end if 367 !In some specific cases, XC has to be handled as non-magnetic 368 non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 369 370 !Local ionic potential and/or pseudo core charge by method 1 371 if (vloc_method==1.or.coredens_method==1) then 372 call timab(552,1,tsec) 373 optv=0;if (vloc_method==1) optv=1 374 optn=0;if (coredens_method==1) optn=n3xccc/nfft 375 optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=1 376 call atm2fft(atindx1,xccc3d,vpsp,dummy_out1,dummy_out2,dummy_out3,dummy_in,& 377 & gmet,gprimd,dummy_out4,dummy_out5,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,& 378 & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,& 379 & dtset%rcut,dummy_in,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,& 380 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 381 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 382 call timab(552,2,tsec) 383 end if 384 if (coretau_method==1) then 385 call timab(552,1,tsec) 386 optv=0;optn=1 387 optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=4 388 call atm2fft(atindx1,xcctau3d,dummy_out6,dummy_out1,dummy_out2,dummy_out3,dummy_in,& 389 & gmet,gprimd,dummy_out4,dummy_out5,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,& 390 & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,& 391 & dtset%rcut,dummy_in,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,& 392 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 393 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 394 call timab(552,2,tsec) 395 end if 396 397 !Local ionic potential by method 2 398 if (vloc_method==2) then 399 option=1 400 ABI_MALLOC(gr_dum,(3,dtset%natom)) 401 ABI_MALLOC(dyfr_dum,(3,3,dtset%natom)) 402 ABI_MALLOC(rhog_dum,(2,nfft)) 403 call mklocl(dtset,dyfr_dum,energies%e_localpsp,gmet,gprimd,& 404 & gr_dum,gsqcut,dummy6,mgfft,mpi_enreg,dtset%natom,nattyp,& 405 & nfft,ngfft,dtset%nspden,ntypat,option,pawtab,ph1d,psps,& 406 & dtset%qprtrb,rhog_dum,rhor,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred) 407 ABI_FREE(gr_dum) 408 ABI_FREE(dyfr_dum) 409 ABI_FREE(rhog_dum) 410 end if 411 412 !3D pseudo core electron density xccc3d by method 2 413 if (coredens_method==2.and.n1xccc/=0) then 414 call timab(91,2,tsec) 415 call timab(92,1,tsec) 416 option=1 417 ABI_MALLOC(gr_dum,(3,dtset%natom)) 418 ABI_MALLOC(dyfr_dum,(3,3,dtset%natom)) 419 if (psps%usewvl==0.and.psps%usepaw==0.and.dtset%icoulomb==0) then 420 call mkcore(dummy6,dyfr_dum,gr_dum,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,& 421 & ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,dtset%typat,ucvol,& 422 & vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred) 423 else if (psps%usewvl==0.and.(psps%usepaw==1.or.dtset%icoulomb==1)) then 424 call mkcore_alt(atindx1,dummy6,dyfr_dum,gr_dum,dtset%icoulomb,mpi_enreg,dtset%natom,& 425 & nfft,dtset%nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,& 426 & ucvol,vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred,pawrad,pawtab,psps%usepaw) 427 else if (psps%usewvl==1.and.psps%usepaw==1) then 428 #if defined HAVE_BIGDFT 429 ! call mkcore_wvl_old(atindx1,dummy6,dyfr_dum,wvl%descr%atoms%astruct%geocode,gr_dum,wvl%descr%h,& 430 ! & dtset%natom,nattyp,nfft,wvl%den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl,:),& 431 ! & dtset%nspden,ntypat,wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,wvl%descr%Glr%d%n2,& 432 ! & wvl%descr%Glr%d%n2i,wvl%descr%Glr%d%n3,wvl%den%denspot%dpbox%n3pi,n3xccc,option,& 433 ! & pawrad,pawtab,psps%gth_params%psppar,rprimd,ucvol_local,vxc,xccc3d,xred,& 434 ! & mpi_comm_wvl=mpi_enreg%comm_wvl) 435 call mkcore_wvl(atindx1,dummy6,gr_dum,dtset%natom,nattyp,nfft,dtset%nspden,ntypat,& 436 & n1xccc,n3xccc,option,pawrad,pawtab,rprimd,vxc,psps%xccc1d,xccc3d,& 437 & psps%xcccrc,xred,wvl%den,wvl%descr,mpi_comm_wvl=mpi_enreg%comm_wvl) 438 #endif 439 end if 440 ABI_FREE(gr_dum) 441 ABI_FREE(dyfr_dum) 442 call timab(92,2,tsec) 443 call timab(91,1,tsec) 444 end if 445 if (coretau_method==2) then 446 call timab(91,2,tsec) 447 call timab(92,1,tsec) 448 option=1 449 ABI_MALLOC(gr_dum,(3,dtset%natom)) 450 ABI_MALLOC(dyfr_dum,(3,3,dtset%natom)) 451 call mkcore_alt(atindx1,dummy6,dyfr_dum,gr_dum,dtset%icoulomb,mpi_enreg,dtset%natom,& 452 & nfft,dtset%nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,& 453 & ucvol,vxc,psps%xcccrc,psps%xccc1d,xcctau3d,xred,pawrad,pawtab,psps%usepaw,& 454 & usekden=.true.) 455 ABI_FREE(gr_dum) 456 ABI_FREE(dyfr_dum) 457 call timab(92,2,tsec) 458 call timab(91,1,tsec) 459 end if 460 461 !Adds the jellium potential to the local part of ionic potential 462 if (dtset%jellslab/=0) then 463 ABI_MALLOC(vjell,(nfft)) 464 ABI_MALLOC(rhojellg,(2,nfft)) 465 ABI_MALLOC(rhojellr,(nfft)) 466 option=1 467 call jellium(gmet,gsqcut,mpi_enreg,nfft,ngfft,dtset%nspden,option,& 468 & dtset%slabwsrad,rhojellg,rhojellr,rprimd,vjell,dtset%slabzbeg,dtset%slabzend) 469 ! Compute background-background energy 470 call dotprod_vn(1,rhojellr,ebb,doti,nfft,nfftot,1,1,vjell,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid) 471 ebb=half*ebb 472 ! Compute electrostatic energy between background and nuclei before adding vjell to vpsp 473 call dotprod_vn(1,rhojellr,ebn,doti,nfft,nfftot,1,1,vpsp,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid) 474 ! Update e_ewald with ebb and ebn 475 energies%e_ewald=energies%e_ewald+ebb+ebn 476 ! Compute gradient of ebn wrt tn 477 ! This is not yet coded for usewvl or icoulomb=1 478 if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then 479 write(message,'(3a)')& 480 & 'The computation of forces due to jellium background',ch10,& 481 & 'has to be verified in the PAW formalism.' 482 ABI_WARNING(message) 483 484 ABI_MALLOC(grtn,(3,dtset%natom)) 485 optatm=0;optdyfr=0;opteltfr=0;optgr=1;optstr=0;optv=1;optn=0;optn2=1 486 call atm2fft(atindx1,dummy_out1,vpsp,dummy_out2,dummy_out3,dummy_out4,dummy_in,& 487 & gmet,gprimd,dummy_out5,grtn,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,& 488 & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,& 489 & dtset%rcut,rhojellg,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,& 490 & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,& 491 & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft) 492 493 ! Update grewtn with gradient of ebn wrt tn 494 do iatom=1,dtset%natom 495 grewtn(1:3,iatom)=grewtn(1:3,iatom)+grtn(1:3,iatom) 496 end do 497 ABI_FREE(grtn) 498 else ! of usepaw==1 499 option=2 500 ABI_MALLOC(dyfr_dum,(3,3,dtset%natom)) 501 ABI_MALLOC(grtn,(3,dtset%natom)) 502 call mklocl(dtset,dyfr_dum,energies%e_localpsp,gmet,gprimd,& 503 & grtn,gsqcut,dummy6,mgfft,mpi_enreg,dtset%natom,nattyp,& 504 & nfft,ngfft,1,ntypat,option,pawtab,ph1d,psps,dtset%qprtrb,rhojellg,& 505 & rhojellr,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred) 506 ! Update grewtn with gradient of ebn wrt tn (reestablish order of atoms) 507 do iatom=1,dtset%natom 508 grewtn(1:3,atindx1(iatom))=grewtn(1:3,atindx1(iatom))+grtn(1:3,iatom) 509 end do 510 ABI_FREE(dyfr_dum) 511 ABI_FREE(grtn) 512 end if ! of usepaw==1 513 vpsp(:)=vpsp(:)+vjell(:) 514 ABI_FREE(vjell) 515 ABI_FREE(rhojellg) 516 ABI_FREE(rhojellr) 517 end if 518 519 !Additional stuff for electron-positron calculation 520 !Compute the electronic/positronic local (Hartree) potential 521 if (ipositron==1) vpsp=-vpsp 522 523 !If we are at the initialisation, or 524 !if the atom positions has changed and the non-linear core correction 525 !is included, or the rhor has changed, one needs to compute the xc stuff. 526 !One needs also to compute the Hartree stuff if the density changed, 527 !or at initialisation. 528 !-------------------------------------------------------------- 529 530 if(istep==1 .or. n1xccc/=0 .or. moved_rhor==1 .or. dtset%positron<0 .or. mod(dtset%fockoptmix,100)==11) then 531 532 option=0 533 if(istep==1 .or. moved_rhor==1 .or. dtset%positron<0 .or. mod(dtset%fockoptmix,100)==11) option=1 534 if (nkxc>0) option=2 535 if (dtset%iscf==-1) option=-2 536 537 if (ipositron/=1) then 538 if (dtset%icoulomb == 0 .and. dtset%usewvl == 0) then 539 540 ! >>>> Hartree potential 541 if(option/=0 .and. option/=10)then 542 call hartre(1,gsqcut,dtset%icutcoul,psps%usepaw,mpi_enreg,nfft,ngfft,& 543 &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr) 544 end if 545 546 ! >>>> Exchange-correlation potential 547 call xcdata_init(xcdata,dtset=dtset) 548 if(mod(dtset%fockoptmix,100)==11)then 549 xcdatahyb=xcdata 550 ! Setup the auxiliary xc functional information 551 call xcdata_init(xcdata,dtset=dtset,auxc_ixc=0,ixc=dtset%auxc_ixc) 552 end if 553 ! Not yet able to deal fully with the full XC kernel in case of GGA + spin 554 option_eff=option;if (option==2.and.xcdata%xclevel==2.and.(nkxc==3-2*mod(xcdata%nspden,2))) option_eff=12 555 nk3xc=1 556 557 ! If we use the XC Tran-Blaha 2009 (modified BJ) functional, update the c value 558 if (dtset%xc_tb09_c>99._dp) then 559 call xc_tb09_update_c(dtset%intxc,dtset%ixc,mpi_enreg,dtset%natom, & 560 & nfft,ngfft,nhat,psps%usepaw,nhatgr,nhatgrdim,dtset%nspden,dtset%ntypat,n3xccc, & 561 & pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,rhor,rprimd,psps%usepaw, & 562 & xccc3d,dtset%xc_denpos,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, & 563 & computation_type='all') 564 end if 565 566 if (ipositron==0) then 567 568 ! Compute energies%e_xc and associated quantities 569 if(.not.is_hybrid_ncpp .or. mod(dtset%fockoptmix,100)==11)then 570 call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,& 571 & nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,& 572 & option_eff,rhor,rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,& 573 & taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_,xcctau3d=xcctau3d,grho1_over_rho1=rpnrm) 574 else 575 ! Only when is_hybrid_ncpp, and moreover, the xc functional is not the auxiliary xc functional, then call xchybrid_ncpp_cc 576 call xchybrid_ncpp_cc(dtset,energies%e_xc,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,& 577 & strsxc,vxcavg,xccc3d,vxc=vxc) 578 end if 579 580 ! Possibly compute energies%e_hybcomp_E0 581 if(mod(dtset%fockoptmix,100)==11)then 582 ! This call to rhotoxc uses the hybrid xc functional 583 if(.not.is_hybrid_ncpp)then 584 call rhotoxc(energies%e_hybcomp_E0,kxc,mpi_enreg,nfft,ngfft,& 585 & nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,& 586 & option_eff,rhor,rprimd,strsxc,usexcnhat,vxc_hybcomp,vxcavg,xccc3d,xcdatahyb,& 587 & taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_) 588 else 589 call xchybrid_ncpp_cc(dtset,energies%e_hybcomp_E0,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,& 590 & strsxc,vxcavg,xccc3d,vxc=vxc_hybcomp) 591 end if 592 593 ! Combine hybrid and auxiliary quantities 594 energies%e_xc=energies%e_xc*dtset%auxc_scal 595 energies%e_hybcomp_E0=energies%e_hybcomp_E0-energies%e_xc 596 vxc(:,:)=vxc(:,:)*dtset%auxc_scal 597 vxc_hybcomp(:,:)=vxc_hybcomp(:,:)-vxc(:,:) 598 599 end if 600 601 else if (ipositron==2) then 602 call rhotoxc(energies%e_xc,kxc,mpi_enreg,nfft,ngfft,& 603 & nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,& 604 & option_eff,rhor,rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,& 605 & taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_,& 606 & electronpositron=electronpositron) 607 end if 608 609 elseif(.not. wvlbigdft) then 610 ! Use the free boundary solver 611 call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, & 612 & dtset%icoulomb, dtset%ixc, & 613 & mpi_enreg, nfft, ngfft,& 614 & nhat,psps%usepaw,& 615 & dtset%nscforder,dtset%nspden,n3xccc,rhor,rprimd, & 616 & usexcnhat,psps%usepaw,dtset%usewvl,vhartr, vxc, & 617 & vxcavg,wvl%descr,wvl%den,wvl%e,& 618 & xccc3d,dtset%xclevel,dtset%xc_denpos) 619 end if 620 else 621 energies%e_xc=zero 622 call rhohxcpositron(electronpositron,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat,nkxc,dtset%nspden,n3xccc,& 623 & dtset%paral_kgb,rhor,strsxc,ucvol,usexcnhat,psps%usepaw,vhartr,vxc,vxcavg,xccc3d,dtset%xc_denpos) 624 end if 625 if (ipositron/=0) then 626 if (optene>=1) then 627 call dotprod_vn(1,rhor,electronpositron%e_hartree,doti,& 628 & nfft,nfftot,1,1,electronpositron%vha_ep,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid) 629 end if 630 vhartr=vhartr+electronpositron%vha_ep 631 end if 632 end if 633 634 !Compute the trial potential 635 !------------------------------------------------------------- 636 if (.not. wvlbigdft) then 637 ! Now, compute trial Hxc potential. Local psp potential will be added later. 638 if(moved_atm_inside==0 .or.dtset%iscf>=10) then 639 640 ! Compute starting Hxc potential. 641 ! Multiply by identity, should not change anything if nspden /= 4 642 do ispden=1,dtset%nspden 643 vtrial(:,ispden)=vhartr(:)*identity(ispden)+vxc(:,ispden) 644 end do 645 646 else 647 648 ! One should be here only when moved_atm_inside==1 649 ! The (H)xc now added corrects the previous one. 650 if(dtset%densfor_pred==1)then 651 ! xc was substracted off. This should be rationalized later 652 do ispden=1,dtset%nspden 653 vtrial(:,ispden)=vtrial(:,ispden)+vxc(:,ispden) 654 end do 655 else if(abs(dtset%densfor_pred)==2.or.abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)then 656 ! Hxc was substracted off. This should be rationalized later 657 do ispden=1,dtset%nspden 658 vtrial(:,ispden)=vtrial(:,ispden)+vhartr(:)*identity(ispden)+vxc(:,ispden) 659 end do 660 end if 661 end if 662 663 ! Adds the local part of the potential 664 if ((moved_atm_inside==0).or.(dtset%densfor_pred/=3)) then 665 do ispden=1,min(2,dtset%nspden) 666 do ifft=1,nfft 667 vtrial(ifft,ispden)=vtrial(ifft,ispden)+vpsp(ifft) 668 end do 669 end do 670 end if 671 672 ! Adds the compensating vxc for hybrids 673 if(mod(dtset%fockoptmix,100)==11)then 674 vtrial(:,:)=vtrial(:,:)+vxc_hybcomp(:,:) 675 end if 676 677 if(dtset%usewvl==1) then 678 call wvl_vtrial_abi2big(1,vtrial,wvl%den) 679 end if 680 681 else 682 683 ! Compute with covering comms the different part of the potential. 684 #if defined HAVE_BIGDFT 685 ! Copy e_ewald. 686 wvl%e%energs%eion = energies%e_ewald 687 ! Setup the mixing, if necessary 688 call denspot_set_history(wvl%den%denspot,dtset%iscf,dtset%nsppol, & 689 & wvl%den%denspot%dpbox%ndims(1),wvl%den%denspot%dpbox%ndims(2)) 690 #endif 691 ABI_MALLOC(xcart,(3, dtset%natom)) 692 call xred2xcart(dtset%natom, rprimd, xcart, xred) 693 call wvl_psitohpsi(dtset%diemix,energies%e_exactX, energies%e_xc, energies%e_hartree, & 694 & energies%e_kinetic, energies%e_localpsp, energies%e_nlpsp_vfock, energies%e_sicdc, & 695 & istep, 1, dtset%iscf, mpi_enreg%me_wvl, dtset%natom, dtset%nfft, mpi_enreg%nproc_wvl, dtset%nspden, & 696 & rpnrm, .true.,evxc, wvl,.true., xcart, strsxc,& 697 & vtrial, vxc) 698 if (optene==3.or.optene==4) energies%e_xcdc=evxc 699 ABI_FREE(xcart) 700 701 end if 702 703 !Add the zeeman field to vtrial 704 if (any(abs(dtset%zeemanfield(:))>tol8)) then 705 vzeeman(:) = zero ! vzeeman_ij = -1/2*sigma_ij^alpha*B_alpha 706 if(dtset%nspden==2)then 707 vzeeman(1) = -half*dtset%zeemanfield(3) ! v_dwndwn = -1/2*B_z 708 vzeeman(2) = half*dtset%zeemanfield(3) ! v_upup = 1/2*B_z 709 do ifft=1,nfft 710 vtrial(ifft,1) = vtrial(ifft,1) + vzeeman(1) !SPr: added 1st component 711 vtrial(ifft,2) = vtrial(ifft,2) + vzeeman(2) 712 end do !ifft 713 end if 714 if(dtset%nspden==4)then 715 vzeeman(1)=-half*dtset%zeemanfield(3) ! v_dwndwn => v_11 716 vzeeman(2)= half*dtset%zeemanfield(3) ! v_upup => v_22 717 vzeeman(3)=-half*dtset%zeemanfield(1) ! Re(v_dwnup) = Re(v_updwn) => Re(v_12) 718 vzeeman(4)= half*dtset%zeemanfield(2) ! Im(v_dwnup) =-Im(v_dwnup) => Im(v_12) 719 do ispden=1,dtset%nspden 720 do ifft=1,nfft 721 vtrial(ifft,ispden) = vtrial(ifft,ispden) + vzeeman(ispden) 722 end do 723 end do 724 end if 725 end if 726 727 !Compute the constrained potential for the magnetic moments 728 if (dtset%magconon==1.or.dtset%magconon==2) then 729 ! Initialize the datastructure constrained_dft, for penalty function constrained magnetization 730 call constrained_dft_ini(dtset%chrgat,constrained_dft,dtset%constraint_kind,dtset%magconon,dtset%magcon_lambda,& 731 & mpi_enreg,dtset%natom,dtset%nfft,dtset%ngfft,dtset%nspden,dtset%ntypat,dtset%ratsm,& 732 & dtset%ratsph,rprimd,dtset%spinat,dtset%typat,xred,dtset%ziontypat) 733 ABI_MALLOC(v_constr_dft_r, (nfft,dtset%nspden)) 734 v_constr_dft_r = zero 735 call mag_penalty(constrained_dft,mpi_enreg,rhor,v_constr_dft_r,xred) 736 if(dtset%nspden==4)then 737 do ispden=1,dtset%nspden ! (SPr: both components should be used? EB: Yes it should be the case, corrected now) 738 do ifft=1,nfft 739 vtrial(ifft,ispden) = vtrial(ifft,ispden) + v_constr_dft_r(ifft,ispden) 740 end do !ifft 741 end do !ispden 742 else if(dtset%nspden==2)then 743 do ifft=1,nfft 744 ! TODO : MJV: check that magnetic constraint works also for nspden 2 or add input variable condition 745 ! EB: ispden=2 is rho_up only: to be tested 746 ! SPr: for ispden=2, both components should be used (e.g. see definition for vzeeman)? 747 vtrial(ifft,1) = vtrial(ifft,1) + v_constr_dft_r(ifft,1) !SPr: added the first component here 748 vtrial(ifft,2) = vtrial(ifft,2) + v_constr_dft_r(ifft,2) 749 end do !ifft 750 end if 751 ABI_FREE(v_constr_dft_r) 752 call constrained_dft_free(constrained_dft) 753 end if 754 755 !Compute parts of total energy depending on potentials 756 !-------------------------------------------------------------- 757 758 !For icoulomb==0 or usewvl Ehartree is calculated in psolver_rhohxc(). 759 !For PAW we recalculate this since nhat was not taken into account 760 !in psolver_rhohxc: E_H= int v_H (n+nhat) dr 761 762 if (optene>=1 .and. .not. wvlbigdft .and. (dtset%icoulomb==0 .or. dtset%usepaw==1 ) ) then 763 ! Compute Hartree energy ehart 764 ! Already available in the Psolver case through psolver_rhohxc(). 765 if (ipositron/=1) then 766 call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol_local,& 767 & mpi_comm_sphgrid=mpi_comm_sphgrid) 768 if (ipositron==0) energies%e_hartree = half * energies%e_hartree 769 if (ipositron==2) energies%e_hartree = half * (energies%e_hartree-electronpositron%e_hartree) 770 else 771 energies%e_hartree=zero 772 end if 773 end if 774 775 if(mod(dtset%fockoptmix,100)==11)then 776 if (.not. wvlbigdft) then 777 call dotprod_vn(1,rhor,energies%e_hybcomp_v0,doti,nfft,nfftot,1,1,vxc_hybcomp,ucvol_local,& 778 & mpi_comm_sphgrid=mpi_comm_sphgrid) 779 energies%e_hybcomp_v=energies%e_hybcomp_v0 780 end if 781 end if 782 783 if (optene==2.or.optene==4 .and. .not. wvlbigdft) then 784 ! Compute local psp energy eei 785 call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfft,nfftot,1,1,vpsp,ucvol_local,& 786 & mpi_comm_sphgrid=mpi_comm_sphgrid) 787 end if 788 789 if (optene==3.or.optene==4 .and. .not. wvlbigdft) then 790 ! Compute double-counting XC energy enxcdc 791 if (ipositron/=1) then 792 if (dtset%usepaw==0.or.usexcnhat/=0) then 793 call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,& 794 & mpi_comm_sphgrid=mpi_comm_sphgrid) 795 else 796 ABI_MALLOC(rhowk,(nfft,dtset%nspden)) 797 rhowk=rhor-nhat 798 call dotprod_vn(1,rhowk,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,& 799 & mpi_comm_sphgrid=mpi_comm_sphgrid) 800 ABI_FREE(rhowk) 801 end if 802 if (with_vxctau) then 803 call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfft,nfftot,dtset%nspden,1,vxctau(:,:,1),& 804 & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid) 805 energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau 806 end if 807 if (ipositron==2) energies%e_xcdc=energies%e_xcdc-electronpositron%e_xcdc 808 else 809 energies%e_xcdc=zero 810 end if 811 end if 812 813 !-------------------------------------------------------------- 814 815 !The initialisation for the new atomic positions has been done 816 moved_atm_inside=0 817 818 call timab(91,2,tsec) 819 820 end subroutine setvtr
m_setvtr/spatialchempot [ Modules ]
[ Top ] [ m_setvtr ] [ Modules ]
NAME
spatialchempot
FUNCTION
Treat spatially varying chemical potential. Compute energy and derivative with respect to dimensionless reduced atom coordinates of the spatially varying chemical potential. No contribution to stresses.
INPUTS
chempot(3,nzchempot,ntypat)=input array with information about the chemical potential (see input variable description) natom=number of atoms in unit cell ntypat=number of type of atoms nzchempot=number of limiting planes for chemical potential typat(natom)=integer label of each type of atom (1,2,...) xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
OUTPUT
e_chempot=chemical potential energy in hartrees grchempottn(3,natom)=grads of e_chempot wrt xred(3,natom), hartrees.
SOURCE
845 subroutine spatialchempot(e_chempot,chempot,grchempottn,natom,ntypat,nzchempot,typat,xred) 846 847 !Arguments ------------------------------------ 848 !scalars 849 integer,intent(in) :: natom,ntypat,nzchempot 850 real(dp),intent(out) :: e_chempot 851 !arrays 852 integer,intent(in) :: typat(natom) 853 real(dp),intent(in) :: chempot(3,nzchempot,ntypat),xred(3,natom) 854 real(dp),intent(out) :: grchempottn(3,natom) 855 856 !Local variables------------------------------- 857 !scalars 858 integer :: iatom,itypat,iz 859 real(dp) :: a_2,a_3,cp0,cp1,dcp0,dcp1,ddz,deltaz,deltaziz 860 real(dp) :: dqz,dz1,qz,zred,z0 861 !character(len=500) :: message 862 863 ! ************************************************************************* 864 865 !DEBUG 866 !write(std_out,'(a)')' spatialchempot : enter ' 867 !write(std_out,'(a,2i6)')' nzchempot,ntypat',nzchempot,ntypat 868 !write(std_out,'(a,6es13.3)') ' chempot(1:3,1:2,1)=',chempot(1:3,1:2,1) 869 !write(std_out,'(a,6es13.3)') ' chempot(1:3,1:2,2)=',chempot(1:3,1:2,2) 870 !ENDDEBUG 871 872 e_chempot=zero 873 grchempottn(:,:)=zero 874 875 !Loop on the different atoms 876 do iatom=1,natom 877 878 itypat=typat(iatom) 879 zred=xred(3,iatom) 880 881 ! Determine the delimiting plane just lower to zred 882 ! First compute the relative zred with respect to the first delimiting plane 883 ! Take into account a tolerance : 884 deltaz=zred-chempot(1,1,itypat) 885 deltaz=modulo(deltaz+tol12,1.0d0)-tol12 886 ! deltaz is positive (or higher than -tol12), and lower than one-tol12. 887 do iz=2,nzchempot+1 888 if(iz/=nzchempot+1)then 889 deltaziz=chempot(1,iz,itypat)-chempot(1,1,itypat) 890 else 891 deltaziz=one 892 end if 893 if(deltaziz>deltaz)exit 894 end do 895 896 ! Defines coordinates and values inside the delimiting interval, 897 ! with respect to the lower delimiting plane 898 z0=chempot(1,iz-1,itypat)-chempot(1,1,itypat) ; cp0=chempot(2,iz-1,itypat) ; dcp0=chempot(3,iz-1,itypat) 899 if(iz/=nzchempot+1)then 900 dz1=chempot(1,iz,itypat)-chempot(1,iz-1,itypat) ; cp1=chempot(2,iz,itypat) ; dcp1=chempot(3,iz,itypat) 901 else 902 dz1=(chempot(1,1,itypat)+one)-chempot(1,nzchempot,itypat) ; cp1=chempot(2,1,itypat) ; dcp1=chempot(3,1,itypat) 903 end if 904 ddz=deltaz-z0 905 906 !DEBUG 907 ! write(std_out,'(a,2i5)')' Delimiting planes, iz-1 and iz=',iz-1,iz 908 ! write(std_out,'(a,2es13.3)')' z0, dz1= :',z0,dz1 909 ! write(std_out,'(a,2es13.3)')' cp0, cp1= :',cp0,cp1 910 ! write(std_out,'(a,2es13.3)')' dcp0, dcp1= :',dcp0,dcp1 911 ! write(std_out,'(a,2es13.3)')' deltaz,ddz=',deltaz,ddz 912 !ENDDEBUG 913 914 ! Determine the coefficient of the third-order polynomial taking z0 as origin 915 ! P(dz=z-z0)= a_3*dz**3 + a_2*dz**2 + a_1*dz + a_0 ; obviously a_0=cp0 and a_1=dcp0 916 ! Define qz=a_3*dz + a_2 and dqz=3*a_3*dz + 2*a_2 917 qz=((cp1-cp0)-dcp0*dz1)/dz1**2 918 dqz=(dcp1-dcp0)/dz1 919 a_3=(dqz-two*qz)/dz1 920 a_2=three*qz-dqz 921 922 ! Compute value and gradient of the chemical potential, at ddz wrt to lower delimiting plane 923 e_chempot=e_chempot+(a_3*ddz**3 + a_2*ddz**2 + dcp0*ddz + cp0) 924 grchempottn(3,iatom)=three*a_3*ddz**2 + two*a_2*ddz + dcp0 925 926 !DEBUG 927 ! write(std_out,'(a,4es16.6)')' qz,dqz=',qz,dqz 928 ! write(std_out,'(a,4es16.6)')' cp0,dcp0,a_2,a_3=',cp0,dcp0,a_2,a_3 929 ! write(std_out,'(a,2es13.3)')' dcp0*ddz + cp0=',dcp0*ddz + cp0 930 ! write(std_out,'(a,2es13.3)')' a_2*ddz**2=',a_2*ddz**2 931 ! write(std_out,'(a,2es13.3)')' a_3*ddz**3=',a_3*ddz**3 932 ! write(std_out,'(a,2es13.3)')' contrib=',a_3*ddz**3 + a_2*ddz**2 + dcp0*ddz + cp0 933 ! write(std_out,'(a,2es13.3)')' e_chempot=',e_chempot 934 ! write(std_out,'(a,3es20.10)')' grchempottn=',grchempottn(:,iatom) 935 !ENDDEBUG 936 937 end do 938 939 !DEBUG 940 !write(std_out,'(a)')' spatialchempot : exit ' 941 !write(std_out,'(a,es16.6)') ' e_chempot=',e_chempot 942 !ENDDEBUG 943 944 end subroutine spatialchempot