Branch data Line data Source code
1 : : !!****m* ABINIT/m_rhotov
2 : : !! NAME
3 : : !! m_rhotov
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 1998-2026 ABINIT group (XG, GMR, MT, EB)
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_rhotov
22 : :
23 : : use defs_basis
24 : : use defs_wvltypes
25 : : use m_errors
26 : : use m_abicore
27 : : use m_abi_mixing
28 : : use m_abi2big
29 : : use m_xmpi
30 : : use m_xomp
31 : : use m_cgtools
32 : : use m_xcdata
33 : : use m_dtset
34 : :
35 : : use defs_abitypes, only : MPI_type
36 : : use m_time, only : timab
37 : : use m_geometry, only : xred2xcart
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_spacepar, only : hartre
42 : : use m_dens, only : constrained_dft_t,mag_penalty,constrained_residual
43 : : use m_rhotoxc, only : rhotoxc
44 : : use m_xchybrid, only : xchybrid_ncpp_cc
45 : : use m_psolver, only : psolver_rhohxc
46 : : use m_wvl_psi, only : wvl_psitohpsi
47 : : use m_pawang, only : pawang_type
48 : : use m_pawrad, only : pawrad_type
49 : : use m_pawrhoij, only : pawrhoij_type
50 : : use m_pawtab, only : pawtab_type
51 : : use m_xc_tb09, only : xc_tb09_update_c
52 : :
53 : : implicit none
54 : :
55 : : private
56 : : !!***
57 : :
58 : : public :: rhotov
59 : : !!***
60 : :
61 : : contains
62 : : !!***
63 : :
64 : : !!****f* ABINIT/rhotov
65 : : !! NAME
66 : : !! rhotov
67 : : !!
68 : : !! FUNCTION
69 : : !! This routine is called to compute, from a given total density
70 : : !! the trial (local) potential and the residual potential.
71 : : !!
72 : : !! INPUTS
73 : : !! [add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy
74 : : !! constrained_dft <type(constrained_dft_t>=data for constrained dft calculations
75 : : !! dtset <type(dataset_type)>=all input variables in this dataset
76 : : !! | spinmagntarget=input variable that governs fixed moment calculation
77 : : !! | natom=number of atoms in cell.
78 : : !! | nspden=number of spin-density components
79 : : !! | ntypat=number of types of atoms in unit cell.
80 : : !! | occopt=option for occupancies
81 : : !! | typat(natom)=type (integer) for each atom
82 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
83 : : !! gsqcut=cutoff on (k+G)^2 (bohr^-2)
84 : : !! mpi_enreg=information about MPI parallelization
85 : : !! nfft=(effective) number of FFT grid points (for this processor)
86 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
87 : : !! nhat(nfft,nspden*usepaw)= -PAW only- compensation density
88 : : !! nhatgr(nfft,nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density
89 : : !! nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise
90 : : !! nkxc=second dimension of the array kxc, see rhotoxc.F90 for a description
91 : : !! n3xccc=dimension of the xccc3d array (0 or nfft).
92 : : !! optene=option for the computation of additional energies
93 : : !! optres=0: the trial potential residual is computed ; the input potential value is kept
94 : : !! 1: the new value of the trial potential is computed in place of the input value
95 : : !! optxc=option to be used for the call to rhotoxc
96 : : !! pawang <type(pawang_type)> =paw angular mesh and related data
97 : : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
98 : : !! pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom)
99 : : !! pawtab(ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
100 : : !! rhog(2,nfft)=array for Fourier transform of electron density
101 : : !! rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
102 : : !! | definition for spin components:
103 : : !! | case of nspden = 2
104 : : !! | rhor(:,1) => rho_up + rho_dwn
105 : : !! | rhor(:,2) => rho_up
106 : : !! | case of nspden = 4
107 : : !! | rhor(:,1) => rho_upup + rho_dwndwn
108 : : !! | rhor(:,2:4) => {m_x,m_y,m_z}
109 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
110 : : !! taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density
111 : : !! ucvol = unit cell volume (Bohr**3)
112 : : !! usepaw= 0 for non paw calculation; =1 for paw calculation
113 : : !! usexcnhat= -PAW only- flag controling use of compensation density in Vxc
114 : : !! usevxctau=1 if if XC functional depends on kinetic energy density
115 : : !! vpsp(nfft)=array for holding local psp
116 : : !! [vxc_hybcomp(nfft,nspden)= compensation xc potential (Hartree) in case of hybrids] Optional output
117 : : !! i.e. difference between the hybrid Vxc at fixed density and the auxiliary Vxc at fixed density
118 : : !! xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
119 : : !! xccctau3d(n3xccc*dtset%usekden)=3D core kinetic energy density for XC core correction (bohr^-3)
120 : : !! ==== if optres==0
121 : : !! vtrial(nfft,nspden)= old value of trial potential
122 : : !!
123 : : !! OUTPUT
124 : : !! energies <type(energies_type)>=all part of total energy.
125 : : !! | e_hartree=Hartree part of total energy (hartree units)
126 : : !! | e_xc=exchange-correlation energy (hartree)
127 : : !! | In case of hybrid compensation algorithm:
128 : : !! | e_hybcomp_v=self-consistent potential compensation term for the exchange-correlation energy (hartree)
129 : : !! ==== if optene==0.or.2
130 : : !! | e_localpsp=local psp energy (hartree)
131 : : !! ==== if optene==1.or.2
132 : : !! | e_xcdc=exchange-correlation double-counting energy (hartree)
133 : : !! grcondft(3,natom)=d(E_constrained_DFT)/d(xred) (hartree)
134 : : !! intgres(nspden,ngrcondft)=integrated residuals from constrained DFT.
135 : : !! They are also Lagrange parameters, or gradients with respect to constraints.
136 : : !! kxc(nfft,nkxc)=exchange-correlation kernel, needed only if optxc==2.
137 : : !! strscondft(6)=constrained DFT contribution to stress tensor (hartree/bohr^3)
138 : : !! strsxc(6)=xc contribution to stress tensor (hartree/bohr^3)
139 : : !! vxc(nfft,nspden)=Vxc(r) (already computed above; gets recomputed below too)
140 : : !! vxcavg=mean of the vxc potential
141 : : !! ==== if optres==0
142 : : !! vresidnew(nfft,nspden)=potential residual
143 : : !! vnew_mean(nspden)=mean of the potential formed from vpsp, vhartr and vxc, might be spin-dependent
144 : : !! vres_mean(nspden)=mean of the potential residual, might be spin-dependent
145 : : !! vres2=square of the norm of the residual
146 : : !! vxctau(nfft,dtset%nspden,4*usevxctau)]=derivative of XC energy density with respect to
147 : : !! kinetic energy density (metaGGA cases)
148 : : !! vtauresid(nfft,nspden*usevxctau)=array for vxctau residue (see vxctau)
149 : : !!
150 : : !! SIDE EFFECTS
151 : : !! Input/Output:
152 : : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
153 : : !! vhartr(nfft)=array for holding Hartree potential
154 : : !! ==== if optres==1
155 : : !! vtrial(nfft,nspden)= new value of trial potential
156 : : !!
157 : : !! NOTES
158 : : !! In case of PAW calculations:
159 : : !! All computations are done on the fine FFT grid.
160 : : !! All variables (nfft,ngfft,mgfft) refer to this fine FFT grid.
161 : : !! All arrays (densities/potentials...) are computed on this fine FFT grid.
162 : : !! ! Developpers have to be careful when introducing others arrays:
163 : : !! they have to be stored on the fine FFT grid.
164 : : !! In case of norm-conserving calculations the FFT grid is the usual FFT grid.
165 : : !!
166 : : !! SOURCE
167 : :
168 [ - + ]: 45879 : subroutine rhotov(constrained_dft,dtset,energies,gprimd,grcondft,gsqcut,intgres,istep,kxc,mpi_enreg,nfft,ngfft,&
169 : 91758 : & nhat,nhatgr,nhatgrdim,nkxc,vresidnew,n3xccc,optene,optres,optxc,&
170 : 91758 : & pawang,pawrad,pawrhoij,pawtab,rhog,rhor,rprimd,strscondft,strsxc,taur,ucvol,usepaw,usexcnhat,usevxctau,&
171 [ - + - + ]: 91758 : & vhartr,vnew_mean,vpsp,vres_mean,vres2,vtrial,vxcavg,vxc,vxctau,vtauresid,wvl,xccc3d,xcctau3d,xred,&
172 [ + - + - : 45879 : & electronpositron,vxc_hybcomp,add_tfw) ! optional arguments
- + ]
173 : :
174 : : !Arguments ------------------------------------
175 : : !scalars
176 : : integer,intent(in) :: n3xccc,nfft,nhatgrdim,nkxc,optene,optres,optxc,usepaw,istep
177 : : integer,intent(in) :: usexcnhat,usevxctau
178 : : logical,intent(in),optional :: add_tfw
179 : : real(dp),intent(in) :: gsqcut,ucvol
180 : : real(dp),intent(out) :: vres2,vxcavg
181 : : type(MPI_type),intent(inout) :: mpi_enreg
182 : : type(constrained_dft_t),intent(inout) :: constrained_dft
183 : : type(dataset_type),intent(in) :: dtset
184 : : type(electronpositron_type),pointer,optional :: electronpositron
185 : : type(energies_type),intent(inout) :: energies
186 : : type(pawang_type),intent(in) :: pawang
187 : : type(wvl_data), intent(inout) :: wvl
188 : : !arrays
189 : : integer,intent(in) :: ngfft(18)
190 : : real(dp),intent(in) :: gprimd(3,3),nhat(nfft,dtset%nspden*usepaw)
191 : : real(dp),intent(in) :: nhatgr(nfft,dtset%nspden,3*nhatgrdim),rhog(2,nfft)
192 : : real(dp),intent(in) :: rprimd(3,3)
193 : : real(dp),intent(in) :: taur(nfft,dtset%nspden*dtset%usekden)
194 : : real(dp),intent(inout) :: rhor(nfft,dtset%nspden),vhartr(nfft),vpsp(nfft)
195 : : real(dp),intent(inout) :: vtrial(nfft,dtset%nspden),vxc(nfft,dtset%nspden)
196 : : real(dp),intent(inout) :: xccc3d(n3xccc),xred(3,dtset%natom)
197 : : real(dp),intent(out) :: grcondft(:,:) ! (3,ngrcondft) ngrcondft=natom when condft is activated
198 : : real(dp),intent(out) :: intgres(:,:) ! (nspden,ngrcondft) ngrcondft=natom when condft is activated
199 : : real(dp),intent(out) :: kxc(nfft,nkxc),strsxc(6),vnew_mean(dtset%nspden)
200 : : real(dp),intent(out) :: strscondft(6)
201 : : real(dp),intent(out) :: vres_mean(dtset%nspden),vresidnew(nfft,dtset%nspden)
202 : : real(dp),intent(inout) :: vtauresid(nfft,dtset%nspden*usevxctau)
203 : : real(dp),intent(out),target :: vxctau(nfft,dtset%nspden,4*usevxctau)
204 : : real(dp),intent(out),optional :: vxc_hybcomp(:,:) ! (nfft,nspden)
205 : : real(dp),intent(out) :: xcctau3d(n3xccc)
206 : : type(pawrhoij_type),intent(in) :: pawrhoij(:)
207 : : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*dtset%usepaw)
208 : : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
209 : :
210 : : !Local variables-------------------------------
211 : : !scalars
212 : : integer :: nk3xc,ifft,ipositron,ispden,nfftot,offset
213 : : integer :: mpi_comm_sphgrid,ixc_current
214 : : integer :: mpicomm,nthreads,nmpi,ierr
215 : : !integer :: ii,jj,kk,ipt,nx,ny,nz !SPr: debug
216 : : !real(dp):: rx,ry,rz !SPr: debug
217 : : real(dp) :: doti,e_xcdc_vxctau
218 : : logical :: add_tfw_,calc_xcdc,non_magnetic_xc,with_vxctau
219 : : logical :: is_hybrid_ncpp,wvlbigdft=.false.
220 : : type(xcdata_type) :: xcdata
221 : : !arrays
222 : 91758 : real(dp) :: evxc,tsec(2),vmean(dtset%nspden),vhspinfield(dtset%nspden)
223 : : real(dp),target :: vxctau_dum(0,0,0)
224 : 45879 : real(dp),allocatable :: rhowk(:,:),v_constr_dft_r(:,:),vnew(:,:),xcart(:,:),tmp(:,:)
225 : : real(dp),pointer :: vxctau_(:,:,:)
226 : : !real(dp),allocatable :: vzeemanHarm(:,:) !SPr: debug Zeeman field q/=0 real space
227 : :
228 : : ! *********************************************************************
229 : :
230 : : DBG_ENTER("COLL")
231 : :
232 : 45879 : call timab(940,1,tsec)
233 : :
234 : : !Check if want to use vxctau (mGGA)
235 : 45879 : with_vxctau = (usevxctau==1)
236 [ + + ]: 45879 : vxctau_ => vxctau_dum ; if (with_vxctau) vxctau_ => vxctau
237 : :
238 : : !Check if we're in hybrid norm conserving pseudopotential with a core correction
239 : : is_hybrid_ncpp=(usepaw==0 .and. n3xccc/=0 .and. &
240 [ + + + + : 45879 : & (dtset%ixc==41.or.dtset%ixc==42.or.libxc_functionals_is_hybrid()))
+ + + + ]
241 : :
242 : : !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed
243 [ - + - - ]: 45879 : wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1)
244 : :
245 : : !mpi communicator for spherical grid
246 : 45879 : mpi_comm_sphgrid=mpi_enreg%comm_fft
247 [ - + ]: 45879 : if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl
248 : :
249 : : !Get size of FFT grid
250 [ + + ]: 183516 : nfftot=PRODUCT(ngfft(1:3))
251 : :
252 [ + - ]: 45879 : ipositron=0;if (present(electronpositron)) ipositron=electronpositron_calctype(electronpositron)
253 [ + - ]: 45879 : add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw
254 : :
255 : : !------Compute Hartree and xc potentials----------------------------------
256 : :
257 : : !allocate vnew here.
258 : : !In wvl: vnew is used at call to wvl_psitohpsi
259 [ + + ]: 45879 : if (optres==0) then
260 [ - + + - : 130312 : ABI_MALLOC(vnew,(nfft,dtset%nspden))
- + + - +
- + - - +
- + ]
261 [ + + + + ]: 113078 : vmean(:)=zero ; vnew_mean(:)=zero
262 : : end if
263 : :
264 [ + + ]: 45879 : if (ipositron/=1) then
265 : : ! if metaGGA, save current value of vxctau potential
266 [ + + + + : 35274283 : if (with_vxctau) vtauresid(:,:)=vxctau(:,:,1)
+ + ]
267 : : ! Compute xc potential (separate up and down if spin-polarized)
268 [ + - + - ]: 45655 : if (dtset%icoulomb == 0 .and. dtset%usewvl == 0) then
269 : :
270 : : ! >>>> Hartree potential
271 : : call hartre(1,gsqcut,dtset%icutcoul,usepaw,mpi_enreg,nfft,ngfft,&
272 : 45655 : &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr)
273 : :
274 : : ! >>>> Exchange-correlation potential
275 : : !Use the proper exchange_correlation energy : either the origin one, or the auxiliary one
276 : 45655 : ixc_current=dtset%ixc
277 [ + + ]: 45655 : if(mod(dtset%fockoptmix,100)==11)ixc_current=dtset%auxc_ixc
278 : 45655 : call xcdata_init(xcdata,dtset=dtset,ixc=ixc_current)
279 [ + + + + ]: 45655 : non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
280 : 45655 : nk3xc=1
281 : :
282 : : ! If we use the XC Tran-Blaha 2009 (modified BJ) functional, update the c value
283 [ + + ]: 45655 : if (dtset%xc_tb09_c>99._dp) then
284 : : call xc_tb09_update_c(dtset%intxc,dtset%ixc,mpi_enreg,dtset%natom, &
285 : : & nfft,ngfft,nhat,usepaw,nhatgr,nhatgrdim,dtset%nspden,dtset%ntypat,n3xccc, &
286 : : & pawang,pawrad,pawrhoij,pawtab,dtset%pawxcdev,rhor,rprimd,usepaw, &
287 : : & xccc3d,dtset%xc_denpos,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab, &
288 : 45645 : & computation_type='all')
289 : : end if
290 : :
291 : : ! Use the periodic solver to compute Hxc.
292 : 45655 : call timab(941,1,tsec)
293 [ + + ]: 45655 : if (ipositron==0) then
294 [ + + - + ]: 45403 : if(.not.is_hybrid_ncpp .or. mod(dtset%fockoptmix,100)==11)then
295 : : call rhotoxc(energies%e_xc,energies%entropy_xc,kxc,mpi_enreg,nfft,ngfft,&
296 : : & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,optxc,&
297 : : & rhor,rprimd,usexcnhat,vxc,vxcavg,xccc3d,xcdata,strsxc=strsxc,&
298 : 45073 : & taur=taur,vhartr=vhartr,vxctau=vxctau_,add_tfw=add_tfw_,xcctau3d=xcctau3d)
299 [ + + ]: 45073 : if(mod(dtset%fockoptmix,100)==11)then
300 : 45 : energies%e_xc=energies%e_xc*dtset%auxc_scal
301 : 45 : energies%entropy_xc=energies%entropy_xc*dtset%auxc_scal
302 [ + + + + ]: 184410 : vxc(:,:)=vxc(:,:)*dtset%auxc_scal
303 : : end if
304 : : else
305 : : call xchybrid_ncpp_cc(dtset,energies%e_xc,energies%entropy_xc,mpi_enreg,nfft,ngfft,n3xccc,rhor,rprimd,&
306 : 330 : & strsxc,vxcavg,xccc3d,vxc=vxc)
307 : : end if
308 : : else
309 : : call rhotoxc(energies%e_xc,energies%entropy_xc,kxc,mpi_enreg,nfft,ngfft,&
310 : : & nhat,usepaw,nhatgr,nhatgrdim,nkxc,nk3xc,non_magnetic_xc,n3xccc,optxc,&
311 : : & rhor,rprimd,usexcnhat,vxc,vxcavg,xccc3d,xcdata,&
312 : : & strsxc=strsxc,taur=taur,vhartr=vhartr,vxctau=vxctau_,add_tfw=add_tfw_,&
313 : 252 : & electronpositron=electronpositron,xcctau3d=xcctau3d)
314 : : end if
315 : :
316 : 45655 : call timab(941,2,tsec)
317 [ # # ]: 0 : elseif (.not. wvlbigdft) then
318 : : ! Use the free boundary solver.
319 : 0 : call timab(943,1,tsec)
320 : : call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, &
321 : : & dtset%icoulomb, dtset%ixc, &
322 : : & mpi_enreg, nfft, &
323 : : & ngfft, nhat,usepaw,&
324 : : & dtset%nscforder, dtset%nspden, n3xccc, rhor,rprimd,&
325 : : & usexcnhat,dtset%usepaw,dtset%usewvl,vhartr, vxc, vxcavg,&
326 : : & wvl%descr,wvl%den,wvl%e,&
327 : 0 : & xccc3d,dtset%xclevel,dtset%xc_denpos)
328 : 0 : call timab(943,2,tsec)
329 : : end if
330 : : ! For icoulomb==0 and usewvl Ehartree is calculated in psolver_rhohxc().
331 : : ! For PAW we recalculate this since nhat was not taken into account
332 : : ! in psolver_rhohxc: E_H= int v_H (n+nhat) dr
333 [ + - - + : 45655 : if(.not. wvlbigdft .and. (dtset%icoulomb==0 .or. dtset%usepaw==1 ) ) then
- - ]
334 : 45655 : call timab(942,1,tsec)
335 : 45655 : call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
336 : 45655 : energies%e_hartree=half*energies%e_hartree
337 : 45655 : call timab(942,2,tsec)
338 : : end if
339 : : else
340 : 224 : call timab(944,1,tsec)
341 : 224 : energies%e_hartree=zero;energies%e_xc=zero;energies%entropy_xc=zero
342 : : call rhohxcpositron(electronpositron,gprimd,kxc,mpi_enreg,nfft,ngfft,nhat,nkxc,dtset%nspden,n3xccc,&
343 : 224 : & dtset%paral_kgb,rhor,strsxc,ucvol,usexcnhat,usepaw,vhartr,vxc,vxcavg,xccc3d,dtset%xc_denpos)
344 : 224 : call timab(944,2,tsec)
345 : : end if
346 : :
347 : 45879 : call timab(945,1,tsec)
348 [ + + ]: 45879 : if (ipositron/=0) then
349 : : call dotprod_vn(1,rhor,electronpositron%e_hartree,doti,&
350 : 476 : & nfft,nfftot,1,1,electronpositron%vha_ep,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
351 [ + + ]: 59085467 : vhartr=vhartr+electronpositron%vha_ep
352 : : end if
353 : :
354 : : !------Compute parts of total energy depending on potentials--------
355 : :
356 [ + - + - ]: 45879 : if ( (optene==0.or.optene==2 ).and. .not. wvlbigdft) then
357 : : ! Compute local psp energy energies%e_localpsp
358 : : call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfft,nfftot,1,1,vpsp,ucvol,&
359 : 45879 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
360 : : end if
361 : :
362 [ + + ]: 45879 : if(mod(dtset%fockoptmix,100)==11)then
363 [ + - ]: 45 : if (.not. wvlbigdft) then
364 : : ! Compute second compensation energy for hybrid functionals
365 : : call dotprod_vn(1,rhor,energies%e_hybcomp_v,doti,nfft,nfftot,1,1,vxc_hybcomp,ucvol,&
366 [ - + - - : 45 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
- - - - -
+ ]
367 : : end if
368 : : end if
369 : :
370 : 45879 : calc_xcdc=.false.
371 [ + + ]: 45879 : if (optene==1.or.optene==2) calc_xcdc=.true.
372 [ + - - - ]: 45879 : if (dtset%usewvl==1.and.dtset%nnsclo>0) calc_xcdc=.true.
373 [ - + ]: 45879 : if (wvlbigdft) calc_xcdc=.false.
374 [ + + ]: 45879 : if (dtset%usefock==1) calc_xcdc=.true.
375 : :
376 [ + + ]: 44529 : if (calc_xcdc) then
377 : :
378 : : ! Compute double-counting XC energy energies%e_xcdc
379 [ + + ]: 15286 : if (ipositron/=1) then
380 [ + + + + ]: 15062 : if (usepaw==0.or.usexcnhat/=0) then
381 : : call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol,&
382 : 6566 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
383 : : else
384 [ - + + - : 33984 : ABI_MALLOC(rhowk,(nfft,dtset%nspden))
- + + - +
- + - - +
- + ]
385 [ - + + + : 375451975 : rhowk=rhor-nhat
+ + ]
386 : : call dotprod_vn(1,rhowk,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol,&
387 : 8496 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
388 [ - + ]: 8496 : ABI_FREE(rhowk)
389 : : end if
390 [ + + ]: 15062 : if (with_vxctau)then
391 : : call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfft,nfftot,dtset%nspden,1,vxctau(:,:,1),&
392 : 492 : & ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
393 : 492 : energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau
394 : : end if
395 [ + + ]: 15062 : if (ipositron==2) energies%e_xcdc=energies%e_xcdc-electronpositron%e_xcdc
396 : : else
397 : 224 : energies%e_xcdc=zero
398 : : end if
399 : :
400 : : end if
401 : :
402 : : !------Produce residual vector and square norm of it-------------
403 : : !(only if requested ; if optres==0)
404 : :
405 : : !Set up array for hspinfield
406 : : !EB vhspinfield(:) = factor*( Hz, Hx+iHy; Hx-iHy, -Hz)
407 : : !EB factor = -g/2 * mu_B * mu_0 = -1/2*B in a.u.
408 : : !EB <-- vhspinfield might have to be allocated correctly --> to be checked
409 : : ! vhspinfield = 1/2 ( -B_z, -B_x + iB_y ; -B_x - iB_y , B_z)
410 [ + + ]: 104096 : vhspinfield(:) = zero
411 : : ! ABI_MALLOC(vzeemanHarm,(nfft,dtset%nspden)) ! SPr: debug stuff
412 : : ! vzeemanHarm(:,:) = zero !
413 [ + + + + : 183434 : if (any(abs(dtset%hspinfield(:))>tol8)) then
+ + ]
414 [ + + ]: 58 : if(dtset%nspden==2)then
415 : : ! EB The collinear case has to be checked :
416 : : ! EB Is it vhspinfield(1) or (2) that has to be added here? to be checked in setvtr and energy as well
417 : : ! SPr: the density components are: rhor(1) => n_upup + n_dwndwn
418 : : ! rhor(2) => n_upup
419 : : ! the convention for the potential components is different:
420 : : ! v(1) => v_upup
421 : : ! v(2) => v_dndn
422 : : ! verified by comparing collinear and non-collinear calculations
423 : :
424 : 12 : vhspinfield(1) =-half*dtset%hspinfield(3) ! v_upup
425 : 12 : vhspinfield(2) = half*dtset%hspinfield(3) ! v_dndn
426 : :
427 : : !vhspinfield(1) = zero ! v_upup
428 : : !vhspinfield(2) = zero ! v_dndn
429 : :
430 : : !nx=ngfft(1); ny=ngfft(2); nz=ngfft(3)
431 : : !do kk=0,nz-1
432 : : ! do jj=0,ny-1
433 : : ! do ii=0,nx-1
434 : : ! ipt=1+ii+nx*(jj+ny*kk)
435 : : ! !rx=(dble(ii)/nx)*rprimd(1,1)+(dble(jj)/ny)*rprimd(1,2)+(dble(kk)/nz)*rprimd(1,3)
436 : : ! !ry=(dble(ii)/nx)*rprimd(2,1)+(dble(jj)/ny)*rprimd(2,2)+(dble(kk)/nz)*rprimd(2,3)
437 : : ! !rz=(dble(ii)/nx)*rprimd(3,1)+(dble(jj)/ny)*rprimd(3,2)+(dble(kk)/nz)*rprimd(3,3)
438 : : ! vzeemanHarm(ipt,1)= -half*dtset%hspinfield(3)*cos(2*PI*(dble(ii)/dble(nx)))
439 : : ! vzeemanHarm(ipt,2)= half*dtset%hspinfield(3)*cos(2*PI*(dble(ii)/dble(nx)))
440 : : ! end do
441 : : ! end do
442 : : !end do
443 : :
444 [ + - ]: 46 : else if(dtset%nspden==4)then
445 : :
446 : 46 : vhspinfield(1)=-half*dtset%hspinfield(3) ! v_upup
447 : 46 : vhspinfield(2)= half*dtset%hspinfield(3) ! v_dndn
448 : 46 : vhspinfield(3)=-half*dtset%hspinfield(1) ! Re(v_updn)
449 : 46 : vhspinfield(4)= half*dtset%hspinfield(2) ! Im(v_updn)
450 : :
451 : : !vhspinfield(1)=0.0
452 : : !vhspinfield(2)=0.0
453 : : !vhspinfield(3)=0.0
454 : : !vhspinfield(4)=0.0
455 : :
456 : : !nx=ngfft(1); ny=ngfft(2); nz=ngfft(3)
457 : : !do kk=0,nz-1
458 : : ! do jj=0,ny-1
459 : : ! do ii=0,nx-1
460 : : ! ipt=1+ii+nx*(jj+ny*kk)
461 : : ! !rx=(dble(ii)/nx)*rprimd(1,1)+(dble(jj)/ny)*rprimd(1,2)+(dble(kk)/nz)*rprimd(1,3)
462 : : ! !ry=(dble(ii)/nx)*rprimd(2,1)+(dble(jj)/ny)*rprimd(2,2)+(dble(kk)/nz)*rprimd(2,3)
463 : : ! !rz=(dble(ii)/nx)*rprimd(3,1)+(dble(jj)/ny)*rprimd(3,2)+(dble(kk)/nz)*rprimd(3,3)
464 : : ! vzeemanHarm(ipt,1)= -half*dtset%hspinfield(3)*cos(2*PI*(dble(ii)/dble(nx)))
465 : : ! vzeemanHarm(ipt,2)= half*dtset%hspinfield(3)*cos(2*PI*(dble(ii)/dble(nx)))
466 : : ! vzeemanHarm(ipt,3)= -half*dtset%hspinfield(1)*cos(2*PI*(dble(ii)/dble(nx)))
467 : : ! vzeemanHarm(ipt,4)= half*dtset%hspinfield(2)*cos(2*PI*(dble(ii)/dble(nx)))
468 : : ! end do
469 : : ! end do
470 : : !end do
471 : :
472 : : end if
473 : : end if
474 : :
475 : : !Compute the constrained potential for the magnetic moments
476 [ - + + - : 183516 : ABI_MALLOC(v_constr_dft_r, (nfft,dtset%nspden))
- + + - +
- + - - +
- + ]
477 [ + + + + ]: 1394711871 : v_constr_dft_r = zero
478 [ + + ]: 45879 : if (dtset%magconon==1.or.dtset%magconon==2) then
479 : 104 : call mag_penalty(constrained_dft,mpi_enreg,rhor,v_constr_dft_r,xred,dtset%qgbt,dtset%use_gbt)
480 : : end if
481 : :
482 [ + + ]: 45879 : if (optres==0) then
483 : :
484 : :
485 : : ! ------ Compute potential residual -------------
486 : :
487 [ + - ]: 32578 : if (.not. wvlbigdft) then
488 : : !$OMP PARALLEL DO COLLAPSE(2)
489 [ + + ]: 69866 : do ispden=1,min(dtset%nspden,2)
490 [ + + ]: 611322374 : do ifft=1,nfft
491 : 611252508 : vnew(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)+vhspinfield(ispden)+v_constr_dft_r(ifft,ispden)
492 : : !vnew(ifft,ispden)=vnew(ifft,ispden)+vzeemanHarm(ifft,ispden)
493 [ + + ]: 611252508 : if(mod(dtset%fockoptmix,100)==11)vnew(ifft,ispden)=vnew(ifft,ispden)+vxc_hybcomp(ifft,ispden)
494 : 611289796 : vresidnew(ifft,ispden)=vnew(ifft,ispden)-vtrial(ifft,ispden)
495 : : end do
496 : : end do
497 [ + + ]: 32578 : if(dtset%nspden==4)then
498 : : !$OMP PARALLEL DO COLLAPSE(2)
499 [ + + ]: 4443 : do ispden=3,4
500 [ + + ]: 36690443 : do ifft=1,nfft
501 : 36686000 : vnew(ifft,ispden)=vxc(ifft,ispden)+vhspinfield(ispden)+v_constr_dft_r(ifft,ispden)
502 : : !vnew(ifft,ispden)=vnew(ifft,ispden)+vzeemanHarm(ifft,ispden)
503 [ - + ]: 36686000 : if(mod(dtset%fockoptmix,100)==11)vnew(ifft,ispden)=vnew(ifft,ispden)+vxc_hybcomp(ifft,ispden)
504 : 36688962 : vresidnew(ifft,ispden)=vnew(ifft,ispden)-vtrial(ifft,ispden)
505 : : end do
506 : : end do
507 : : end if
508 : :
509 : : !If constrained_dft, must take into account the constraints, and recompute the residual and the new potential
510 [ + + + + : 76432 : if( any(dtset%constraint_kind(:)/=0))then
+ + ]
511 : : call constrained_residual(constrained_dft,energies%e_constrained_dft,&
512 : 1304 : & grcondft,intgres,mpi_enreg,rhor,strscondft,vresidnew,xred,dtset%qgbt,dtset%use_gbt)
513 [ + + + + ]: 61818215 : vnew(:,1:dtset%nspden)=vtrial(:,1:dtset%nspden)+vresidnew(:,1:dtset%nspden)
514 : : endif
515 : :
516 : : ! /!\ ---- DO NOT CHANGE THESE LINES WITHOUT CORE DEVELOPERS PERMISSION ---- /!\
517 : : ! LB-03/2026:
518 : : ! A noise can accumulate in nvresid after each SCF cycle,
519 : : ! resulting in different densities/potentials for different MPI processes.
520 : : ! This has been observed using threads, but could happen in other contexts.
521 : : ! This slowly worsens the SCF cycle, leading to wrong results after many iterations.
522 : : ! So here we compute the mean of nvresid over all MPI processes to reduce the noise.
523 : : ! This error is difficult to test as it is observed in long runs only, so BE VERY CAREFUL.
524 : 32578 : mpicomm = mpi_enreg%comm_kptband
525 : 32578 : nmpi = xmpi_comm_size(mpicomm)
526 : 32578 : nthreads = xomp_get_num_threads(open_parallel=.true.)
527 [ - + ]: 32578 : if (nmpi>1.and.nthreads>1) then
528 [ # # # # : 0 : ABI_MALLOC(tmp,(size(vresidnew,1),size(vresidnew,2)))
# # # # #
# # # # #
# # ]
529 [ # # # # ]: 0 : tmp(:,:) = vresidnew(:,:) / nmpi
530 : 0 : call xmpi_sum(tmp,mpicomm,ierr)
531 [ # # ]: 0 : if (ierr/=0) then
532 : 0 : ABI_ERROR("Error in mpi sum (vresidnew)")
533 : : end if
534 [ # # # # ]: 0 : vresidnew(:,:) = tmp(:,:)
535 [ # # ]: 0 : ABI_FREE(tmp)
536 : : end if
537 : : ! /!\--------------------/!\
538 : :
539 : 32578 : offset = 0
540 : :
541 [ - + - - : 32578 : if (dtset%iscf==0) vtrial=vnew
- - ]
542 : :
543 : : ! Pass vtrial to BigDFT object
544 [ - + ]: 32578 : if(dtset%usewvl==1) then
545 : 0 : call wvl_vtrial_abi2big(1,vnew,wvl%den)
546 : : ! call wvl_vtrial_abi2big(1,vtrial,wvl%den)
547 : : end if
548 : :
549 : : else
550 : : ! Compute with covering comms the different part of the potential.
551 : : ! only for wvlbigdft
552 [ # # # # : 0 : ABI_MALLOC(xcart,(3, dtset%natom))
# # # # #
# # # ]
553 : 0 : call xred2xcart(dtset%natom, rprimd, xcart, xred)
554 : : call wvl_psitohpsi(dtset%diemix,energies%e_exactX, energies%e_xc, energies%e_hartree, &
555 : : & energies%e_kinetic, energies%e_localpsp, energies%e_nlpsp_vfock, energies%e_sicdc, &
556 : : & istep + 1, 1, dtset%iscf, mpi_enreg%me_wvl, dtset%natom, dtset%nfft,&
557 : : & mpi_enreg%nproc_wvl, dtset%nspden, &
558 : : & vres2, .true., energies%e_xcdc, wvl,&
559 : 0 : & wvlbigdft, xcart, strsxc,vtrial=vnew,vxc=vxc)
560 [ # # ]: 0 : ABI_FREE(xcart)
561 : :
562 [ # # # # ]: 0 : vresidnew = vnew - vtrial
563 [ # # # # ]: 0 : vtrial = vnew
564 : :
565 : : call mean_fftr(vxc, vmean(1:1), nfft, nfftot, dtset%nspden,&
566 : 0 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
567 : 0 : vxcavg = vmean(1)
568 : 0 : offset = 0
569 : : end if
570 : :
571 : : ! Compute mean values of potential and residual
572 : : call mean_fftr(vnew(1+offset, 1),vnew_mean,nfft,nfftot,dtset%nspden,&
573 : 32578 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
574 : : call mean_fftr(vresidnew(1+offset, 1),vmean,nfft,nfftot,dtset%nspden,&
575 : 32578 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
576 : :
577 [ - + ]: 32578 : ABI_FREE(vnew)
578 : :
579 : : ! Subtract the mean of the residual
580 : : ! Must take into account fixed occupation number in case of spin-polarized
581 [ + + ]: 72828 : do ispden=1,dtset%nspden
582 [ + + + + : 72828 : if (dtset%nspden==2.and.dtset%occopt>=3.and. abs(dtset%spinmagntarget+99.99_dp)<1.0d-10)then
+ + ]
583 : 1398 : vres_mean(ispden)=(vmean(1)+vmean(2))*half
584 : : else
585 : 38852 : vres_mean(ispden)=vmean(ispden)
586 : : end if
587 : : end do
588 : :
589 : : !$OMP PARALLEL DO COLLAPSE(2)
590 [ + + ]: 72828 : do ispden=1,dtset%nspden
591 [ + + ]: 648011336 : do ifft=1,nfft
592 : 647978758 : vresidnew(ifft,ispden)=vresidnew(ifft,ispden)-vres_mean(ispden)
593 : : end do
594 : : end do
595 : :
596 : : ! Compute square norm vres2 of potential residual vresid
597 : 32578 : call sqnorm_v(1,nfft,vres2,dtset%nspden,optres,vresidnew(1+offset, 1),mpi_comm_sphgrid=mpi_comm_sphgrid)
598 : :
599 : : ! Now take care of Vxctau residual (metaGGA)
600 [ + + ]: 32578 : if (with_vxctau) then
601 [ + - + + : 8447065 : if (ipositron/=1) vtauresid(:,:)=vxctau(:,:,1)-vtauresid(:,:)
+ + ]
602 [ - + - - : 149 : if (ipositron==1) vtauresid(:,:)=zero
- - ]
603 : : end if
604 : :
605 : : else ! optres/=0
606 : :
607 : : ! ------Produce new value of trial potential-------------
608 : :
609 [ + - ]: 13301 : if (.not. wvlbigdft) then
610 : : !$OMP PARALLEL DO COLLAPSE(2)
611 [ + + ]: 29286 : do ispden=1,min(dtset%nspden,2)
612 [ + + ]: 720953421 : do ifft=1,nfft
613 : 720924135 : vtrial(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)+vhspinfield(ispden)+v_constr_dft_r(ifft,ispden)
614 : : !vtrial(ifft,ispden)=vtrial(ifft,ispden)+vzeemanHarm(ifft,ispden)
615 [ - + ]: 720940120 : if(mod(dtset%fockoptmix,100)==11)vtrial(ifft,ispden)=vtrial(ifft,ispden)+vxc_hybcomp(ifft,ispden)
616 : : end do
617 : : end do
618 [ + + ]: 13301 : if(dtset%nspden==4) then
619 : : !$OMP PARALLEL DO
620 [ + + ]: 12873557 : do ifft=1,nfft
621 [ + + ]: 38617698 : vtrial(ifft,3:4)=vxc(ifft,3:4)+vhspinfield(3:4)+v_constr_dft_r(ifft,3:4)
622 : : !vtrial(ifft,3:4)=vtrial(ifft,3:4)+vzeemanHarm(ifft,3:4)
623 [ - + - - ]: 12873557 : if(mod(dtset%fockoptmix,100)==11)vtrial(ifft,3:4)=vtrial(ifft,3:4)+vxc_hybcomp(ifft,3:4)
624 : : end do
625 : : end if
626 : : ! Pass vtrial to BigDFT object
627 [ - + ]: 13301 : if(dtset%usewvl==1) then
628 : 0 : call wvl_vtrial_abi2big(1,vtrial,wvl%den)
629 : : end if
630 : : else
631 : : ! Compute with covering comms the different part of the potential.
632 [ # # # # : 0 : ABI_MALLOC(xcart,(3, dtset%natom))
# # # # #
# # # ]
633 : 0 : call xred2xcart(dtset%natom, rprimd, xcart, xred)
634 : : call wvl_psitohpsi(dtset%diemix,energies%e_exactX, energies%e_xc, energies%e_hartree, &
635 : : & energies%e_kinetic, energies%e_localpsp, energies%e_nlpsp_vfock, energies%e_sicdc, &
636 : : & istep + 1, 1, dtset%iscf, mpi_enreg%me_wvl, &
637 : : & dtset%natom, dtset%nfft, mpi_enreg%nproc_wvl,&
638 : : & dtset%nspden,vres2, .true.,energies%e_xcdc, wvl,&
639 : 0 : & wvlbigdft, xcart, strsxc, vtrial, vxc)
640 [ # # ]: 0 : ABI_FREE(xcart)
641 : : ! Compute vxcavg
642 : : call mean_fftr(vxc, vmean(1:1), nfft, nfftot, dtset%nspden,&
643 : 0 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
644 : 0 : vxcavg = vmean(1)
645 : : end if
646 : :
647 : : end if
648 : :
649 [ - + ]: 45879 : ABI_FREE(v_constr_dft_r)
650 : : !ABI_FREE(vzeemanHarm) !SPr: debug for q/=0 magnetic field
651 : :
652 : 45879 : call timab(945,2,tsec)
653 : 45879 : call timab(940,2,tsec)
654 : :
655 : : DBG_EXIT("COLL")
656 : :
657 [ - + - + ]: 45879 : end subroutine rhotov
658 : : !!***
659 : :
660 : : end module m_rhotov
661 : : !!***
|
