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 45405 : subroutine rhotov(constrained_dft,dtset,energies,gprimd,grcondft,gsqcut,intgres,istep,kxc,mpi_enreg,nfft,ngfft,&
169 90810 : & nhat,nhatgr,nhatgrdim,nkxc,vresidnew,n3xccc,optene,optres,optxc,&
170 90810 : & pawang,pawrad,pawrhoij,pawtab,rhog,rhor,rprimd,strscondft,strsxc,taur,ucvol,usepaw,usexcnhat,usevxctau,&
171 90810 : & vhartr,vnew_mean,vpsp,vres_mean,vres2,vtrial,vxcavg,vxc,vxctau,vtauresid,wvl,xccc3d,xcctau3d,xred,&
172 45405 : & 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 90810 : real(dp) :: evxc,tsec(2),vmean(dtset%nspden),vhspinfield(dtset%nspden)
223 : real(dp),target :: vxctau_dum(0,0,0)
224 45405 : 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 45405 : call timab(940,1,tsec)
233 :
234 : !Check if want to use vxctau (mGGA)
235 45405 : with_vxctau = (usevxctau==1)
236 45405 : 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 45405 : & (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 45405 : wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1)
244 :
245 : !mpi communicator for spherical grid
246 45405 : mpi_comm_sphgrid=mpi_enreg%comm_fft
247 45405 : if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl
248 :
249 : !Get size of FFT grid
250 181620 : nfftot=PRODUCT(ngfft(1:3))
251 :
252 45405 : ipositron=0;if (present(electronpositron)) ipositron=electronpositron_calctype(electronpositron)
253 45405 : 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 45405 : if (optres==0) then
260 129932 : ABI_MALLOC(vnew,(nfft,dtset%nspden))
261 112421 : vmean(:)=zero ; vnew_mean(:)=zero
262 : end if
263 :
264 45405 : if (ipositron/=1) then
265 : ! if metaGGA, save current value of vxctau potential
266 35273809 : if (with_vxctau) vtauresid(:,:)=vxctau(:,:,1)
267 : ! Compute xc potential (separate up and down if spin-polarized)
268 45181 : 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 45181 : &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 45181 : ixc_current=dtset%ixc
277 45181 : if(mod(dtset%fockoptmix,100)==11)ixc_current=dtset%auxc_ixc
278 45181 : call xcdata_init(xcdata,dtset=dtset,ixc=ixc_current)
279 45181 : non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
280 45181 : nk3xc=1
281 :
282 : ! If we use the XC Tran-Blaha 2009 (modified BJ) functional, update the c value
283 45181 : 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 45171 : & computation_type='all')
289 : end if
290 :
291 : ! Use the periodic solver to compute Hxc.
292 45181 : call timab(941,1,tsec)
293 45181 : if (ipositron==0) then
294 44929 : 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 44599 : & taur=taur,vhartr=vhartr,vxctau=vxctau_,add_tfw=add_tfw_,xcctau3d=xcctau3d)
299 44599 : 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 45181 : 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 45181 : if(.not. wvlbigdft .and. (dtset%icoulomb==0 .or. dtset%usepaw==1 ) ) then
334 45181 : call timab(942,1,tsec)
335 45181 : call dotprod_vn(1,rhor,energies%e_hartree,doti,nfft,nfftot,1,1,vhartr,ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
336 45181 : energies%e_hartree=half*energies%e_hartree
337 45181 : 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 45405 : call timab(945,1,tsec)
348 45405 : 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 45405 : 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 45405 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
360 : end if
361 :
362 45405 : 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 45405 : calc_xcdc=.false.
371 45405 : if (optene==1.or.optene==2) calc_xcdc=.true.
372 45405 : if (dtset%usewvl==1.and.dtset%nnsclo>0) calc_xcdc=.true.
373 45405 : if (wvlbigdft) calc_xcdc=.false.
374 45405 : if (dtset%usefock==1) calc_xcdc=.true.
375 :
376 44055 : if (calc_xcdc) then
377 :
378 : ! Compute double-counting XC energy energies%e_xcdc
379 14907 : if (ipositron/=1) then
380 14683 : 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 6527 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
383 : else
384 32624 : ABI_MALLOC(rhowk,(nfft,dtset%nspden))
385 323867083 : rhowk=rhor-nhat
386 : call dotprod_vn(1,rhowk,energies%e_xcdc,doti,nfft,nfftot,dtset%nspden,1,vxc,ucvol,&
387 8156 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
388 8156 : ABI_FREE(rhowk)
389 : end if
390 14683 : 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 14683 : 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 102923 : vhspinfield(:) = zero
411 : ! ABI_MALLOC(vzeemanHarm,(nfft,dtset%nspden)) ! SPr: debug stuff
412 : ! vzeemanHarm(:,:) = zero !
413 181538 : 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 181620 : ABI_MALLOC(v_constr_dft_r, (nfft,dtset%nspden))
477 1329396891 : v_constr_dft_r = zero
478 45405 : 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 45405 : if (optres==0) then
483 :
484 :
485 : ! ------ Compute potential residual -------------
486 :
487 32483 : if (.not. wvlbigdft) then
488 : !$OMP PARALLEL DO COLLAPSE(2)
489 69614 : do ispden=1,min(dtset%nspden,2)
490 608659055 : do ifft=1,nfft
491 608589441 : 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 608589441 : if(mod(dtset%fockoptmix,100)==11)vnew(ifft,ispden)=vnew(ifft,ispden)+vxc_hybcomp(ifft,ispden)
494 608626572 : vresidnew(ifft,ispden)=vnew(ifft,ispden)-vtrial(ifft,ispden)
495 : end do
496 : end do
497 32483 : if(dtset%nspden==4)then
498 : !$OMP PARALLEL DO COLLAPSE(2)
499 4257 : do ispden=3,4
500 34249565 : do ifft=1,nfft
501 34245308 : 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 34245308 : if(mod(dtset%fockoptmix,100)==11)vnew(ifft,ispden)=vnew(ifft,ispden)+vxc_hybcomp(ifft,ispden)
504 34248146 : 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 76233 : 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 32483 : mpicomm = mpi_enreg%comm_kptband
525 32483 : nmpi = xmpi_comm_size(mpicomm)
526 32483 : nthreads = xomp_get_num_threads(open_parallel=.true.)
527 32483 : 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 32483 : offset = 0
540 :
541 32483 : if (dtset%iscf==0) vtrial=vnew
542 :
543 : ! Pass vtrial to BigDFT object
544 32483 : 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 32483 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
574 : call mean_fftr(vresidnew(1+offset, 1),vmean,nfft,nfftot,dtset%nspden,&
575 32483 : & mpi_comm_sphgrid=mpi_comm_sphgrid)
576 :
577 32483 : 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 72452 : do ispden=1,dtset%nspden
582 72452 : 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 38571 : vres_mean(ispden)=vmean(ispden)
586 : end if
587 : end do
588 :
589 : !$OMP PARALLEL DO COLLAPSE(2)
590 72452 : do ispden=1,dtset%nspden
591 642907201 : do ifft=1,nfft
592 642874718 : 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 32483 : 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 32483 : 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 12922 : if (.not. wvlbigdft) then
610 : !$OMP PARALLEL DO COLLAPSE(2)
611 28489 : do ispden=1,min(dtset%nspden,2)
612 660742576 : do ifft=1,nfft
613 660714087 : 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 660729654 : if(mod(dtset%fockoptmix,100)==11)vtrial(ifft,ispden)=vtrial(ifft,ispden)+vxc_hybcomp(ifft,ispden)
616 : end do
617 : end do
618 12922 : 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 12922 : 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 45405 : ABI_FREE(v_constr_dft_r)
650 : !ABI_FREE(vzeemanHarm) !SPr: debug for q/=0 magnetic field
651 :
652 45405 : call timab(945,2,tsec)
653 45405 : call timab(940,2,tsec)
654 :
655 : DBG_EXIT("COLL")
656 :
657 45405 : end subroutine rhotov
658 : !!***
659 :
660 : end module m_rhotov
661 : !!***
|
