Line data Source code
1 : !!****m* ABINIT/m_dft_energy
2 : !! NAME
3 : !! m_dft_energy
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, AR, MB, MT, EB)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_dft_energy
23 :
24 : use defs_basis
25 : use defs_wvltypes
26 : use m_abicore
27 : use m_hamiltonian
28 : use m_errors
29 : use m_xmpi
30 : use m_xcdata
31 : use m_cgtools
32 : use m_dtset
33 : use m_extfpmd
34 : use m_ompgpu_utils
35 :
36 : use defs_datatypes, only : pseudopotential_type
37 : use defs_abitypes, only : MPI_type
38 : use m_time, only : timab
39 : use m_geometry, only : metric
40 : use m_kg, only : mkkin
41 : use m_energies, only : energies_type
42 : use m_electronpositron, only : electronpositron_type, electronpositron_calctype, rhohxcpositron
43 : use m_bandfft_kpt, only : bandfft_kpt, bandfft_kpt_type, prep_bandfft_tabs, &
44 : bandfft_kpt_savetabs, bandfft_kpt_restoretabs
45 : use m_pawang, only : pawang_type
46 : use m_pawtab, only : pawtab_type
47 : use m_paw_ij, only : paw_ij_type
48 : use m_pawfgrtab, only : pawfgrtab_type
49 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_init_unpacked, &
50 : pawrhoij_mpisum_unpacked, pawrhoij_free_unpacked, pawrhoij_inquire_dim, &
51 : & pawrhoij_symrhoij
52 : use m_pawcprj, only : pawcprj_type,pawcprj_alloc,pawcprj_free,pawcprj_gather_spin
53 : use m_pawfgr, only : pawfgr_type
54 : use m_paw_dmft, only : paw_dmft_type
55 : use m_paw_nhat, only : pawmknhat
56 : use m_paw_mkrho, only : pawmkrho
57 : use m_paw_occupancies, only : pawaccrhoij
58 : use m_rcpaw, only : rcpaw_type
59 : use m_fft, only : fftpac, fourdp
60 : use m_spacepar, only : meanvalue_g, hartre
61 : use m_dens, only : constrained_dft_t,mag_penalty
62 : use m_mkrho, only : mkrho
63 : use m_mkffnl, only : mkffnl
64 : use m_getghc, only : getghc
65 : use m_rhotoxc, only : rhotoxc
66 : use m_mpinfo, only : proc_distrb_cycle
67 : use m_nonlop, only : nonlop
68 : use m_fourier_interpol, only : transgrid
69 : use m_prep_kgb, only : prep_getghc, prep_nonlop
70 : use m_psolver, only : psolver_rhohxc
71 : use m_gemm_nonlop_projectors, only : set_gemm_nonlop_ikpt, gemm_nonlop_use_gemm
72 :
73 : #ifdef HAVE_FC_ISO_C_BINDING
74 : use, intrinsic :: iso_c_binding, only : c_int64_t
75 : #endif
76 :
77 : #if defined HAVE_GPU_CUDA
78 : use m_manage_cuda
79 : #endif
80 :
81 : implicit none
82 :
83 : private
84 : !!***
85 :
86 : public :: energy,entropy
87 : !!***
88 :
89 : contains
90 : !!***
91 :
92 : !!****f* ABINIT/energy
93 : !! NAME
94 : !! energy
95 : !!
96 : !! FUNCTION
97 : !! Compute electronic energy terms
98 : !! energies%e_eigenvalues, ek and enl from arbitrary (orthonormal) provided wf,
99 : !! ehart, enxc, and eei from provided density and potential,
100 : !! energies%e_eigenvalues=Sum of the eigenvalues - Band energy (Hartree)
101 : !! energies%e_hspinfield=Spin magnetic energy from applied magnetic field -m.B
102 : !! ek=kinetic energy, ehart=Hartree electron-electron energy,
103 : !! enxc,enxcdc=exchange-correlation energies, eei=local pseudopotential energy,
104 : !! enl=nonlocal pseudopotential energy
105 : !! Also, compute new density from provided wfs, after the evaluation
106 : !! of ehart, enxc, and eei.
107 : !! WARNING XG180913 : At present, Fock energy not computed !
108 : !!
109 : !! NOTE that this routine is callned in m_scfcv_core only when nstep == 0
110 : !!
111 : !! INPUTS
112 : !! [add_tfw]=flag controling the addition of Weiszacker gradient correction to Thomas-Fermi kin energy
113 : !! cg(2,mpw*nspinor*mband*mkmem*nsppol)=<G|Cnk>=Fourier coefficients of wavefunction
114 : !! dtset <type(dataset_type)>=all input variables for this dataset
115 : !! | mband=maximum number of bands
116 : !! | mgfft=maximum size of 1D FFTs
117 : !! | mkmem=number of k points treated by this node.
118 : !! | mpw=maximum dimension for number of planewaves
119 : !! | natom=number of atoms in unit cell
120 : !! | nfft=(effective) number of FFT grid points (for this processor)
121 : !! | nkpt=number of k points
122 : !! | nspden=number of spin-density components
123 : !! | nsppol=1 for unpolarized, 2 for polarized
124 : !! | nspinor=number of spinorial components
125 : !! | nsym=number of symmetry elements in space group (at least 1)
126 : !! | occopt=option for occupancies
127 : !! | tsmear=smearing energy or temperature (if metal)
128 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
129 : !! extfpmd <type(extfpmd_type)>=extended first-principles molecular dynamics type
130 : !! gsqcut=G^2 cutoff from gsqcut=ecut/(2 Pi^2)
131 : !! indsym(4,nsym,natom)=indirect indexing array for atom labels
132 : !! irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
133 : !! kg(3,mpw*mkmem)=work array for coordinates of G vectors in basis
134 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
135 : !! mpi_enreg=information about MPI parallelization
136 : !! my_natom=number of atoms treated by current processor
137 : !! nfftf= -PAW ONLY- number of FFT grid points for the fine grid
138 : !! (nfftf=nfft for norm-conserving potential runs)
139 : !! ngfftf(18)= -PAW ONLY- contain all needed information about 3D FFT for the fine grid
140 : !! (ngfftf=ngfft for norm-conserving potential runs)
141 : !! nhatgr(nfft,nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density
142 : !! nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise
143 : !! npwarr(nkpt)=number of planewaves at each k point, and boundary
144 : !! n3xccc=dimension of the xccc3d array (0 or nfftf).
145 : !! occ(mband*nkpt*nsppol)=occupation numbers of bands (usually 2) at each k point
146 : !! optene=option for the computation of total energy (direct scheme or double-counting scheme)
147 : !! paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
148 : !! pawang <type(pawang_type)>=paw angular mesh and related data
149 : !! pawfgr(natom) <type(pawfgr_type)>=fine grid parameters and related data
150 : !! pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
151 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
152 : !! phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
153 : !! ph1d(2,3*(2*mgfft+1)*natom)=phase information related to structure factor
154 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
155 : !! | mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
156 : !! | ntypat=number of types of atoms in cell
157 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
158 : !! symrec(3,3,nsym)=symmetry operations in reciprocal space
159 : !! usexcnhat= -PAW only- flag controling use of compensation density in Vxc
160 : !! vpsp(nfftf)=local pseudopotential in real space (hartree)
161 : !! wfs <type(wvl_projector_type)>=wavefunctions information for wavelets.
162 : !! wvl <type(wvl_internal_type)>=wavelets internal data
163 : !! xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
164 : !! xcctau3d(n3xccc)=3D core electron kinetic energy density for XC core correction (bohr^-3)
165 : !! xred(3,natom)=reduced coordinates of atoms (dimensionless)
166 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
167 : !!
168 : !! OUTPUT
169 : !! compch_fft=-PAW only- compensation charge inside spheres computed over fine fft grid
170 : !! etotal=total energy (hartree):
171 : !! - computed by direct scheme if optene=0 or 2
172 : !! - computed by double-counting scheme if optene=1 or 3
173 : !! resid(mband*nkpt*nsppol)=residuals for each band over all k points (hartree^2)
174 : !! strsxc(6)=exchange-correlation contribution to stress tensor
175 : !! vhartr(nfftf)=work space to hold Hartree potential in real space (hartree)
176 : !! vtrial(nfftf,nspden)=total local potential (hartree)
177 : !! vxc(nfftf,nspden)=work space to hold Vxc(r) in real space (hartree)
178 : !! [vxctau(nfftf,nspden,4*usevxctau)]=(only for meta-GGA): derivative of XC energy density
179 : !! with respect to kinetic energy density (depsxcdtau). The arrays vxctau contains also
180 : !! the gradient of vxctau (gvxctau) in vxctau(:,:,2:4)
181 : !!
182 : !! SIDE EFFECTS
183 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
184 : !! energies <type(energies_type)>=all part of total energy.
185 : !! | entropy(IN)=entropy due to the occupation number smearing (if metal)
186 : !! | e_chempot(IN)=energy from spatially varying chemical potential (hartree)
187 : !! | e_ewald(IN)=Ewald energy (hartree)
188 : !! | e_vdw_dftd(IN)=VdW DFT-D energy
189 : !! | e_corepsp(IN)=psp core-core energy
190 : !! | paw%epaw(IN)=PAW spherical part energy
191 : !! | paw%epaw_dc(IN)=PAW spherical part double-counting energy
192 : !! | e_eigenvalues(OUT)=Sum of the eigenvalues - Band energy (Hartree)
193 : !! | e_hartree(OUT)=Hartree part of total energy (hartree units)
194 : !! | e_kinetic(OUT)=kinetic energy part of total energy.
195 : !! | e_nlpsp_vfock(OUT)=nonlocal psp + potential Fock ACE part of total energy.
196 : !! | e_xc(OUT)=exchange-correlation energy (hartree)
197 : !! ==== if optene==0, 2 or 3
198 : !! | e_localpsp(OUT)=local psp energy (hartree)
199 : !! ==== if optene==1, 2 or 3
200 : !! | e_xcdc(OUT)=exchange-correlation double-counting energy (hartree)
201 : !! rhog(2,nfftf)=work space for rho(G); save intact on return (? MT 08-12-2008: is that true now ?)
202 : !! rhor(nfftf,nspden)=work space for rho(r); save intact on return (? MT 08-12-2008: is that true now ?)
203 : !! taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic energy density
204 : !! taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density
205 : !! nspinor should not be modified in the call of rdnpw
206 : !! === if psps%usepaw==1 ===
207 : !! nhat(nfftf,nspden*usepaw)= compensation charge density
208 : !! pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
209 : !!
210 : !! NOTES
211 : !! Be careful to the meaning of nfft (size of FFT grids):
212 : !! - In case of norm-conserving calculations the FFT grid is the usual FFT grid.
213 : !! - In case of PAW calculations:
214 : !! Two FFT grids are used; one with nfft points (coarse grid) for
215 : !! the computation of wave functions ; one with nfftf points
216 : !! (fine grid) for the computation of total density.
217 : !!
218 : !! There is a large amount of overhead in the way this routine do the computation of the energy !
219 : !! For example, the density has already been precomputed, so why to compute it again here ??
220 : !!
221 : !! SOURCE
222 :
223 70 : subroutine energy(cg,compch_fft,constrained_dft,dtset,electronpositron,&
224 35 : & energies,eigen,etotal,gsqcut,extfpmd,indsym,irrzon,kg,mcg,mpi_enreg,my_natom,nfftf,ngfftf,nhat,&
225 35 : & nhatgr,nhatgrdim,npwarr,n3xccc,occ,optene,paw_dmft,paw_ij,pawang,pawfgr,&
226 35 : & pawfgrtab,pawrhoij,pawtab,phnons,ph1d,psps,resid,rhog,rhor,rprimd,strsxc,symrec,&
227 35 : & taug,taur,usexcnhat,vhartr,vtrial,vpsp,vxc,wfs,wvl,wvl_den,wvl_e,xccc3d,xred,ylm,&
228 35 : & add_tfw,vxctau,xcctau3d,rcpaw) ! optional argument
229 :
230 : !Arguments ------------------------------------
231 : !scalars
232 : integer,intent(in) :: mcg,my_natom,n3xccc,nfftf,nhatgrdim,optene,usexcnhat
233 : logical,intent(in),optional :: add_tfw
234 : real(dp),intent(in) :: gsqcut
235 : real(dp),intent(out) :: compch_fft,etotal
236 : type(MPI_type),intent(inout) :: mpi_enreg
237 : type(constrained_dft_t),intent(in) :: constrained_dft
238 : type(dataset_type),intent(in) :: dtset
239 : type(electronpositron_type),pointer :: electronpositron
240 : type(energies_type),intent(inout) :: energies
241 : type(extfpmd_type),pointer,intent(inout) :: extfpmd
242 : type(paw_dmft_type), intent(inout) :: paw_dmft
243 : type(pawang_type),intent(in) :: pawang
244 : type(pawfgr_type),intent(in) :: pawfgr
245 : type(pseudopotential_type),intent(in) :: psps
246 : type(wvl_internal_type), intent(in) :: wvl
247 : type(wvl_wf_type),intent(inout) :: wfs
248 : type(wvl_denspot_type), intent(inout) :: wvl_den
249 : type(wvl_energy_terms),intent(inout) ::wvl_e
250 : type(rcpaw_type),pointer,intent(in),optional :: rcpaw
251 : !arrays
252 : !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
253 : integer, intent(in) :: indsym(4,dtset%nsym,dtset%natom)
254 : integer :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)),kg(3,dtset%mpw*dtset%mkmem)
255 : integer, intent(in) :: ngfftf(18),npwarr(dtset%nkpt),symrec(3,3,dtset%nsym)
256 : real(dp), intent(in) :: cg(2,mcg),eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
257 : real(dp), intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol),ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)
258 : real(dp), intent(inout) :: nhat(nfftf,dtset%nspden*psps%usepaw)
259 : real(dp),intent(in) :: nhatgr(nfftf,dtset%nspden,3*nhatgrdim)
260 : !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
261 : real(dp), intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
262 : real(dp), intent(inout) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol)
263 : real(dp), intent(inout) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden)
264 : real(dp), intent(inout) :: taug(2,nfftf*dtset%usekden),taur(nfftf,dtset%nspden*dtset%usekden)
265 : real(dp), intent(out) :: strsxc(6)
266 : real(dp), intent(in) :: rprimd(3,3),vpsp(nfftf),xccc3d(n3xccc),xred(3,dtset%natom)
267 : real(dp), intent(in) :: xcctau3d(nfftf*dtset%usekden)
268 : real(dp), intent(out) :: vhartr(nfftf),vtrial(nfftf,dtset%nspden),vxc(nfftf,dtset%nspden)
269 : real(dp),intent(out),optional,target :: vxctau(:,:,:) ! vxctau(nfftf,dtset%nspden,4*usevxctau)
270 : real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
271 : type(paw_ij_type), intent(in) :: paw_ij(my_natom*psps%usepaw)
272 : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
273 : type(pawrhoij_type),target,intent(inout) :: pawrhoij(my_natom*psps%usepaw)
274 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
275 :
276 : !Local variables-------------------------------
277 : !scalars
278 : integer :: bdtot_index,blocksize,choice,cplex,cplex_rhoij,cpopt,dimffnl
279 : integer :: iband,iband_last,iblock,iblocksize,icg,ider,idir,ierr,ifft,ikg,ikpt,ilm
280 : integer :: ipert,ipositron,iresid,ispden,isppol,istwf_k,itypat,izero
281 : integer :: me_distrb,mpi_comm_sphgrid,my_ikpt,my_nspinor,n1,n2,n3,n4,n5,n6
282 : integer :: nband_k,nblockbd,nfftotf,nkpg,nkxc,nk3xc,nnlout,npw_k,nspden_rhoij,option
283 : integer :: option_rhoij,paw_opt,signs,spaceComm,tim_mkrho,tim_nonlop
284 : logical :: add_tfw_,paral_atom,use_timerev,use_zeromag,with_vxctau
285 : logical :: non_magnetic_xc,wvlbigdft=.false.
286 : real(dp) :: dotr,doti,eeigk,ekk,enlk,evxc,e_xcdc_vxctau,ucvol,ucvol_local,vxcavg
287 : !character(len=500) :: message
288 35 : type(gs_hamiltonian_type) :: gs_hamk
289 : type(xcdata_type) :: xcdata
290 : !arrays
291 35 : integer,allocatable :: kg_k(:,:)
292 : real(dp) :: gmet(3,3),gprimd(3,3),kpg_dum(0,0),kpoint(3),nonlop_out(1,1)
293 : real(dp) :: qpt(3),rhodum(1),rmet(3,3),tsec(2),ylmgr_dum(1,1,1),vhspinfield(4)
294 70 : real(dp) :: magvec(dtset%nspden)
295 : real(dp),target :: vxctau_dum(0,0,0)
296 35 : real(dp),allocatable :: buffer(:)
297 35 : real(dp),allocatable :: cwavef(:,:),eig_k(:),enlout(:),ffnl(:,:,:,:),ffnl_sav(:,:,:,:)
298 35 : real(dp),allocatable :: kinpw(:),kinpw_sav(:),kxc(:,:),occ_k(:),occblock(:)
299 35 : real(dp),allocatable :: ph3d(:,:,:),ph3d_sav(:,:,:)
300 35 : real(dp),allocatable :: resid_k(:),rhowfg(:,:),rhowfr(:,:),vlocal(:,:,:,:)
301 35 : real(dp),allocatable :: vxctaulocal(:,:,:,:,:),ylm_k(:,:),v_constr_dft_r(:,:)
302 35 : real(dp),pointer :: vxctau_(:,:,:)
303 : type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null()
304 35 : type(pawcprj_type),target,allocatable :: cwaveprj(:,:)
305 35 : type(pawcprj_type),pointer :: cwaveprj_gat(:,:)
306 35 : type(pawrhoij_type),pointer :: pawrhoij_unsym(:)
307 : ! *************************************************************************
308 :
309 : DBG_ENTER("COLL")
310 :
311 : !Test size of kinetic energy potential Vxctau
312 35 : with_vxctau = (present(vxctau))
313 140 : if (with_vxctau) with_vxctau = (size(vxctau)>0.and.dtset%usekden/=0)
314 : if (with_vxctau) then
315 0 : if (size(vxctau)/=nfftf*dtset%nspden*4) then
316 0 : ABI_BUG("Wrong size for vxctau!")
317 : end if
318 : end if
319 35 : vxctau_ => vxctau_dum ; if (with_vxctau) vxctau_ => vxctau
320 :
321 : !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW)
322 140 : nfftotf=PRODUCT(ngfftf(1:3))
323 35 : if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then
324 0 : ABI_BUG('wrong values for nfft, nfftf!')
325 : end if
326 :
327 : !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed
328 35 : wvlbigdft=(dtset%usewvl==1 .and. dtset%wvl_bigdft_comp==1)
329 :
330 35 : call timab(59,1,tsec)
331 :
332 : !Data for parallelism
333 35 : spaceComm=mpi_enreg%comm_cell
334 35 : if(mpi_enreg%paral_kgb==1) spaceComm=mpi_enreg%comm_kpt
335 35 : if(mpi_enreg%paral_hf==1) spaceComm=mpi_enreg%comm_kpt
336 35 : mpi_comm_sphgrid=mpi_enreg%comm_fft
337 35 : if(dtset%usewvl==1) then
338 0 : spaceComm=mpi_enreg%comm_wvl
339 0 : mpi_comm_sphgrid=mpi_enreg%comm_wvl
340 : end if
341 35 : me_distrb=mpi_enreg%me_kpt
342 35 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
343 35 : paral_atom=(my_natom/=dtset%natom)
344 :
345 : !Compute gmet, gprimd and ucvol from rprimd
346 35 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
347 35 : if (dtset%usewvl == 0) then
348 35 : ucvol_local = ucvol
349 : #if defined HAVE_BIGDFT
350 : else
351 : ! We need to tune the volume when wavelets are used because, not
352 : ! all FFT points are used.
353 : ! ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3)
354 : ucvol_local = product(wvl_den%denspot%dpbox%hgrids) * real(product(wvl_den%denspot%dpbox%ndims), dp)
355 : #endif
356 : end if
357 :
358 : !Compute Hxc potential from density
359 35 : option=1;nkxc=0
360 35 : ipositron=electronpositron_calctype(electronpositron)
361 35 : add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw
362 35 : non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4)
363 35 : if (ipositron/=1) then
364 :
365 35 : if (dtset%icoulomb == 0) then
366 : ! Use the periodic solver to compute Hxc.
367 : call hartre(1,gsqcut,dtset%icutcoul,psps%usepaw,mpi_enreg,nfftf,ngfftf,&
368 35 : &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr)
369 35 : call xcdata_init(xcdata,dtset=dtset)
370 35 : ABI_MALLOC(kxc,(1,nkxc))
371 : ! to be adjusted for the call to rhotoxc
372 35 : nk3xc=1
373 35 : if (ipositron==0) then
374 : call rhotoxc(energies%e_xc,energies%entropy_xc,kxc, &
375 : & mpi_enreg,nfftf,ngfftf,nhat,psps%usepaw,nhatgr,nhatgrdim, &
376 : & nkxc,nk3xc,non_magnetic_xc,n3xccc,option,rhor,rprimd, &
377 : & usexcnhat,vxc,vxcavg,xccc3d,xcdata,taur=taur,vhartr=vhartr, &
378 : & vxctau=vxctau_,exc_vdw_out=energies%e_xc_vdw,add_tfw=add_tfw_, &
379 35 : & xcctau3d=xcctau3d,strsxc=strsxc)
380 : else
381 : call rhotoxc(energies%e_xc,energies%entropy_xc,kxc, &
382 : & mpi_enreg,nfftf,ngfftf,nhat,psps%usepaw,nhatgr,nhatgrdim, &
383 : & nkxc,nk3xc,non_magnetic_xc,n3xccc,option,rhor,rprimd, &
384 : & usexcnhat,vxc,vxcavg,xccc3d,xcdata, &
385 : & electronpositron=electronpositron,taur=taur,vhartr=vhartr, &
386 : & vxctau=vxctau_,exc_vdw_out=energies%e_xc_vdw,add_tfw=add_tfw_, &
387 0 : & xcctau3d=xcctau3d,strsxc=strsxc)
388 : end if
389 70 : ABI_FREE(kxc)
390 0 : else if (dtset%usewvl == 0) then
391 : ! Use the free boundary solver.
392 : call psolver_rhohxc(energies%e_hartree, energies%e_xc, evxc, &
393 : & dtset%icoulomb, dtset%ixc, mpi_enreg, nfftf, &
394 : & ngfftf,nhat,psps%usepaw,&
395 : & dtset%nscforder,dtset%nspden,n3xccc,rhor,rprimd, &
396 : & usexcnhat,psps%usepaw,dtset%usewvl,&
397 : & vhartr, vxc, vxcavg, wvl,wvl_den,wvl_e,&
398 0 : & xccc3d,dtset%xclevel,dtset%xc_denpos)
399 : end if
400 : else
401 0 : energies%e_xc=zero
402 : call rhohxcpositron(electronpositron,gprimd,kxc,mpi_enreg,nfftf,ngfftf,nhat,nkxc,dtset%nspden,n3xccc,&
403 0 : & dtset%paral_kgb,rhor,strsxc,ucvol,usexcnhat,psps%usepaw,vhartr,vxc,vxcavg,xccc3d,dtset%xc_denpos)
404 : end if
405 35 : if (ipositron/=0) then
406 : call dotprod_vn(1,rhor,electronpositron%e_hartree,doti,nfftf,nfftotf,1,1,electronpositron%vha_ep,&
407 0 : & ucvol,mpi_comm_sphgrid=mpi_comm_sphgrid)
408 0 : vhartr=vhartr+electronpositron%vha_ep
409 : end if
410 :
411 : !Total local potential (for either spin channel) is
412 : !Hartree + local psp + Vxc(spin), minus its mean
413 : !(Note : this potential should agree with the input vtrial)
414 70 : do ispden=1,min(dtset%nspden,2)
415 535080 : do ifft=1,nfftf
416 535045 : vtrial(ifft,ispden)=vhartr(ifft)+vpsp(ifft)+vxc(ifft,ispden)
417 : end do
418 : end do
419 35 : if (dtset%nspden==4) then
420 0 : do ifft=1,nfftf
421 0 : vtrial(ifft,3:4)=vxc(ifft,3:4)
422 : end do
423 : end if
424 :
425 : !Add the vhspinfield pot in the trial pot
426 : !Vhspinfield might have to be allocated correctly --> to be checked
427 140 : if (any(abs(dtset%hspinfield(:))>tol8)) then
428 0 : vhspinfield(:) = zero
429 0 : if(dtset%nspden==2)then
430 : !TODO: check this against rhotov and setvtr, where the potential is -1/2 and +1/2 for the 2 spin components.
431 : ! see comment by SPr in rhotov
432 : ! TODO: check this 1/2 factor is for the electron spin magnetic moment.
433 0 : vhspinfield(1) = -half*dtset%hspinfield(3) ! For collinear ispden=1 potential is v_upup
434 0 : vhspinfield(2) = +half*dtset%hspinfield(3) ! For collinear ispden=2 potential is v_dndn
435 : end if
436 0 : if(dtset%nspden==4)then
437 0 : vhspinfield(1)=-half*dtset%hspinfield(3)
438 0 : vhspinfield(2)= half*dtset%hspinfield(3)
439 0 : vhspinfield(3)=-half*dtset%hspinfield(1)
440 0 : vhspinfield(4)= half*dtset%hspinfield(2)
441 : end if
442 0 : magvec = zero
443 0 : do ispden=1,dtset%nspden
444 0 : do ifft=1,nfftf
445 : !TODO: the full cell magnetization will need extra PAW terms, and is certainly calculated elsewhere.
446 : !The calculation of the spin magnetic energy can be moved there
447 0 : magvec(ispden) = magvec(ispden) + rhor(ifft,ispden)
448 0 : vtrial(ifft,ispden)=vtrial(ifft,ispden)+vhspinfield(ispden)
449 : end do
450 : end do
451 0 : if(dtset%nspden==2)then
452 0 : energies%e_hspinfield = -half*dtset%hspinfield(3)*(two*magvec(2)-magvec(1)) ! diff rho = rhoup-rhodown = 2 rhoup - rho
453 0 : else if(dtset%nspden==4)then
454 : energies%e_hspinfield = -half * (dtset%hspinfield(1)*magvec(2)& ! x
455 : & +dtset%hspinfield(2)*magvec(3)& ! y
456 0 : & +dtset%hspinfield(3)*magvec(4)) ! z
457 : end if
458 : end if
459 :
460 : !Compute the constrained potential for the magnetic moments
461 : !NOTE: here in energy.F90 rhor and vtrial are given on nfftf grid
462 : !the values coming from mag_penalty may be different from those calculated
463 : !calling mag_penalty with nfft in setvtr and rhotov
464 35 : if (dtset%magconon==1.or.dtset%magconon==2) then
465 0 : ABI_MALLOC(v_constr_dft_r, (nfftf,dtset%nspden))
466 0 : v_constr_dft_r = zero
467 0 : call mag_penalty(constrained_dft,mpi_enreg,rhor,v_constr_dft_r,xred,dtset%qgbt,dtset%use_gbt)
468 : ! call mag_penalty(dtset%natom, dtset%spinat, dtset%nspden, dtset%magconon, dtset%magcon_lambda, rprimd, &
469 : !& mpi_enreg, nfftf, dtset%ngfft, dtset%ntypat, dtset%ratsph, rhor, &
470 : !& dtset%typat, v_constr_dft_r, xred)
471 0 : do ispden=1,dtset%nspden
472 0 : do ifft=1,nfftf
473 0 : vtrial(ifft,ispden)=vtrial(ifft,ispden)+v_constr_dft_r(ifft,ispden)
474 : end do
475 : end do
476 0 : ABI_FREE(v_constr_dft_r)
477 : end if
478 :
479 : !Compute Hartree energy - use up+down rhor
480 35 : if (ipositron/=1) then
481 : call dotprod_vn(1,rhor,energies%e_hartree ,doti,nfftf,nfftotf,1,1,vhartr,&
482 35 : & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid)
483 35 : if (ipositron==0) energies%e_hartree=half*energies%e_hartree
484 35 : if (ipositron==2) energies%e_hartree = half *(energies%e_hartree-electronpositron%e_hartree)
485 : else
486 0 : energies%e_hartree=zero
487 : end if
488 :
489 : !Compute local psp energy - use up+down rhor
490 35 : if (optene/=1) then
491 : call dotprod_vn(1,rhor,energies%e_localpsp,doti,nfftf,nfftotf,1,1,vpsp,&
492 35 : & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid)
493 : end if
494 :
495 : !Compute DC-xc energy - use up+down rhor
496 35 : if (optene>0) then
497 3 : if (ipositron/=1) then
498 3 : if (psps%usepaw==0.or.usexcnhat/=0) then
499 : call dotprod_vn(1,rhor,energies%e_xcdc,doti,nfftf,nfftotf,dtset%nspden,1,vxc,&
500 1 : & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid)
501 1 : if (with_vxctau)then
502 : call dotprod_vn(1,taur,e_xcdc_vxctau,doti,nfftf,nfftotf,dtset%nspden,1,vxctau(:,:,1),&
503 0 : & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid)
504 0 : energies%e_xcdc=energies%e_xcdc+e_xcdc_vxctau
505 : end if
506 : else
507 8 : ABI_MALLOC(rhowfr,(nfftf,dtset%nspden))
508 54006 : rhowfr=rhor-nhat
509 : call dotprod_vn(1,rhowfr,energies%e_xcdc,doti,nfftf,nfftotf,dtset%nspden,1,vxc,&
510 2 : & ucvol_local,mpi_comm_sphgrid=mpi_comm_sphgrid)
511 2 : ABI_FREE(rhowfr)
512 : end if
513 3 : if (ipositron==2) energies%e_xcdc=energies%e_xcdc-electronpositron%e_xcdc
514 : else
515 0 : energies%e_xcdc=zero
516 : end if
517 : end if
518 :
519 35 : energies%e_eigenvalues=zero
520 35 : energies%e_kinetic=zero
521 35 : energies%e_nlpsp_vfock=zero
522 35 : energies%e_fock0=zero
523 35 : bdtot_index=0
524 35 : icg=0
525 :
526 35 : n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
527 35 : n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
528 :
529 : !============================================
530 : !==== Initialize most of the Hamiltonian ====
531 : !============================================
532 : !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
533 : !2) Perform the setup needed for the non-local factors:
534 : !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
535 : !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
536 :
537 : call gs_hamk%init(psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,&
538 : & dtset%natom,dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,&
539 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
540 : & paw_ij=paw_ij,ph1d=ph1d,electronpositron=electronpositron,&
541 35 : & nucdipmom=dtset%nucdipmom,gpu_option=dtset%gpu_option)
542 :
543 210 : ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamk%nvloc))
544 35 : if (with_vxctau) then
545 0 : ABI_MALLOC(vxctaulocal,(n4,n5,n6,gs_hamk%nvloc,4))
546 : end if
547 :
548 : !PAW: additional initializations
549 35 : if (psps%usepaw==1) then
550 59 : ABI_MALLOC(cwaveprj,(dtset%natom,my_nspinor))
551 3 : call pawcprj_alloc(cwaveprj,0,gs_hamk%dimcprj)
552 3 : if (mpi_enreg%paral_spinor==1) then
553 0 : ABI_MALLOC(cwaveprj_gat,(dtset%natom,dtset%nspinor))
554 0 : call pawcprj_alloc(cwaveprj_gat,0,gs_hamk%dimcprj)
555 : else
556 3 : cwaveprj_gat => cwaveprj
557 : end if
558 3 : if (paral_atom) then
559 0 : ABI_MALLOC(pawrhoij_unsym,(dtset%natom))
560 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
561 0 : & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cpxocc=dtset%pawcpxocc)
562 : call pawrhoij_alloc(pawrhoij_unsym,cplex_rhoij,nspden_rhoij,dtset%nspinor,&
563 0 : & dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1)
564 : else
565 3 : pawrhoij_unsym => pawrhoij
566 3 : call pawrhoij_init_unpacked(pawrhoij_unsym)
567 : end if
568 3 : option_rhoij=1
569 3 : use_timerev=(dtset%kptopt>0.and.dtset%kptopt<3)
570 6 : use_zeromag=(pawrhoij_unsym(1)%nspden==4.and.dtset%nspden==1)
571 : else
572 32 : ABI_MALLOC(cwaveprj,(0,0))
573 : end if
574 :
575 : !LOOP OVER SPINS
576 70 : do isppol=1,dtset%nsppol
577 35 : ikg=0
578 :
579 : ! Set up local potential vlocal on the coarse FFT mesh from vtrial taking into account the spin.
580 : ! Also take into account the spin.
581 :
582 : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
583 35 : dtset%nspden, gs_hamk%nvloc, 1, pawfgr, mpi_enreg, vtrial, vlocal)
584 35 : call gs_hamk%load_spin(isppol,vlocal=vlocal,with_nonlocal=.true.)
585 :
586 35 : if (with_vxctau) then
587 : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
588 0 : dtset%nspden, gs_hamk%nvloc, 4, pawfgr, mpi_enreg, vxctau, vxctaulocal)
589 0 : call gs_hamk%load_spin(isppol, vxctaulocal=vxctaulocal)
590 : end if
591 :
592 : ! Loop over k points
593 280 : do ikpt=1,dtset%nkpt
594 210 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
595 210 : istwf_k=dtset%istwfk(ikpt)
596 210 : npw_k=npwarr(ikpt)
597 :
598 : ! Skip this k-point if not the proper processor
599 210 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
600 0 : resid(1+bdtot_index : nband_k+bdtot_index) = zero
601 : bdtot_index=bdtot_index+nband_k
602 : cycle
603 : end if
604 :
605 : ! Parallelism over FFT and/or bands: define sizes and tabs
606 210 : if (mpi_enreg%paral_kgb==1) then
607 0 : my_ikpt=mpi_enreg%my_kpttab(ikpt)
608 0 : nblockbd=nband_k/(mpi_enreg%nproc_band*mpi_enreg%bandpp)
609 0 : my_bandfft_kpt => bandfft_kpt(my_ikpt)
610 : else
611 210 : my_ikpt=ikpt
612 210 : nblockbd=nband_k
613 : !if (nband_k/=nblockbd*mpi_enreg%nproc_fft) nblockbd=nblockbd+1
614 : end if
615 210 : blocksize=nband_k/nblockbd
616 :
617 630 : ABI_MALLOC(eig_k,(nband_k))
618 420 : ABI_MALLOC(occ_k,(nband_k))
619 420 : ABI_MALLOC(resid_k,(nband_k))
620 630 : ABI_MALLOC(cwavef,(2,npw_k*my_nspinor*blocksize))
621 2198 : resid_k(:)=zero
622 840 : kpoint(:)=dtset%kptns(:,ikpt)
623 2198 : occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
624 2198 : eig_k(:)=eigen(1+bdtot_index:nband_k+bdtot_index)
625 2577 : if (minval(eig_k)>1.d100) eig_k=zero
626 210 : eeigk=zero ; ekk=zero ; enlk=zero
627 :
628 630 : ABI_MALLOC(kg_k,(3,npw_k))
629 346246 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
630 :
631 840 : ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
632 210 : if (psps%useylm==1) then
633 370 : do ilm=1,psps%mpsang*psps%mpsang
634 61934 : ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
635 : end do
636 : end if
637 :
638 : ! Compute kinetic energy
639 630 : ABI_MALLOC(kinpw,(npw_k))
640 210 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw,kpoint,npw_k,0,0)
641 :
642 : ! Compute kinetic energy of each band
643 2198 : do iblock=1,nblockbd
644 4186 : do iblocksize=1,blocksize
645 1988 : iband=(iblock-1)*blocksize+iblocksize
646 3976 : if (abs(occ_k(iband))>tol8) then
647 : cwavef(1:2,1:npw_k*my_nspinor)= &
648 2115366 : & cg(:,1+(iband-1)*npw_k*my_nspinor+icg:iband*npw_k*my_nspinor+icg)
649 906 : call meanvalue_g(dotr,kinpw,0,istwf_k,mpi_enreg,npw_k,my_nspinor,cwavef,cwavef,0)
650 906 : energies%e_kinetic=energies%e_kinetic+dtset%wtk(ikpt)*occ_k(iband)*dotr
651 : end if
652 : end do
653 : end do
654 :
655 : ! Compute nonlocal form factors ffnl at all (k+G):
656 210 : ider=0;dimffnl=1;nkpg=0
657 1050 : ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
658 : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,&
659 : & gmet,gprimd,ider,ider,psps%indlmn,kg_k,kpg_dum,kpoint,psps%lmnmax,&
660 : & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
661 : & npw_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,&
662 210 : & psps%usepaw,psps%useylm,ylm_k,ylmgr_dum)
663 :
664 : ! Load k-dependent part in the Hamiltonian datastructure
665 : ! - Compute 3D phase factors
666 : ! - Prepare various tabs in case of band-FFT parallelism
667 : ! - Load k-dependent quantities in the Hamiltonian
668 840 : ABI_MALLOC(ph3d,(2,npw_k,gs_hamk%matblk))
669 : call gs_hamk%load_k(kpt_k=dtset%kptns(:,ikpt),istwf_k=istwf_k,npw_k=npw_k,&
670 : & kinpw_k=kinpw,kg_k=kg_k,ffnl_k=ffnl,ph3d_k=ph3d,&
671 210 : & compute_ph3d=.true.,compute_gbound=(mpi_enreg%paral_kgb/=1))
672 :
673 : ! Load band-FFT tabs (transposed k-dependent arrays)
674 210 : if (mpi_enreg%paral_kgb==1) then
675 0 : call bandfft_kpt_savetabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kinpw=kinpw_sav)
676 0 : call prep_bandfft_tabs(gs_hamk,ikpt,dtset%mkmem,mpi_enreg)
677 : call gs_hamk%load_k(npw_fft_k=my_bandfft_kpt%ndatarecv, &
678 : & gbound_k =my_bandfft_kpt%gbound, &
679 : & kinpw_k =my_bandfft_kpt%kinpw_gather, &
680 : & kg_k =my_bandfft_kpt%kg_k_gather, &
681 : & ffnl_k =my_bandfft_kpt%ffnl_gather, &
682 0 : & ph3d_k =my_bandfft_kpt%ph3d_gather)
683 : end if
684 :
685 : ! If OpenMP GPU, load "hamiltonian" on GPU device
686 210 : if (dtset%gpu_option == ABI_GPU_OPENMP) then
687 0 : if(dtset%paral_kgb==0) then
688 0 : call ompgpu_load_hamilt_buffers(gs_hamk%kg_k,gs_hamk%kg_kp,gs_hamk%ffnl_k,gs_hamk%ph3d_k)
689 0 : else if(istwf_k==1) then
690 0 : call ompgpu_load_hamilt_buffers(gs_hamk%kg_k,gs_hamk%kg_kp,gs_hamk%ffnl_k,gs_hamk%ph3d_k,kg_k_gather=my_bandfft_kpt%kg_k_gather)
691 0 : else if(istwf_k==2) then
692 0 : call ompgpu_load_hamilt_buffers(gs_hamk%kg_k,gs_hamk%kg_kp,gs_hamk%ffnl_k,gs_hamk%ph3d_k,kg_k_gather=my_bandfft_kpt%kg_k_gather_sym)
693 : else
694 0 : ABI_ERROR("istwfk > 2 is not handled with OpenMP GPU offload mode !")
695 : end if
696 : end if
697 :
698 210 : choice=1-gs_hamk%usepaw ; signs=1 ; idir=0 ; nnlout=blocksize
699 :
700 : ! Setup gemm_nonlop
701 210 : if (gemm_nonlop_use_gemm) then
702 : call set_gemm_nonlop_ikpt(my_ikpt,gs_hamk%npw_fft_k,gs_hamk%istwf_k,gs_hamk%indlmn,&
703 0 : & gs_hamk%ntypat,gs_hamk%nattyp,gs_hamk%gpu_option)
704 : end if
705 :
706 : #if defined HAVE_GPU_CUDA
707 : if (dtset%gpu_option==ABI_GPU_LEGACY .or. dtset%gpu_option==ABI_GPU_KOKKOS) then
708 : call gpu_update_ffnl_ph3d( &
709 : & ph3d, INT(size(ph3d,dim=1),c_int64_t)*size(ph3d,dim=2)*size(ph3d,dim=3), &
710 : & ffnl, INT(size(ffnl),c_int64_t) )
711 : end if
712 : #endif
713 :
714 : ! Compute nonlocal psp energy (NCPP) or Rhoij (PAW)
715 630 : ABI_MALLOC(enlout,(blocksize))
716 420 : ABI_MALLOC(occblock,(blocksize))
717 2198 : do iblock=1,nblockbd
718 1988 : iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k)
719 1988 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle
720 :
721 : ! Select occupied bands
722 3976 : occblock(:)=occ(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
723 6174 : if(abs(maxval(occblock))>=tol8 ) then
724 : cwavef(:,1:npw_k*my_nspinor*blocksize)=&
725 2115366 : & cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg)
726 :
727 906 : paw_opt=gs_hamk%usepaw;cpopt=gs_hamk%usepaw-1
728 :
729 906 : if (mpi_enreg%paral_kgb/=1) then
730 906 : tim_nonlop=3
731 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,(/zero/),mpi_enreg,blocksize,nnlout,&
732 906 : & paw_opt,signs,nonlop_out,tim_nonlop,cwavef,cwavef)
733 : else
734 0 : tim_nonlop=14
735 : call prep_nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,(/zero/),blocksize,&
736 0 : & mpi_enreg,nnlout,paw_opt,signs,nonlop_out,tim_nonlop,cwavef,cwavef)
737 : end if
738 :
739 1812 : do iblocksize=1,blocksize
740 906 : iband=(iblock-1)*blocksize+iblocksize
741 906 : energies%e_eigenvalues=energies%e_eigenvalues+dtset%wtk(ikpt)*occ_k(iband)*eig_k(iband)
742 : ! WARNING : the Fock contribution is NOT computed !!!
743 1812 : energies%e_nlpsp_vfock=energies%e_nlpsp_vfock+dtset%wtk(ikpt)*occ_k(iband)*enlout(iblocksize)
744 : end do
745 :
746 : ! PAW: accumulate rhoij
747 906 : if (psps%usepaw==1) then
748 452 : cplex=merge(1,2,istwf_k>1)
749 452 : if (mpi_enreg%paral_spinor==1) then
750 : call pawcprj_gather_spin(cwaveprj,cwaveprj_gat,dtset%natom,1,my_nspinor,dtset%nspinor,&
751 0 : & mpi_enreg%comm_spinor,ierr)
752 : call pawaccrhoij(gs_hamk%atindx,cplex,cwaveprj_gat,cwaveprj_gat,0,isppol,dtset%natom,dtset%natom,&
753 0 : & dtset%nspinor,occ_k(iband),option_rhoij,pawrhoij_unsym,use_timerev,use_zeromag,dtset%wtk(ikpt))
754 : else
755 : call pawaccrhoij(gs_hamk%atindx,cplex,cwaveprj,cwaveprj,0,isppol,dtset%natom,dtset%natom,&
756 452 : & dtset%nspinor,occ_k(iband),option_rhoij,pawrhoij_unsym,use_timerev,use_zeromag,dtset%wtk(ikpt))
757 : end if
758 : end if
759 :
760 : ! End loop on bands
761 : end if
762 : end do
763 :
764 : ! Compute residual of each band (for informative purposes)
765 210 : call mkresi(cg,eig_k,gs_hamk,icg,ikpt,isppol,mcg,mpi_enreg,nband_k,dtset%prtvol,resid_k)
766 2198 : resid(1+bdtot_index : nband_k+bdtot_index) = resid_k(:)
767 :
768 : ! Restore the bandfft tabs
769 210 : if (mpi_enreg%paral_kgb==1) then
770 0 : call bandfft_kpt_restoretabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kinpw=kinpw_sav)
771 : end if
772 :
773 : ! Incremente indexes
774 210 : bdtot_index=bdtot_index+nband_k
775 210 : if (dtset%mkmem/=0) then
776 210 : icg=icg+npw_k*my_nspinor*nband_k
777 210 : ikg=ikg+npw_k
778 : end if
779 :
780 : #if defined HAVE_GPU_CUDA
781 : if(dtset%gpu_option==ABI_GPU_LEGACY .or. dtset%gpu_option==ABI_GPU_KOKKOS) then
782 : call gpu_finalize_ffnl_ph3d()
783 : end if
784 : #endif
785 :
786 210 : ABI_FREE(eig_k)
787 210 : ABI_FREE(occ_k)
788 210 : ABI_FREE(resid_k)
789 210 : ABI_FREE(enlout)
790 210 : ABI_FREE(occblock)
791 210 : ABI_FREE(ffnl)
792 210 : ABI_FREE(kinpw)
793 210 : ABI_FREE(ph3d)
794 210 : ABI_FREE(cwavef)
795 210 : ABI_FREE(kg_k)
796 245 : ABI_FREE(ylm_k)
797 :
798 : ! End loops on isppol and ikpt
799 : end do
800 : end do
801 :
802 35 : call gs_hamk%free()
803 35 : if ( dtset%gpu_option == ABI_GPU_OPENMP) then
804 0 : call ompgpu_free_hamilt_buffers()
805 : end if
806 :
807 : if(xmpi_paral==1)then
808 : ! Accumulate enl eeig and ek on all proc.
809 105 : ABI_MALLOC(buffer,(3+dtset%mband*dtset%nkpt*dtset%nsppol))
810 35 : buffer(1)=energies%e_nlpsp_vfock ; buffer(2)=energies%e_kinetic ; buffer(3)=energies%e_eigenvalues
811 2023 : do iresid=1,dtset%mband*dtset%nkpt*dtset%nsppol
812 2023 : buffer(iresid+3)=resid(iresid)
813 : end do
814 35 : call timab(48,1,tsec)
815 35 : call xmpi_sum(buffer,spaceComm,ierr)
816 35 : call timab(48,2,tsec)
817 35 : energies%e_nlpsp_vfock=buffer(1) ; energies%e_kinetic=buffer(2) ; energies%e_eigenvalues=buffer(3)
818 2023 : do iresid=1,dtset%mband*dtset%nkpt*dtset%nsppol
819 2023 : resid(iresid)=buffer(iresid+3)
820 : end do
821 35 : ABI_FREE(buffer)
822 : ! Accumulate rhoij_
823 70 : if (psps%usepaw==1) then
824 3 : call pawrhoij_mpisum_unpacked(pawrhoij_unsym,spaceComm,comm2=mpi_enreg%comm_band)
825 : end if
826 : end if
827 :
828 : !Compute total (free) energy
829 35 : if (optene==0.or.optene==2) then
830 : etotal = energies%e_kinetic + energies%e_hartree + energies%e_xc + &
831 : !& energies%e_nlpsp_vfock - energies%e_fock0 +
832 : ! Should compute the e_fock0 energy !! Also, the Fock contribution to e_nlpsp_vfock
833 34 : & energies%e_nlpsp_vfock + energies%e_localpsp + energies%e_corepsp
834 34 : if (psps%usepaw==1) etotal=etotal + energies%paw%epaw
835 1 : else if (optene==1.or.optene==3) then
836 : etotal = energies%e_eigenvalues - energies%e_hartree + energies%e_xc - &
837 1 : & energies%e_xcdc + energies%e_corepsp - energies%e_corepspdc
838 1 : if (psps%usepaw==1) etotal=etotal + energies%paw%epaw_dc
839 : end if
840 35 : etotal = etotal + energies%e_ewald + energies%e_chempot + energies%e_vdw_dftd
841 :
842 : !Add the contribution of extfpmd to the entropy
843 35 : if(associated(extfpmd)) then
844 0 : energies%e_extfpmd=extfpmd%e_kinetic
845 0 : energies%edc_extfpmd=extfpmd%edc_kinetic
846 0 : if(optene==0.or.optene==2) etotal=etotal+energies%e_extfpmd
847 0 : if(optene==1.or.optene==3) etotal=etotal+energies%edc_extfpmd
848 : end if
849 :
850 : ! Add the contribution from cores
851 35 : if(present(rcpaw)) then
852 35 : if(associated(rcpaw)) then
853 0 : energies%paw%epaw_core=rcpaw%ehnzc+rcpaw%ekinc
854 0 : energies%paw%epaw_core_dc=rcpaw%eeigc-rcpaw%edcc+rcpaw%ehnzc
855 0 : if (ipositron/=1) then
856 0 : do itypat=1,dtset%ntypat
857 0 : energies%paw%epaw_core=energies%paw%epaw_core+pawtab(itypat)%exccore*rcpaw%atm(itypat)%mult
858 0 : energies%paw%epaw_core_dc=energies%paw%epaw_core_dc+pawtab(itypat)%exccore*rcpaw%atm(itypat)%mult
859 : enddo
860 : endif
861 0 : if(optene==0.or.optene==2) etotal=etotal+energies%paw%epaw_core
862 0 : if(optene==1.or.optene==3) etotal=etotal+energies%paw%epaw_core_dc
863 : endif
864 : endif
865 :
866 35 : call entropy(dtset,energies)
867 35 : etotal=etotal+energies%e_entropy
868 :
869 : !Additional stuff for electron-positron
870 35 : if (dtset%positron/=0) then
871 0 : if (ipositron==0) then
872 0 : energies%e_electronpositron =zero
873 0 : energies%edc_electronpositron=zero
874 : else
875 0 : energies%e_electronpositron =electronpositron%e_hartree+electronpositron%e_xc
876 0 : energies%edc_electronpositron=electronpositron%e_hartree+electronpositron%e_xcdc
877 0 : if (psps%usepaw==1) then
878 0 : energies%e_electronpositron =energies%e_electronpositron +electronpositron%e_paw
879 0 : energies%edc_electronpositron=energies%edc_electronpositron+electronpositron%e_pawdc
880 : end if
881 : end if
882 0 : if (optene==0.or.optene==2) electronpositron%e0=etotal
883 0 : if (optene==1.or.optene==3) electronpositron%e0=etotal-energies%edc_electronpositron
884 0 : etotal=electronpositron%e0+energies%e0_electronpositron+energies%e_electronpositron
885 : end if
886 :
887 : !Compute new charge density based on incoming wf
888 : !Keep rhor and rhog intact for later use e.g. in stress. (? MT 08-12-2008: is that true now ?)
889 : !=== Norm-conserving psps: simply compute rho from WFs
890 : !paw_dmft%use_dmft=0 ! dmft not used here
891 : !paw_dmft%use_sc_dmft=0 ! dmft not used here
892 35 : if (psps%usepaw==0) then
893 32 : tim_mkrho=3
894 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,&
895 : & npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl_den,wfs,&
896 32 : & extfpmd=extfpmd)
897 32 : if(dtset%usekden==1)then
898 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,&
899 0 : & npwarr,occ,paw_dmft,phnons,taug,taur,rprimd,tim_mkrho,ucvol,wvl_den,wfs,option=1)
900 : end if
901 : else
902 : ! === PAW case: symmetrize rhoij and add compensation charge density
903 3 : tim_mkrho=3;option=1;choice=1
904 : call pawrhoij_symrhoij(pawrhoij,pawrhoij_unsym,choice,gprimd,indsym,0,dtset%natom,dtset%nsym,&
905 : & dtset%ntypat,option,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,dtset%typat,&
906 3 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
907 3 : ider=0;izero=0;cplex=1;ipert=0;idir=0;qpt(:)=zero
908 : call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,&
909 : & my_natom,dtset%natom,nfftf,ngfftf,&
910 : & 0,dtset%nspden,dtset%ntypat,pawang,pawfgrtab,rhodum,nhat,pawrhoij,pawrhoij,&
911 : & pawtab,qpt,rprimd,ucvol_local,dtset%usewvl,xred,&
912 : & comm_fft=mpi_enreg%comm_fft,paral_kgb=dtset%paral_kgb,me_g0=mpi_enreg%me_g0,&
913 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
914 3 : & distribfft=mpi_enreg%distribfft,mpi_comm_wvl=mpi_enreg%comm_wvl)
915 :
916 12 : ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden))
917 9 : ABI_MALLOC(rhowfg,(2,dtset%nfft))
918 98602 : rhowfr(:,:)=zero
919 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,&
920 3 : & npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,rprimd,tim_mkrho,ucvol_local,wvl_den,wfs)
921 :
922 : call pawmkrho(1,compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,&
923 : & my_natom,dtset%natom,dtset%nspden,dtset%nsym,dtset%ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,&
924 : & dtset%pawprtvol,pawrhoij,pawrhoij_unsym,pawtab,qpt,rhowfg,rhowfr,rhor,rprimd,dtset%symafm,&
925 3 : & symrec,dtset%typat,ucvol,dtset%usewvl,xred,pawnhat=nhat,rhog=rhog,extfpmd=extfpmd)
926 3 : call pawrhoij_free_unpacked(pawrhoij_unsym)
927 3 : if (paral_atom) then
928 0 : call pawrhoij_free(pawrhoij_unsym)
929 0 : ABI_FREE(pawrhoij_unsym)
930 : end if
931 3 : if(dtset%usekden==1)then
932 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
933 0 : & rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl_den,wfs,option=1)
934 0 : call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,taug,rhowfr,taur)
935 : end if
936 3 : ABI_FREE(rhowfr)
937 3 : ABI_FREE(rhowfg)
938 : end if
939 :
940 35 : ABI_COMMENT('New density rho(r) made from input wfs')
941 :
942 35 : call timab(59,2,tsec)
943 :
944 35 : ABI_FREE(vlocal)
945 35 : if (with_vxctau) then
946 0 : ABI_FREE(vxctaulocal)
947 : end if
948 :
949 35 : if (psps%usepaw==1) then
950 3 : call pawcprj_free(cwaveprj)
951 47 : ABI_FREE(cwaveprj)
952 3 : if (mpi_enreg%paral_spinor==1) then
953 0 : call pawcprj_free(cwaveprj_gat)
954 0 : ABI_FREE(cwaveprj_gat)
955 : else
956 3 : nullify(cwaveprj_gat)
957 : end if
958 : end if
959 :
960 : DBG_EXIT("COLL")
961 :
962 70 : end subroutine energy
963 : !!***
964 :
965 : !!****f* ABINIT/mkresi
966 : !! NAME
967 : !! mkresi
968 : !!
969 : !! FUNCTION
970 : !! Make residuals from knowledge of wf in G space and application of Hamiltonian.
971 : !!
972 : !! INPUTS
973 : !! cg(2,mcg)=<G|Cnk>=Fourier coefficients of wavefunction
974 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
975 : !! icg=shift to be applied on the location of data in the array cg
976 : !! ikpt=index of k-point
977 : !! isppol=index of spin
978 : !! mcg=second dimension of the cg array
979 : !! mpi_enreg=information about MPI parallelization
980 : !! nband=number of bands involved in subspace matrix.
981 : !! npw=number of planewaves in basis sphere at this k point.
982 : !! prtvol=control print volume and debugging output
983 : !! usepaw= 0 for non paw calculation; =1 for paw calculation
984 : !!
985 : !! OUTPUT
986 : !! eig_k(nband)$= \langle C_n \mid H \mid C_n \rangle $ for each band.
987 : !! resid_k(nband)=residual for each band
988 : !! $= \langle C_n \mid H H \mid C_n \rangle- \langle C_n \mid H \mid C_n \rangle^2 $.
989 : !!
990 : !! SOURCE
991 :
992 210 : subroutine mkresi(cg,eig_k,gs_hamk,icg,ikpt,isppol,mcg,mpi_enreg,nband,prtvol,resid_k)
993 :
994 : !Arguments ------------------------------------
995 : !scalars
996 : integer,intent(in) :: icg,ikpt,isppol,mcg,nband,prtvol
997 : type(MPI_type),intent(inout) :: mpi_enreg
998 : type(gs_hamiltonian_type),intent(inout) :: gs_hamk
999 : !arrays
1000 : real(dp),intent(in) :: cg(2,mcg)
1001 : real(dp),intent(out) :: eig_k(nband),resid_k(nband)
1002 :
1003 : !Local variables-------------------------------
1004 : !scalars
1005 : integer,parameter :: tim_getghc=3
1006 : integer :: blocksize,cpopt,iband,iband_last,iblock,iblocksize,ipw,ipw_shift
1007 : integer :: my_nspinor,nblockbd,npw_k
1008 : real(dp) :: doti,dotr
1009 : !arrays
1010 : real(dp) :: tsec(2)
1011 420 : real(dp),allocatable,target :: cwavef(:,:),ghc(:,:),gsc(:,:),gvnlxc(:,:)
1012 420 : real(dp), contiguous, pointer :: cwavef_ptr(:,:),ghc_ptr(:,:),gsc_ptr(:,:)
1013 1470 : type(pawcprj_type) :: cwaveprj(1,1)
1014 : ! *************************************************************************
1015 :
1016 : !Keep track of total time spent in mkresi
1017 210 : call timab(13,1,tsec)
1018 :
1019 : !Parallelism over FFT and/or bands: define sizes and tabs
1020 210 : my_nspinor=max(1,gs_hamk%nspinor/mpi_enreg%nproc_spinor)
1021 210 : if (mpi_enreg%paral_kgb==1) then
1022 0 : nblockbd=nband/(mpi_enreg%nproc_band*mpi_enreg%bandpp)
1023 : else
1024 210 : nblockbd=nband/mpi_enreg%nproc_fft
1025 210 : if (nband/=nblockbd*mpi_enreg%nproc_fft) nblockbd=nblockbd+1
1026 : end if
1027 210 : blocksize=nband/nblockbd
1028 :
1029 210 : npw_k=gs_hamk%npw_k
1030 630 : ABI_MALLOC(cwavef,(2,npw_k*my_nspinor))
1031 420 : ABI_MALLOC(ghc,(2,npw_k*my_nspinor))
1032 420 : ABI_MALLOC(gvnlxc,(2,npw_k*my_nspinor))
1033 210 : if (gs_hamk%usepaw==1) then
1034 146 : ABI_MALLOC(gsc,(2,npw_k*my_nspinor))
1035 : else
1036 137 : ABI_MALLOC(gsc,(0,0))
1037 : end if
1038 :
1039 : !Loop over (blocks of) bands
1040 2198 : do iblock=1,nblockbd
1041 1988 : iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband)
1042 1988 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,mpi_enreg%me_kpt)) cycle
1043 :
1044 : ! Load |Cn>
1045 1988 : ipw_shift=(iblock-1)*npw_k*my_nspinor*blocksize+icg
1046 : !$OMP PARALLEL DO
1047 892572 : do ipw=1,npw_k*my_nspinor*blocksize
1048 890584 : cwavef(1,ipw)=cg(1,ipw+ipw_shift)
1049 892572 : cwavef(2,ipw)=cg(2,ipw+ipw_shift)
1050 : end do
1051 :
1052 : ! Compute H|Cn>
1053 1988 : cpopt=-1
1054 1988 : if (mpi_enreg%paral_kgb==0) then
1055 : call getghc(cpopt,cwavef,cwaveprj,ghc,gsc,gs_hamk,gvnlxc,zero,mpi_enreg,1,&
1056 1988 : & prtvol,gs_hamk%usepaw,tim_getghc,0)
1057 : else
1058 : call prep_getghc(cwavef,gs_hamk,gvnlxc,ghc,gsc,zero,nband,mpi_enreg,&
1059 : & prtvol,gs_hamk%usepaw,cpopt,cwaveprj,&
1060 0 : & already_transposed=.false.)
1061 : end if
1062 :
1063 : !call cg_get_eigens(usepaw, istwf_k, npwsp, nband, cg, ghc, gsc, eig, me_g0, comm_bsf)
1064 : !call cg_get_residvecs(usepaw, npwsp, nband, eig, cg, ghc, gsc, gwork)
1065 : !call cg_norm2g(istwf_k, npwsp, nband, gwork, resid, me_g0, comm_bsf)
1066 : ! MG: Communicators are wrongi if paral_kgb. One should use mpi_enreg%comm_bandspinorfft
1067 :
1068 : ! Compute the residual, <Cn|(H-<Cn|H|Cn>)**2|Cn>:
1069 4186 : do iblocksize=1,blocksize
1070 1988 : iband=(iblock-1)*blocksize+iblocksize
1071 1988 : ipw_shift=(iblocksize-1)*npw_k*my_nspinor
1072 1988 : cwavef_ptr => cwavef(:,1+ipw_shift:npw_k*my_nspinor+ipw_shift)
1073 1988 : ghc_ptr => ghc (:,1+ipw_shift:npw_k*my_nspinor+ipw_shift)
1074 :
1075 : ! First get eigenvalue <Cn|H|Cn>:
1076 : call dotprod_g(dotr,doti,gs_hamk%istwf_k,npw_k*my_nspinor,1,cwavef_ptr,ghc_ptr,&
1077 1988 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
1078 1988 : eig_k(iband)=dotr
1079 :
1080 : ! Next need <G|(H-S<Cn|H|Cn>)|Cn> (in ghc):
1081 1988 : if (gs_hamk%usepaw==0) then
1082 : !$OMP PARALLEL DO PRIVATE(ipw) SHARED(cwavef_ptr,ghc_ptr,eig_k,iband,npw_k,my_nspinor)
1083 362941 : do ipw=1,npw_k*my_nspinor
1084 362050 : ghc_ptr(1,ipw)=ghc_ptr(1,ipw)-eig_k(iband)*cwavef_ptr(1,ipw)
1085 362941 : ghc_ptr(2,ipw)=ghc_ptr(2,ipw)-eig_k(iband)*cwavef_ptr(2,ipw)
1086 : end do
1087 : else
1088 1097 : gsc_ptr => gsc(:,1+ipw_shift:npw_k*my_nspinor+ipw_shift)
1089 : !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gsc_ptr,ghc_ptr,eig_k,iband,npw_k,my_nspinor)
1090 529631 : do ipw=1,npw_k*my_nspinor
1091 528534 : ghc_ptr(1,ipw)=ghc_ptr(1,ipw)-eig_k(iband)*gsc_ptr(1,ipw)
1092 529631 : ghc_ptr(2,ipw)=ghc_ptr(2,ipw)-eig_k(iband)*gsc_ptr(2,ipw)
1093 : end do
1094 : end if
1095 :
1096 : ! Then simply square the result:
1097 : call sqnorm_g(dotr,gs_hamk%istwf_k,npw_k*my_nspinor,ghc_ptr,&
1098 1988 : & mpi_enreg%me_g0,mpi_enreg%comm_fft)
1099 3976 : resid_k(iband)=dotr
1100 :
1101 : end do ! iblocksize
1102 :
1103 : end do ! iblock
1104 :
1105 210 : ABI_FREE(cwavef)
1106 210 : ABI_FREE(ghc)
1107 210 : ABI_FREE(gvnlxc)
1108 210 : ABI_FREE(gsc)
1109 :
1110 210 : call timab(13,2,tsec)
1111 :
1112 840 : end subroutine mkresi
1113 : !!***
1114 :
1115 : !!****f* ABINIT/entropy
1116 : !! NAME
1117 : !! entropy
1118 : !!
1119 : !! FUNCTION
1120 : !! Compute electronic entropy terms
1121 : !! This subroutine returns the total entropy and entropy energy. In the most
1122 : !! common case, at finite temperature, the electronic entropy is mainly constitued
1123 : !! of the non-interacting entropy (entropy_ks). Finite-temperature exchange-correlation
1124 : !! functionals or other methods may introduce additional entropy terms.
1125 : !!
1126 : !! NOTE
1127 : !!
1128 : !! INPUTS
1129 : !! dtset <type(dataset_type)>=all input variables for this dataset
1130 : !! | occopt=option for occupancies
1131 : !! | tsmear=smearing energy or temperature (if metal)
1132 : !! | tphysel=electornic temperature for particular values of occopt
1133 : !! energies <type(energies_type)>=all part of total energy.
1134 : !! | entropy_ks(IN)=non-interacting entropy of the kohn-sham states
1135 : !! | entropy_paw(IN)=entropy due to paw corrections (for finite-temperature xc functionals)
1136 : !! | entropy_xc(IN)=exchange-correlation entropy (for finite-temperature xc functionals)
1137 : !! | entropy_extfpmd(IN)=entropy of extfpmd model
1138 : !! | entropy_imp(IN)=entropy of impurity electrons (for DFT+DMFT)
1139 : !!
1140 : !! OUTPUT
1141 : !! energies <type(energies_type)>=all part of total energy.
1142 : !! | entropy(OUT)=total entropy
1143 : !! | e_entropy(OUT)=total entropy energy (hartree units)
1144 : !!
1145 : !! SOURCE
1146 53849 : subroutine entropy(dtset,energies)
1147 : !Arguments ------------------------------------
1148 : !scalars
1149 : type(dataset_type),intent(in) :: dtset
1150 : type(energies_type),intent(inout) :: energies
1151 : ! *************************************************************************
1152 :
1153 : !In case we have other sources of entropy than kohn-sham states occupation,
1154 : !we sum all entropy terms. %entropy is now total entropy.
1155 : !Examples of other sources of entropy: finite-temperature xc functionals, extfpmd, ...
1156 53849 : energies%entropy=energies%entropy_ks
1157 53849 : if(abs(energies%paw%entropy_paw)>tiny(zero)) energies%entropy=energies%entropy+energies%paw%entropy_paw
1158 53849 : if(abs(energies%entropy_xc)>tiny(zero)) energies%entropy=energies%entropy+energies%entropy_xc
1159 53849 : if(abs(energies%entropy_extfpmd)>tiny(zero)) energies%entropy=energies%entropy+energies%entropy_extfpmd
1160 53849 : if(abs(energies%entropy_imp)>tiny(zero)) energies%entropy=energies%entropy+energies%entropy_imp
1161 :
1162 : !When the finite-temperature VG broadening scheme is used,
1163 : !the total entropy contribution "tsmear*entropy" has a meaning,
1164 : !and gather the two last terms of Eq.8 of VG paper
1165 : !Warning : might have to be changed for fixed moment calculations
1166 53849 : if(dtset%occopt>=3 .and. dtset%occopt<=8) then
1167 14965 : if (abs(dtset%tphysel) < tol10) then
1168 14919 : energies%e_entropy = - dtset%tsmear * energies%entropy
1169 : else
1170 46 : energies%e_entropy = - dtset%tphysel * energies%entropy
1171 : end if
1172 : else
1173 38884 : energies%e_entropy = zero
1174 : end if
1175 :
1176 53849 : end subroutine entropy
1177 :
1178 : end module m_dft_energy
1179 : !!***
|
