Line data Source code
1 : !!****m* ABINIT/m_setvtr
2 : !! NAME
3 : !! m_setvtr
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1998-2026 ABINIT group (XG, GMR, FJ, MT, EB, SPr)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
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
67 : !!***
68 :
69 : public :: setvtr
70 : !!***
71 :
72 : contains
73 : !!***
74 :
75 : !!****f* m_setvtr/setvtr
76 : !! NAME
77 : !! setvtr
78 : !!
79 : !! FUNCTION
80 : !! Set up the trial potential and some energy terms
81 : !!
82 : !! INPUTS
83 : !! [add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy
84 : !! atindx1(dtset%natom)=index table for atoms, inverse of atindx
85 : !! dtset <type(dataset_type)>=all input variables in this dataset
86 : !! | if =0,1 no xc kernel, =2 spin-averaged (LDA) kernel
87 : !! | densfor_pred=govern the choice of preconditioner for the SCF cycle
88 : !! | iscf=determines the way the SCF cycle is handled
89 : !! | natom=number of atoms in cell.
90 : !! | nspden=number of spin-density components
91 : !! | qprtrb(3)= integer wavevector of possible perturbing potential
92 : !! | in basis of reciprocal lattice translations
93 : !! | typat(natom)=type integer for each atom in cell
94 : !! | vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero,
95 : !! | perturbing potential is added of the form
96 : !! | V(G)=(vprtrb(1)+I*vprtrb(2))/2 at the values G=qprtrb and
97 : !! | (vprtrb(1)-I*vprtrb(2))/2 at G=-qprtrb (integers)
98 : !! | for each type of atom, from psp (used in norm-conserving only)
99 : !! gmet(3,3)=metric tensor for G vecs (in bohr**-2)
100 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
101 : !! gsqcut=cutoff on (k+G)^2 (bohr^-2) (sphere for density and potential)
102 : !! istep=step number in the main loop of scfcv
103 : !! mgfft=maximum size of 1D FFTs
104 : !! moved_rhor=1 if the density was moved just before
105 : !! mpi_enreg=information about MPI parallelization
106 : !! nattyp(ntypat)=number of atoms of each type in cell.
107 : !! nfft=(effective) number of FFT grid points (for this processor)
108 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
109 : !! nhat(nfft,nspden*usepaw)= -PAW only- compensation density
110 : !! nhatgr(nfft,nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density
111 : !! nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise
112 : !! nkxc=second dimension of the array kxc
113 : !! ntypat=number of types of atoms in unit cell.
114 : !! n1xccc=dimension of xccc1d; 0 if no XC core correction is used
115 : !! n3xccc=dimension of the xccc3d array (0 or nfft).
116 : !! optene=>0 if some additional energies have to be computed
117 : !! pawang <type(pawang_type)> =paw angular mesh and related data
118 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
119 : !! pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom)
120 : !! pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
121 : !! ph1d(2,3*(2*mgfft+1)*natom)=phase (structure factor) information.
122 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
123 : !! rhog(2,nfft)=Fourier transform of electron density
124 : !! rhor(nfft,nspden)=electron density in electrons/bohr**3.
125 : !! | definition for spin components:
126 : !! | case of nspden = 2
127 : !! | rhor(:,1) => rho_up + rho_dwn
128 : !! | rhor(:,2) => rho_up
129 : !! | case of nspden = 4
130 : !! | rhor(:,1) => rho_upup + rho_dwndwn
131 : !! | rhor(:,2:4) => {m_x,m_y,m_z}
132 : !! rmet(3,3)=real space metric (bohr**2)
133 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
134 : !! ucvol = unit cell volume (bohr^3)
135 : !! usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
136 : !! xred(3,natom)=reduced dimensionless atomic coordinates
137 : !! [taur(nfftf,nspden*dtset%usekden)]=array for kinetic energy density
138 : !!
139 : !! OUTPUT
140 : !! energies <type(energies_type)>=all part of total energy.
141 : !! | e_xc=exchange-correlation energy (hartree)
142 : !! | In case of hybrid compensation algorithm:
143 : !! | e_hybcomp_E0=energy compensation term for hybrid exchange-correlation energy (hartree) at fixed density
144 : !! | e_hybcomp_v0=potential compensation term for hybrid exchange-correlation energy (hartree) at fixed density
145 : !! | e_hybcomp_v=potential compensation term for hybrid exchange-correlation energy (hartree) at self-consistent density
146 : !! ==== if optene==2 or 4
147 : !! | e_localpsp=local psp energy (hartree)
148 : !! ==== if dtset%icoulomb == 0
149 : !! | e_ewald=Ewald energy (hartree)
150 : !! ==== if optene>=1
151 : !! | e_hartree=Hartree part of total energy (hartree)
152 : !! ==== if optene==3 or 4
153 : !! | e_xcdc=exchange-correlation double-counting energy (hartree)
154 : !! ==== if dtset%vdw_xc == 5 or 6 or 7
155 : !! | e_vdw_dftd=Dispersion energy from DFT-D Van der Waals correction (hartree)
156 : !! grchempottn(3,natom)=grads of spatially-varying chemical energy (hartree)
157 : !! grewtn(3,natom)=grads of Ewald energy (hartree)
158 : !! grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D2 dispersion (hartree)
159 : !! kxc(nfft,nkxc)=exchange-correlation kernel, will be computed if nkxc/=0 .
160 : !! see routine rhotoxc for a more complete description
161 : !! strsxc(6)=xc contribution to stress tensor (hartree/bohr^3)
162 : !! vxcavg=mean of the vxc potential
163 : !!
164 : !! SIDE EFFECTS
165 : !! moved_atm_inside=1 if the atomic positions were moved inside the SCF loop.
166 : !! vhartr(nfft)=Hartree potential (Hartree)
167 : !! vpsp(nfft)=local psp (Hartree)
168 : !! vtrial(nfft,nspden)= trial potential (Hartree)
169 : !! vxc(nfft,nspden)= xc potential (Hartree)
170 : !! [electronpositron <type(electronpositron_type)>]=quantities for the electron-positron annihilation (optional argument)
171 : !! [vxc_hybcomp(nfft,nspden)= compensation xc potential (Hartree) in case of hybrids] Optional output
172 : !! i.e. difference between the hybrid Vxc at fixed density and the auxiliary Vxc at fixed density
173 : !! [vxctau(nfftf,dtset%nspden,4*usevxctau)]=derivative of XC energy density with respect to
174 : !! kinetic energy density (metaGGA cases) (optional output)
175 : !! xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
176 : !! [xcctau3d(n3xccc*usekden)]=3D core electron kinetic energy density for XC core correction, bohr^-3
177 : !!
178 : !! NOTES
179 : !! In case of PAW calculations:
180 : !! All computations are done on the fine FFT grid.
181 : !! All variables (nfft,ngfft,mgfft) refer to this fine FFT grid.
182 : !! All arrays (densities/potentials...) are computed on this fine FFT grid.
183 : !! Developers have to be careful when introducing others arrays: they have to be stored on the fine FFT grid.
184 : !!
185 : !! In case of norm-conserving calculations the FFT grid is the usual FFT grid.
186 : !!
187 : !! SOURCE
188 :
189 18638 : subroutine setvtr(atindx1,dtset,energies,gmet,gprimd,grchempottn,grewtn,grvdw,gsqcut,&
190 9319 : & istep,kxc,mgfft,moved_atm_inside,moved_rhor,mpi_enreg,&
191 27957 : & nattyp,nfft,ngfft,ngrvdw,nhat,nhatgr,nhatgrdim,nkxc,ntypat,n1xccc,n3xccc,&
192 9319 : & optene,pawang,pawrad,pawrhoij,pawtab,ph1d,psps,rhog,rhor,rmet,rprimd,strsxc,&
193 9319 : & ucvol,usexcnhat,vhartr,vpsp,vtrial,vxc,vxcavg,wvl,xccc3d,xred,&
194 36583 : & 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(:,:,:) !vxctau(nfft,dtset%nspden,4*usevxctau)
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,dum
240 : logical :: add_tfw_,is_hybrid_ncpp,non_magnetic_xc,with_vxctau,wvlbigdft,lewald
241 : character(len=500) :: message
242 9319 : type(constrained_dft_t) :: constrained_dft
243 : type(xcdata_type) :: xcdata,xcdatahyb
244 : !arrays
245 9319 : real(dp),allocatable :: xcart(:,:)
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) :: vhspinfield(4)
253 9319 : real(dp),allocatable :: grtn(:,:),dyfr_dum(:,:,:),gr_dum(:,:)
254 9319 : real(dp),allocatable :: rhojellg(:,:),rhojellr(:),rhowk(:,:),vjell(:)
255 9319 : real(dp),allocatable :: v_constr_dft_r(:,:),rhog_dum(:,:)
256 :
257 : ! *********************************************************************
258 :
259 9319 : call timab(91,1,tsec)
260 :
261 : !Test size of kinetic energy potential Vxctau
262 9319 : with_vxctau = (present(vxctau).and.present(taur))
263 35914 : if (with_vxctau) with_vxctau = (size(vxctau)>0.and.dtset%usekden/=0)
264 : if (with_vxctau) then
265 204 : if (size(vxctau)/=nfft*dtset%nspden*4) then
266 0 : ABI_BUG("Wrong size for vxctau!")
267 : end if
268 : end if
269 :
270 : !Check if we're in hybrid norm conserving pseudopotential with a core correction
271 : is_hybrid_ncpp=(dtset%usepaw==0 .and. n3xccc/=0 .and. &
272 9319 : & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid()))
273 :
274 : !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed
275 18638 : wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1)
276 :
277 : !Get size of FFT grid
278 37276 : nfftot=PRODUCT(ngfft(1:3))
279 :
280 : !mpi
281 9319 : mpi_comm_sphgrid=mpi_enreg%comm_fft
282 9319 : if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl
283 :
284 : !Test electron-positron case
285 9319 : ipositron=0;if (present(electronpositron)) ipositron=electronpositron_calctype(electronpositron)
286 :
287 : !Test addition of Weiszacker gradient correction to Thomas-Fermi kin energy
288 9319 : add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw
289 :
290 : !Get Ewald energy and Ewald forces, as well as vdW-DFTD energy and forces, and chemical potential energy and forces.
291 : !-------------------------------------------------------------------------------------------------------------------
292 9319 : call timab(5,1,tsec)
293 9319 : if (ipositron/=1) then
294 9276 : lewald = .true.
295 9276 : if (present(calc_ewald)) lewald=calc_ewald
296 1755 : if (lewald) then
297 7574 : if (dtset%icoulomb == 0 .or. (dtset%usewvl == 0 .and. dtset%icoulomb == 2)) then
298 : ! Periodic system, need to compute energy and forces due to replica and
299 : ! to correct the shift in potential calculation.
300 : call ewald(energies%e_ewald,gmet,grewtn,gsqcut,dtset%icutcoul,dtset%natom,ngfft,dtset%nkpt,ntypat,&
301 7574 : &dtset%rcut,rmet,rprimd,dtset%typat,ucvol,dtset%vcutgeo,xred,psps%ziontypat)
302 : ! For a periodic system bearing a finite charge, the monopole correction to the
303 : ! energy is relevant.
304 : ! See Leslie and Gillan, JOURNAL OF PHYSICS C-SOLID STATE PHYSICS 18, 973 (1985)
305 7574 : if(abs(dtset%cellcharge(1))>tol8) then
306 : call ewald(energies%e_monopole,gmet,grewtn_fake,gsqcut,dtset%icutcoul,1,ngfft,dtset%nkpt,1,&
307 104 : &dtset%rcut,rmet,rprimd,(/1/),ucvol,dtset%vcutgeo,(/0.0_dp,0.0_dp,0.0_dp/),(/dtset%cellcharge(1)/))
308 52 : energies%e_monopole=-energies%e_monopole
309 : end if
310 0 : else if (dtset%icoulomb == 1) then
311 : ! In a non periodic system (real space computation), the G=0 divergence
312 : ! doesn't occur and ewald is not needed. Only the ion/ion interaction
313 : ! energy is relevant and used as Ewald energy and gradient.
314 0 : call ionion_realSpace(dtset, energies%e_ewald, grewtn, rprimd, xred, psps%ziontypat)
315 0 : else if (dtset%icoulomb == 2) then
316 : call ionion_surface(dtset, energies%e_ewald, grewtn, mpi_enreg%me_wvl, mpi_enreg%nproc_wvl, rprimd, &
317 0 : & wvl%descr, wvl%den, xred)
318 : end if
319 : end if
320 9276 : if (dtset%nzchempot>0) then
321 23 : call spatialchempot(energies%e_chempot,dtset%chempot,grchempottn,dtset%natom,ntypat,dtset%nzchempot,dtset%typat,xred)
322 : end if
323 9276 : if (dtset%vdw_xc==5.and.ngrvdw==dtset%natom) then
324 : call vdw_dftd2(energies%e_vdw_dftd,dtset%ixc,dtset%natom,ntypat,1,dtset%typat,rprimd,&
325 10 : & dtset%vdw_tol,xred,psps%znucltypat,gred_vdw_dftd2=grvdw)
326 : end if
327 9276 : if ((dtset%vdw_xc==6.or.dtset%vdw_xc==7).and.ngrvdw==dtset%natom) then
328 : call vdw_dftd3(energies%e_vdw_dftd,dtset%ixc,dtset%natom,&
329 : & ntypat,1,dtset%typat,rprimd,dtset%vdw_xc,dtset%vdw_tol,dtset%vdw_tol_3bt,&
330 8 : & xred,psps%znucltypat,gred_vdw_dftd3=grvdw)
331 : end if
332 : else
333 43 : energies%e_ewald=zero
334 43 : energies%e_chempot=zero
335 1975 : grchempottn=zero
336 1975 : grewtn=zero
337 43 : energies%e_vdw_dftd=zero
338 43 : if (ngrvdw>0) grvdw=zero
339 : end if
340 9319 : call timab(5,2,tsec)
341 :
342 : !Compute parts of total energy depending on potentials
343 : !--------------------------------------------------------------
344 9319 : if (dtset%usewvl == 0) then
345 9319 : ucvol_local = ucvol
346 : #if defined HAVE_BIGDFT
347 : else
348 : ! We need to tune the volume when wavelets are used because, not all FFT points are used.
349 : ! ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3)
350 : ucvol_local = product(wvl%den%denspot%dpbox%hgrids) * real(product(wvl%den%denspot%dpbox%ndims), dp)
351 : #endif
352 : end if
353 :
354 : !Determine by which method the local ionic potential and/or the pseudo core charge density
355 : ! have to be computed
356 : !Local ionic potential:
357 : ! Method 1: PAW
358 : ! Method 2: Norm-conserving PP, icoulomb>0, wavelets
359 9319 : vloc_method=1;if (psps%usepaw==0) vloc_method=2
360 9319 : if (dtset%icoulomb>0) vloc_method=2
361 9319 : if (psps%usewvl==1) vloc_method=2
362 : !Pseudo core charge density:
363 : ! Method 1: PAW, nc_xccc_gspace
364 : ! Method 2: Norm-conserving PP, wavelets
365 : coredens_method=1;if (psps%usepaw==0) coredens_method=2
366 9319 : if (psps%nc_xccc_gspace==1) coredens_method=1
367 9319 : if (psps%nc_xccc_gspace==0) coredens_method=2
368 9319 : if (psps%usewvl==1) coredens_method=2
369 9319 : coretau_method=0
370 9319 : if (dtset%usekden==1) then
371 74 : coretau_method=1;if (psps%nc_xccc_gspace==0) coretau_method=2
372 : end if
373 : !In some specific cases, XC has to be handled as non-magnetic
374 9319 : non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
375 :
376 : !Local ionic potential and/or pseudo core charge by method 1
377 9319 : if (vloc_method==1.or.coredens_method==1) then
378 3130 : call timab(552,1,tsec)
379 3130 : optv=0;if (vloc_method==1) optv=1
380 3130 : optn=0;if (coredens_method==1) optn=n3xccc/nfft
381 3130 : optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=1
382 : call atm2fft(atindx1,xccc3d,vpsp,dummy_out1,dummy_out2,dummy_out3,dummy_in,&
383 : & gmet,gprimd,dummy_out4,dummy_out5,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,&
384 : & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,&
385 : & dtset%rcut,dummy_in,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,&
386 : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
387 3130 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
388 3130 : call timab(552,2,tsec)
389 : end if
390 9319 : if (coretau_method==1) then
391 34 : call timab(552,1,tsec)
392 34 : optv=0;optn=1
393 34 : optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optn2=4
394 : call atm2fft(atindx1,xcctau3d,dummy_out6,dummy_out1,dummy_out2,dummy_out3,dummy_in,&
395 : & gmet,gprimd,dummy_out4,dummy_out5,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,&
396 : & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,&
397 : & dtset%rcut,dummy_in,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,&
398 : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
399 34 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
400 34 : call timab(552,2,tsec)
401 : end if
402 :
403 : !Local ionic potential by method 2
404 9319 : if (vloc_method==2) then
405 6192 : option=1
406 18576 : ABI_MALLOC(gr_dum,(3,dtset%natom))
407 18576 : ABI_MALLOC(dyfr_dum,(3,3,dtset%natom))
408 18576 : ABI_MALLOC(rhog_dum,(2,nfft))
409 : call mklocl(dtset,dyfr_dum,energies%e_localpsp,gmet,gprimd,&
410 : & gr_dum,gsqcut,dummy6,mgfft,mpi_enreg,dtset%natom,nattyp,&
411 : & nfft,ngfft,dtset%nspden,ntypat,option,pawtab,ph1d,psps,&
412 6192 : & dtset%qprtrb,rhog_dum,rhor,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred)
413 6192 : ABI_FREE(gr_dum)
414 6192 : ABI_FREE(dyfr_dum)
415 6192 : ABI_FREE(rhog_dum)
416 : end if
417 :
418 : !3D pseudo core electron density xccc3d by method 2
419 9319 : if (coredens_method==2.and.n1xccc/=0) then
420 2853 : call timab(91,2,tsec)
421 2853 : call timab(92,1,tsec)
422 2853 : option=1
423 8559 : ABI_MALLOC(gr_dum,(3,dtset%natom))
424 8559 : ABI_MALLOC(dyfr_dum,(3,3,dtset%natom))
425 2853 : if (psps%usewvl==0.and.psps%usepaw==0.and.dtset%icoulomb==0) then
426 : call mkcore(dummy6,dyfr_dum,gr_dum,mpi_enreg,dtset%natom,nfft,dtset%nspden,ntypat,&
427 : & ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,dtset%typat,ucvol,&
428 2853 : & vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred)
429 0 : else if (psps%usewvl==0.and.(psps%usepaw==1.or.dtset%icoulomb==1)) then
430 : call mkcore_alt(atindx1,dummy6,dyfr_dum,gr_dum,dtset%icoulomb,mpi_enreg,dtset%natom,&
431 : & nfft,dtset%nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,&
432 0 : & ucvol,vxc,psps%xcccrc,psps%xccc1d,xccc3d,xred,pawrad,pawtab,psps%usepaw)
433 : else if (psps%usewvl==1.and.psps%usepaw==1) then
434 : #if defined HAVE_BIGDFT
435 : ! call mkcore_wvl_old(atindx1,dummy6,dyfr_dum,wvl%descr%atoms%astruct%geocode,gr_dum,wvl%descr%h,&
436 : ! & dtset%natom,nattyp,nfft,wvl%den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl,:),&
437 : ! & dtset%nspden,ntypat,wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,wvl%descr%Glr%d%n2,&
438 : ! & wvl%descr%Glr%d%n2i,wvl%descr%Glr%d%n3,wvl%den%denspot%dpbox%n3pi,n3xccc,option,&
439 : ! & pawrad,pawtab,psps%gth_params%psppar,rprimd,ucvol_local,vxc,xccc3d,xred,&
440 : ! & mpi_comm_wvl=mpi_enreg%comm_wvl)
441 : call mkcore_wvl(atindx1,dummy6,gr_dum,dtset%natom,nattyp,nfft,dtset%nspden,ntypat,&
442 : & n1xccc,n3xccc,option,pawrad,pawtab,rprimd,vxc,psps%xccc1d,xccc3d,&
443 : & psps%xcccrc,xred,wvl%den,wvl%descr,mpi_comm_wvl=mpi_enreg%comm_wvl)
444 : #endif
445 : end if
446 2853 : ABI_FREE(gr_dum)
447 2853 : ABI_FREE(dyfr_dum)
448 2853 : call timab(92,2,tsec)
449 2853 : call timab(91,1,tsec)
450 : end if
451 9319 : if (coretau_method==2) then
452 40 : call timab(91,2,tsec)
453 40 : call timab(92,1,tsec)
454 40 : option=1
455 120 : ABI_MALLOC(gr_dum,(3,dtset%natom))
456 120 : ABI_MALLOC(dyfr_dum,(3,3,dtset%natom))
457 : call mkcore_alt(atindx1,dummy6,dyfr_dum,gr_dum,dtset%icoulomb,mpi_enreg,dtset%natom,&
458 : & nfft,dtset%nspden,nattyp,ntypat,ngfft(1),n1xccc,ngfft(2),ngfft(3),option,rprimd,&
459 : & ucvol,vxc,psps%xcccrc,psps%xccc1d,xcctau3d,xred,pawrad,pawtab,psps%usepaw,&
460 40 : & usekden=.true.)
461 40 : ABI_FREE(gr_dum)
462 40 : ABI_FREE(dyfr_dum)
463 40 : call timab(92,2,tsec)
464 40 : call timab(91,1,tsec)
465 : end if
466 :
467 : !Adds the jellium potential to the local part of ionic potential
468 9319 : if (dtset%jellslab/=0) then
469 18 : ABI_MALLOC(vjell,(nfft))
470 18 : ABI_MALLOC(rhojellg,(2,nfft))
471 12 : ABI_MALLOC(rhojellr,(nfft))
472 6 : option=1
473 : call jellium(gmet,gsqcut,mpi_enreg,nfft,ngfft,dtset%nspden,option,&
474 6 : & dtset%slabwsrad,rhojellg,rhojellr,rprimd,vjell,dtset%slabzbeg,dtset%slabzend)
475 : ! Compute background-background energy
476 6 : call dotprod_vn(1,rhojellr,ebb,doti,nfft,nfftot,1,1,vjell,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
477 6 : ebb=half*ebb
478 : ! Compute electrostatic energy between background and nuclei before adding vjell to vpsp
479 6 : call dotprod_vn(1,rhojellr,ebn,doti,nfft,nfftot,1,1,vpsp,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
480 : ! Update e_ewald with ebb and ebn
481 6 : energies%e_ewald=energies%e_ewald+ebb+ebn
482 : ! Compute gradient of ebn wrt tn
483 : ! This is not yet coded for usewvl or icoulomb=1
484 6 : if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
485 : write(message,'(3a)')&
486 0 : & 'The computation of forces due to jellium background',ch10,&
487 0 : & 'has to be verified in the PAW formalism.'
488 0 : ABI_WARNING(message)
489 :
490 0 : ABI_MALLOC(grtn,(3,dtset%natom))
491 0 : optatm=0;optdyfr=0;opteltfr=0;optgr=1;optstr=0;optv=1;optn=0;optn2=1
492 : call atm2fft(atindx1,dummy_out1,vpsp,dummy_out2,dummy_out3,dummy_out4,dummy_in,&
493 : & gmet,gprimd,dummy_out5,grtn,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,ntypat,&
494 : & optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,&
495 : & dtset%rcut,rhojellg,rprimd,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,dummy_in,dummy_in,dummy_in,dtset%vprtrb,psps%vlspl,&
496 : & comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
497 0 : & paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
498 :
499 : ! Update grewtn with gradient of ebn wrt tn
500 0 : do iatom=1,dtset%natom
501 0 : grewtn(1:3,iatom)=grewtn(1:3,iatom)+grtn(1:3,iatom)
502 : end do
503 0 : ABI_FREE(grtn)
504 : else ! of usepaw==1
505 6 : option=2
506 18 : ABI_MALLOC(dyfr_dum,(3,3,dtset%natom))
507 18 : ABI_MALLOC(grtn,(3,dtset%natom))
508 : call mklocl(dtset,dyfr_dum,energies%e_localpsp,gmet,gprimd,&
509 : & grtn,gsqcut,dummy6,mgfft,mpi_enreg,dtset%natom,nattyp,&
510 : & nfft,ngfft,1,ntypat,option,pawtab,ph1d,psps,dtset%qprtrb,rhojellg,&
511 6 : & rhojellr,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred)
512 : ! Update grewtn with gradient of ebn wrt tn (reestablish order of atoms)
513 12 : do iatom=1,dtset%natom
514 30 : grewtn(1:3,atindx1(iatom))=grewtn(1:3,atindx1(iatom))+grtn(1:3,iatom)
515 : end do
516 6 : ABI_FREE(dyfr_dum)
517 6 : ABI_FREE(grtn)
518 : end if ! of usepaw==1
519 104886 : vpsp(:)=vpsp(:)+vjell(:)
520 6 : ABI_FREE(vjell)
521 6 : ABI_FREE(rhojellg)
522 6 : ABI_FREE(rhojellr)
523 : end if
524 :
525 : !Additional stuff for electron-positron calculation
526 : !Compute the electronic/positronic local (Hartree) potential
527 4917190 : if (ipositron==1) vpsp=-vpsp
528 :
529 : !If we are at the initialisation, or
530 : !if the atom positions has changed and the non-linear core correction
531 : !is included, or the rhor has changed, one needs to compute the xc stuff.
532 : !One needs also to compute the Hartree stuff if the density changed,
533 : !or at initialisation.
534 : !--------------------------------------------------------------
535 :
536 9319 : if(istep==1 .or. n1xccc/=0 .or. moved_rhor==1 .or. dtset%positron<0 .or. mod(dtset%fockoptmix,100)==11) then
537 :
538 9319 : option=0
539 9319 : if(istep==1 .or. moved_rhor==1 .or. dtset%positron<0 .or. mod(dtset%fockoptmix,100)==11) option=1
540 9319 : if (nkxc>0) option=2
541 9319 : if (dtset%iscf==-1) option=-2
542 :
543 9319 : if (ipositron/=1) then
544 9276 : if (dtset%icoulomb == 0 .and. dtset%usewvl == 0) then
545 :
546 : ! >>>> Hartree potential
547 9276 : if(option/=0 .and. option/=10)then
548 : call hartre(1,gsqcut,dtset%icutcoul,psps%usepaw,mpi_enreg,nfft,ngfft,&
549 9202 : &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr)
550 : end if
551 :
552 : ! >>>> Exchange-correlation potential
553 9276 : call xcdata_init(xcdata,dtset=dtset)
554 9276 : if(mod(dtset%fockoptmix,100)==11)then
555 18 : xcdatahyb=xcdata
556 : ! Setup the auxiliary xc functional information
557 18 : call xcdata_init(xcdata,dtset=dtset,auxc_ixc=0,ixc=dtset%auxc_ixc)
558 : end if
559 : ! Not yet able to deal fully with the full XC kernel in case of GGA + spin
560 9276 : option_eff=option;if (option==2.and.xcdata%xclevel==2.and.(nkxc==3-2*mod(xcdata%nspden,2))) option_eff=12
561 9276 : nk3xc=1
562 :
563 : ! If we use the XC Tran-Blaha 2009 (modified BJ) functional, update the c value
564 9276 : if (dtset%xc_tb09_c>99._dp) then
565 : call xc_tb09_update_c(dtset%intxc,dtset%ixc,mpi_enreg,dtset%natom, &
566 : & nfft,ngfft,nhat,psps%usepaw,nhatgr,nhatgrdim,dtset%nspden,dtset%ntypat,n3xccc, &
567 : & pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,rhor,rprimd,psps%usepaw, &
568 : & xccc3d,dtset%xc_denpos,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, &
569 9275 : & computation_type='all')
570 : end if
571 :
572 9276 : if (ipositron==0) then
573 :
574 : ! Compute energies%e_xc and associated quantities
575 9240 : if(.not.is_hybrid_ncpp .or. mod(dtset%fockoptmix,100)==11)then
576 : call rhotoxc(energies%e_xc,energies%entropy_xc,kxc,mpi_enreg,nfft,ngfft,&
577 : & nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,&
578 : & option_eff,rhor,rprimd,usexcnhat,vxc,vxcavg,xccc3d,xcdata,strsxc=strsxc,&
579 10350 : & taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_,xcctau3d=xcctau3d,grho1_over_rho1=rpnrm)
580 : else
581 : ! Only when is_hybrid_ncpp, and moreover, the xc functional is not the auxiliary xc functional, then call xchybrid_ncpp_cc
582 : call xchybrid_ncpp_cc(dtset,energies%e_xc,energies%entropy_xc,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,&
583 33 : & strsxc,vxcavg,xccc3d,vxc=vxc)
584 : end if
585 :
586 : ! Possibly compute energies%e_hybcomp_E0
587 9240 : if(mod(dtset%fockoptmix,100)==11)then
588 : ! This call to rhotoxc uses the hybrid xc functional
589 18 : if(.not.is_hybrid_ncpp)then
590 : call rhotoxc(energies%e_hybcomp_E0,dum,kxc,mpi_enreg,nfft,ngfft,&
591 : & nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,&
592 : & option_eff,rhor,rprimd,usexcnhat,vxc_hybcomp,vxcavg,xccc3d,xcdatahyb,&
593 18 : & strsxc=strsxc,taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_)
594 : else
595 : call xchybrid_ncpp_cc(dtset,energies%e_hybcomp_E0,dum,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,&
596 0 : & strsxc,vxcavg,xccc3d,vxc=vxc_hybcomp)
597 : end if
598 :
599 : ! Combine hybrid and auxiliary quantities
600 18 : energies%e_xc=energies%e_xc*dtset%auxc_scal
601 18 : energies%entropy_xc=energies%entropy_xc*dtset%auxc_scal
602 18 : energies%e_hybcomp_E0=energies%e_hybcomp_E0-energies%e_xc
603 73764 : vxc(:,:)=vxc(:,:)*dtset%auxc_scal
604 73764 : vxc_hybcomp(:,:)=vxc_hybcomp(:,:)-vxc(:,:)
605 : end if
606 :
607 36 : else if (ipositron==2) then
608 : call rhotoxc(energies%e_xc,energies%entropy_xc,kxc,mpi_enreg,nfft,ngfft,&
609 : & nhat,psps%usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,&
610 : & option_eff,rhor,rprimd,usexcnhat,vxc,vxcavg,xccc3d,xcdata,strsxc=strsxc,&
611 : & taur=taur,vhartr=vhartr,vxctau=vxctau,add_tfw=add_tfw_,&
612 36 : & electronpositron=electronpositron)
613 : end if
614 :
615 0 : elseif(.not. wvlbigdft) then
616 : ! Use the free boundary solver
617 : call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, &
618 : & dtset%icoulomb, dtset%ixc, &
619 : & mpi_enreg, nfft, ngfft,&
620 : & nhat,psps%usepaw,&
621 : & dtset%nscforder,dtset%nspden,n3xccc,rhor,rprimd, &
622 : & usexcnhat,psps%usepaw,dtset%usewvl,vhartr, vxc, &
623 : & vxcavg,wvl%descr,wvl%den,wvl%e,&
624 0 : & xccc3d,dtset%xclevel,dtset%xc_denpos)
625 : end if
626 : else
627 43 : energies%e_xc=zero
628 : call rhohxcpositron(electronpositron,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat,nkxc,dtset%nspden,n3xccc,&
629 43 : & dtset%paral_kgb,rhor,strsxc,ucvol,usexcnhat,psps%usepaw,vhartr,vxc,vxcavg,xccc3d,dtset%xc_denpos)
630 : end if
631 9319 : if (ipositron/=0) then
632 79 : if (optene>=1) then
633 : call dotprod_vn(1,rhor,electronpositron%e_hartree,doti,&
634 79 : & nfft,nfftot,1,1,electronpositron%vha_ep,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
635 : end if
636 9924334 : vhartr=vhartr+electronpositron%vha_ep
637 : end if
638 : end if
639 :
640 : !Compute the trial potential
641 : !-------------------------------------------------------------
642 9319 : if (.not. wvlbigdft) then
643 : ! Now, compute trial Hxc potential. Local psp potential will be added later.
644 9319 : if(moved_atm_inside==0 .or.dtset%iscf>=10) then
645 :
646 : ! Compute starting Hxc potential.
647 : ! Multiply by identity, should not change anything if nspden /= 4
648 19549 : do ispden=1,dtset%nspden
649 229968800 : vtrial(:,ispden)=vhartr(:)*identity(ispden)+vxc(:,ispden)
650 : end do
651 :
652 : else
653 :
654 : ! One should be here only when moved_atm_inside==1
655 : ! The (H)xc now added corrects the previous one.
656 106 : if(dtset%densfor_pred==1)then
657 : ! xc was substracted off. This should be rationalized later
658 86 : do ispden=1,dtset%nspden
659 461910 : vtrial(:,ispden)=vtrial(:,ispden)+vxc(:,ispden)
660 : end do
661 71 : else if(abs(dtset%densfor_pred)==2.or.abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)then
662 : ! Hxc was substracted off. This should be rationalized later
663 84 : do ispden=1,dtset%nspden
664 1573623 : vtrial(:,ispden)=vtrial(:,ispden)+vhartr(:)*identity(ispden)+vxc(:,ispden)
665 : end do
666 : end if
667 : end if
668 :
669 : ! Adds the local part of the potential
670 9319 : if ((moved_atm_inside==0).or.(dtset%densfor_pred/=3)) then
671 19423 : do ispden=1,min(2,dtset%nspden)
672 228287743 : do ifft=1,nfft
673 228278424 : vtrial(ifft,ispden)=vtrial(ifft,ispden)+vpsp(ifft)
674 : end do
675 : end do
676 : end if
677 :
678 : ! Adds the compensating vxc for hybrids
679 9319 : if(mod(dtset%fockoptmix,100)==11)then
680 73764 : vtrial(:,:)=vtrial(:,:)+vxc_hybcomp(:,:)
681 : end if
682 :
683 9319 : if(dtset%usewvl==1) then
684 0 : call wvl_vtrial_abi2big(1,vtrial,wvl%den)
685 : end if
686 :
687 : else
688 :
689 : ! Compute with covering comms the different part of the potential.
690 : #if defined HAVE_BIGDFT
691 : ! Copy e_ewald.
692 : wvl%e%energs%eion = energies%e_ewald
693 : ! Setup the mixing, if necessary
694 : call denspot_set_history(wvl%den%denspot,dtset%iscf,dtset%nsppol, &
695 : & wvl%den%denspot%dpbox%ndims(1),wvl%den%denspot%dpbox%ndims(2))
696 : #endif
697 0 : ABI_MALLOC(xcart,(3, dtset%natom))
698 0 : call xred2xcart(dtset%natom, rprimd, xcart, xred)
699 : call wvl_psitohpsi(dtset%diemix,energies%e_exactX, energies%e_xc, energies%e_hartree, &
700 : & energies%e_kinetic, energies%e_localpsp, energies%e_nlpsp_vfock, energies%e_sicdc, &
701 : & istep, 1, dtset%iscf, mpi_enreg%me_wvl, dtset%natom, dtset%nfft, mpi_enreg%nproc_wvl, dtset%nspden, &
702 : & rpnrm, .true.,evxc, wvl,.true., xcart, strsxc,&
703 0 : & vtrial, vxc)
704 0 : if (optene==3.or.optene==4) energies%e_xcdc=evxc
705 0 : ABI_FREE(xcart)
706 :
707 : end if
708 :
709 : !Add the hspinfield to vtrial
710 37258 : if (any(abs(dtset%hspinfield(:))>tol8)) then
711 12 : vhspinfield(:) = zero ! vhspinfield_ij = -1/2*sigma_ij^alpha*B_alpha
712 12 : if(dtset%nspden==2)then
713 2 : vhspinfield(1) = -half*dtset%hspinfield(3) ! v_dwndwn = -1/2*B_z
714 2 : vhspinfield(2) = half*dtset%hspinfield(3) ! v_upup = 1/2*B_z
715 39368 : do ifft=1,nfft
716 39366 : vtrial(ifft,1) = vtrial(ifft,1) + vhspinfield(1) !SPr: added 1st component
717 39368 : vtrial(ifft,2) = vtrial(ifft,2) + vhspinfield(2)
718 : end do !ifft
719 : end if
720 12 : if(dtset%nspden==4)then
721 10 : vhspinfield(1)=-half*dtset%hspinfield(3) ! v_dwndwn => v_11
722 10 : vhspinfield(2)= half*dtset%hspinfield(3) ! v_upup => v_22
723 10 : vhspinfield(3)=-half*dtset%hspinfield(1) ! Re(v_dwnup) = Re(v_updwn) => Re(v_12)
724 10 : vhspinfield(4)= half*dtset%hspinfield(2) ! Im(v_dwnup) =-Im(v_dwnup) => Im(v_12)
725 50 : do ispden=1,dtset%nspden
726 787370 : do ifft=1,nfft
727 787360 : vtrial(ifft,ispden) = vtrial(ifft,ispden) + vhspinfield(ispden)
728 : end do
729 : end do
730 : end if
731 : end if
732 :
733 : !Compute the constrained potential for the magnetic moments
734 9319 : if (dtset%magconon==1.or.dtset%magconon==2) then
735 : ! Initialize the datastructure constrained_dft, for penalty function constrained magnetization
736 : call constrained_dft_ini(dtset%chrgat,constrained_dft,dtset%constraint_kind,dtset%magconon,dtset%magcon_lambda,&
737 : & mpi_enreg,dtset%natom,dtset%nfft,dtset%ngfft,dtset%nspden,dtset%ntypat,dtset%ratsm,&
738 5 : & dtset%ratsph,rprimd,dtset%spinat,dtset%typat,xred,dtset%ziontypat,dtset%znucl,dtset%qgbt,dtset%use_gbt,dtset%spinaxis)
739 20 : ABI_MALLOC(v_constr_dft_r, (nfft,dtset%nspden))
740 65557 : v_constr_dft_r = zero
741 5 : call mag_penalty(constrained_dft,mpi_enreg,rhor,v_constr_dft_r,xred,dtset%qgbt,dtset%use_gbt)
742 5 : if(dtset%nspden==4)then
743 15 : do ispden=1,dtset%nspden ! (SPr: both components should be used? EB: Yes it should be the case, corrected now)
744 49167 : do ifft=1,nfft
745 49164 : vtrial(ifft,ispden) = vtrial(ifft,ispden) + v_constr_dft_r(ifft,ispden)
746 : end do !ifft
747 : end do !ispden
748 2 : else if(dtset%nspden==2)then
749 8194 : do ifft=1,nfft
750 : ! TODO : MJV: check that magnetic constraint works also for nspden 2 or add input variable condition
751 : ! EB: ispden=2 is rho_up only: to be tested
752 : ! SPr: for ispden=2, both components should be used (e.g. see definition for vhspinfield)?
753 8192 : vtrial(ifft,1) = vtrial(ifft,1) + v_constr_dft_r(ifft,1) !SPr: added the first component here
754 8194 : vtrial(ifft,2) = vtrial(ifft,2) + v_constr_dft_r(ifft,2)
755 : end do !ifft
756 : end if
757 5 : ABI_FREE(v_constr_dft_r)
758 5 : call constrained_dft_free(constrained_dft)
759 : end if
760 :
761 : !Compute parts of total energy depending on potentials
762 : !--------------------------------------------------------------
763 :
764 : !For icoulomb==0 or usewvl Ehartree is calculated in psolver_rhohxc().
765 : !For PAW we recalculate this since nhat was not taken into account
766 : !in psolver_rhohxc: E_H= int v_H (n+nhat) dr
767 :
768 9319 : if (optene>=1 .and. .not. wvlbigdft .and. (dtset%icoulomb==0 .or. dtset%usepaw==1 ) ) then
769 : ! Compute Hartree energy ehart
770 : ! Already available in the Psolver case through psolver_rhohxc().
771 3714 : if (ipositron/=1) then
772 : call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol_local,&
773 3671 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
774 3671 : if (ipositron==0) energies%e_hartree = half * energies%e_hartree
775 3671 : if (ipositron==2) energies%e_hartree = half * (energies%e_hartree-electronpositron%e_hartree)
776 : else
777 43 : energies%e_hartree=zero
778 : end if
779 : end if
780 :
781 9319 : if(mod(dtset%fockoptmix,100)==11)then
782 18 : if (.not. wvlbigdft) then
783 : call dotprod_vn(1,rhor,energies%e_hybcomp_v0,doti,nfft,nfftot,1,1,vxc_hybcomp,ucvol_local,&
784 18 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
785 18 : energies%e_hybcomp_v=energies%e_hybcomp_v0
786 : end if
787 : end if
788 :
789 9319 : if (optene==2.or.optene==4 .and. .not. wvlbigdft) then
790 : ! Compute local psp energy eei
791 : call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfft,nfftot,1,1,vpsp,ucvol_local,&
792 3714 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
793 : end if
794 :
795 9319 : if (optene==3.or.optene==4 .and. .not. wvlbigdft) then
796 : ! Compute double-counting XC energy enxcdc
797 3714 : if (ipositron/=1) then
798 3671 : if (dtset%usepaw==0.or.usexcnhat/=0) then
799 : call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,&
800 1314 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
801 : else
802 9428 : ABI_MALLOC(rhowk,(nfft,dtset%nspden))
803 53785705 : rhowk=rhor-nhat
804 : call dotprod_vn(1,rhowk,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol_local,&
805 2357 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
806 2357 : ABI_FREE(rhowk)
807 : end if
808 3671 : if (with_vxctau) then
809 : call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfft,nfftot,dtset%nspden,1,vxctau(:,:,1),&
810 34 : & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid)
811 34 : energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau
812 : end if
813 3671 : if (ipositron==2) energies%e_xcdc=energies%e_xcdc-electronpositron%e_xcdc
814 : else
815 43 : energies%e_xcdc=zero
816 : end if
817 : end if
818 :
819 : !--------------------------------------------------------------
820 :
821 : !The initialisation for the new atomic positions has been done
822 9319 : moved_atm_inside=0
823 :
824 9319 : call timab(91,2,tsec)
825 :
826 27928 : end subroutine setvtr
827 : !!***
828 :
829 : !!****m* m_setvtr/spatialchempot
830 : !! NAME
831 : !! spatialchempot
832 : !! FUNCTION
833 : !! Treat spatially varying chemical potential.
834 : !! Compute energy and derivative with respect to dimensionless reduced atom coordinates of the
835 : !! spatially varying chemical potential. No contribution to stresses.
836 : !!
837 : !! INPUTS
838 : !! chempot(3,nzchempot,ntypat)=input array with information about the chemical potential (see input variable description)
839 : !! natom=number of atoms in unit cell
840 : !! ntypat=number of type of atoms
841 : !! nzchempot=number of limiting planes for chemical potential
842 : !! typat(natom)=integer label of each type of atom (1,2,...)
843 : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
844 : !!
845 : !! OUTPUT
846 : !! e_chempot=chemical potential energy in hartrees
847 : !! grchempottn(3,natom)=grads of e_chempot wrt xred(3,natom), hartrees.
848 : !!
849 : !! SOURCE
850 :
851 23 : subroutine spatialchempot(e_chempot,chempot,grchempottn,natom,ntypat,nzchempot,typat,xred)
852 :
853 : !Arguments ------------------------------------
854 : !scalars
855 : integer,intent(in) :: natom,ntypat,nzchempot
856 : real(dp),intent(out) :: e_chempot
857 : !arrays
858 : integer,intent(in) :: typat(natom)
859 : real(dp),intent(in) :: chempot(3,nzchempot,ntypat),xred(3,natom)
860 : real(dp),intent(out) :: grchempottn(3,natom)
861 :
862 : !Local variables-------------------------------
863 : !scalars
864 : integer :: iatom,itypat,iz
865 : real(dp) :: a_2,a_3,cp0,cp1,dcp0,dcp1,ddz,deltaz,deltaziz
866 : real(dp) :: dqz,dz1,qz,zred,z0
867 : !character(len=500) :: message
868 :
869 : ! *************************************************************************
870 :
871 : !DEBUG
872 : !write(std_out,'(a)')' spatialchempot : enter '
873 : !write(std_out,'(a,2i6)')' nzchempot,ntypat',nzchempot,ntypat
874 : !write(std_out,'(a,6es13.3)') ' chempot(1:3,1:2,1)=',chempot(1:3,1:2,1)
875 : !write(std_out,'(a,6es13.3)') ' chempot(1:3,1:2,2)=',chempot(1:3,1:2,2)
876 : !ENDDEBUG
877 :
878 23 : e_chempot=zero
879 207 : grchempottn(:,:)=zero
880 :
881 : !Loop on the different atoms
882 69 : do iatom=1,natom
883 :
884 46 : itypat=typat(iatom)
885 46 : zred=xred(3,iatom)
886 :
887 : ! Determine the delimiting plane just lower to zred
888 : ! First compute the relative zred with respect to the first delimiting plane
889 : ! Take into account a tolerance :
890 46 : deltaz=zred-chempot(1,1,itypat)
891 46 : deltaz=modulo(deltaz+tol12,1.0d0)-tol12
892 : ! deltaz is positive (or higher than -tol12), and lower than one-tol12.
893 56 : do iz=2,nzchempot+1
894 56 : if(iz/=nzchempot+1)then
895 46 : deltaziz=chempot(1,iz,itypat)-chempot(1,1,itypat)
896 : else
897 : deltaziz=one
898 : end if
899 56 : if(deltaziz>deltaz)exit
900 : end do
901 :
902 : ! Defines coordinates and values inside the delimiting interval,
903 : ! with respect to the lower delimiting plane
904 46 : z0=chempot(1,iz-1,itypat)-chempot(1,1,itypat) ; cp0=chempot(2,iz-1,itypat) ; dcp0=chempot(3,iz-1,itypat)
905 46 : if(iz/=nzchempot+1)then
906 36 : dz1=chempot(1,iz,itypat)-chempot(1,iz-1,itypat) ; cp1=chempot(2,iz,itypat) ; dcp1=chempot(3,iz,itypat)
907 : else
908 10 : dz1=(chempot(1,1,itypat)+one)-chempot(1,nzchempot,itypat) ; cp1=chempot(2,1,itypat) ; dcp1=chempot(3,1,itypat)
909 : end if
910 46 : ddz=deltaz-z0
911 :
912 : !DEBUG
913 : ! write(std_out,'(a,2i5)')' Delimiting planes, iz-1 and iz=',iz-1,iz
914 : ! write(std_out,'(a,2es13.3)')' z0, dz1= :',z0,dz1
915 : ! write(std_out,'(a,2es13.3)')' cp0, cp1= :',cp0,cp1
916 : ! write(std_out,'(a,2es13.3)')' dcp0, dcp1= :',dcp0,dcp1
917 : ! write(std_out,'(a,2es13.3)')' deltaz,ddz=',deltaz,ddz
918 : !ENDDEBUG
919 :
920 : ! Determine the coefficient of the third-order polynomial taking z0 as origin
921 : ! 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
922 : ! Define qz=a_3*dz + a_2 and dqz=3*a_3*dz + 2*a_2
923 46 : qz=((cp1-cp0)-dcp0*dz1)/dz1**2
924 46 : dqz=(dcp1-dcp0)/dz1
925 46 : a_3=(dqz-two*qz)/dz1
926 46 : a_2=three*qz-dqz
927 :
928 : ! Compute value and gradient of the chemical potential, at ddz wrt to lower delimiting plane
929 46 : e_chempot=e_chempot+(a_3*ddz**3 + a_2*ddz**2 + dcp0*ddz + cp0)
930 69 : grchempottn(3,iatom)=three*a_3*ddz**2 + two*a_2*ddz + dcp0
931 :
932 : !DEBUG
933 : ! write(std_out,'(a,4es16.6)')' qz,dqz=',qz,dqz
934 : ! write(std_out,'(a,4es16.6)')' cp0,dcp0,a_2,a_3=',cp0,dcp0,a_2,a_3
935 : ! write(std_out,'(a,2es13.3)')' dcp0*ddz + cp0=',dcp0*ddz + cp0
936 : ! write(std_out,'(a,2es13.3)')' a_2*ddz**2=',a_2*ddz**2
937 : ! write(std_out,'(a,2es13.3)')' a_3*ddz**3=',a_3*ddz**3
938 : ! write(std_out,'(a,2es13.3)')' contrib=',a_3*ddz**3 + a_2*ddz**2 + dcp0*ddz + cp0
939 : ! write(std_out,'(a,2es13.3)')' e_chempot=',e_chempot
940 : ! write(std_out,'(a,3es20.10)')' grchempottn=',grchempottn(:,iatom)
941 : !ENDDEBUG
942 :
943 : end do
944 :
945 : !DEBUG
946 : !write(std_out,'(a)')' spatialchempot : exit '
947 : !write(std_out,'(a,es16.6)') ' e_chempot=',e_chempot
948 : !ENDDEBUG
949 :
950 23 : end subroutine spatialchempot
951 : !!***
952 :
953 : !!****f* ABINIT/ionion_realspace
954 : !!
955 : !! NAME
956 : !! ionion_realspace
957 : !!
958 : !! FUNCTION
959 : !! Compute the ion/ion interaction energies and forces in real space
960 : !! case. Use ewald() instead if computations are done in reciprocal
961 : !! space since it also includes the correction for the shift done in
962 : !! potentials calculations and includes replica interactions.
963 : !!
964 : !! INPUTS
965 : !! dtset <type(dataset_type)>=all input variables in this dataset
966 : !! rmet(3,3)=metric tensor in real space (bohr^2)
967 : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
968 : !! zion(ntypat)=charge on each type of atom (real number)
969 : !!
970 : !! OUTPUT
971 : !! eew=final ion/ion energy in hartrees
972 : !! grewtn(3,natom)=grads of ion/ion wrt xred(3,natom), hartrees.
973 : !!
974 : !! SOURCE
975 :
976 0 : subroutine ionion_realSpace(dtset, eew, grewtn, rprimd, xred, zion)
977 :
978 : !Arguments ------------------------------------
979 : !scalars
980 : real(dp),intent(out) :: eew
981 : type(dataset_type),intent(in) :: dtset
982 : !arrays
983 : real(dp),intent(in) :: rprimd(3,3),zion(dtset%ntypat)
984 : real(dp),intent(in) :: xred(3,dtset%natom)
985 : real(dp),intent(out) :: grewtn(3,dtset%natom)
986 :
987 : !Local variables-------------------------------
988 : !scalars
989 : integer :: ia1,ia2,iatom,igeo
990 : real(dp) :: r
991 : !arrays
992 0 : real(dp),allocatable :: grew_cart(:,:),xcart(:,:)
993 :
994 : ! *************************************************************************
995 :
996 : !Store xcart for each atom
997 0 : ABI_MALLOC(xcart,(3, dtset%natom))
998 0 : call xred2xcart(dtset%natom, rprimd, xcart, xred)
999 :
1000 : !Summing the interaction between ions.
1001 0 : eew = 0._dp
1002 0 : do ia1 = 1, dtset%natom, 1
1003 0 : do ia2 = ia1 + 1, dtset%natom, 1
1004 : r = sqrt((xcart(1, ia1) - xcart(1, ia2)) ** 2 + &
1005 : & (xcart(2, ia1) - xcart(2, ia2)) ** 2 + &
1006 0 : & (xcart(3, ia1) - xcart(3, ia2)) ** 2)
1007 0 : eew = eew + zion(dtset%typat(ia1)) * zion(dtset%typat(ia2)) / r
1008 : end do
1009 : end do
1010 :
1011 : !Allocate temporary array to store cartesian gradients.
1012 0 : ABI_MALLOC(grew_cart,(3, dtset%natom))
1013 :
1014 : !Summing the forces for each atom
1015 0 : do ia1 = 1, dtset%natom, 1
1016 0 : grew_cart(:, ia1) = 0._dp
1017 0 : do ia2 = 1, dtset%natom, 1
1018 0 : if (ia1 /= ia2) then
1019 : r = (xcart(1, ia1) - xcart(1, ia2)) ** 2 + &
1020 : & (xcart(2, ia1) - xcart(2, ia2)) ** 2 + &
1021 0 : & (xcart(3, ia1) - xcart(3, ia2)) ** 2
1022 0 : do igeo = 1, 3, 1
1023 : grew_cart(igeo, ia1) = grew_cart(igeo, ia1) - (xcart(igeo, ia1) - xcart(igeo, ia2)) * &
1024 0 : & zion(dtset%typat(ia1)) * zion(dtset%typat(ia2)) / (r ** 1.5_dp)
1025 : end do
1026 : end if
1027 : end do
1028 : end do
1029 :
1030 0 : ABI_FREE(xcart)
1031 :
1032 : !Transform cartesian gradients to reduced gradients.
1033 0 : do iatom = 1, dtset%natom, 1
1034 0 : do igeo = 1, 3, 1
1035 : grewtn(igeo, iatom) = rprimd(1, igeo) * grew_cart(1, iatom) + &
1036 : & rprimd(2, igeo) * grew_cart(2, iatom) + &
1037 0 : & rprimd(3, igeo) * grew_cart(3, iatom)
1038 : end do
1039 : end do
1040 0 : ABI_FREE(grew_cart)
1041 :
1042 0 : end subroutine ionion_realSpace
1043 : !!***
1044 :
1045 : !!****f* ABINIT/ionion_surface
1046 : !!
1047 : !! NAME
1048 : !! ionion_surface
1049 : !!
1050 : !! FUNCTION
1051 : !! Compute the ion/ion interaction energies and forces in real space
1052 : !! case. Use ewald() instead if computations are done in reciprocal
1053 : !! space since it also includes the correction for the shift done in
1054 : !! potentials calculations and includes replica interactions.
1055 : !!
1056 : !! INPUTS
1057 : !! dtset <type(dataset_type)>=all input variables in this dataset
1058 : !! rmet(3,3)=metric tensor in real space (bohr^2)
1059 : !! xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
1060 : !! zion(ntypat)=charge on each type of atom (real number)
1061 : !!
1062 : !! OUTPUT
1063 : !! eew=final ion/ion energy in hartrees
1064 : !! grewtn(3,natom)=grads of ion/ion wrt xred(3,natom), hartrees.
1065 : !!
1066 : !! SOURCE
1067 :
1068 0 : subroutine ionion_surface(dtset, eew, grewtn, me, nproc, rprimd, wvl, wvl_den, xred)
1069 :
1070 : #if defined HAVE_BIGDFT
1071 : use BigDFT_API, only: IonicEnergyandForces
1072 : #endif
1073 :
1074 : !Arguments ------------------------------------
1075 : !scalars
1076 : integer, intent(in) :: me, nproc
1077 : real(dp),intent(out) :: eew
1078 : type(dataset_type),intent(in) :: dtset
1079 : type(wvl_internal_type), intent(in) :: wvl
1080 : type(wvl_denspot_type), intent(inout) :: wvl_den
1081 : !arrays
1082 : real(dp),intent(in) :: rprimd(3,3)
1083 : real(dp),intent(in) :: xred(3,dtset%natom)
1084 : real(dp),intent(out) :: grewtn(3,dtset%natom)
1085 :
1086 : !Local variables-------------------------------
1087 : !scalars
1088 : integer :: dispersion, iatom, igeo
1089 : real(dp) :: psoffset
1090 : !arrays
1091 0 : real(dp),allocatable :: xcart(:,:)
1092 0 : real(dp),pointer :: grew_cart(:,:),fdisp(:,:)
1093 : #if defined HAVE_BIGDFT
1094 : real(dp) :: edisp
1095 : real(dp) :: ewaldstr(6)
1096 : #endif
1097 :
1098 : ! *************************************************************************
1099 :
1100 : !Store xcart for each atom
1101 0 : ABI_MALLOC(xcart,(3, dtset%natom))
1102 0 : call xred2xcart(dtset%natom, rprimd, xcart, xred)
1103 :
1104 0 : nullify(fdisp)
1105 0 : nullify(grew_cart)
1106 0 : dispersion = 0
1107 0 : psoffset = 0._dp
1108 : #if defined HAVE_BIGDFT
1109 : call IonicEnergyandForces(me, nproc, wvl_den%denspot%dpbox,&
1110 : & wvl%atoms, dtset%efield, xcart, &
1111 : & eew, grew_cart, dispersion, edisp, fdisp,&
1112 : & ewaldstr,wvl%Glr%d%n1,wvl%Glr%d%n2,wvl%Glr%d%n3,&
1113 : & wvl_den%denspot%V_ext, wvl_den%denspot%pkernel,psoffset)
1114 :
1115 : if (associated(fdisp)) then
1116 : ABI_FREE(fdisp)
1117 : end if
1118 : #endif
1119 :
1120 0 : ABI_FREE(xcart)
1121 :
1122 : !Transform cartesian gradients to reduced gradients.
1123 0 : do iatom = 1, dtset%natom, 1
1124 0 : do igeo = 1, 3, 1
1125 : grewtn(igeo, iatom) = -rprimd(1, igeo) * grew_cart(1, iatom) - &
1126 : & rprimd(2, igeo) * grew_cart(2, iatom) - &
1127 0 : & rprimd(3, igeo) * grew_cart(3, iatom)
1128 : end do
1129 : end do
1130 : if (associated(grew_cart)) then
1131 : ABI_FREE(grew_cart)
1132 : end if
1133 :
1134 : #if !defined HAVE_BIGDFT
1135 : if (.false.) write(std_out,*) me,nproc,wvl%h(1),wvl_den%symObj
1136 : #endif
1137 :
1138 0 : end subroutine ionion_surface
1139 : !!***
1140 :
1141 : end module m_setvtr
1142 : !!***
|
