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