Branch data 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 : 54215 : 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 : 54215 : energies%entropy=energies%entropy_ks
1157 [ + + ]: 54215 : if(abs(energies%paw%entropy_paw)>tiny(zero)) energies%entropy=energies%entropy+energies%paw%entropy_paw
1158 [ + + ]: 54215 : if(abs(energies%entropy_xc)>tiny(zero)) energies%entropy=energies%entropy+energies%entropy_xc
1159 [ + + ]: 54215 : if(abs(energies%entropy_extfpmd)>tiny(zero)) energies%entropy=energies%entropy+energies%entropy_extfpmd
1160 [ - + ]: 54215 : 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 [ + + ]: 54215 : 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 : 39250 : energies%e_entropy = zero
1174 : : end if
1175 : :
1176 : 54215 : end subroutine entropy
1177 : :
1178 : : end module m_dft_energy
1179 : : !!***
|
