Line data Source code
1 : !!****m* ABINIT/m_vtorho
2 : !! NAME
3 : !! m_vtorho
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, MF, AR, MM, MT, FJ, MB, MT, TR)
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 : ! nvtx related macro definition
23 : #include "nvtx_macros.h"
24 :
25 : module m_vtorho
26 :
27 : use iso_fortran_env, only : int32,int64,real32,real64
28 : use defs_basis
29 : use defs_wvltypes
30 : use m_abicore
31 : use m_xmpi
32 : use m_xomp
33 : use m_abi_mixing
34 : use m_errors
35 : use m_wffile
36 : use m_efield
37 : use m_cgtools
38 : use m_hdr
39 : use m_dtset
40 : use m_dtfil
41 : use m_extfpmd
42 : use m_ompgpu_utils
43 : use m_xg_nonlop
44 : use m_ebands
45 : use m_crystal
46 :
47 : use defs_datatypes, only : pseudopotential_type
48 : use defs_abitypes, only : MPI_type
49 : use m_fstrings, only : sjoin, itoa
50 : use m_time, only : timab
51 : use m_geometry, only : xred2xcart
52 : use m_occ, only : newocc
53 : use m_pawang, only : pawang_type
54 : use m_pawrad, only : pawrad_type
55 : use m_pawtab, only : pawtab_type
56 : use m_paw_ij, only : paw_ij_type
57 : use m_pawfgrtab, only : pawfgrtab_type
58 : use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_io, pawrhoij_inquire_dim
59 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_getdim
60 : use m_pawfgr, only : pawfgr_type
61 : use m_rcpaw, only : rcpaw_type
62 : use m_energies, only : energies_type
63 : use m_hamiltonian, only : gs_hamiltonian_type, gspot_transgrid_and_pack
64 : use m_bandfft_kpt, only : bandfft_kpt, bandfft_kpt_type, bandfft_kpt_set_ikpt, &
65 : bandfft_kpt_savetabs, bandfft_kpt_restoretabs, prep_bandfft_tabs
66 : use m_electronpositron, only : electronpositron_type,electronpositron_calctype
67 : use m_paw_dmft, only : paw_dmft_type,init_dmft,destroy_dmft,print_dmft,saveocc_dmft
68 : use m_paw_correlations, only : setnoccmmp
69 : use m_paw_occupancies, only : pawmkrhoij
70 : use m_paw_mkrho, only : pawmkrho
71 : use m_results_gs, only : results_gs_type
72 : use m_oper, only : oper_type,init_oper,destroy_oper
73 : use m_io_tools, only : flush_unit
74 : use m_abi2big, only : wvl_occ_abi2big, wvl_rho_abi2big, wvl_occopt_abi2big, wvl_eigen_abi2big
75 : use m_fock, only : fock_type, fock_ACE_type, fock_updateikpt, fock_calc_ene
76 : use m_invovl, only : make_invovl
77 : use m_tddft, only : tddft
78 : use m_kg, only : mkkin, mkkpg
79 : use m_suscep_stat, only : suscep_stat
80 : use m_fft, only : fftpac
81 : use m_spacepar, only : symrhg
82 : use m_vtowfk, only : vtowfk
83 : use m_mkrho, only : mkrho, prtrhomxmn
84 : use m_mkffnl, only : mkffnl
85 : use m_mpinfo, only : proc_distrb_cycle
86 : use m_common, only : prteigrs,get_gemm_nonlop_ompgpu_blocksize
87 : use m_dmft, only : dmft_solve
88 : use m_datafordmft, only : datafordmft
89 : use m_fourier_interpol, only : transgrid
90 : use m_cgprj, only : ctocprj
91 : use m_wvl_rho, only : wvl_mkrho
92 : use m_wvl_psi, only : wvl_hpsitopsi, wvl_psitohpsi, wvl_nl_gradient
93 : use m_inwffil, only : cg_from_atoms
94 : use m_gemm_nonlop_projectors, only : set_gemm_nonlop_ikpt, gemm_nonlop_use_gemm, &
95 : gemm_nonlop_block_size, gemm_nonlop_is_distributed
96 : use m_alloc_hamilt_gpu, only : hamilt_gpu_nfft_blocks
97 :
98 : use m_abstract_wf, only : abstract_wf, init_mywfc
99 : use m_mlwfovlp, only : mlwfovlp
100 : #if defined HAVE_PYTHON_INVOCATION
101 : use m_invoke_python
102 : #endif
103 : use ISO_C_BINDING
104 :
105 : #if defined HAVE_GPU_CUDA
106 : use m_manage_cuda
107 : #endif
108 :
109 : #if defined HAVE_YAKL
110 : use gator_mod
111 : #endif
112 :
113 : #if defined HAVE_BIGDFT
114 : use BigDFT_API, only : last_orthon, evaltoocc, write_energies, eigensystem_info
115 : #endif
116 :
117 : #if defined(HAVE_GPU_MARKERS)
118 : use m_nvtx_data
119 : #endif
120 :
121 : !#ifdef HAVE_FC_ISO_C_BINDING
122 : ! use, intrinsic :: iso_c_binding, only : c_int64_t
123 : !#endif
124 :
125 : implicit none
126 :
127 : private
128 : !!***
129 :
130 : public :: vtorho
131 : !!***
132 :
133 : contains
134 : !!***
135 :
136 : !!****f* ABINIT/vtorho
137 : !! NAME
138 : !! vtorho
139 : !!
140 : !! FUNCTION
141 : !! This routine compute the new density from a fixed potential (vtrial)
142 : !! but might also simply compute eigenvectors and eigenvalues.
143 : !! The main part of it is a wf update over all k points.
144 : !!
145 : !! INPUTS
146 : !! itimes(2)=itime array, contain itime=itimes(1) and itimimage_gstate=itimes(2) from outer loops
147 : !! afford=used to dimension susmat
148 : !! atindx(natom)=index table for atoms (see gstate.f)
149 : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
150 : !! cpus= cpu time limit in seconds
151 : !! dbl_nnsclo=if 1, will double the value of dtset%nnsclo
152 : !! dielop= if positive, the dielectric matrix must be computed.
153 : !! dielstrt=number of the step at which the dielectric preconditioning begins.
154 : !! dmatpawu= fixed occupation matrix of correlated orbitals (DFT+U or DMFT only)
155 : !! dtefield <type(efield_type)> = variables related to Berry phase
156 : !! calculations (see initberry.f)
157 : !! dtfil <type(datafiles_type)>=variables related to files
158 : !! dtset <type(dataset_type)>=all input variables for this dataset
159 : !! | mband=maximum number of bands
160 : !! | mgfft=maximum size of 1D FFTs
161 : !! | mkmem =number of k points treated by this node.
162 : !! | mpw=maximum dimensioned size of npw
163 : !! | nfft=(effective) number of FFT grid points (for this processor)
164 : !! | nkpt=number of k points.
165 : !! | nspden=number of spin-density components
166 : !! | nsppol=1 for unpolarized, 2 for spin-polarized
167 : !! | nsym=number of symmetry elements in space group
168 : !! | typat= array of types of the natoms
169 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
170 : !! etotal=total energy (Ha) - only needed for tddft
171 : !! extfpmd <type(extfpmd_type)>=extended first-principles molecular dynamics type
172 : !! fock <type(fock_type)>= quantities to calculate Fock exact exchange
173 : !! gbound_diel(2*mgfftdiel+8,2)=G sphere boundary for the dielectric matrix
174 : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2.
175 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
176 : !! (3x3 tensor) and grads wrt atomic coordinates (3*natom)
177 : !! gsqcut=cutoff on (k+G)^2 (bohr^-2)
178 : !! hdr <type(hdr_type)>=the header of wf, den and pot files
179 : !! indsym(4,nsym,natom)=indirect indexing array for atom labels
180 : !! irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
181 : !! irrzondiel(nfftdiel**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data for diel matrix
182 : !! nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
183 : !! istep=index of the number of steps in the routine scfcv
184 : !! istep_mix=index of the number of steps for the SCF mixing (can be <istep)
185 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
186 : !! kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
187 : !! kxc(nfftf,nkxc)=exchange-correlation kernel, needed only if nkxc/=0 .
188 : !! lmax_diel=1+max. value of l angular momentum used for dielectric matrix
189 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
190 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
191 : !! mgfftdiel=maximum size of 1D FFTs, for the computation of the dielectric matrix
192 : !! mpi_enreg=information about MPI parallelization
193 : !! my_natom=number of atoms treated by current processor
194 : !! natom=number of atoms in cell.
195 : !! nattyp(ntypat)= # atoms of each type.
196 : !! nfftf= -PAW ONLY- number of FFT grid points for the fine grid
197 : !! (nfftf=nfft for norm-conserving potential runs)
198 : !! nfftdiel=number of fft grid points for the computation of the diel matrix
199 : !! ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix,
200 : !! see ~abinit/doc/variables/vargs.htm#ngfft
201 : !! nkxc=second dimension of the array kxc, see rhotoxc.f for a description
202 : !! npwarr(nkpt)=number of planewaves in basis at this k point
203 : !! npwdiel=size of the susmat array.
204 : !! ntypat=number of types of atoms in unit cell.
205 : !! optforces=option for the computation of forces (0: no force;1: forces)
206 : !! optres=0: the new value of the density is computed in place of the input value
207 : !! 1: only the density residual is computed ; the input density is kept
208 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
209 : !! paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
210 : !! pawang <type(pawang)>=paw angular mesh and related data
211 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
212 : !! pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
213 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and
214 : !! related data
215 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
216 : !! phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
217 : !! nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
218 : !! phnonsdiel(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases,
219 : !! for diel matr
220 : !! nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
221 : !! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
222 : !! ph1ddiel(2,3*(2*mgfftdiel+1)*natom*usepaw)=one-dimensional structure factor information
223 : !! for the dielectric matrix
224 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
225 : !! pwind(pwind_alloc,2,3) = array used to compute
226 : !! the overlap matrix smat between k-points (see initberry.f)
227 : !! pwind_alloc = first dimension of pwind
228 : !! pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations (see initberry.f)
229 : !! results_gs <type(results_gs_type)>=results (energy and its components,
230 : !! forces and its components, the stress tensor) of a ground-state
231 : !! computation (should be made a pure output quantity)
232 : !! rmet(3,3)=real space metric (bohr**2)
233 : !! rprimd(3,3)=dimensional primitive vectors
234 : !! symrec(3,3,nsym)=symmetry operations in reciprocal space
235 : !! ucvol=unit cell volume in bohr**3.
236 : !! usecprj=1 if cprj datastructure is stored in memory
237 : !! usevxctau=1 if if XC functional depends on kinetic energy density
238 : !! wffnew,unit numbers for wf disk files.
239 : !! with_vectornd = 1 if vectornd allocated
240 : !! vectornd(with_vectornd*nfftf,nspden,3)=nuclear dipole moment vector potential
241 : !! vtrial(nfftf,nspden)=INPUT potential Vtrial(r).
242 : !! [vxctau(nfftf,nspden,4*usevxctau)]=(only for meta-GGA): derivative of XC energy density
243 : !! with respect to kinetic energy density (depsxcdtau). The arrays vxctau contains also
244 : !! the gradient of vxctau (gvxctau) in vxctau(:,:,2:4)
245 : !! xred(3,natom)=reduced dimensionless atomic coordinates
246 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
247 : !! ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
248 : !! ylmdiel(npwdiel,lmax_diel**2)= real spherical harmonics for each G and k point
249 : !! for the dielectric matrix
250 : !!
251 : !! OUTPUT
252 : !! compch_fft=-PAW only- compensation charge inside spheres computed over fine fft grid
253 : !! dphase(3) : dphase(idir) = accumulated change in the string-averaged
254 : !! Zak phase along the idir-th direction caused by the update of all
255 : !! the occupied Bloch states at all the k-points (only if finite electric field)
256 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
257 : !! resid(mband*nkpt*nsppol)=residuals for each band over all k points.
258 : !! residm=maximum value from resid array (except for nbdbuf highest bands)
259 : !! susmat(2,npwdiel*afford,nspden,npwdiel,nspden)=
260 : !! the susceptibility (or density-density response) matrix in reciprocal space
261 : !! === if optforces>0 ===
262 : !! grnl(3*natom)=stores grads of nonlocal energy wrt length scales
263 : !! ==== if optres==1
264 : !! nres2=square of the norm of the residual
265 : !! nvresid(nfftf,nspden)=density residual
266 : !! ==== if psps%usepaw==1
267 : !! cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors:
268 : !! cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector.
269 : !! nhat(nfftf,nspden*psps%usepaw)=compensation charge density on rectangular grid in real space
270 : !!
271 : !! SIDE EFFECTS
272 : !! cg(2,mpw*dtset%nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions.
273 : !! At output contains updated wavefunctions coefficients;
274 : !! if nkpt>1, these are kept in a disk file.
275 : !! energies <type(energies_type)>=storage for energies computed here :
276 : !! | e_eigenvalues=Sum of the eigenvalues - Band energy (Hartree)
277 : !! | e_kinetic=kinetic energy part of total energy
278 : !! | e_nlpsp_vfock=nonlocal psp + potential Fock ACE part of total energy
279 : !! | e_fermie=fermi energy (Hartree)
280 : !! occ(mband*nkpt*nsppol)=occupation number for each band for each k.
281 : !! (input if insulator - occopt<3 - ; output if metallic)
282 : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
283 : !! rhog(2,nfftf)=Fourier transform of total electron density
284 : !! rhor(nfftf,nspden)=total electron density (el/bohr**3)
285 : !! taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic energy density
286 : !! taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density
287 : !! tauresid(nfftf,nspden*dtset%usekden)=array for kinetic energy density residual
288 : !! wvl <type(wvl_data)>=wavelets structures in case of wavelets basis.
289 : !! ==== if (usepaw==1) ====
290 : !! cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors:
291 : !! cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector.
292 : !! rmm_diis_status= Status of the RMM-DIIS eigensolver. See m_rmm_diis
293 : !!
294 : !! NOTES
295 : !! Be careful to the meaning of nfft (size of FFT grids):
296 : !! - In case of norm-conserving calculations the FFT grid is the usual FFT grid.
297 : !! - In case of PAW calculations:
298 : !! Two FFT grids are used; one with nfft points (coarse grid) for
299 : !! the computation of wave functions ; one with nfftf points
300 : !! (fine grid) for the computation of total density.
301 : !!
302 : !! The total electronic density (rhor,rhog) is divided into two terms:
303 : !! - The density related to WFs =Sum[Psi**2]
304 : !! - The compensation density (nhat) - only in PAW
305 : !!
306 : !! The parallelisation needed for the electric field should be
307 : !! made an independent subroutine, so that this routine could be put
308 : !! back in the 95_drive directory.
309 : !!
310 : !! SOURCE
311 :
312 93986 : subroutine vtorho(afford,atindx,atindx1,cg,compch_fft,cprj,cpus,dbl_nnsclo,&
313 93986 : & dielop,dielstrt,dmatpawu,dphase,dtefield,dtfil,dtset,&
314 93986 : & eigen,electronpositron,energies,etotal,gbound_diel,&
315 93986 : & gmet,gprimd,grnl,gsqcut,hdr,extfpmd,indsym,irrzon,irrzondiel,&
316 93986 : & istep,istep_mix,itimes,kg,kg_diel,kxc,lmax_diel,mcg,mcprj,mgfftdiel,mpi_enreg,&
317 46993 : & my_natom,natom,nattyp,nfftf,nfftdiel,ngfftdiel,nhat,nkxc,&
318 46993 : & npwarr,npwdiel,nres2,ntypat,nvresid,occ,optforces,&
319 93986 : & optres,paw_dmft,paw_ij,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,&
320 93986 : & phnons,phnonsdiel,ph1d,ph1ddiel,psps,fock,&
321 93986 : & pwind,pwind_alloc,pwnsfac,results_gs,resid,residm,rhog,rhor,&
322 93986 : & rmet,rprimd,susmat,symrec,taug,taur,tauresid,&
323 46993 : & ucvol,usecprj,usevxctau,wffnew,with_vectornd,vectornd,vtrial,vxctau,wvl,&
324 46993 : & xg_nonlop,xred,ylm,ylmgr,ylmdiel,rmm_diis_status,rcpaw)
325 :
326 : !Arguments -------------------------------
327 : !scalars
328 : integer, intent(in) :: afford,dbl_nnsclo,dielop,dielstrt,istep,istep_mix,lmax_diel,mcg,mcprj,mgfftdiel
329 : integer, intent(in) :: my_natom,natom,nfftf,nfftdiel,nkxc,npwdiel
330 : integer, intent(in) :: ntypat,optforces,optres,pwind_alloc,usecprj,usevxctau,with_vectornd
331 : real(dp), intent(in) :: cpus,etotal,gsqcut,ucvol
332 : real(dp), intent(out) :: compch_fft,nres2,residm
333 : type(MPI_type), intent(inout) :: mpi_enreg
334 : type(datafiles_type), intent(in) :: dtfil
335 : type(dataset_type), intent(inout) :: dtset
336 : type(efield_type), intent(inout) :: dtefield
337 : type(electronpositron_type),pointer :: electronpositron
338 : type(energies_type), intent(inout) :: energies
339 : type(hdr_type), intent(inout) :: hdr
340 : type(extfpmd_type), pointer, intent(inout) :: extfpmd
341 : type(paw_dmft_type), intent(inout) :: paw_dmft
342 : type(pawang_type), intent(in) :: pawang
343 : type(pawfgr_type), intent(in) :: pawfgr
344 : type(pseudopotential_type), intent(in) :: psps
345 : type(results_gs_type),intent(inout) :: results_gs
346 : type(fock_type),pointer, intent(inout) :: fock
347 : type(wffile_type), intent(inout) :: wffnew
348 : type(wvl_data), intent(inout) :: wvl
349 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
350 : type(rcpaw_type),pointer,intent(inout) :: rcpaw
351 : !arrays
352 : integer, intent(in) :: atindx(natom),atindx1(natom),gbound_diel(2*mgfftdiel+8,2)
353 : integer, intent(in) :: indsym(4,dtset%nsym,natom)
354 : integer, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
355 : integer, intent(in) :: irrzondiel(nfftdiel**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
356 : integer, intent(in) :: itimes(2)
357 : integer, intent(in) :: kg(3,dtset%mpw*dtset%mkmem),kg_diel(3,npwdiel),nattyp(ntypat),ngfftdiel(18),npwarr(dtset%nkpt)
358 : integer, intent(in) :: pwind(pwind_alloc,2,3),symrec(3,3,dtset%nsym)
359 : integer, intent(inout) :: rmm_diis_status(2, dtset%nkpt, dtset%nsppol)
360 : real(dp), intent(in) :: dmatpawu(:,:,:,:),gmet(3,3),gprimd(3,3),ph1d(2,3*(2*dtset%mgfft+1)*natom)
361 : real(dp), intent(in) :: ph1ddiel(2,(3*(2*mgfftdiel+1)*natom)*psps%usepaw)
362 : real(dp), intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
363 : real(dp), intent(in) :: phnonsdiel(2,nfftdiel**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
364 : real(dp), intent(in) :: pwnsfac(2,pwind_alloc),rmet(3,3),rprimd(3,3)
365 : real(dp), intent(inout) :: vectornd(with_vectornd*nfftf,dtset%nspden,3),vtrial(nfftf,dtset%nspden)
366 : real(dp), intent(inout) :: xred(3,natom)
367 : real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
368 : real(dp), intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
369 : real(dp), intent(in) :: ylmdiel(npwdiel,lmax_diel**2)
370 : real(dp), intent(out) :: dphase(3),grnl(3*natom)
371 : real(dp), intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
372 : real(dp), intent(out) :: nhat(nfftf,dtset%nspden*psps%usepaw)
373 : real(dp), intent(inout) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol)
374 : real(dp), intent(out) :: nvresid(nfftf,dtset%nspden)
375 : real(dp), intent(out) :: susmat(2,npwdiel*afford,dtset%nspden,npwdiel,dtset%nspden)
376 : real(dp), intent(inout) :: cg(2,mcg)
377 : real(dp), intent(inout) :: kxc(nfftf,nkxc),occ(dtset%mband*dtset%nkpt*dtset%nsppol)
378 : real(dp), intent(inout) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden)
379 : real(dp), intent(inout) :: taug(2,nfftf*dtset%usekden),taur(nfftf,dtset%nspden*dtset%usekden)
380 : real(dp), intent(inout) :: tauresid(nfftf,dtset%nspden*dtset%usekden)
381 : real(dp), intent(inout),optional :: vxctau(nfftf,dtset%nspden,4*usevxctau)
382 : type(pawcprj_type),pointer,intent(inout) :: cprj(:,:)
383 : type(paw_ij_type),intent(inout) :: paw_ij(my_natom*psps%usepaw)
384 : type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
385 : type(pawrad_type), intent(in) :: pawrad(psps%ntypat*psps%usepaw)
386 : type(pawrhoij_type),target,intent(inout) :: pawrhoij(my_natom*psps%usepaw)
387 : type(pawtab_type),intent(inout) :: pawtab(ntypat*psps%usepaw)
388 :
389 : !Local variables-------------------------------
390 : !scalars
391 : ! integer,parameter :: level=111
392 : integer,parameter :: tim_mkrho=2
393 : !integer,save :: nwarning=0
394 : integer :: bdtot_index,counter,cplex,cplex_rhoij,dimffnl,enunit,iband,iband1,ibdkpt
395 : integer :: ibg,icg,ider,idir,ierr,ifft,ifor,ifor1,ii,ikg,ikpt,itypat
396 : integer :: ikpt_loc,ikpt1,my_ikpt,ikxc,ilm,imagn,index1,iorder_cprj,ipert
397 : integer :: iscf,ispden,isppol,istwf_k,mband_cprj,mbdkpsp,mb2dkpsp
398 : integer :: mcgq,mcprj_local,mcprj_tmp,me_distrb,mkgq,mpi_comm_sphgrid
399 : integer :: my_nspinor,n1,n2,n3,n4,n5,n6,nband_eff,nbdbuf_eff !mwarning,
400 : integer :: nband_k,nband_cprj_k,nbuf,neglect_pawhat,nfftot,nkpg,nkpt1,nnsclo_now
401 : integer :: nproc_distrb,npw_k,nspden_rhoij,option,prtvol,quit,nblk_gemm_nonlop
402 : integer :: spaceComm_distrb,usecprj_local,usefock_ACE,usetimerev
403 : #if defined HAVE_GPU_CUDA
404 : integer(c_int64_t) :: ph3d_size
405 : #endif
406 : integer :: nthreads,nmpi,mpicomm
407 :
408 : logical :: berryflag,computesusmat,fixed_occ,has_vectornd,step_cond
409 : logical :: locc_test,paral_atom,remove_inv,usefock,with_vxctau,transfer_cg_once
410 : logical :: do_last_ortho,wvlbigdft=.false.,do_invS,calc_ffnl_ph3d,gpu_mem_estimated
411 : integer :: dmft_dftocc
412 : real(dp) :: nelect,min_eigv
413 : real(dp) :: edmft,ebandlda,ebanddmft,ebandldatot,ekindmft,ekindmft2,ekinlda
414 : real(dp) :: min_occ,vxcavg_dum,strsxc(6)
415 : character(len=500) :: msg
416 : type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null()
417 46993 : type(gs_hamiltonian_type) :: gs_hamk
418 : !arrays
419 : integer(int32), contiguous, pointer :: kg_k(:,:) => null()
420 : real(dp) :: dielar(7),dphase_k(3),kpoint(3),qpt(3),rhodum(1),tsec(2),ylmgr_dum(0,0,0), kphq(3)
421 93986 : real(dp),allocatable :: EigMin(:,:),buffer1(:),cgq(:,:)
422 93986 : real(dp),allocatable :: cgrkxc(:,:),doccde(:)
423 46993 : real(dp),allocatable :: dphasek(:,:),ek_k(:),ek_k_nd(:,:,:),eknk(:),eknk_nd(:,:,:,:,:),end_k(:)
424 46993 : real(dp),allocatable :: enlx_k(:),enlxnk(:),focknk(:),fockfornk(:,:,:),ffnl(:,:,:,:), ffnl_kphq(:,:,:,:)
425 93986 : real(dp),allocatable :: grnlnk(:,:), grnl_k(:,:), xcart(:,:)
426 46993 : real(dp),allocatable :: nvresid_tmp(:,:)
427 :
428 : #if defined HAVE_GPU && defined HAVE_YAKL
429 : real(c_double), contiguous, pointer :: kinpw(:) => null()
430 : real(c_double), contiguous, pointer :: eig_k(:) => null()
431 : #else
432 46993 : real(dp),allocatable :: kinpw(:), eig_k(:)
433 : #endif
434 :
435 46993 : real(dp),allocatable :: kpg_k(:,:),kpg_kphq(:,:),occ_k(:),ph3d(:,:,:), ph3d_kphq(:,:,:)
436 46993 : real(dp),allocatable :: pwnsfacq(:,:), kinpw_kphq(:)
437 :
438 : #if defined HAVE_GPU && defined HAVE_YAKL
439 : real(c_double), contiguous, pointer :: resid_k(:) => null()
440 : real(c_double), contiguous, pointer :: rhoaug(:,:,:,:) => null()
441 : #else
442 46993 : real(dp),allocatable :: resid_k(:), rhoaug(:,:,:,:)
443 : #endif
444 :
445 46993 : real(dp),allocatable :: rhowfg(:,:),rhowfr(:,:),tauwfg(:,:),tauwfr(:,:), vectornd_pac(:,:,:,:,:)
446 :
447 : #if defined HAVE_GPU && defined HAVE_YAKL
448 : real(real64), contiguous, pointer :: vlocal(:,:,:,:) => null()
449 : #else
450 46993 : real(dp), allocatable :: vlocal(:,:,:,:)
451 : #endif
452 :
453 46993 : real(dp),allocatable :: vxctaulocal(:,:,:,:,:),ylm_k(:,:),zshift(:)
454 46993 : type(pawcprj_type),allocatable :: cprj_tmp(:,:)
455 46993 : type(pawcprj_type),allocatable,target:: cprj_local(:,:)
456 46993 : type(oper_type) :: dft_occup
457 46993 : type(pawrhoij_type),pointer :: pawrhoij_unsym(:)
458 2443636 : type(crystal_t) :: cryst_struc
459 : integer :: idum1(0),idum3(0,0,0)
460 : real(dp) :: rdum2(0,0),rdum4(0,0,0,0)
461 : #if defined HAVE_BIGDFT
462 : integer :: occopt_bigdft
463 : #endif
464 :
465 : #if defined(HAVE_PYTHON_INVOCATION)
466 : integer :: bantot
467 : logical :: exists
468 : character(len=500) :: filename
469 : class(abstract_wf), pointer :: mywfc
470 : integer :: exclude_bands(hdr%mband, hdr%nsppol)
471 : logical :: exclude_bands_ind(hdr%mband, hdr%nsppol)
472 : type(ebands_t) :: ebands
473 : real(dp), allocatable :: occnd_tmp(:)
474 : #endif
475 : ! *********************************************************************
476 :
477 : DBG_ENTER("COLL")
478 :
479 : !Keep track of total time spent in vtorho
480 46993 : call timab(980,1,tsec)
481 46993 : call timab(981,1,tsec)
482 :
483 : !Structured debugging if prtvol==-level
484 46993 : prtvol=dtset%prtvol
485 :
486 : ! Electric fields: set flag to turn on various behaviors
487 327901 : berryflag = any(dtset%berryopt == [4, 14, 6, 16, 7, 17])
488 :
489 : !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed
490 46993 : wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1)
491 :
492 : !Several inits
493 46993 : n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
494 46993 : n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
495 46993 : usecprj_local=0;if (psps%usepaw==1) usecprj_local=1
496 46993 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
497 46993 : paral_atom=(my_natom/=natom)
498 46993 : compch_fft=-1.d5
499 46993 : gpu_mem_estimated=.false.; if(istep>1) gpu_mem_estimated=.true.
500 :
501 : !Check that usekden is not 0 if want to use vxctau
502 46993 : with_vxctau = (present(vxctau).and.dtset%usekden/=0.and.usevxctau/=0)
503 :
504 : !Check that fock is present if want to use fock option
505 46993 : usefock = (dtset%usefock==1 .and. associated(fock))
506 1375 : usefock_ACE=0
507 1375 : if (usefock) usefock_ACE=fock%fock_common%use_ACE
508 :
509 : !Init MPI
510 46993 : spaceComm_distrb=mpi_enreg%comm_cell
511 46993 : if (mpi_enreg%paral_kgb==1) spaceComm_distrb=mpi_enreg%comm_kpt
512 46993 : if (mpi_enreg%paral_hf ==1) spaceComm_distrb=mpi_enreg%comm_kpt
513 46993 : nproc_distrb=xmpi_comm_size(spaceComm_distrb)
514 46993 : me_distrb=xmpi_comm_rank(spaceComm_distrb)
515 46993 : mpi_comm_sphgrid=mpi_enreg%comm_fft
516 46993 : if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl
517 : !if (mpi_enreg%me_img/=0) nwarning=nwarning+1
518 :
519 : !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW)
520 46993 : if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then
521 0 : ABI_BUG('wrong values for nfft, nfftf!')
522 : end if
523 :
524 : !Test optforces (to prevent memory overflow)
525 46993 : if (optforces/=0.and.optforces/=1) then
526 0 : ABI_BUG(sjoin('wrong value for optforces: ',itoa(optforces)))
527 : end if
528 :
529 46993 : iscf=dtset%iscf
530 46993 : fixed_occ=(dtset%occopt<3.or.electronpositron_calctype(electronpositron)==1)
531 46993 : if(.not. wvlbigdft) then
532 46993 : energies%e_eigenvalues = zero
533 46993 : energies%e_kinetic = zero
534 46993 : energies%e_nucdip = zero
535 46993 : energies%e_nlpsp_vfock = zero
536 46993 : if (usefock) then
537 1375 : energies%e_fock=zero
538 1375 : energies%e_fockdc=zero
539 : end if
540 425005 : grnl(:)=zero
541 46993 : if (berryflag) then
542 72799 : resid(:) = zero ! JWZ 13 May 2010. resid and eigen need to be fully zeroed each time before use
543 : end if
544 : ! MG: The previous line is not compatible with the RMM-DIIS since rmm_diis recevies the previous resid_k
545 : ! to select the accuracy level.
546 : ! For the time being, we set resid to zero if berryflag to avoid breaking the CG solver with E-field
547 : ! but it's clear that the treatment of resid should be rationalized and that the previous values should be passed to vtowfk
548 3066385 : eigen(:) = zero
549 46993 : bdtot_index=0
550 46993 : ibg=0;icg=0
551 46993 : mbdkpsp=dtset%mband*dtset%nkpt*dtset%nsppol
552 : if(paw_dmft%use_dmft==1) mb2dkpsp=2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol
553 : end if
554 :
555 46993 : if(dtset%usewvl==0) then
556 140979 : ABI_MALLOC(eknk,(mbdkpsp))
557 93986 : ABI_MALLOC(enlxnk,(mbdkpsp))
558 281958 : ABI_MALLOC(eknk_nd,(2,dtset%mband,dtset%mband*paw_dmft%use_dmft,dtset%nkpt,dtset%nsppol))
559 140979 : ABI_MALLOC(EigMin,(2,dtset%mband))
560 234965 : ABI_MALLOC(grnlnk,(3*natom,mbdkpsp*optforces))
561 46993 : if (usefock) then
562 2750 : ABI_MALLOC(focknk,(mbdkpsp))
563 44143 : focknk=zero
564 1375 : if (optforces>0)then
565 880 : ABI_MALLOC(fockfornk,(3,natom,mbdkpsp))
566 13216 : fockfornk=zero
567 : end if
568 : end if
569 6085777 : eknk(:)=zero;enlxnk(:)=zero
570 6592766 : if (optforces>0) grnlnk(:,:)=zero
571 2256446 : if (paw_dmft%use_dmft==1) eknk_nd(:,:,:,:,:)=zero
572 : end if !usewvl==0
573 :
574 : !Initialize rhor if needed; store old rhor
575 46993 : if(iscf>=0 .or. iscf==-3) then
576 46517 : if (optres==1) then
577 779365295 : nvresid=rhor ; tauresid=taur
578 : end if
579 : ! NC and plane waves
580 46517 : if (psps%usepaw==0 .and. dtset%usewvl==0) then
581 645002257 : rhor=zero ; taur=zero
582 15106 : else if(psps%usepaw==1) then
583 : ! PAW
584 60424 : ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden))
585 45318 : ABI_MALLOC(rhowfg,(2,dtset%nfft))
586 60424 : ABI_MALLOC(tauwfr,(dtset%nfft,dtset%nspden*dtset%usekden))
587 45318 : ABI_MALLOC(tauwfg,(2,dtset%nfft*dtset%usekden))
588 259435641 : rhowfr(:,:)=zero ; tauwfr(:,:)=zero
589 : end if
590 : end if
591 :
592 : ! Here we set the max number of non-self-consistent loops nnsclo_now used in vtowfk
593 46993 : if (iscf<0) then
594 : ! Non self-consistent case
595 496 : nnsclo_now=dtset%nstep
596 : else
597 : ! Self-consistent case
598 46497 : if (dtset%nnsclo>0) then
599 : ! Use input variable if specified and > 0
600 446 : nnsclo_now=dtset%nnsclo
601 46051 : else if (dtset%nnsclo < 0) then
602 : ! imposed during abs(nnsclo) steps
603 50 : nnsclo_now=1
604 50 : if (istep<=abs(dtset%nnsclo)) nnsclo_now=merge(5,dtset%useria,dtset%useria==0)
605 : else
606 : ! Default branch for self-consistent case.
607 : ! Perform 2 NSCF loops for the first two iterations. This is important especially wfs have
608 : ! been initialized with random numbers.
609 46001 : nnsclo_now = 1
610 46001 : if (dtset%usewvl == 0) then
611 : ! Plane waves
612 46001 : if (istep <= 2 .and. iscf /= 0) nnsclo_now = 2
613 : ! MG: I don't understand why we need to perform 2 NSCF loops after the first SCF cycle
614 : ! when we are relaxing the structure as the initial density and wavefunctions should be already good enough.
615 : ! Here I change the default behavior to avoid the extra loop but only if RMM-DIIS is used.
616 : ! XG 20210312 : I prefectly agree with you. This is historical, and should be changed, after testing and update of reference files.
617 46001 : if ((itimes(1) > 1 .or. (itimes(2)>1)) .and. dtset%rmm_diis /= 0) nnsclo_now = 1
618 : else
619 : ! Wavelets
620 0 : if (iscf==0) then
621 0 : nnsclo_now=0
622 0 : else if (istep<=2) then
623 0 : nnsclo_now=3
624 0 : else if (istep<=4) then
625 0 : nnsclo_now=2
626 : end if
627 : end if
628 : end if
629 : ! Double the value if required
630 46497 : if (dbl_nnsclo==1) nnsclo_now=nnsclo_now*2
631 : end if
632 :
633 46993 : if(dtset%wfoptalg==2)nnsclo_now=40 ! UNDER DEVELOPMENT
634 :
635 : transfer_cg_once = dtset%gpu_option==ABI_GPU_OPENMP &
636 : & .and. (dtset%wfoptalg == 114 .or. dtset%wfoptalg == 112 .or. dtset%wfoptalg == 111) &
637 : & .and. dtset%mkmem == 1 .and. dtset%nsppol == 1
638 : #ifdef HAVE_OPENMP_OFFLOAD
639 : ! cg is allocated once on GPU if we're using "XG" algos
640 : ! and only one k-point and spin are in use
641 : !$OMP TARGET ENTER DATA MAP(alloc:cg) IF(transfer_cg_once)
642 : #endif
643 :
644 46993 : if (dtset%prtvol > 0) then
645 4392 : write(msg, '(a,i0,a,3(i0,1x))' ) ' vtorho: nnsclo_now = ',nnsclo_now,&
646 8784 : ', note that nnsclo, dbl_nnsclo, istep= ',dtset%nnsclo,dbl_nnsclo,istep
647 4392 : call wrtout(std_out,msg)
648 : else
649 42601 : if (nnsclo_now > 1) call wrtout(std_out, sjoin(" Max number of non-self-consistent loops:", itoa(nnsclo_now)))
650 : end if
651 :
652 : !==== Initialize most of the Hamiltonian ====
653 : ! Allocate all arrays and initialize quantities that do not depend on k and spin.
654 : call gs_hamk%init(psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,natom,&
655 : dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,&
656 : paw_ij=paw_ij,ph1d=ph1d,usecprj=usecprj_local,electronpositron=electronpositron,fock=fock,&
657 : comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
658 : nucdipmom=dtset%nucdipmom,gpu_option=dtset%gpu_option,spinaxis=dtset%spinaxis,&
659 46993 : use_gbt=dtset%use_gbt,nfft_blocks=hamilt_gpu_nfft_blocks)
660 :
661 46993 : if (dtset%cprj_in_memory==1) then
662 2406 : call xg_nonlop_update_weight(xg_nonlop,ucvol) ! ucvol could have changed in mover
663 2406 : if (xg_nonlop%paw) call xg_nonlop_make_Dij(xg_nonlop,paw_ij,dtset%nsppol,atindx)
664 : end if
665 :
666 : !Initializations for PAW (projected wave functions)
667 46993 : mcprj_local=0 ; mband_cprj=0
668 46993 : if (psps%usepaw==1) then
669 15126 : mband_cprj=dtset%mband
670 15126 : if (dtset%paral_kgb/=0) mband_cprj=mband_cprj/mpi_enreg%nproc_band
671 15126 : iorder_cprj=0 ; mcprj_local=mcprj
672 15126 : if (usecprj==0) then
673 14486 : mcprj_local=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol
674 : !This is a check but should always be true since scfcv allocated cprj anyway
675 : if (allocated(cprj_local)) then
676 : !Was allocated in scfcv so we just destroy and reconstruct it as desired
677 : call pawcprj_free(cprj_local)
678 : ABI_FREE(cprj_local)
679 : end if
680 4706603 : ABI_MALLOC(cprj_local,(dtset%natom,mcprj_local))
681 14486 : call pawcprj_alloc(cprj_local,0,gs_hamk%dimcprj)
682 : cprj => null()
683 14486 : cprj => cprj_local
684 : end if
685 : end if
686 :
687 46993 : call timab(981,2,tsec)
688 :
689 : !===================================================================
690 : ! WAVELETS - Branching with a separate VTORHO procedure
691 : !===================================================================
692 :
693 46993 : if (dtset%usewvl == 1) then
694 : #ifndef HAVE_BIGDFT
695 0 : BIGDFT_NOTENABLED_ERROR()
696 : #else
697 :
698 : ! do_last_ortho in case of diagonalization scheme
699 : if ( wvlbigdft) do_last_ortho=(dtset%iscf/=0)
700 : if (.not.wvlbigdft) do_last_ortho=(.true.)
701 :
702 : ABI_MALLOC(xcart,(3, dtset%natom))
703 : call xred2xcart(dtset%natom, rprimd, xcart, xred)
704 :
705 : if(wvlbigdft) then
706 : ! NSCF loop for wvlbigdt:
707 : call wvl_nscf_loop_bigdft()
708 : else
709 : ! NSCF loop for WVL: (not wvlbigdft)
710 : call wvl_nscf_loop()
711 : end if
712 :
713 : ! Eventually orthogonalize WFs now
714 : if (do_last_ortho) then
715 : call write_energies(ii,0,wvl%e%energs,0.d0,0.d0,"final")
716 : call last_orthon(me_distrb, nproc_distrb, ii, wvl%wfs%ks, wvl%e%energs%evsum, .true.)
717 : if(wvlbigdft) energies%e_xcdc = wvl%e%energs%evxc
718 : ! If occupation numbers are not changed...
719 : if (fixed_occ .or. (iscf<0 .and. iscf/=-3)) then
720 : call wvl_comm_eigen()
721 : end if
722 : ! ... or update occupations:
723 : if( ( .not.fixed_occ) .and. (iscf>0.or.iscf==-3)) then
724 : if(wvlbigdft) then
725 : call wvl_occ_bigdft()
726 : else
727 : ! Communicate eigenvalues:
728 : call wvl_comm_eigen()
729 : ! Update occ and Fermi level
730 : call wvl_occ()
731 : end if
732 : end if
733 : ! This might accelerate convergence:
734 : wvl%wfs%ks%diis%energy_min=one
735 : wvl%wfs%ks%diis%alpha=two
736 : end if !do_last_ortho
737 :
738 : ! Compute eigenvalues energy
739 : if(.not. wvlbigdft .and. nnsclo_now>0) then
740 : call e_eigen(eigen,energies%e_eigenvalues,dtset%mband,dtset%nband,dtset%nkpt,&
741 : & dtset%nsppol,occ,dtset%wtk)
742 : else
743 : energies%e_eigenvalues = energies%e_kinetic + energies%e_localpsp &
744 : & + energies%e_xcdc + two*energies%e_hartree +energies%e_nlpsp_vfock
745 : end if
746 :
747 : if (optforces == 1) then ! not compatible with iscf=0 and wvlbigdftcomp=1 + PAW
748 : call wvl_nl_gradient(grnl, mpi_enreg, dtset%natom, rprimd, wvl, xcart)
749 : end if
750 :
751 : ! For iscf<0 we do not update the density
752 : if (dtset%iscf>=0) then !(dtset%iscf>=0 .and. .not. wvlbigdft ) then
753 : call wvl_mkrho(dtset,irrzon,mpi_enreg,phnons,rhor,wvl%wfs,wvl%den)
754 : end if
755 : ABI_FREE(xcart)
756 :
757 : ! Note in WVL+NC: the rest will be skipped.
758 : ! For PAW: we will compute Rho_ij at the end.
759 : !if(wvlbigdft) return
760 : #endif
761 : else
762 :
763 : !===================================================================
764 : ! PLANE WAVES - Standard VTORHO procedure
765 : !===================================================================
766 :
767 : ! Electric field: allocate dphasek
768 46993 : nkpt1 = dtset%nkpt
769 46993 : if ( berryflag ) then
770 525 : ABI_MALLOC(dphasek,(3,dtset%nkpt*dtset%nsppol))
771 50519 : dphasek(:,:) = zero
772 175 : nkpt1 = dtefield%mkmem_max
773 : end if
774 :
775 46993 : if(dtset%gpu_option==ABI_GPU_KOKKOS) then
776 : #if defined HAVE_GPU && defined HAVE_YAKL
777 : ABI_MALLOC_MANAGED(rhoaug, (/n4,n5,n6,gs_hamk%nvloc/))
778 : ABI_MALLOC_MANAGED(vlocal, (/n4,n5,n6,gs_hamk%nvloc/))
779 : #endif
780 : else
781 281958 : ABI_MALLOC(rhoaug,(n4,n5,n6,gs_hamk%nvloc))
782 234965 : ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamk%nvloc))
783 : end if
784 :
785 46993 : if(with_vxctau) then
786 4536 : ABI_MALLOC(vxctaulocal,(n4,n5,n6,gs_hamk%nvloc,4))
787 : end if
788 :
789 46993 : has_vectornd = (with_vectornd .EQ. 1)
790 46993 : if(has_vectornd) then
791 2891 : ABI_MALLOC(vectornd_pac,(n4,n5,n6,gs_hamk%nvloc,3))
792 99602741 : vectornd_pac=zero
793 : end if
794 :
795 46993 : nbdbuf_eff = dtset%nbdbuf
796 : ! In metallic case, at first iteration, occupations could be 0. So residm should be computed as usual
797 3113378 : if (dtset%nbdbuf==-101.and..not.fixed_occ.and.istep==1.and.minval(occ)<tol10) then
798 1 : write(msg,*) 'vtorho: nbdbuf is set to 0 for this step'
799 1 : call wrtout(std_out,msg,'COLL')
800 1 : nbdbuf_eff = 0
801 : end if
802 :
803 : ! LOOP OVER SPINS
804 98669 : do isppol=1,dtset%nsppol
805 51676 : call timab(982,1,tsec)
806 :
807 51676 : ikpt_loc = 0
808 51676 : ikg=0
809 :
810 51676 : if (dtset%cprj_in_memory==1) then
811 2466 : if (xg_nonlop%paw) call xg_nonlop_set_Dij_spin(xg_nonlop,isppol)
812 : end if
813 : ! Set up local potential vlocal on the coarse FFT mesh from vtrial taking into account the spin.
814 : ! Also, continue to initialize the Hamiltonian.
815 :
816 : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
817 51676 : dtset%nspden, gs_hamk%nvloc, 1, pawfgr, mpi_enreg, vtrial, vlocal)
818 : !vlocal = zero
819 51676 : call gs_hamk%load_spin(isppol, vlocal=vlocal, with_nonlocal=.true.)
820 :
821 51676 : if (with_vxctau) then
822 : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
823 671 : dtset%nspden, gs_hamk%nvloc, 4, pawfgr, mpi_enreg, vxctau, vxctaulocal)
824 671 : call gs_hamk%load_spin(isppol, vxctaulocal=vxctaulocal)
825 : end if
826 :
827 963632980 : rhoaug(:,:,:,:)=zero
828 :
829 : ! if vectornd is present, set it up for addition to gs_hamk similarly to how it's done for
830 : ! vtrial. Note that it must be done for the three Cartesian directions. Also, the following
831 : ! code assumes explicitly and implicitly that nvloc = 1. This should eventually be generalized.
832 51676 : if (has_vectornd) then
833 : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, &
834 423 : & dtset%nspden, gs_hamk%nvloc, 3, pawfgr, mpi_enreg, vectornd, vectornd_pac)
835 423 : call gs_hamk%load_spin(isppol, vectornd=vectornd_pac)
836 : end if
837 :
838 51676 : call timab(982,2,tsec)
839 :
840 : ! BIG FAT k POINT LOOP
841 : ! MVeithen: I had to modify the structure of this loop in order to implement MPI // of the electric field
842 : ! Note that the loop here differs from the similar one in berryphase_new.F90.
843 : ! here, ikpt_loc numbers the kpts treated by the current processor.
844 : ! in berryphase_new.F90, ikpt_loc ALSO includes info about value of isppol.
845 :
846 51676 : ikpt = 0
847 490388 : do while (ikpt_loc < nkpt1)
848 :
849 438712 : call timab(997,1,tsec)
850 :
851 438712 : if ( .not.berryflag ) then
852 429182 : ikpt_loc = ikpt_loc + 1
853 429182 : ikpt = ikpt_loc
854 429182 : my_ikpt = mpi_enreg%my_kpttab(ikpt)
855 : else
856 9530 : if (ikpt_loc < dtset%mkmem) ikpt = ikpt + 1
857 9530 : if ((ikpt > dtset%nkpt).and.(ikpt_loc < dtset%mkmem)) exit
858 9530 : my_ikpt=ikpt
859 : end if
860 :
861 438712 : dphase_k(:) = zero
862 438712 : counter=100*ikpt+isppol
863 438712 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
864 438712 : nband_cprj_k=nband_k/mpi_enreg%nproc_band
865 438712 : istwf_k=dtset%istwfk(ikpt)
866 438712 : npw_k=npwarr(ikpt)
867 :
868 438712 : mcgq=1 ; mkgq=1
869 438712 : if (.not.berryflag) then
870 429182 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
871 460275 : eigen(1+bdtot_index : nband_k+bdtot_index) = zero
872 460275 : resid(1+bdtot_index : nband_k+bdtot_index) = zero
873 : bdtot_index=bdtot_index+nband_k
874 : cycle
875 : end if
876 : else
877 9530 : if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) .and.(ikpt_loc <= dtset%mkmem)) then
878 21280 : eigen(1+bdtot_index : nband_k+bdtot_index) = zero
879 21280 : resid(1+bdtot_index : nband_k+bdtot_index) = zero
880 : bdtot_index = bdtot_index + nband_k
881 : cycle
882 : end if
883 6474 : ikpt_loc = ikpt_loc + 1
884 6474 : mcgq = dtset%mpw*my_nspinor*nband_k*dtefield%nneigh(ikpt)
885 6474 : ikg = dtefield%kgindex(ikpt)
886 6474 : mkgq = 6*dtset%mpw
887 : end if ! berryflag
888 :
889 343291 : call timab(997,2,tsec)
890 :
891 : ! In case of MPI // of a finite field calculation
892 : ! build the cgq array that stores the wavefunctions for the
893 : ! neighbours of ikpt, and the pwnsfacq array that stores the
894 : ! corresponding phase factors (in case of tnons)
895 1029873 : ABI_MALLOC(cgq,(2,mcgq))
896 1029873 : ABI_MALLOC(pwnsfacq,(2,mkgq))
897 343291 : if ( berryflag ) then
898 : call cgq_builder(berryflag,cg,cgq,dtefield,dtset,ikpt,ikpt_loc,isppol,mcg,mcgq,&
899 : me_distrb,mkgq,mpi_enreg,my_nspinor,nband_k,nproc_distrb,&
900 6474 : npwarr,pwnsfac,pwnsfacq,pwind_alloc,spaceComm_distrb)
901 6474 : if (ikpt_loc > dtset%mkmem) then
902 0 : ABI_FREE(cgq)
903 0 : ABI_FREE(pwnsfacq)
904 0 : cycle
905 : end if
906 : end if !berryopt
907 :
908 343291 : call timab(984,1,tsec)
909 :
910 343291 : if (mpi_enreg%paral_kgb==1) my_bandfft_kpt => bandfft_kpt(my_ikpt)
911 343291 : call bandfft_kpt_set_ikpt(ikpt,mpi_enreg)
912 : !my_ikpt = ikpt
913 : !if (mpi_enreg%paral_kgb==1) then
914 : ! my_ikpt = mpi_enreg%my_kpttab(ikpt)
915 : ! my_bandfft_kpt => bandfft_kpt(my_ikpt)
916 : ! call bandfft_kpt_set_ikpt(ikpt,mpi_enreg)
917 : !end if
918 :
919 1029873 : ABI_MALLOC(ek_k,(nband_k))
920 1373164 : ABI_MALLOC(ek_k_nd,(2,nband_k,nband_k*paw_dmft%use_dmft))
921 686582 : ABI_MALLOC(end_k,(nband_k))
922 686582 : ABI_MALLOC(enlx_k,(nband_k))
923 686582 : ABI_MALLOC(occ_k,(nband_k))
924 686582 : ABI_MALLOC(zshift,(nband_k))
925 1373164 : ABI_MALLOC(grnl_k,(3*natom,nband_k*optforces))
926 :
927 343291 : if(dtset%gpu_option==ABI_GPU_KOKKOS) then
928 : #if defined HAVE_GPU && defined HAVE_YAKL
929 : ABI_MALLOC_MANAGED(eig_k,(/nband_k/))
930 : ABI_MALLOC_MANAGED(resid_k,(/nband_k/))
931 : #endif
932 : else
933 686582 : ABI_MALLOC(eig_k,(nband_k))
934 686582 : ABI_MALLOC(resid_k,(nband_k))
935 : end if
936 :
937 2956768 : eig_k(:)=zero
938 2956768 : ek_k(:)=zero
939 2956768 : end_k(:)=zero
940 2956768 : enlx_k(:)=zero
941 1892755 : if(paw_dmft%use_dmft==1) ek_k_nd(:,:,:)=zero
942 6144784 : if (optforces>0) grnl_k(:,:)=zero
943 1373164 : kpoint(:)=dtset%kptns(:,ikpt)
944 :
945 343291 : if (dtset%use_gbt /= 0) then
946 : ! If GBT is activate, kpoint becomes k-q/2 so that we can reuse all the calls
947 : ! to mkkin and mkffnl, and we only have to deal with k+q/2.
948 3104 : kphq = kpoint + half * dtset%qgbt
949 3104 : kpoint(:) = dtset%kptns(:,ikpt) - half * dtset%qgbt
950 : end if
951 :
952 2956768 : occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
953 2956768 : resid_k(:) = resid(1+bdtot_index : nband_k+bdtot_index)
954 : !resid_k(:)=zero
955 2956768 : zshift(:)=dtset%eshift
956 :
957 343291 : if(dtset%gpu_option==ABI_GPU_KOKKOS) then
958 : #if defined HAVE_GPU && defined HAVE_YAKL
959 : ABI_MALLOC_MANAGED(kg_k, (/3,npw_k/))
960 : #endif
961 : else
962 1029873 : ABI_MALLOC(kg_k,(3,npw_k))
963 : end if
964 :
965 1373164 : ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
966 299702359 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
967 343291 : if (psps%useylm==1) then
968 751758 : do ilm=1,psps%mpsang*psps%mpsang
969 143085941 : ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
970 : end do
971 : end if
972 :
973 : ! Set up remaining of the Hamiltonian
974 : ! Compute (1/2) (2 Pi)**2 (k+G)**2:
975 343291 : if(dtset%gpu_option==ABI_GPU_KOKKOS) then
976 : #if defined HAVE_GPU && defined HAVE_YAKL
977 : ABI_MALLOC_MANAGED(kinpw,(/npw_k/))
978 : #endif
979 : else
980 1029873 : ABI_MALLOC(kinpw,(npw_k))
981 : end if
982 :
983 343291 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw,kpoint,npw_k,0,0)
984 :
985 : ! Compute (k+G) vectors (only if useylm=1)
986 343291 : if (dtset%cprj_in_memory/=1) then
987 326207 : nkpg=3*optforces*dtset%nloalg(3)
988 1304828 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
989 326207 : if ((mpi_enreg%paral_kgb/=1.or.istep<=1).and.nkpg>0) then
990 9334 : call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
991 9334 : if (dtset%use_gbt /= 0) then
992 0 : ABI_MALLOC(kpg_kphq,(npw_k,nkpg))
993 0 : call mkkpg(kg_k,kpg_kphq,kphq,nkpg,npw_k)
994 : end if
995 : end if
996 : else ! cprj_in_memory = 1
997 17084 : nkpg=3*optforces
998 68336 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
999 17084 : if (optforces/=0) then
1000 9664 : call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
1001 9664 : if (dtset%use_gbt /= 0) then
1002 0 : ABI_MALLOC(kpg_kphq,(npw_k,nkpg))
1003 0 : call mkkpg(kg_k,kpg_kphq,kphq,nkpg,npw_k)
1004 : end if
1005 : end if
1006 : end if
1007 :
1008 : ! Compute nonlocal form factors ffnl at all (k+G):
1009 343291 : ider=0;idir=0;dimffnl=1
1010 :
1011 1716455 : ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat))
1012 343291 : calc_ffnl_ph3d=.false.
1013 343291 : if(mpi_enreg%paral_kgb/=1.or.istep<=1) calc_ffnl_ph3d=.true.
1014 343291 : if(associated(rcpaw)) then
1015 21 : if(rcpaw%istep>=rcpaw%updatepaw(1)+1.and.rcpaw%istep<=rcpaw%updatepaw(2)+1.and.&
1016 : (dtset%wfoptalg==111.or.dtset%wfoptalg == 1)) then
1017 0 : calc_ffnl_ph3d=.true.
1018 : endif
1019 : endif
1020 343291 : if (calc_ffnl_ph3d) then
1021 : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,&
1022 : gmet,gprimd,ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,&
1023 : psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
1024 : npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,&
1025 310427 : psps%usepaw,psps%useylm,ylm_k,ylmgr,kinpw=kinpw)
1026 : end if
1027 :
1028 343291 : if (dtset%use_gbt /= 0) then
1029 : ! Compute (1/2) (2 Pi)**2 (k+q/2+G)**2:
1030 2328 : ABI_MALLOC(kinpw_kphq, (npw_k))
1031 776 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw_kphq,kphq,npw_k,0,0)
1032 :
1033 : ! Compute nonlocal form factors ffnl at all (k+q/2+G):
1034 : ! TODO: useylm = 1 requires ylm_kphq, ylmgr_kphq
1035 : ! and ylm should be computed with k-q/2 in the caller.
1036 3104 : ABI_MALLOC(ffnl_kphq,(npw_k,dimffnl,psps%lmnmax,ntypat))
1037 : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl_kphq,psps%ffspl,&
1038 : gmet,gprimd,ider,idir,psps%indlmn,kg_k,kpg_kphq,kphq,psps%lmnmax,&
1039 : psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
1040 : npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,&
1041 776 : psps%usepaw,psps%useylm,ylm_k,ylmgr,kinpw=kinpw_kphq)
1042 : end if
1043 :
1044 : ! Load k-dependent part in the Hamiltonian datastructure
1045 : ! - Compute 3D phase factors
1046 : ! - Prepare various tabs in case of band-FFT parallelism
1047 : ! - Load k-dependent quantities in the Hamiltonian
1048 1373164 : ABI_MALLOC(ph3d,(2,npw_k,gs_hamk%matblk))
1049 343291 : if (dtset%use_gbt /= 0) then
1050 2328 : ABI_MALLOC(ph3d_kphq,(2,npw_k,gs_hamk%matblk))
1051 : end if
1052 :
1053 343291 : if (usefock_ACE/=0) then
1054 : call gs_hamk%load_k(kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,&
1055 : kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,fockACE_k=fock%fockACE(ikpt,isppol),ph3d_k=ph3d,&
1056 6063 : compute_ph3d=calc_ffnl_ph3d, compute_gbound=(mpi_enreg%paral_kgb/=1))
1057 :
1058 6063 : if (dtset%use_gbt /= 0) then
1059 0 : ABI_ERROR("GBT with fock_ace not implemented")
1060 : end if
1061 : else
1062 : call gs_hamk%load_k(kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,&
1063 : kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,&
1064 337228 : compute_ph3d=calc_ffnl_ph3d, compute_gbound=(mpi_enreg%paral_kgb/=1))
1065 :
1066 337228 : if (dtset%use_gbt /= 0) then
1067 : call gs_hamk%load_kprime(kpt_kp=kphq,&
1068 : kinpw_kp=kinpw_kphq,kpg_kp=kpg_kphq,ffnl_kp=ffnl_kphq,ph3d_kp=ph3d_kphq,&
1069 776 : compute_ph3d=calc_ffnl_ph3d, compute_gbound=(mpi_enreg%paral_kgb/=1))
1070 : end if
1071 : end if
1072 :
1073 : ! Load band-FFT tabs (transposed k-dependent arrays)
1074 343291 : if (mpi_enreg%paral_kgb==1) then
1075 36640 : if (istep<=1) call prep_bandfft_tabs(gs_hamk,ikpt,dtset%mkmem,mpi_enreg)
1076 : call gs_hamk%load_k(npw_fft_k=my_bandfft_kpt%ndatarecv, &
1077 : gbound_k =my_bandfft_kpt%gbound, &
1078 : kinpw_k =my_bandfft_kpt%kinpw_gather, &
1079 : kg_k =my_bandfft_kpt%kg_k_gather, &
1080 : kpg_k =my_bandfft_kpt%kpg_k_gather, &
1081 : ffnl_k =my_bandfft_kpt%ffnl_gather, &
1082 36640 : ph3d_k =my_bandfft_kpt%ph3d_gather)
1083 : end if
1084 :
1085 : ! If OpenMP GPU, load "hamiltonian" on GPU device
1086 343291 : if (gs_hamk%gpu_option == ABI_GPU_OPENMP) then
1087 0 : if(dtset%paral_kgb==0) then
1088 0 : call ompgpu_load_hamilt_buffers(gs_hamk%kg_k,gs_hamk%kg_kp,gs_hamk%ffnl_k,gs_hamk%ph3d_k)
1089 0 : else if(istwf_k==1) then
1090 : call ompgpu_load_hamilt_buffers(gs_hamk%kg_k,gs_hamk%kg_kp,gs_hamk%ffnl_k,gs_hamk%ph3d_k,&
1091 0 : kg_k_gather=bandfft_kpt(my_ikpt)%kg_k_gather)
1092 : else
1093 : call ompgpu_load_hamilt_buffers(gs_hamk%kg_k,gs_hamk%kg_kp,gs_hamk%ffnl_k,gs_hamk%ph3d_k,&
1094 0 : kg_k_gather=bandfft_kpt(my_ikpt)%kg_k_gather_sym)
1095 : end if
1096 : end if
1097 :
1098 343291 : if(gemm_nonlop_use_gemm .and. istep <= 1 .and. dtset%gpu_option==ABI_GPU_OPENMP) then
1099 0 : gemm_nonlop_block_size = dtset%gpu_nl_splitsize
1100 : call get_gemm_nonlop_ompgpu_blocksize(ikpt,gs_hamk,mpi_enreg%bandpp,nband_k,&
1101 : & dtset%nspinor,dtset%nspden,mpi_enreg%paral_kgb,mpi_enreg%nproc_band,&
1102 : & 0,0,dtset%wfoptalg,gs_hamk%gpu_option,(dtset%gpu_nl_distrib/=0),&
1103 : & gemm_nonlop_block_size,nblk_gemm_nonlop,hamilt_gpu_nfft_blocks,&
1104 0 : & disable_output=gpu_mem_estimated)
1105 0 : gs_hamk%nfft_blocks = hamilt_gpu_nfft_blocks
1106 0 : gemm_nonlop_is_distributed = (dtset%gpu_nl_distrib/=0 .and. nblk_gemm_nonlop > 0)
1107 0 : gpu_mem_estimated=.true.
1108 : end if
1109 :
1110 : ! Build inverse of overlap matrix for chebfi or slice
1111 343291 : if(associated(rcpaw)) then
1112 : step_cond=istep<=1.or.(rcpaw%istep>=rcpaw%updatepaw(1)+1.and.rcpaw%istep<=rcpaw%updatepaw(2)+1.and.&
1113 24 : (dtset%wfoptalg==111.or.dtset%wfoptalg == 1))
1114 : else
1115 343270 : step_cond=istep <= 1
1116 : endif
1117 :
1118 343291 : if (dtset%cprj_in_memory==0) then
1119 : if(psps%usepaw == 1 .and. (dtset%wfoptalg == 1 .or. dtset%wfoptalg == 111 .or. dtset%wfoptalg == 112) &
1120 326207 : & .and. step_cond) then
1121 475 : call make_invovl(gs_hamk, dimffnl, ffnl, ph3d, mpi_enreg)
1122 : end if
1123 : end if
1124 :
1125 : ! Setup gemm_nonlop
1126 343291 : if (gemm_nonlop_use_gemm) then
1127 : call set_gemm_nonlop_ikpt(my_ikpt,gs_hamk%npw_fft_k,gs_hamk%istwf_k,gs_hamk%indlmn,&
1128 660 : & gs_hamk%ntypat,gs_hamk%nattyp,gs_hamk%gpu_option)
1129 : end if
1130 :
1131 : #if defined HAVE_GPU_CUDA
1132 : if (gs_hamk%gpu_option==ABI_GPU_LEGACY .or. gs_hamk%gpu_option==ABI_GPU_KOKKOS) then
1133 : if (mpi_enreg%paral_kgb==1) then
1134 : ph3d_size=INT(size(my_bandfft_kpt%ph3d_gather,dim=1),c_int64_t) &
1135 : & * size(my_bandfft_kpt%ph3d_gather,dim=2) * size(my_bandfft_kpt%ph3d_gather,dim=3)
1136 : call gpu_update_ffnl_ph3d( &
1137 : & my_bandfft_kpt%ph3d_gather, ph3d_size, &
1138 : & my_bandfft_kpt%ffnl_gather, INT(size(my_bandfft_kpt%ffnl_gather),c_int64_t) )
1139 : else
1140 : ph3d_size=INT(size(ph3d,dim=1),c_int64_t)*size(ph3d,dim=2)*size(ph3d,dim=3)
1141 : call gpu_update_ffnl_ph3d( &
1142 : & ph3d, ph3d_size, &
1143 : & ffnl, INT(size(ffnl),c_int64_t) )
1144 : end if
1145 : end if
1146 : #endif
1147 :
1148 343291 : call timab(984,2,tsec)
1149 :
1150 : ! Update the value of ikpt,isppol in fock_exchange and allocate the memory space to perform HF calculation.
1151 343291 : if (usefock) call fock_updateikpt(fock%fock_common,ikpt,isppol)
1152 343291 : if (psps%usepaw==1 .and. usefock) then
1153 596 : if ((fock%fock_common%optfor).and.(usefock_ACE==0)) fock%fock_common%forces_ikpt=zero
1154 : end if
1155 :
1156 343291 : if (dtset%cprj_in_memory==1) then
1157 17084 : do_invS=xg_nonlop%paw.and.(dtset%wfoptalg==111.or.dtset%wfoptalg==112)
1158 : call xg_nonlop_make_k(xg_nonlop,my_ikpt,istwf_k,mpi_enreg%me_g0,mpi_enreg%me_g0_fft,npw_k,ffnl,ph3d,kpg_k,&
1159 17084 : & step_cond,compute_invS_approx=do_invS,compute_gram=do_invS)
1160 : end if
1161 :
1162 : ! Here we initialize the wavefunctions with atomic orbitals at the first GS iteration of the first
1163 : ! relaxation step (if any).
1164 : ! NB: Not all the cases are presently supported.
1165 : !print *, "istep, itimes(1), wfinit", istep, itimes(1), dtset%wfinit
1166 : ! FIXME: This check is not enough as I need to check whether cg have been read from WFK file
1167 343291 : if (istep == 1 .and. itimes(1) == 0 .and. dtset%wfinit /= 0) then
1168 : call cg_from_atoms(ikpt, isppol, rprimd, xred, kg_k, cg(:,icg+1:), dtset, psps, eig_k, gs_hamk, &
1169 4 : mpi_enreg, nband_k, npw_k, my_nspinor)
1170 : end if
1171 :
1172 : #ifdef HAVE_OPENMP_OFFLOAD
1173 : !$OMP TARGET UPDATE TO(cg) IF(transfer_cg_once)
1174 : #endif
1175 :
1176 : ABI_NVTX_START_RANGE(NVTX_VTOWFK)
1177 : ! Compute the eigenvalues, wavefunction, residuals,
1178 : ! contributions to kinetic energy, nuclear dipole energy, nonlocal energy, forces,
1179 : ! and update of rhor to this k-point and this spin polarization.
1180 :
1181 : call vtowfk(cg,cgq,cprj,cpus,dphase_k,dtefield,dtfil,&
1182 : dtset,eig_k,ek_k,ek_k_nd,end_k,enlx_k,fixed_occ,grnl_k,gs_hamk,&
1183 : ibg,icg,ikpt,iscf,isppol,kg_k,kinpw,mband_cprj,mcg,mcgq,mcprj_local,mkgq,&
1184 : mpi_enreg,dtset%mpw,natom,nband_k,nbdbuf_eff,dtset%nkpt,istep,nnsclo_now,npw_k,npwarr,&
1185 : occ_k,optforces,prtvol,pwind,pwind_alloc,pwnsfac,pwnsfacq,resid_k,&
1186 343291 : rhoaug,paw_dmft,dtset%wtk(ikpt),xg_nonlop,zshift, rmm_diis_status(:,ikpt,isppol))
1187 : ABI_NVTX_END_RANGE()
1188 :
1189 : ! LB-01/03/2024: Very weird compiler error on eos-nvhpc23.1 if the second call of timab(985,...) is included...
1190 : ! Drastic short-term solution : disable this timing for nvhpc... In fact this part is not important unless fock is activated
1191 : ! Note: Should we keep nvhpc-23.1 in eos?
1192 : #ifndef FC_NVHPC
1193 343291 : call timab(985,1,tsec)
1194 : #endif
1195 :
1196 : #if defined HAVE_GPU_CUDA
1197 : if(gs_hamk%gpu_option==ABI_GPU_LEGACY .or. gs_hamk%gpu_option==ABI_GPU_KOKKOS) call gpu_finalize_ffnl_ph3d()
1198 : #endif
1199 343291 : ABI_FREE(ffnl)
1200 :
1201 343291 : if(dtset%gpu_option==ABI_GPU_KOKKOS) then
1202 : #if defined HAVE_GPU && defined HAVE_YAKL
1203 : ABI_FREE_MANAGED(kinpw)
1204 : ABI_FREE_MANAGED(kg_k)
1205 : #endif
1206 : else
1207 343291 : ABI_FREE(kinpw)
1208 343291 : ABI_FREE(kg_k)
1209 343291 : ABI_SFREE(kinpw_kphq)
1210 343291 : ABI_SFREE(ffnl_kphq)
1211 : end if
1212 :
1213 343291 : ABI_FREE(kpg_k)
1214 343291 : ABI_SFREE(kpg_kphq)
1215 343291 : ABI_FREE(ylm_k)
1216 343291 : ABI_FREE(ph3d)
1217 343291 : ABI_SFREE(ph3d_kphq)
1218 343291 : ABI_FREE(cgq)
1219 343291 : ABI_FREE(pwnsfacq)
1220 :
1221 343291 : if (berryflag) then
1222 : ! electric field
1223 25896 : dphasek(:,ikpt + (isppol - 1)*dtset%nkpt) = dphase_k(:)
1224 :
1225 : ! The overlap matrices for all first neighbours of ikpt are no more up to date
1226 25896 : do idir = 1, 3
1227 64740 : do ifor = 1, 2
1228 38844 : ikpt1 = dtefield%ikpt_dk(dtefield%i2fbz(ikpt),ifor,idir)
1229 38844 : ikpt1 = dtefield%indkk_f2ibz(ikpt1,1)
1230 38844 : ifor1 = -1*ifor + 3 ! ifor = 1 -> ifor1 = 2 & ifor = 2 -> ifor1 = 1
1231 275322 : dtefield%sflag(:,ikpt1+(isppol-1)*dtset%nkpt,ifor1,idir) = 0
1232 : end do
1233 : end do
1234 : end if ! berryflag
1235 :
1236 : ! Save eigenvalues (hartree), residuals (hartree**2)
1237 2956768 : eigen(1+bdtot_index : nband_k+bdtot_index) = eig_k(:)
1238 2956768 : eknk (1+bdtot_index : nband_k+bdtot_index) = ek_k (:)
1239 343291 : if(usefock) then
1240 35071 : focknk (1+bdtot_index : nband_k+bdtot_index) = fock%fock_common%eigen_ikpt (:)
1241 19059 : if (optforces>0) fockfornk(:,:,1+bdtot_index : nband_k+bdtot_index) = fock%fock_common%forces_ikpt(:,:,:)
1242 : end if
1243 1892755 : if(paw_dmft%use_dmft==1) eknk_nd(:,:,:,ikpt,isppol) = ek_k_nd(:,:,:)
1244 2956768 : resid(1+bdtot_index : nband_k+bdtot_index) = resid_k(:)
1245 6144784 : if (optforces>0) grnlnk(:,1+bdtot_index : nband_k+bdtot_index) = grnl_k(:,:)
1246 2956768 : enlxnk(1+bdtot_index : nband_k+bdtot_index) = enlx_k(:)
1247 :
1248 343291 : if(iscf>0 .or. iscf==-3)then
1249 : ! Accumulate sum over k points for band, nonlocal and kinetic energies,
1250 : ! also accumulate gradients of Enonlocal:
1251 2820239 : do iband=1,nband_k
1252 2820239 : if (abs(occ_k(iband))>tol8) then
1253 1590689 : energies%e_kinetic = energies%e_kinetic + dtset%wtk(ikpt)*occ_k(iband)*ek_k(iband)
1254 1590689 : energies%e_nucdip = energies%e_nucdip + dtset%wtk(ikpt)*occ_k(iband)*end_k(iband)
1255 1590689 : energies%e_eigenvalues = energies%e_eigenvalues + dtset%wtk(ikpt)*occ_k(iband)*eig_k(iband)
1256 1590689 : energies%e_nlpsp_vfock = energies%e_nlpsp_vfock + dtset%wtk(ikpt)*occ_k(iband)*enlx_k(iband)
1257 4922987 : if (optforces>0) grnl(:)=grnl(:)+dtset%wtk(ikpt)*occ_k(iband)*grnl_k(:,iband)
1258 1590689 : if (usefock) then
1259 22651 : energies%e_fock=energies%e_fock + half*fock%fock_common%eigen_ikpt(iband)*occ_k(iband)*dtset%wtk(ikpt)
1260 22651 : if (usefock_ACE==0) energies%e_fock0=energies%e_fock
1261 : endif
1262 : end if
1263 : end do
1264 :
1265 : ! Calculate Fock contribution to the total energy if required
1266 329209 : if ((psps%usepaw==1).and.(usefock)) then
1267 596 : if ((fock%fock_common%optfor).and.(usefock_ACE==0)) then
1268 : !WARNING : this routine actually does NOT compute the Fock contrib to total energy, but modifies the force ONLY.
1269 0 : call fock_calc_ene(dtset,fock%fock_common,energies%e_exactX,ikpt,nband_k,occ_k)
1270 : end if
1271 : end if
1272 : end if
1273 :
1274 343291 : if ( gs_hamk%gpu_option == ABI_GPU_OPENMP) then
1275 0 : call ompgpu_free_hamilt_buffers()
1276 : end if
1277 : ! LB-01/03/2024: Very weird compiler error on eos-nvhpc23.1 if the second call of timab(985,...) is included...
1278 : ! Drastic short-term solution : disable this timing for nvhpc... In fact this part is not important unless fock is activated
1279 : #ifndef FC_NVHPC
1280 343291 : call timab(985,2,tsec)
1281 : #endif
1282 :
1283 343291 : ABI_FREE(ek_k)
1284 343291 : ABI_FREE(ek_k_nd)
1285 343291 : ABI_FREE(end_k)
1286 343291 : ABI_FREE(grnl_k)
1287 343291 : ABI_FREE(occ_k)
1288 343291 : ABI_FREE(zshift)
1289 343291 : ABI_FREE(enlx_k)
1290 :
1291 343291 : if(dtset%gpu_option==ABI_GPU_KOKKOS) then
1292 : #if defined HAVE_GPU && defined HAVE_YAKL
1293 : ABI_FREE_MANAGED(eig_k)
1294 : ABI_FREE_MANAGED(resid_k)
1295 : #endif
1296 : else
1297 343291 : ABI_FREE(eig_k)
1298 343291 : ABI_FREE(resid_k)
1299 : end if
1300 :
1301 : ! Keep track of total number of bands (all k points so far, even for k points not treated by me)
1302 343291 : bdtot_index=bdtot_index+nband_k
1303 :
1304 : ! Also shift array memory if dtset%mkmem/=0
1305 394967 : if (dtset%mkmem/=0) then
1306 343291 : ibg=ibg+my_nspinor*nband_cprj_k
1307 343291 : icg=icg+npw_k*my_nspinor*nband_k
1308 343291 : ikg=ikg+npw_k
1309 : end if
1310 :
1311 : end do ! End big k point loop
1312 :
1313 51676 : call timab(986,1,tsec)
1314 :
1315 51676 : if (fixed_occ .and. mpi_enreg%paral_kgb==1) then
1316 5554 : call xmpi_sum(rhoaug,mpi_enreg%comm_bandspinorfft,ierr) !Sum the contributions over bands/FFT/spinors
1317 : end if
1318 :
1319 : ! Transfer density on augmented fft grid to normal fft grid in real space
1320 : ! Also take into account the spin.
1321 51676 : if(iscf>0.or.iscf==-3)then
1322 51173 : if (psps%usepaw==0) then
1323 34151 : call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug(:,:,:,1),1)
1324 34151 : if(dtset%nspden==4)then
1325 5528 : do imagn=2,4
1326 5528 : call fftpac(imagn,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug(:,:,:,imagn),1)
1327 : end do
1328 : end if
1329 : else
1330 17022 : call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhowfr,rhoaug(:,:,:,1),1)
1331 17022 : if(dtset%nspden==4)then
1332 4324 : do imagn=2,4
1333 4324 : call fftpac(imagn,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhowfr,rhoaug(:,:,:,imagn),1)
1334 : end do
1335 : end if
1336 : end if
1337 : end if
1338 :
1339 98669 : call timab(986,2,tsec)
1340 : end do ! End loop over spins
1341 :
1342 46993 : call timab(988,1,tsec)
1343 :
1344 46993 : if (dtset%cprj_in_memory==1) then
1345 2406 : if (xg_nonlop%paw) call xg_nonlop_destroy_Dij(xg_nonlop)
1346 : end if
1347 :
1348 46993 : if (usefock) then
1349 1375 : if (usefock_ACE==0) then
1350 0 : call xmpi_sum(energies%e_fock0,mpi_enreg%comm_kpt,ierr)
1351 : end if
1352 1375 : if(fock%fock_common%optfor) call xmpi_sum(fock%fock_common%forces,mpi_enreg%comm_kpt,ierr)
1353 : end if
1354 : ! Electric field: compute string-averaged change in Zak phase
1355 : ! along each direction, store it in dphase(idir)
1356 : ! ji: it is not convenient to do this anymore. Remove. Set dphase(idir)=0.0_dp.
1357 : ! eventually, dphase(idir) will have to go...
1358 46993 : if (berryflag) then
1359 175 : dphase(:) = zero
1360 : ! In case of MPI // of a finite field calculation, send dphasek to all cpus
1361 175 : call xmpi_sum(dphasek,spaceComm_distrb,ierr)
1362 175 : ABI_FREE(dphasek)
1363 : end if ! berryflag
1364 :
1365 46993 : if(dtset%gpu_option==ABI_GPU_KOKKOS) then
1366 : #if defined HAVE_GPU && defined HAVE_YAKL
1367 : ABI_FREE_MANAGED(rhoaug)
1368 : ABI_FREE_MANAGED(vlocal)
1369 : #endif
1370 : else
1371 46993 : ABI_FREE(rhoaug)
1372 46993 : ABI_FREE(vlocal)
1373 : end if
1374 :
1375 46993 : if(with_vxctau) then
1376 648 : ABI_FREE(vxctaulocal)
1377 : end if
1378 46993 : if(has_vectornd) then
1379 413 : ABI_FREE(vectornd_pac)
1380 : end if
1381 :
1382 46993 : call timab(988,2,tsec)
1383 :
1384 140979 : ABI_MALLOC(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
1385 3066385 : doccde(:)=zero
1386 :
1387 : ! Treat now varying occupation numbers, in the self-consistent case
1388 46993 : if((.not.fixed_occ) .and. (iscf>0.or.iscf==-3)) then
1389 :
1390 : ! Parallel case
1391 13358 : if (mpi_enreg%nproc_spkpt>1) then
1392 :
1393 1576 : call timab(989,1,tsec)
1394 :
1395 : ! If needed, exchange the values of eigen,resid,eknk,enlxnk,grnlnk
1396 4728 : ABI_MALLOC(buffer1,((4+3*natom*optforces+dtset%usefock+3*natom*dtset%usefock*optforces)*mbdkpsp))
1397 : ! Pack eigen,resid,eknk,enlxnk,grnlnk in buffer1
1398 277696 : buffer1(1 : mbdkpsp)=eigen(:)
1399 277696 : buffer1(1+ mbdkpsp:2*mbdkpsp)=resid(:)
1400 277696 : buffer1(1+2*mbdkpsp:3*mbdkpsp)=eknk(:)
1401 277696 : buffer1(1+3*mbdkpsp:4*mbdkpsp)=enlxnk(:)
1402 1576 : index1=4*mbdkpsp
1403 1576 : if (optforces>0) then
1404 720 : buffer1(index1+1:index1+3*natom*mbdkpsp)=reshape(grnlnk,(/(3*natom)*mbdkpsp/) )
1405 360 : index1=index1+3*natom*mbdkpsp
1406 : end if
1407 1576 : if (usefock) then
1408 20898 : buffer1(1+index1:index1+mbdkpsp)=focknk(:)
1409 258 : if (optforces>0) then
1410 0 : index1=index1+mbdkpsp
1411 0 : buffer1(index1+1:index1+3*natom*mbdkpsp)=reshape(fockfornk,(/(3*natom)*mbdkpsp/) )
1412 : end if
1413 : end if
1414 : ! Build sum of everything
1415 1576 : call timab(48,1,tsec)
1416 1576 : call xmpi_sum(buffer1,mpi_enreg%comm_kpt,ierr)
1417 : ! if(mpi_enreg%paral_kgb/=1.and.paw_dmft%use_dmft==1) then
1418 1576 : if(paw_dmft%use_dmft==1) then
1419 28 : call xmpi_sum(eknk_nd(:,:,:,:,:),mpi_enreg%comm_kpt,ierr)
1420 : end if
1421 1576 : call timab(48,2,tsec)
1422 :
1423 : ! Unpack eigen,resid,eknk,enlxnk,grnlnk in buffer1
1424 277696 : eigen(:) =buffer1(1 : mbdkpsp)
1425 277696 : resid(:) =buffer1(1+ mbdkpsp:2*mbdkpsp)
1426 277696 : eknk(:) =buffer1(1+2*mbdkpsp:3*mbdkpsp)
1427 277696 : enlxnk(:) =buffer1(1+3*mbdkpsp:4*mbdkpsp)
1428 1576 : index1=4*mbdkpsp
1429 1576 : if (optforces>0) then
1430 1080 : grnlnk(:,:)=reshape(buffer1(index1+1:index1+3*natom*mbdkpsp),(/3*natom,mbdkpsp/) )
1431 : end if
1432 1576 : if (usefock) then
1433 20898 : focknk(:)=buffer1(1+index1:index1+mbdkpsp)
1434 258 : if (optforces>0) then
1435 0 : index1=index1+mbdkpsp
1436 0 : fockfornk(:,:,:)=reshape(buffer1(index1+1:index1+3*natom*mbdkpsp),(/3,natom,mbdkpsp/) )
1437 : end if
1438 : end if
1439 1576 : ABI_FREE(buffer1)
1440 3152 : call timab(989,2,tsec)
1441 :
1442 : end if ! nproc_spkpt>1
1443 :
1444 : ! Compute extfpmd energy shift
1445 13358 : if(associated(extfpmd)) then
1446 : call extfpmd%compute_eshift(eigen,eknk,dtset%mband,dtset%nband,&
1447 90 : nfftf,dtset%nkpt,dtset%nsppol,dtset%nspden,dtset%wtk,vtrial)
1448 : end if
1449 :
1450 : ! RCPAW
1451 13358 : nelect=dtset%nelect
1452 13358 : if(associated(rcpaw)) then
1453 21 : if(rcpaw%istep==1) then
1454 102 : min_eigv=minval(eigen)
1455 6 : do itypat=1,dtset%ntypat
1456 6 : if(allocated(rcpaw%atm(itypat)%eig)) then
1457 15 : rcpaw%atm(itypat)%eig=rcpaw%atm(itypat)%eig-rcpaw%atm(itypat)%min_eigv+min_eigv
1458 : endif
1459 : enddo
1460 : endif
1461 21 : nelect=nelect+rcpaw%nelect_core_orig
1462 21 : if(rcpaw%frocc) then
1463 0 : if(rcpaw%istep>rcpaw%updateocc) then
1464 0 : nelect=nelect-rcpaw%nelect_core
1465 : endif
1466 : endif
1467 : endif
1468 :
1469 : ! Compute occupations
1470 13358 : call timab(990,1,tsec)
1471 : call newocc(doccde,eigen,energies%entropy_ks,energies%e_fermie,energies%e_fermih,dtset%ivalence,&
1472 : & dtset%spinmagntarget,dtset%mband,dtset%nband,nelect,dtset%ne_qFD,dtset%nh_qFD,&
1473 : & dtset%nkpt,dtset%nspinor,dtset%nsppol,occ,dtset%occopt,prtvol,dtset%tphysel,&
1474 13358 : & dtset%tsmear,dtset%wtk,prtstm=dtset%prtstm,stmbias=dtset%stmbias,extfpmd=extfpmd,rcpaw=rcpaw)
1475 13358 : call timab(990,2,tsec)
1476 :
1477 : ! !========= DMFT call begin ============================================
1478 : ! ! Also not sure what to do for Wannier90 DMFT
1479 13358 : dmft_dftocc=0
1480 13358 : if(paw_dmft%use_dmft==1.and.psps%usepaw==1.and.dtset%nbandkss==0) then
1481 113 : call timab(991,1,tsec)
1482 :
1483 : ! energies%entropy is the non-interacting entropy. This is obviously
1484 : ! wrong in DFT+DMFT (except if U=J=0), so we set it to 0.
1485 : if (dtset%dmftcheck>=0.and.dtset%usedmft>=1.and.(sum(pawtab(:)%upawu)>=tol8.or. &
1486 599 : & sum(pawtab(:)%jpawu)>tol8).and.dtset%dmft_entropy==0) energies%entropy_ks=zero
1487 :
1488 : ! == 0 to a dmft calculation and do not use lda occupations
1489 : ! == 1 to a lda calculation with the dmft loop
1490 113 : if(dtset%dmftcheck==-1) dmft_dftocc=1
1491 :
1492 : ! == initialise occnd
1493 2155838 : paw_dmft%occnd(:,:,:,:,:) = zero
1494 :
1495 113 : bdtot_index = 1
1496 284 : do isppol=1,dtset%nsppol
1497 1482 : do ikpt=1,dtset%nkpt
1498 29769 : do iband=1,dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1499 28400 : paw_dmft%occnd(1,iband,iband,ikpt,isppol) = occ(bdtot_index)
1500 29598 : bdtot_index = bdtot_index + 1
1501 : end do
1502 : end do
1503 : end do
1504 :
1505 :
1506 113 : if(dmft_dftocc==0) then
1507 93 : if(dtset%occopt/=3) then
1508 0 : ABI_ERROR('occopt should be equal to 3 in dmft')
1509 : end if
1510 : ! == initialize edmft
1511 93 : if(paw_dmft%use_dmft>=1) edmft = zero
1512 :
1513 : ! Compute residm to check the value
1514 93 : ibdkpt=1
1515 93 : residm=zero
1516 224 : do isppol=1,dtset%nsppol
1517 1302 : do ikpt=1,dtset%nkpt
1518 1078 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1519 1078 : if (nbdbuf_eff>=0) then
1520 1078 : nband_eff=max(1,nband_k-nbdbuf_eff)
1521 26796 : residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_eff-1)))
1522 0 : else if (nbdbuf_eff==-101) then
1523 0 : residm=max(residm,maxval(occ(ibdkpt:ibdkpt+nband_k-1)*resid(ibdkpt:ibdkpt+nband_k-1)))
1524 : else
1525 0 : ABI_ERROR('Bad value of nbdbuf_eff')
1526 : end if
1527 1209 : ibdkpt=ibdkpt+nband_k
1528 : end do
1529 : end do
1530 :
1531 : ! Test residm
1532 93 : if (paw_dmft%use_dmft>0 .and. residm>tol4 .and. dtset%dmftcheck>=0) then
1533 1 : if(dtset%dmft_entropy>0) then
1534 : write(msg,'(a,e12.3)')&
1535 1 : ' WARNING: Wavefunctions not converged: DFT+DMFT calculation cannot be carried out safely ',residm
1536 1 : call wrtout(std_out,msg)
1537 : else
1538 : write(msg,'(a,e12.3)')&
1539 0 : ' ERROR: Wavefunctions not converged: DFT+DMFT calculation cannot be carried out safely ',residm
1540 0 : call wrtout(std_out,msg)
1541 0 : write(msg,'(a,i0)')' Action: increase nline (or mdeg_filter) and nnsclo',dtset%nstep
1542 0 : ABI_ERROR(msg)
1543 : end if
1544 :
1545 92 : else if (paw_dmft%use_dmft>0 .and. residm>tol10.and. dtset%dmftcheck>=0) then
1546 37 : write(msg,'(3a)')ch10,&
1547 74 : ' Wavefunctions not converged: DFT+DMFT calculation might not be carried out safely ',ch10
1548 37 : ABI_WARNING(msg)
1549 : end if
1550 :
1551 : ! == gather crystal structure date into data "cryst_struc"
1552 93 : remove_inv=.false.
1553 93 : if(dtset%nspden==4) remove_inv=.true.
1554 : call cryst_struc%init(dtset%amu_orig(:,1),dtset%spgroup,natom,dtset%npsp,ntypat, &
1555 : dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,&
1556 : dtset%nspden==2.and.dtset%nsppol==1,remove_inv,hdr%title,&
1557 186 : dtset%symrel,dtset%tnons,dtset%symafm)
1558 :
1559 : ! == allocate paw_dmft%chipsi and paw_dmft%eigen_dft
1560 : call init_dmft(cryst_struc,dmatpawu(:,:,:,:),dtset,energies%e_fermie,dtfil%filctqmcdatain, &
1561 : & dtfil%filselfin,dtfil%filnam_ds(3),dtfil%fnameabo_app,dtfil%ireadctqmcdata, &
1562 93 : & dtfil%ireadself,paw_dmft,pawtab(:))
1563 93 : call print_dmft(paw_dmft,dtset%pawprtvol)
1564 :
1565 : ! == compute chipsi
1566 93 : call xmpi_barrier(spaceComm_distrb)
1567 93 : call init_oper(paw_dmft,dft_occup,opt_ksloc=2)
1568 93 : call flush_unit(std_out)
1569 93 : call timab(620,1,tsec)
1570 :
1571 : call datafordmft(cg(:,:),cprj(:,:),cryst_struc,dft_occup,gs_hamk%dimcprj(:),dtset,eigen(:),mband_cprj,mcg,&
1572 93 : & mpi_enreg,my_nspinor,occ(:),paw_dmft,paw_ij(:),pawtab(:),usecprj_local,dtset%nbandkss)
1573 :
1574 93 : call timab(620,2,tsec)
1575 93 : call flush_unit(std_out)
1576 :
1577 : ! == solve dmft loop
1578 93 : call xmpi_barrier(spaceComm_distrb)
1579 :
1580 93 : call dmft_solve(cryst_struc,istep,dft_occup,mpi_enreg,paw_dmft,pawang,pawtab(:),dtset%pawprtvol)
1581 93 : edmft=paw_dmft%e_hu-paw_dmft%e_dc
1582 93 : energies%e_dc=paw_dmft%e_dc
1583 93 : energies%e_hu=paw_dmft%e_hu
1584 : if (dtset%dmft_triqs_entropy == 1 .and. dtset%dmft_triqs_compute_integral > 0 &
1585 93 : & .and. (dtset%dmft_solv == 6 .or. dtset%dmft_solv == 7)) then
1586 0 : energies%entropy_ks = paw_dmft%sdmft - paw_dmft%simp
1587 0 : energies%entropy_imp = paw_dmft%simp
1588 : end if
1589 93 : call flush_unit(std_out)
1590 : ! paw_dmft%occnd(:,:,:,:,:)=0.5_dp
1591 :
1592 : ! For compatibility with old test, do not use for calculation
1593 206913 : if(dtset%dmft_occnd_imag==0) paw_dmft%occnd(2,:,:,:,:)=zero
1594 :
1595 : ! call print_dmft(paw_dmft,dtset%pawprtvol)
1596 : ! if(dtset%paral_kgb==1) then
1597 : ! write(msg,'(5a)')ch10,&
1598 : !& ' Parallelization over bands is not yet compatible with self-consistency in DMFT ',ch10,&
1599 : !& ' Calculation of density does not taken into account non diagonal occupations',ch10
1600 : ! call wrtout(std_out,msg)
1601 : ! call wrtout(ab_out,msg)
1602 : !! ABI_ERROR(msg)
1603 : ! if(dtset%nstep>1) then
1604 : ! write(msg,'(a,i0)')' Action: use nstep=1 instead of nstep=',dtset%nstep
1605 : ! ABI_ERROR(msg)
1606 : ! end if
1607 : ! residm=zero
1608 : ! end if
1609 : ! if(dtset%nspinor==2) then
1610 : ! call flush_unit(ab_out)
1611 : ! write(msg,'(3a)')&
1612 : ! & ' Self consistent DFT+DMFT with nspinor==2 is not possible yet ',ch10,&
1613 : ! & ' Calculation are restricted to nstep =1'
1614 : ! ! ABI_ERROR(msg)
1615 : ! if(dtset%nstep>1) then
1616 : ! write(msg,'(a,i0)')' Action: use nstep=1 instead of nstep=',dtset%nstep
1617 : ! ! ABI_ERROR(msg)
1618 : ! endif
1619 : ! end if
1620 :
1621 93 : if(paw_dmft%myproc==0) then
1622 : ABI_NVTX_START_RANGE(NVTX_DMFT_SAVEOCC)
1623 42 : call timab(628,1,tsec)
1624 42 : call saveocc_dmft(paw_dmft)
1625 42 : call timab(628,2,tsec)
1626 : ABI_NVTX_END_RANGE()
1627 : end if
1628 93 : call destroy_dmft(paw_dmft)
1629 :
1630 : ! == destroy crystal_t cryst_struc
1631 93 : call cryst_struc%free()
1632 93 : call destroy_oper(dft_occup)
1633 : end if ! dmft_dftocc
1634 113 : call timab(991,2,tsec)
1635 :
1636 13245 : else if (dtset%usedmft == 10) then
1637 0 : write(std_out, *) "pawrad%int_meshsz: ", pawrad%int_meshsz
1638 0 : write(std_out, *) "dtfil%fnameabo_w90: ", dtfil%fnameabo_w90
1639 0 : write(std_out, *) "results_gs%energies%e_fermie: ", results_gs%energies%e_fermie
1640 : #if defined HAVE_PYTHON_INVOCATION
1641 : ! xcryst_struct
1642 : remove_inv=.false.
1643 : call cryst_struc%init(dtset%amu_orig(:,1),dtset%spgroup,natom,dtset%npsp,ntypat,&
1644 : & dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,&
1645 : & dtset%nspden==2.and.dtset%nsppol==1,remove_inv,hdr%title,&
1646 : & dtset%symrel,dtset%tnons,dtset%symafm)
1647 :
1648 : ! ebands
1649 : bantot = dtset%mband*dtset%nkpt*dtset%nsppol
1650 : call ebands%init(bantot,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%ivalence,&
1651 : & doccde,eigen,hdr%istwfk,hdr%kptns,hdr%nband,&
1652 : & hdr%nkpt,hdr%npwarr,hdr%nsppol,hdr%nspinor,hdr%tphysel,hdr%tsmear,hdr%occopt,hdr%occ,hdr%wtk,&
1653 : & hdr%cellcharge, hdr%kptopt, hdr%kptrlatt_orig, hdr%nshiftk_orig, hdr%shiftk_orig, &
1654 : & hdr%kptrlatt, hdr%nshiftk, hdr%shiftk)
1655 : ebands%fermie = results_gs%energies%e_fermie
1656 : ebands%fermih = results_gs%energies%e_fermih
1657 : ebands%entropy = results_gs%energies%entropy
1658 :
1659 : ! hdr : OK, in func args
1660 : ! atindx1 : OK, in func args
1661 : ! cg : OK, in func args
1662 : ! cprj : OK, in func args
1663 : ! dtset : OK, in func args
1664 : ! dtfil : OK, in func args
1665 : ! eigen : OK, in func args
1666 : ! gprimd : OK, in func args
1667 : ! kg : OK, in func args
1668 : ! mband : OK, from dtset
1669 : ! mcg : OK, in func args
1670 : ! mcprj : OK, in func args
1671 : ! mgfftc : OK, from dtset called mgfft
1672 : ! mkmem : OK, from dtset
1673 : ! mpi_enreg : OK, in func args
1674 : ! mpw : OK, from dtset
1675 : ! natom : OK, from dtset
1676 : ! nattyp : OK, in func args
1677 : ! nfft : OK, from dtset
1678 : ! ngfft : OK, from dtset
1679 : ! nkpt : OK, from dtset
1680 : ! npwarr : OK, was added to vtorho func
1681 : ! nsppol : OK, from dtset
1682 : ! ntypat : OK, from dtset
1683 : ! occ : OK, in func args
1684 : ! pawang : OK, in func args
1685 : ! pawrad : OK, was added to func's args
1686 : ! pawtab : OK, in func args
1687 : ! prtvol : OK, from dtset
1688 : ! psps : OK, in func args
1689 : ! rprimd : OK, in func args
1690 : ! ucvol : OK, in func args
1691 : ! xred : OK, in func args
1692 :
1693 : ! Call Wannier90 and print the _hr.dat file with the Hamiltonian
1694 : write(msg,'(6a)') &
1695 : & ch10, ' ======================================================'&
1696 : & ,ch10,' ===== Calling Wannier90 ========'&
1697 : & ,ch10,' ======================================================'
1698 : call wrtout(std_out,msg,'COLL')
1699 :
1700 : call init_mywfc(mywfc=mywfc, ebands=ebands, cg=cg, cprj=cprj, &
1701 : cryst=cryst_struc, dtset=dtset, dtfil=dtfil, hdr=hdr, &
1702 : MPI_enreg=mpi_enreg, nprocs=mpi_enreg%nproc, psps=psps, pawtab=pawtab, &
1703 : rank=mpi_enreg%me, comm=mpi_enreg%comm_world)
1704 :
1705 : call mlwfovlp(mywfc=mywfc, crystal=cryst_struc, ebands=mywfc%ebands, hdr=mywfc%hdr, &
1706 : atindx1=cryst_struc%atindx1, dtset=mywfc%dtset, dtfil=dtfil, eigen=mywfc%ebands%eig, &
1707 : gprimd=cryst_struc%gprimd, kg=kg, mband=dtset%mband, mcg=mcg, mcprj=mcprj, &
1708 : mgfftc=dtset%mgfft, mkmem=dtset%mkmem, mpi_enreg=mpi_enreg, &
1709 : mpw=dtset%mpw, natom=dtset%natom, nattyp=nattyp, nfft=dtset%nfft, ngfft=dtset%ngfft, &
1710 : nkpt=dtset%nkpt, &
1711 : npwarr=npwarr, nsppol=dtset%nsppol, ntypat=dtset%ntypat, occ=occ, pawang=pawang, &
1712 : pawrad=pawrad, pawtab=pawtab, prtvol=dtset%prtvol, &
1713 : psps=psps, rprimd=rprimd, ucvol=ucvol, xred=xred, exclude_bands=exclude_bands)
1714 : call xmpi_barrier(spaceComm_distrb)
1715 :
1716 : ! Need to know if a band is in or not for the calculation, for correcting occupations and density
1717 : do iband=1, hdr%mband
1718 : do isppol=1, hdr%nsppol
1719 : exclude_bands_ind(iband,isppol) = .true.
1720 : enddo
1721 : enddo
1722 : do iband=1, hdr%mband
1723 : do isppol=1, hdr%nsppol
1724 : if (exclude_bands(iband,isppol) /= 0) then
1725 : exclude_bands_ind(exclude_bands(iband,isppol),isppol) = .false.
1726 : endif
1727 : enddo
1728 : enddo
1729 :
1730 : do isppol=1, hdr%nsppol
1731 : do iband=1, paw_dmft%dmftbandi-1
1732 : if (exclude_bands_ind(iband,isppol) .eqv. .true.) then
1733 : write(msg, '(a)') "ERROR: Inconsistency between the excluded bands for Wannier90 and the dmftbandi/dmftbandf ABINIT keywords. TODO: remove this duplication of information."
1734 : ABI_ERROR(msg)
1735 : endif
1736 : enddo
1737 : do iband=paw_dmft%dmftbandi, paw_dmft%dmftbandf
1738 : if (exclude_bands_ind(iband,isppol) .eqv. .false.) then
1739 : write(msg, '(a)') "ERROR: Inconsistency between the excluded bands for Wannier90 and the dmftbandi/dmftbandf ABINIT keywords. TODO: remove this duplication of information."
1740 : ABI_ERROR(msg)
1741 : endif
1742 : enddo
1743 : do iband=paw_dmft%dmftbandf+1, hdr%mband
1744 : if (exclude_bands_ind(iband,isppol) .eqv. .true.) then
1745 : write(msg, '(a)') "ERROR: Inconsistency between the excluded bands for Wannier90 and the dmftbandi/dmftbandf ABINIT keywords. TODO: remove this duplication of information."
1746 : ABI_ERROR(msg)
1747 : endif
1748 : enddo
1749 : enddo
1750 :
1751 : ! Print data needed for full charge self-consistency
1752 : write(filename, '(a, a)') trim(dtfil%filnam_ds(4)), "_w90.abinit"
1753 : write(msg,'(3a)') " ===== Printing data needed for full charge self-consistency in DFT+DMFT in ", trim(filename), " file"
1754 : call wrtout(std_out, msg, 'COLL')
1755 : open(unit=100, file=filename)
1756 : write(100, "(a, f12.6)") "Fermie", energies%e_fermie
1757 : write(100, "(a, i5)") "Nkpt", dtset%nkpt
1758 : write(100, "(a, i5)") "Nband", dtset%mband
1759 : do ikpt=1,dtset%nkpt
1760 : write(100,'(a,3f10.4,a,i4,a)' ) ' k-point ', kpoint, ' number ',ikpt,' :'
1761 : do ii=0,(dtset%mband-1)/12
1762 : write(100,'(12f10.4)') occ(1+ii*12+(ikpt-1)*dtset%mband:min(12+ii*12,dtset%mband)+(ikpt-1)*dtset%mband)
1763 : end do
1764 : end do
1765 : close(100)
1766 :
1767 : call xmpi_barrier(spaceComm_distrb)
1768 :
1769 : ! Perform DMFT. No need for most of the pipeline in dmft_solve.
1770 : ! We directly call the python_invocation.
1771 : write(msg,'(6a)') &
1772 : & ch10, ' ======================================================'&
1773 : & ,ch10,' ===== DMFT starts here ========'&
1774 : & ,ch10,' ======================================================'
1775 : call wrtout(std_out,msg,'COLL')
1776 : call flush_unit(std_out)
1777 :
1778 : ! Invoking python to execute the script
1779 : write(msg, '(a)') trim(dtfil%filnam_ds(3))
1780 : call invoke_python_run_script (istep, paw_dmft%myproc, msg, mpi_enreg%comm_world)
1781 : call xmpi_barrier(paw_dmft%spacecomm)
1782 : call flush_unit(std_out)
1783 :
1784 : ! Need new DMFT occupations
1785 : ! They can be found in the file dtfil%filenam_ds(4)_w90.deltaN
1786 : ! Check if file exists
1787 : write(filename, '(a, i4.4, 3a)') "dft", istep, "_", trim(dtfil%filnam_ds(4)), '_w90.deltaN'
1788 : inquire(file=filename, exist=exists)
1789 : if (.not.exists) then
1790 : write(msg, '(3a)') " ERROR: The file ", trim(filename), " does not exist. It means there was a problem with the DMFT. Abort."
1791 : ABI_ERROR(msg)
1792 : endif
1793 :
1794 : ABI_MALLOC(occnd_tmp, (2*(paw_dmft%dmftbandf+1-paw_dmft%dmftbandi)))
1795 :
1796 : ! Read the new occupations obtained from the DMFT calculation.
1797 : do isppol=1,dtset%nsppol
1798 : open(unit=101, file=filename, status='old')
1799 : rewind(unit=101)
1800 : read(101, *) msg
1801 : do ikpt=1,dtset%nkpt
1802 : read(101, *) msg
1803 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1804 : do iband=1,nband_k
1805 : bdtot_index = iband+dtset%mband*(ikpt-1)+dtset%mband*dtset%nkpt*(isppol-1)
1806 : paw_dmft%occnd(1,iband,iband,ikpt,isppol) = occ(bdtot_index)
1807 : if (paw_dmft%band_in(iband)) then
1808 : read(101, *) occnd_tmp
1809 : ! Works in the presence of spin-symmetry, at this point
1810 : do iband1=1,paw_dmft%dmftbandf+1-paw_dmft%dmftbandi
1811 : paw_dmft%occnd(1,iband,paw_dmft%dmftbandi-1+iband1,ikpt,isppol) = paw_dmft%occnd(1,iband,paw_dmft%dmftbandi-1+iband1,ikpt,isppol) + occnd_tmp(2*(iband1-1)+1)
1812 : ! if (paw_dmft%dmftbandi-1+iband1==iband) then
1813 : ! paw_dmft%occnd(1,iband,paw_dmft%dmftbandi-1+iband1,ikpt,isppol) = paw_dmft%occnd(1,iband,paw_dmft%dmftbandi-1+iband1,ikpt,isppol) + occ(bdtot_index)
1814 : ! end if
1815 : paw_dmft%occnd(2,iband,paw_dmft%dmftbandi-1+iband1,ikpt,isppol) = + occnd_tmp(2*(iband1-1)+2)
1816 : enddo
1817 : endif
1818 : end do
1819 : end do
1820 : close(101)
1821 : end do
1822 :
1823 : ABI_FREE(occnd_tmp)
1824 :
1825 : write(msg, '(6a)') &
1826 : & ch10, ' ================================================== ',&
1827 : & ch10, ' ===== DMFT : END ======== ',&
1828 : & ch10, ' ================================================== '
1829 : call wrtout(std_out,msg,'COLL')
1830 : call flush_unit(std_out)
1831 :
1832 : #else
1833 0 : ABI_ERROR('Cannot use use_dmft == 10 with #HAVE_PYTHON_INVOCATION set to false.')
1834 : #endif
1835 :
1836 : end if ! usedmft
1837 :
1838 13358 : if(dtset%nbandkss/=0) then
1839 : write(msg,'(a,i3,2a,i3,4a)') &
1840 0 : " dtset%nbandkss = ",dtset%nbandkss,ch10,&
1841 0 : " and dtset%usedmft = ",dtset%usedmft,ch10,&
1842 0 : " a DFT loop is carried out without DMFT.",ch10,&
1843 0 : " Only chipsi's will be written at convergence of the DFT loop."
1844 0 : call wrtout(std_out,msg)
1845 : end if
1846 : ! !========= DMFT call end ============================================
1847 :
1848 13358 : call timab(992,1,tsec)
1849 :
1850 : ! Compute eeig, ek,enl and grnl from the new occ, and the shared eknk,enlxnk,grnlnk
1851 13358 : energies%e_eigenvalues = zero
1852 13358 : energies%e_kinetic = zero
1853 13358 : energies%e_nlpsp_vfock = zero
1854 13358 : if (usefock) then
1855 897 : energies%e_fock = zero
1856 1241 : if (optforces>0) fock%fock_common%forces=zero
1857 : end if
1858 34529 : if (optforces>0) grnl(:)=zero
1859 13358 : if(paw_dmft%use_dmft>=1) then
1860 113 : ebandlda = zero
1861 113 : ebanddmft = zero
1862 113 : ebandldatot = zero
1863 113 : ekindmft = zero
1864 113 : ekindmft2 = zero
1865 113 : ekinlda = zero
1866 : end if
1867 :
1868 : ! Compute new energy terms due to non diagonal occupations and DMFT.
1869 : ! It uses the new occupations stored in paw_dmft%occnd.
1870 13358 : bdtot_index=1
1871 28465 : do isppol=1,dtset%nsppol
1872 141353 : do ikpt=1,dtset%nkpt
1873 112888 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1874 1501395 : do iband=1,nband_k
1875 :
1876 1373400 : locc_test = abs(occ(bdtot_index))>tol8
1877 : ! dmft
1878 1373400 : if(paw_dmft%use_dmft>=1.and.dtset%nbandkss==0) then
1879 28400 : if(paw_dmft%band_in(iband).or.(paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7)) then
1880 13212 : if( paw_dmft%use_dmft == 1 .and. dmft_dftocc == 1 ) then ! test of the code
1881 3360 : paw_dmft%occnd(1,iband,iband,ikpt,isppol) = occ(bdtot_index)
1882 : end if
1883 : locc_test = abs(paw_dmft%occnd(1,iband,iband,ikpt,isppol))+&
1884 13212 : & abs(paw_dmft%occnd(2,iband,iband,ikpt,isppol))>tol8
1885 : end if
1886 : end if
1887 :
1888 1373400 : if (locc_test) then
1889 : ! dmft
1890 826735 : if((paw_dmft%use_dmft==1.or.paw_dmft%use_dmft==10).and.dtset%nbandkss==0) then
1891 21291 : ebandldatot=ebandldatot+dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index)
1892 21291 : if(paw_dmft%band_in(iband).or.(paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7)) then
1893 12963 : ebandlda=ebandlda+dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index)
1894 12963 : ekinlda=ekinlda+dtset%wtk(ikpt)*occ(bdtot_index)*eknk(bdtot_index)
1895 12963 : occ(bdtot_index)=paw_dmft%occnd(1,iband,iband,ikpt,isppol)
1896 12963 : ebanddmft=ebanddmft+dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index)
1897 12963 : ekindmft=ekindmft+dtset%wtk(ikpt)*occ(bdtot_index)*eknk(bdtot_index)
1898 : end if
1899 : end if
1900 :
1901 : energies%e_eigenvalues = energies%e_eigenvalues + &
1902 826735 : dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index)
1903 : energies%e_kinetic = energies%e_kinetic + &
1904 826735 : dtset%wtk(ikpt)*occ(bdtot_index)*eknk(bdtot_index)
1905 : energies%e_nlpsp_vfock = energies%e_nlpsp_vfock + &
1906 826735 : dtset%wtk(ikpt)*occ(bdtot_index)*enlxnk(bdtot_index)
1907 :
1908 826735 : if (usefock) then
1909 30595 : energies%e_fock=energies%e_fock + half*focknk(bdtot_index)*occ(bdtot_index)*dtset%wtk(ikpt)
1910 30595 : if (optforces>0) fock%fock_common%forces(:,:)=fock%fock_common%forces(:,:)+&
1911 1161 : dtset%wtk(ikpt)*occ(bdtot_index)*fockfornk(:,:,bdtot_index)
1912 : end if
1913 1699981 : if (optforces>0) grnl(:)=grnl(:)+dtset%wtk(ikpt)*occ(bdtot_index)*grnlnk(:,bdtot_index)
1914 : end if
1915 1373400 : bdtot_index=bdtot_index+1
1916 1486288 : if((paw_dmft%use_dmft==1.or.paw_dmft%use_dmft==10).and.dtset%nbandkss==0) then
1917 737052 : do iband1=1,nband_k
1918 708652 : if((paw_dmft%band_in(iband).and.paw_dmft%band_in(iband1)).or. &
1919 28400 : & ((paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7).and.iband==iband1)) then
1920 : ! write(std_out,*) "II+", isppol,ikpt,iband,iband1
1921 : ekindmft2=ekindmft2 + dtset%wtk(ikpt)*paw_dmft%occnd(1,iband,iband1,ikpt,isppol)*&
1922 219684 : eknk_nd(1,iband,iband1,ikpt,isppol)
1923 : ekindmft2=ekindmft2 - dtset%wtk(ikpt)*paw_dmft%occnd(2,iband,iband1,ikpt,isppol)*&
1924 219684 : eknk_nd(2,iband,iband1,ikpt,isppol)
1925 : ! write(std_out,*) "II", occnd(1,iband,iband1,ikpt,isppol),eknk_nd(isppol,ikpt,iband,iband1)
1926 : end if
1927 : end do
1928 : end if
1929 : end do
1930 : end do
1931 : end do
1932 :
1933 : ! Compute extended plane waves contributions
1934 13358 : if(associated(extfpmd)) then
1935 90 : extfpmd%nelect_res=-extfpmd%nelect
1936 90 : extfpmd%nelect=zero
1937 : call extfpmd%compute_nelect(energies%e_fermie,dtset%nband,extfpmd%nelect,dtset%nkpt,&
1938 90 : dtset%nspinor,dtset%nsppol,dtset%wtk)
1939 90 : extfpmd%nelect_res=extfpmd%nelect_res+extfpmd%nelect
1940 : call extfpmd%compute_e_kinetic(energies%e_fermie,dtset%nkpt,dtset%nspinor,&
1941 90 : dtset%nsppol,dtset%nband,dtset%wtk)
1942 : call extfpmd%compute_entropy(energies%entropy_extfpmd,energies%e_fermie,dtset%nkpt,&
1943 90 : dtset%nsppol,dtset%nspinor,dtset%wtk,dtset%nband)
1944 : end if
1945 :
1946 13358 : if(paw_dmft%use_dmft==1) then
1947 113 : energies%e_kinetic = energies%e_kinetic -ekindmft+ekindmft2
1948 113 : if(abs(dtset%pawprtvol)>=2) then
1949 : write(msg,'(4a,7(2x,a,2x,e14.7,a),a)') &
1950 69 : "-----------------------------------------------",ch10,&
1951 69 : "--- Energy for DMFT and tests (in Ha) ",ch10,&
1952 69 : "--- Ebandldatot (Ha.) = ",ebandldatot,ch10,&
1953 69 : "--- Ebandlda (Ha.) = ",ebandlda,ch10,&
1954 69 : "--- Ebanddmft (Ha.) = ",ebanddmft,ch10,&
1955 69 : "--- Ekinlda (Ha.) = ",ekinlda,ch10, &
1956 69 : "--- Ekindmftdiag (Ha.) = ",ekindmft,ch10,&
1957 69 : "--- Ekindmftnondiag(Ha.) = ",ekindmft2,ch10,&
1958 69 : "--- Edmft= (Ha.) = ",edmft,ch10,&
1959 138 : "-----------------------------------------------"
1960 69 : call wrtout(std_out,msg)
1961 : end if
1962 : ! if(paw_dmft%use_dmft==1.and.mpi_enreg%paral_kgb==1) paw_dmft%use_dmft=0
1963 : end if
1964 :
1965 : ABI_NVTX_START_RANGE(NVTX_MKRHO)
1966 :
1967 13358 : if (psps%usepaw==0) then
1968 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
1969 7408 : rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,extfpmd=extfpmd,nfft_blocks=gs_hamk%nfft_blocks)
1970 : else
1971 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
1972 5950 : rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,nfft_blocks=gs_hamk%nfft_blocks)
1973 : end if
1974 :
1975 : ABI_NVTX_END_RANGE()
1976 13358 : call timab(992,2,tsec)
1977 :
1978 : else
1979 : ! Treat fixed occupation numbers or non-self-consistent case
1980 :
1981 33635 : if (mpi_enreg%nproc_spkpt>1) then
1982 :
1983 6814 : call timab(989,1,tsec)
1984 :
1985 6814 : nbuf=2*mbdkpsp+dtset%nfft*dtset%nspden+4+3*natom*optforces
1986 : ! If Hartree-Fock calculation, the exact exchange energy is k-dependent.
1987 6814 : if(dtset%usefock==1) then
1988 0 : nbuf=nbuf+1
1989 0 : if (optforces>0) nbuf=nbuf+3*natom
1990 : end if
1991 6814 : if(iscf==-1 .or. iscf==-2)nbuf=2*mbdkpsp
1992 20442 : ABI_MALLOC(buffer1,(nbuf))
1993 : ! Pack eigen,resid,rho[wf]r,grnl,enl,ek
1994 336718 : buffer1(1:mbdkpsp)=eigen(:)
1995 336718 : buffer1(1+mbdkpsp:2*mbdkpsp)=resid(:)
1996 6814 : index1=2*mbdkpsp
1997 6814 : if(iscf/=-1 .and. iscf/=-2)then
1998 6728 : if (psps%usepaw==0) then
1999 11744 : buffer1(index1+1:index1+dtset%nfft*dtset%nspden)=reshape(rhor,(/dtset%nfft*dtset%nspden/))
2000 : else
2001 1712 : buffer1(index1+1:index1+dtset%nfft*dtset%nspden)=reshape(rhowfr,(/dtset%nfft*dtset%nspden/))
2002 : end if
2003 6728 : index1=index1+dtset%nfft*dtset%nspden
2004 6728 : buffer1(index1+1) = energies%e_kinetic
2005 6728 : buffer1(index1+2) = energies%e_eigenvalues
2006 6728 : buffer1(index1+3) = energies%e_nlpsp_vfock
2007 6728 : buffer1(index1+4) = energies%e_nucdip
2008 6728 : index1=index1+4
2009 : ! If Hartree-Fock calculation, save e_fock in buffer1
2010 6728 : if (dtset%usefock==1) then
2011 0 : buffer1(index1+1) = energies%e_fock
2012 0 : index1=index1+1
2013 0 : if (optforces>0)then
2014 0 : buffer1(index1+1:index1+3*natom)=reshape(fock%fock_common%forces,(/3*natom/))
2015 0 : index1=index1+3*natom
2016 : end if
2017 : end if
2018 38012 : if (optforces>0) buffer1(index1+1:index1+3*natom)=grnl(1:3*natom)
2019 : end if
2020 :
2021 : ! Build sum of everything
2022 6814 : call timab(48,1,tsec)
2023 6814 : call xmpi_sum(buffer1,nbuf,mpi_enreg%comm_kpt ,ierr)
2024 6814 : call timab(48,2,tsec)
2025 :
2026 : ! Unpack the final result
2027 336718 : eigen(:)=buffer1(1:mbdkpsp)
2028 336718 : resid(:)=buffer1(1+mbdkpsp:2*mbdkpsp)
2029 6814 : index1=2*mbdkpsp
2030 6814 : if(iscf/=-1 .and. iscf/=-2)then
2031 6728 : if (psps%usepaw==0) then
2032 5872 : ii=1
2033 11816 : do ispden=1,dtset%nspden
2034 33722734 : do ifft=1,dtset%nfft
2035 33710918 : rhor(ifft,ispden)=buffer1(index1+ii)
2036 33716862 : ii=ii+1
2037 : end do
2038 : end do
2039 : else
2040 856 : ii=1
2041 1952 : do ispden=1,dtset%nspden
2042 5282016 : do ifft=1,dtset%nfft
2043 5280064 : rhowfr(ifft,ispden)=buffer1(index1+ii)
2044 5281160 : ii=ii+1
2045 : end do
2046 : end do
2047 : end if
2048 6728 : index1=index1+dtset%nfft*dtset%nspden
2049 6728 : energies%e_kinetic = buffer1(index1+1)
2050 6728 : energies%e_eigenvalues = buffer1(index1+2)
2051 6728 : energies%e_nlpsp_vfock = buffer1(index1+3)
2052 6728 : energies%e_nucdip = buffer1(index1+4)
2053 6728 : index1=index1+4
2054 : ! If Hartree-Fock calculation, save e_fock in buffer1
2055 6728 : if (dtset%usefock==1) then
2056 0 : energies%e_fock = buffer1(index1+1)
2057 0 : index1=index1+1
2058 0 : if (optforces>0) then
2059 0 : fock%fock_common%forces(:,:)=reshape(buffer1(index1+1:index1+3*natom),(/3,natom/))
2060 0 : index1=index1+3*natom
2061 : end if
2062 : end if
2063 38012 : if (optforces>0) grnl(1:3*natom)=buffer1(index1+1:index1+3*natom)
2064 : end if
2065 6814 : ABI_FREE(buffer1)
2066 13628 : call timab(989,2,tsec)
2067 :
2068 : end if ! nproc_spkpt>1
2069 :
2070 : ! Compute extended plane waves contributions
2071 33635 : if(associated(extfpmd)) then
2072 : call extfpmd%compute_eshift(eigen,eknk,dtset%mband,dtset%nband,&
2073 0 : nfftf,dtset%nkpt,dtset%nsppol,dtset%nspden,dtset%wtk,vtrial)
2074 0 : extfpmd%nelect=zero
2075 : call extfpmd%compute_nelect(energies%e_fermie,dtset%nband,extfpmd%nelect,dtset%nkpt,&
2076 0 : dtset%nspinor,dtset%nsppol,dtset%wtk)
2077 : call extfpmd%compute_e_kinetic(energies%e_fermie,dtset%nkpt,dtset%nspinor,&
2078 0 : dtset%nsppol,dtset%nband,dtset%wtk)
2079 : call extfpmd%compute_entropy(energies%entropy_extfpmd,energies%e_fermie,dtset%nkpt,&
2080 0 : dtset%nsppol,dtset%nspinor,dtset%wtk,dtset%nband)
2081 : end if
2082 :
2083 : ! Compute the highest occupied eigenenergy
2084 33635 : if(iscf/=-1 .and. iscf/=-2)then
2085 33159 : call timab(993,1,tsec)
2086 33159 : energies%e_fermie = -huge(one)
2087 33159 : if (dtset%occopt==9) then
2088 0 : energies%e_fermih = -huge(one)
2089 : end if
2090 33159 : bdtot_index=1
2091 69225 : do isppol=1,dtset%nsppol
2092 381159 : do ikpt=1,dtset%nkpt
2093 311934 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
2094 348000 : if (dtset%occopt/=9) then
2095 1813772 : do iband=1,nband_k
2096 1501838 : if(abs(occ(bdtot_index))>tol8 .and. eigen(bdtot_index)>energies%e_fermie+tol10) then
2097 178871 : energies%e_fermie=eigen(bdtot_index)
2098 : end if
2099 1813772 : bdtot_index=bdtot_index+1
2100 : end do
2101 : else
2102 : ! In case occopt 9, computing the fermi level for the electrons remaining in the VB = fermi level of holes
2103 0 : do iband=1,dtset%ivalence
2104 0 : if(abs(occ(bdtot_index))>tol8 .and. eigen(bdtot_index)>energies%e_fermih+tol10) then
2105 0 : energies%e_fermih=eigen(bdtot_index)
2106 : end if
2107 0 : bdtot_index=bdtot_index+1
2108 : end do
2109 : ! In case occopt 9, computing the fermi level for the electrons thermalized in the conduction bands
2110 0 : do iband=dtset%ivalence+1,nband_k
2111 0 : if(abs(occ(bdtot_index))>tol8 .and. eigen(bdtot_index)>energies%e_fermie+tol10) then
2112 0 : energies%e_fermie=eigen(bdtot_index)
2113 : end if
2114 0 : bdtot_index=bdtot_index+1
2115 : end do
2116 : end if
2117 : end do
2118 : end do
2119 33159 : call xmpi_max(energies%e_fermie,spaceComm_distrb,ierr)
2120 33159 : if (dtset%occopt == 9) then
2121 0 : call xmpi_max(energies%e_fermih,spaceComm_distrb,ierr)
2122 : end if
2123 33159 : call timab(993,2,tsec)
2124 : end if
2125 :
2126 : ! If needed, compute rhog, and symmetrizes the density
2127 33635 : if (iscf > 0 .or. iscf==-3 ) then
2128 : ! energies%e_fermie=zero ! Actually, should determine the maximum of the valence band XG20020802
2129 33159 : nfftot=dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3)
2130 :
2131 33159 : call timab(994,1,tsec)
2132 33159 : if (psps%usepaw==0) then
2133 : call symrhg(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,&
2134 24003 : & dtset%nsppol,dtset%nsym,phnons,rhog,rhor,rprimd,dtset%symafm,dtset%symrel,dtset%tnons)
2135 : else
2136 : call symrhg(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,&
2137 9156 : & dtset%nsppol,dtset%nsym,phnons,rhowfg,rhowfr,rprimd,dtset%symafm,dtset%symrel,dtset%tnons)
2138 : end if
2139 33159 : call timab(994,2,tsec)
2140 : ! We now have both rho(r) and rho(G), symmetrized, and if dtset%nsppol=2
2141 : ! we also have the spin-up density, symmetrized, in rhor(:,2).
2142 : end if
2143 :
2144 : end if ! End of test on varying or fixed occupation numbers
2145 :
2146 46993 : call timab(994,1,tsec)
2147 :
2148 : ! Compute the kinetic energy density
2149 46993 : if(dtset%usekden==1 .and. (iscf > 0 .or. iscf==-3 ) )then
2150 805 : if (psps%usepaw==0) then
2151 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
2152 323 : taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1,nfft_blocks=gs_hamk%nfft_blocks)
2153 : else
2154 : call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
2155 482 : tauwfg,tauwfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1,nfft_blocks=gs_hamk%nfft_blocks)
2156 : end if
2157 : end if
2158 :
2159 46993 : ABI_FREE(eknk)
2160 46993 : if (usefock) then
2161 1375 : ABI_FREE(focknk)
2162 1375 : if (optforces>0)then
2163 220 : ABI_FREE(fockfornk)
2164 : end if
2165 : end if
2166 46993 : ABI_FREE(eknk_nd)
2167 46993 : ABI_FREE(grnlnk)
2168 46993 : ABI_FREE(enlxnk)
2169 :
2170 : ! In the non-self-consistent case, print eigenvalues and residuals
2171 46993 : if(iscf<=0 .and. me_distrb==0)then
2172 434 : option=2 ; enunit=1 ; vxcavg_dum=zero
2173 : call prteigrs(eigen,enunit,energies%e_fermie,energies%e_fermih,&
2174 : & dtfil%fnameabo_app_eig,ab_out,iscf,dtset%kptns,dtset%kptopt,dtset%mband,&
2175 : & dtset%nband,nbdbuf_eff,dtset%nkpt,nnsclo_now,dtset%nsppol,occ,dtset%occopt,option,&
2176 434 : & dtset%prteig,prtvol,resid,dtset%tolwfr,vxcavg_dum,dtset%wtk)
2177 : end if
2178 :
2179 : ! Find largest residual over bands, k points, and spins, except for nbdbuf highest bands
2180 46993 : ibdkpt=1
2181 46993 : residm=zero
2182 98669 : do isppol=1,dtset%nsppol
2183 540437 : do ikpt=1,dtset%nkpt
2184 441768 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
2185 441768 : if (nbdbuf_eff>=0) then
2186 440046 : nband_eff=max(1,nband_k-nbdbuf_eff)
2187 3643432 : residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_eff-1)))
2188 1722 : else if (nbdbuf_eff==-101) then
2189 20148 : residm=max(residm,maxval(occ(ibdkpt:ibdkpt+nband_k-1)*resid(ibdkpt:ibdkpt+nband_k-1)))
2190 : else
2191 0 : ABI_ERROR('Bad value of nbdbuf_eff')
2192 : end if
2193 493444 : ibdkpt=ibdkpt+nband_k
2194 : end do
2195 : end do
2196 :
2197 : end if !usewvl==0
2198 :
2199 : !===================================================================
2200 : ! End of PLANE WAVES section
2201 : !===================================================================
2202 :
2203 : !In the self-consistent case, diagnose lack of unoccupied state (for each spin and k-point).
2204 :
2205 : !Print a warning if the number of such messages already written does not exceed mwarning.
2206 : ! MG: This is not a good idea as this is a typical mistake done by beginners and we should
2207 : ! keep on spamming this warning message in the log file.
2208 : !mwarning=5
2209 : !if(nwarning<mwarning .and. iscf>=0)then
2210 : !nwarning=nwarning+1
2211 :
2212 46993 : if(iscf>=0)then
2213 46497 : bdtot_index=1
2214 46497 : quit=0
2215 72434 : do isppol=1,dtset%nsppol
2216 209433 : do ikpt=1,dtset%nkpt
2217 183496 : min_occ=two
2218 183496 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
2219 2288180 : do iband=1,nband_k
2220 : ! if ndbbuf_eff<=0, compute min_occ as usual
2221 : ! if nbdbuf_eff>0, compute min_occ only for bands not in the buffer
2222 2104684 : if (nbdbuf_eff<=0.or.iband<=nband_k-nbdbuf_eff) then
2223 1988438 : if(occ(bdtot_index)<min_occ)min_occ=occ(bdtot_index)
2224 : end if
2225 2288180 : bdtot_index=bdtot_index+1
2226 : end do
2227 209433 : if(min_occ>0.01_dp .and. .not. associated(extfpmd))then
2228 23598 : if(dtset%nsppol==1)then
2229 : write(msg, '(a,i0,3a,f7.3,5a)' )&
2230 21933 : 'For k-point number: ',ikpt,',',ch10,&
2231 21933 : 'The minimal occupation factor is: ',min_occ,'.',ch10,&
2232 21933 : 'An adequate monitoring of convergence requires it to be at most 0.01_dp.',ch10,&
2233 43866 : 'Action: increase slightly the number of bands (or decrease nbdbuf).'
2234 : else
2235 : write(msg, '(a,i0,3a,i0,a,f7.3,5a)' )&
2236 1665 : 'For k-point number: ',ikpt,', and',ch10,&
2237 1665 : 'for spin polarization: ',isppol, ' the minimal occupation factor is: ',min_occ,'.',ch10,&
2238 1665 : 'An adequate monitoring of convergence requires it to be at most 0.01_dp.',ch10,&
2239 3330 : 'Action: increase slightly the number of bands (or decrease nbdbuf).'
2240 : end if
2241 23598 : ABI_WARNING(msg)
2242 : quit=1
2243 : exit ! It is enough if one lack of adequate occupation is identified, so exit.
2244 : end if
2245 : end do
2246 22899 : if (quit==1) exit
2247 : end do
2248 : end if
2249 :
2250 : ABI_NVTX_START_RANGE(NVTX_VTORHO_EXTRA)
2251 46993 : if (iscf>0.or.iscf==-3 .or. (dtset%usewvl==1 .and. iscf==0)) then
2252 :
2253 : ! PAW: Build new rhoij quantities from new occ then symetrize them
2254 : ! Compute and add the compensation density to rhowfr to get the total density
2255 46517 : if (psps%usepaw==1) then
2256 15106 : call timab(555,1,tsec)
2257 15106 : if (paral_atom) then
2258 27202 : ABI_MALLOC(pawrhoij_unsym,(natom))
2259 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij,&
2260 5042 : & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cpxocc=dtset%pawcpxocc)
2261 : call pawrhoij_alloc(pawrhoij_unsym,cplex_rhoij,nspden_rhoij,dtset%nspinor,&
2262 5042 : & dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0)
2263 : else
2264 10064 : pawrhoij_unsym => pawrhoij
2265 : end if
2266 15106 : if (usecprj_local==1) then
2267 : call pawmkrhoij(atindx,atindx1,cprj,gs_hamk%dimcprj,dtset%istwfk,dtset%kptopt,dtset%mband,mband_cprj,&
2268 : & mcprj_local,dtset%mkmem,mpi_enreg,natom,dtset%nband,dtset%nkpt,dtset%nspden,dtset%nspinor,&
2269 15106 : & dtset%nsppol,occ,dtset%paral_kgb,paw_dmft,pawrhoij_unsym,dtfil%unpaw,dtset%usewvl,dtset%wtk)
2270 : else
2271 0 : mcprj_tmp=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol
2272 0 : ABI_MALLOC(cprj_tmp,(natom,mcprj_tmp))
2273 0 : call pawcprj_alloc(cprj_tmp,0,gs_hamk%dimcprj)
2274 : call ctocprj(atindx,cg,1,cprj_tmp,gmet,gprimd,0,0,0,dtset%istwfk,kg,dtset%kptns,&
2275 : & mcg,mcprj_tmp,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,&
2276 : & dtset%natom,nattyp,dtset%nband,dtset%natom,dtset%ngfft,dtset%nkpt,dtset%nloalg,&
2277 : & npwarr,dtset%nspinor,dtset%nsppol,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,&
2278 0 : & ucvol,dtfil%unpaw,xred,ylm,ylmgr_dum)
2279 : call pawmkrhoij(atindx,atindx1,cprj_tmp,gs_hamk%dimcprj,dtset%istwfk,dtset%kptopt,&
2280 : & dtset%mband,mband_cprj,mcprj_tmp,dtset%mkmem,mpi_enreg,natom,dtset%nband,dtset%nkpt,&
2281 : & dtset%nspden,dtset%nspinor,dtset%nsppol,occ,dtset%paral_kgb,paw_dmft,pawrhoij_unsym,&
2282 0 : & dtfil%unpaw,dtset%usewvl,dtset%wtk)
2283 0 : call pawcprj_free(cprj_tmp)
2284 0 : ABI_FREE(cprj_tmp)
2285 : end if
2286 15106 : call timab(555,2,tsec)
2287 : ! Build symetrized packed rhoij and compensated pseudo density
2288 15106 : cplex=1;ipert=0;idir=0;qpt(:)=zero
2289 15106 : if(dtset%usewvl==0) then
2290 : call pawmkrho(1,compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,&
2291 : & my_natom,natom,dtset%nspden,dtset%nsym,ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,&
2292 : & dtset%pawprtvol,pawrhoij,pawrhoij_unsym,pawtab,qpt,rhowfg,rhowfr,rhor,rprimd,dtset%symafm,&
2293 15106 : & symrec,dtset%typat,ucvol,dtset%usewvl,xred,pawnhat=nhat,rhog=rhog,extfpmd=extfpmd)
2294 15106 : if (dtset%usekden==1) then
2295 : ! DO WE NEED TAUG?
2296 482 : call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,tauwfg,taug,tauwfr,taur)
2297 : end if
2298 : else
2299 : ! here do not pass rhog, we do not use it
2300 : call pawmkrho(1,compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,&
2301 : & my_natom,natom,dtset%nspden,dtset%nsym,ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,&
2302 : & dtset%pawprtvol,pawrhoij,pawrhoij_unsym,pawtab,qpt,rhowfg,rhowfr,rhor,rprimd,dtset%symafm,&
2303 0 : & symrec,dtset%typat,ucvol,dtset%usewvl,xred,pawnhat=nhat)
2304 : ! In WVL: copy density to BigDFT object:
2305 0 : call wvl_rho_abi2big(1,rhor,wvl%den)
2306 : end if
2307 15106 : if (paral_atom) then
2308 5042 : call pawrhoij_free(pawrhoij_unsym)
2309 5042 : ABI_FREE(pawrhoij_unsym)
2310 : end if
2311 : end if ! psps%usepaw==1
2312 :
2313 46517 : if(paw_dmft%use_dmft==1) then
2314 : ! == check noccmmp
2315 : call setnoccmmp(1,0,rdum4,0,0,idum3,my_natom,natom,0,1,dtset%nsppol,0,ntypat,&
2316 : & paw_ij,pawang,dtset%pawprtvol,pawrhoij,pawtab,rdum2,idum1,dtset%typat,0,dtset%usepawu,&
2317 113 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
2318 : end if
2319 :
2320 : ! Find and print minimum and maximum total electron density and locations
2321 : ! Compute density residual (if required) and its squared norm
2322 46517 : if (iscf>=0) then
2323 46497 : if (psps%usepaw==0) then
2324 31392 : call prtrhomxmn(std_out,mpi_enreg,dtset%nfft,dtset%ngfft,dtset%nspden,1,rhor,ucvol=ucvol)
2325 : else
2326 15105 : call prtrhomxmn(std_out,mpi_enreg,nfftf,pawfgr%ngfft,dtset%nspden,1,rhor,ucvol=ucvol)
2327 : end if
2328 46497 : if (optres==1) then
2329 750913085 : nvresid=rhor-nvresid
2330 : ! /!\ ---- DO NOT CHANGE THESE LINES WITHOUT CORE DEVELOPERS PERMISSION ---- /!\
2331 : ! LB-03/2026:
2332 : ! A noise can accumulate in nvresid after each SCF cycle,
2333 : ! resulting in different densities/potentials for different MPI processes.
2334 : ! This has been observed using threads, but could happen in other contexts.
2335 : ! This slowly worsens the SCF cycle, leading to wrong results after many iterations.
2336 : ! So here we compute the mean of nvresid over all MPI processes to reduce the noise.
2337 : ! This error is difficult to test as it is observed in long runs only, so BE VERY CAREFUL.
2338 14329 : mpicomm = mpi_enreg%comm_kptband
2339 14329 : nmpi = xmpi_comm_size(mpicomm)
2340 14329 : nthreads = xomp_get_num_threads(open_parallel=.true.)
2341 14329 : if (nmpi>1.and.nthreads>1) then
2342 0 : ABI_MALLOC(nvresid_tmp,(nfftf,dtset%nspden))
2343 0 : nvresid_tmp(:,:) = nvresid(:,:) / nmpi
2344 0 : call xmpi_sum(nvresid_tmp,mpicomm,ierr)
2345 0 : if (ierr/=0) then
2346 0 : ABI_ERROR("Error in mpi sum (nvresid)")
2347 : end if
2348 0 : nvresid(:,:) = nvresid_tmp(:,:)
2349 0 : ABI_FREE(nvresid_tmp)
2350 : end if
2351 : ! /!\--------------------/!\
2352 : !
2353 14329 : call sqnorm_v(1,nfftf,nres2,dtset%nspden,optres,nvresid,mpi_comm_sphgrid=mpi_comm_sphgrid)
2354 14329 : if (dtset%usekden==1) then
2355 : if (optres==1) then
2356 28452708 : tauresid=taur-tauresid
2357 : ! /!\ ---- DO NOT CHANGE THESE LINES WITHOUT CORE DEVELOPERS PERMISSION ---- /!\
2358 : if (nmpi>1.and.nthreads>1) then
2359 1992 : ABI_MALLOC(nvresid_tmp,(nfftf,dtset%nspden))
2360 28452708 : nvresid_tmp(:,:) = tauresid(:,:) / nmpi
2361 498 : call xmpi_sum(nvresid_tmp,mpicomm,ierr)
2362 498 : if (ierr/=0) then
2363 0 : ABI_ERROR("Error in mpi sum (tauresid)")
2364 : end if
2365 28452708 : tauresid(:,:) = nvresid_tmp(:,:)
2366 996 : ABI_FREE(nvresid_tmp)
2367 : end if
2368 : ! /!\--------------------/!\
2369 : end if
2370 : end if
2371 : end if
2372 : end if
2373 :
2374 : end if ! iscf>0 or iscf=-3
2375 : ABI_NVTX_END_RANGE()
2376 :
2377 46993 : if(psps%usepaw==1.and.(iscf>=0.or.iscf==-3)) then
2378 15106 : ABI_FREE(rhowfr)
2379 15106 : ABI_FREE(rhowfg)
2380 15106 : if (dtset%usekden==1) then
2381 482 : ABI_FREE(tauwfr)
2382 482 : ABI_FREE(tauwfg)
2383 : end if
2384 : end if
2385 :
2386 46993 : call timab(994,2,tsec)
2387 :
2388 46993 : if (iscf==-1) then
2389 : ! Eventually compute the excited states within tddft
2390 22 : call timab(995,1,tsec)
2391 22 : if (psps%usepaw==1) then
2392 : ! In case of PAW calculation, have to transfer kxc from the fine to the coarse grid:
2393 0 : ABI_MALLOC(cgrkxc,(dtset%nfft,nkxc))
2394 0 : do ikxc=1,nkxc
2395 0 : call transgrid(1,mpi_enreg,1,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrkxc(:,ikxc),kxc(:,ikxc))
2396 : end do
2397 : call tddft(cg,dtfil,dtset,eigen,etotal,gmet,gprimd,gsqcut,&
2398 : & kg,cgrkxc,dtset%mband,mgfftdiel,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%nfft,&
2399 0 : & ngfftdiel,dtset%nkpt,nkxc,npwarr,dtset%nspinor,dtset%nsppol,occ,ucvol,wffnew)
2400 0 : ABI_FREE(cgrkxc)
2401 : else
2402 : call tddft(cg,dtfil,dtset,eigen,etotal,gmet,gprimd,gsqcut,&
2403 : & kg,kxc,dtset%mband,mgfftdiel,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%nfft,&
2404 22 : & ngfftdiel,dtset%nkpt,nkxc,npwarr,dtset%nspinor,dtset%nsppol,occ,ucvol,wffnew)
2405 : end if
2406 22 : call timab(995,2,tsec)
2407 :
2408 : else
2409 : ! Eventually compute the susceptibility matrix and the
2410 : ! dielectric matrix when istep_mix is equal to 1 or dielstrt
2411 46971 : call timab(996,1,tsec)
2412 46971 : computesusmat = dtset%testsusmat(dielop, dielstrt, istep_mix) !test if the matrix is to be computed
2413 46971 : if(computesusmat) then
2414 82 : dielar(1)=dtset%diecut;dielar(2)=dtset%dielng
2415 82 : dielar(3)=dtset%diemac;dielar(4)=dtset%diemix
2416 82 : dielar(5)=dtset%diegap;dielar(6)=dtset%dielam
2417 82 : dielar(7)=dtset%diemix;if (iscf>=10) dielar(7)=dtset%diemixmag
2418 82 : usetimerev=1
2419 82 : if (psps%usepaw==1.and.dtset%pawspnorb>0.and.dtset%kptopt/=1.and.dtset%kptopt/=2) usetimerev=0
2420 82 : neglect_pawhat=1-dtset%pawsushat
2421 : call suscep_stat(atindx,atindx1,cg,cprj,dielar,&
2422 : & gs_hamk%dimcprj,doccde,eigen,gbound_diel,gprimd,&
2423 : & irrzondiel,dtset%istwfk,kg,kg_diel,lmax_diel,&
2424 : & dtset%mband,mcg,mcprj_local,mgfftdiel,dtset%mkmem,mpi_enreg,dtset%mpw,natom,dtset%nband,&
2425 : & neglect_pawhat,nfftdiel,ngfftdiel,&
2426 : & dtset%nkpt,npwarr,npwdiel,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%nsym,ntypat,&
2427 : & occ,dtset%occopt,pawang,pawtab,phnonsdiel,ph1ddiel,rprimd,&
2428 : & susmat,dtset%symafm,dtset%symrel,dtset%tnons,dtset%typat,ucvol,&
2429 82 : & dtfil%unpaw,usecprj_local,psps%usepaw,usetimerev,dtset%wtk,ylmdiel)
2430 : end if
2431 46971 : call timab(996,2,tsec)
2432 :
2433 : end if ! end condition on iscf
2434 :
2435 46993 : call gs_hamk%free()
2436 :
2437 46993 : if (psps%usepaw==1) then
2438 15126 : if (usecprj==0) then
2439 14486 : call pawcprj_free(cprj_local)
2440 3347716 : ABI_FREE(cprj_local)
2441 : end if
2442 : end if
2443 :
2444 : #ifdef HAVE_OPENMP_OFFLOAD
2445 : !$OMP TARGET EXIT DATA MAP(from:cg) IF(transfer_cg_once)
2446 : #endif
2447 :
2448 46993 : if(dtset%usewvl==0) then
2449 46993 : ABI_FREE(EigMin)
2450 46993 : ABI_FREE(doccde)
2451 : #if defined HAVE_GPU_CUDA
2452 : if(dtset%gpu_option==ABI_GPU_LEGACY .or. dtset%gpu_option==ABI_GPU_KOKKOS) call gpu_finalize_ham_data()
2453 : #endif
2454 : end if
2455 :
2456 93986 : call timab(980,2,tsec)
2457 :
2458 : DBG_EXIT("COLL")
2459 :
2460 : contains
2461 : !!***
2462 :
2463 : !!****f* ABINIT/wvl_nscf_loop
2464 : !! NAME
2465 : !! wvl_nscf_loop
2466 : !!
2467 : !! FUNCTION
2468 : !! Non-self-consistent field cycle in Wavelets
2469 : !! See also "wvl_nscf_loop_bigdft"
2470 : !!
2471 : !! INPUTS
2472 : !! nnsclo= number of non-self consistent field iterations
2473 : !!
2474 : !! OUTPUT
2475 : !!
2476 : !! SOURCE
2477 :
2478 : subroutine wvl_nscf_loop()
2479 :
2480 : !Arguments ------------------------------------
2481 : ! integer, intent(in) :: istep,mcprj,nfft,nnsclo
2482 : ! real(dp), intent(inout) :: residm
2483 : ! type(dataset_type), intent(in) :: dtset
2484 : ! type(MPI_type), intent(in) :: mpi_enreg
2485 : ! type(energies_type), intent(inout) :: energies
2486 : ! type(wvl_data), intent(inout) :: wvl
2487 : ! !arrays
2488 : ! real(dp), intent(inout) :: xcart(3, dtset%natom)
2489 : ! real(dp), dimension(6), intent(out) :: strsxc
2490 : ! type(pawcprj_type),dimension(dtset%natom,mcprj),intent(out)::cprj
2491 :
2492 : !Local variables-------------------------------
2493 : integer :: inonsc,ii
2494 : integer,parameter :: iscf_=-1 !do not do a SCF cycle
2495 : logical,parameter :: do_scf=.false. !do not do a SCF cycle
2496 : logical,parameter :: wvlbigdft=.false.
2497 : real(dp) :: dum,eexctx,eh,ekin,eloc,enl,esicdc,evxc,exc
2498 : ! *************************************************************************
2499 :
2500 : DBG_ENTER("COLL")
2501 :
2502 : if(nnsclo_now>0) then
2503 : do inonsc=1,nnsclo_now
2504 : call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,&
2505 : & istep,inonsc,iscf_,mpi_enreg%me_wvl,dtset%natom,&
2506 : & nfftf,mpi_enreg%nproc_wvl,dtset%nspden,&
2507 : & dum,do_scf,evxc,wvl,wvlbigdft,xcart,strsxc)
2508 : call wvl_hpsitopsi(cprj,dtset,energies,inonsc,mcprj_local,mpi_enreg, &
2509 : & residm,wvl,xcart)
2510 : if(residm<dtset%tolwfr) exit !Exit loop if converged
2511 : end do
2512 :
2513 : else
2514 : do ii=1, dtset%nline
2515 : ! Direct minimization technique: no diagonalization
2516 : call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,&
2517 : & istep,ii,iscf_,mpi_enreg%me_wvl,dtset%natom,&
2518 : & nfftf,mpi_enreg%nproc_wvl,dtset%nspden,&
2519 : & dum,do_scf,evxc,wvl,wvlbigdft,xcart,strsxc)
2520 : call wvl_hpsitopsi(cprj,dtset,energies,ii,mcprj_local,mpi_enreg, &
2521 : & residm,wvl,xcart)
2522 : if(residm<dtset%tolwfr) exit !Exit loop if converged
2523 : end do
2524 : end if
2525 :
2526 : ! Update energies depending on new WF
2527 : energies%e_kinetic=ekin
2528 : energies%e_nlpsp_vfock=enl
2529 : energies%e_exactX=eexctx
2530 : energies%e_sicdc=esicdc
2531 :
2532 : ! Eventually update energies depending on density
2533 : if (dtset%iscf<10) then
2534 : energies%e_localpsp=eloc
2535 : energies%e_hartree=eh
2536 : energies%e_xc=exc ; energies%e_xcdc=evxc
2537 : else if (nnsclo_now==0) then
2538 : energies%e_localpsp=eloc
2539 : end if
2540 :
2541 : ! End of nscf iterations
2542 : if (do_last_ortho) then
2543 : ! !Don't update energies (nscf cycle has been done); just recompute potential
2544 : inonsc=nnsclo_now;if (nnsclo_now==0) inonsc=dtset%nline
2545 : call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,&
2546 : & istep,inonsc,iscf_,mpi_enreg%me_wvl,dtset%natom,&
2547 : & nfftf,mpi_enreg%nproc_wvl,dtset%nspden,&
2548 : & dum,do_scf,evxc,wvl,wvlbigdft,xcart,strsxc)
2549 : end if
2550 :
2551 : DBG_EXIT("COLL")
2552 :
2553 : end subroutine wvl_nscf_loop
2554 : !!***
2555 :
2556 : !!****f* ABINIT/wvl_nscf_loop_bigdft
2557 : !! NAME
2558 : !! wvl_nscf_loop_bigdft
2559 : !!
2560 : !! FUNCTION
2561 : !! Non-self-consistent field cycle in Wavelets
2562 : !! It follows the BigDFT scheme.
2563 : !! See also "wvl_nscf_loop"
2564 : !!
2565 : !! INPUTS
2566 : !! nnsclo= number of non-self consistent field iterations
2567 : !!
2568 : !! OUTPUT
2569 : !! argout(sizeout)=description
2570 : !!
2571 : !! SOURCE
2572 :
2573 : subroutine wvl_nscf_loop_bigdft()
2574 :
2575 : !Arguments ------------------------------------
2576 : ! integer, intent(in) :: istep,mcprj,nfft,nnsclo
2577 : ! real(dp), intent(inout) :: residm
2578 : ! real(dp), intent(out) :: nres2
2579 : ! type(dataset_type), intent(in) :: dtset
2580 : ! type(MPI_type), intent(in) :: mpi_enreg
2581 : ! type(energies_type), intent(inout) :: energies
2582 : ! type(wvl_data), intent(inout) :: wvl
2583 : !arrays
2584 : ! real(dp), intent(inout) :: xcart(3, dtset%natom)
2585 : ! real(dp), dimension(6), intent(out) :: strsxc
2586 : ! type(pawcprj_type),dimension(dtset%natom,mcprj),intent(out)::cprj
2587 :
2588 : !Local variables-------------------------------
2589 : integer :: inonsc
2590 : integer,parameter :: iscf_=-1 !do not do a SCF cycle
2591 : logical,parameter :: do_scf=.false. !do not do a SCF cycle
2592 : logical,parameter :: wvlbigdft=.true.
2593 : real(dp) :: eexctx,eh,ekin,eloc,enl,esicdc,evxc,exc
2594 : ! *************************************************************************
2595 :
2596 : DBG_ENTER("COLL")
2597 :
2598 : call wvl_hpsitopsi(cprj,dtset, energies, istep, mcprj_local,mpi_enreg, &
2599 : & residm, wvl,xcart)
2600 :
2601 : if (nnsclo_now>2) then
2602 : do inonsc = 2, nnsclo_now-1
2603 : call wvl_psitohpsi(dtset%diemix, energies%e_exactX, energies%e_xc, &
2604 : & energies%e_hartree, energies%e_kinetic, energies%e_localpsp, &
2605 : & energies%e_nlpsp_vfock, energies%e_sicdc, istep, inonsc, iscf_, &
2606 : & mpi_enreg%me_wvl, dtset%natom, nfftf, mpi_enreg%nproc_wvl,&
2607 : & dtset%nspden, nres2, do_scf,energies%e_xcdc, &
2608 : & wvl, wvlbigdft, xcart, strsxc)
2609 : call wvl_hpsitopsi(cprj,dtset, energies, inonsc, mcprj_local,mpi_enreg, &
2610 : & residm, wvl,xcart)
2611 : end do
2612 : end if
2613 :
2614 : ! End of nscf iterations
2615 : if (do_last_ortho.and.nnsclo_now<=1) then
2616 : ! !Don't update energies (nscf cycle has been done); just recompute potential
2617 : call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc, &
2618 : & istep, 1, iscf_, mpi_enreg%me_wvl, dtset%natom, nfftf, &
2619 : & mpi_enreg%nproc_wvl,dtset%nspden, nres2, do_scf,evxc, &
2620 : & wvl, wvlbigdft, xcart, strsxc)
2621 : else if (do_last_ortho.and.nnsclo_now>1) then
2622 : ! !Update energies and potential (nscf cycles are not finished)
2623 : call wvl_psitohpsi(dtset%diemix, energies%e_exactX, energies%e_xc, &
2624 : & energies%e_hartree,energies%e_kinetic, energies%e_localpsp, &
2625 : & energies%e_nlpsp_vfock, energies%e_sicdc, istep, nnsclo_now, iscf_, &
2626 : & mpi_enreg%me_wvl, dtset%natom, nfftf, mpi_enreg%nproc_wvl,&
2627 : & dtset%nspden, nres2, do_scf,energies%e_xcdc, &
2628 : & wvl, wvlbigdft, xcart, strsxc)
2629 : end if
2630 :
2631 : DBG_EXIT("COLL")
2632 :
2633 : end subroutine wvl_nscf_loop_bigdft
2634 : !!***
2635 :
2636 : !!****f* ABINIT/e_eigen
2637 : !! NAME
2638 : !! e_eigen
2639 : !!
2640 : !! FUNCTION
2641 : !! Computes eigenvalues energy from eigen, occ, kpt, wtk
2642 : !!
2643 : !! INPUTS
2644 : !! eigen(nkpt*nsppol)=eigenvalues
2645 : !! mband= maximum number of bands
2646 : !! nband(nkpt*nsppol)= number of bands for each k-point and spin
2647 : !! nkpt= number of k-points
2648 : !! nsppol= number of spin polarization
2649 : !! occ(mband*nkpt*nsppol)=occupations
2650 : !! wtk(nkpt)= k-point weights
2651 : !!
2652 : !! OUTPUT
2653 : !! e_eigenvalues= eigenvalues energy
2654 : !!
2655 : !! SOURCE
2656 :
2657 : subroutine e_eigen(eigen,e_eigenvalues,mband,nband,nkpt,nsppol,occ,wtk)
2658 :
2659 : !Arguments ------------------------------------
2660 : integer , intent(in) :: mband,nkpt,nsppol
2661 : integer , intent(in) :: nband(nkpt*nsppol)
2662 : real(dp) , intent(in) :: eigen(mband*nkpt*nsppol)
2663 : real(dp) , intent(in) :: occ(mband*nkpt*nsppol)
2664 : real(dp) , intent(in) :: wtk(nkpt)
2665 : real(dp) , intent(out) :: e_eigenvalues
2666 :
2667 : !Local variables-------------------------------
2668 : integer :: ib,iband,ii,ikpt,isppol,nband_k
2669 : real(dp) :: wtk_k
2670 : ! *************************************************************************
2671 :
2672 : DBG_ENTER("COLL")
2673 : ii=0;ib=0
2674 : do isppol=1,nsppol
2675 : do ikpt=1,nkpt
2676 : ii=ii+1
2677 : nband_k=nband(ii) ; wtk_k=wtk(ii)
2678 : do iband=1,nband_k
2679 : ib=ib+1
2680 : if(abs(occ(ib)) > tol8) then
2681 : e_eigenvalues = e_eigenvalues + wtk_k*occ(ib)*eigen(ib)
2682 : end if
2683 : end do
2684 : end do
2685 : end do
2686 :
2687 : DBG_EXIT("COLL")
2688 :
2689 : end subroutine e_eigen
2690 : !!***
2691 :
2692 : !!****f* ABINIT/wvl_occ
2693 : !! NAME
2694 : !! wvl_occ
2695 : !!
2696 : !! FUNCTION
2697 : !! Computes occupations for the wavelet case
2698 : !!
2699 : !! NOTES
2700 : !! for the wvlbigdft case, see the routine 'wvl_occ_bigdft'
2701 : !!
2702 : !! SOURCE
2703 :
2704 : subroutine wvl_occ()
2705 :
2706 : !Local variables-------------------------------
2707 : real(dp):: doccde_(dtset%mband*dtset%nkpt*dtset%nsppol)
2708 : ! *************************************************************************
2709 :
2710 : DBG_ENTER("COLL")
2711 :
2712 : ! Compute the new occupation numbers from eigen
2713 : call newocc(doccde_,eigen,energies%entropy_ks,energies%e_fermie,energies%e_fermih,dtset%ivalence,dtset%spinmagntarget,&
2714 : & dtset%mband,dtset%nband,dtset%nelect,dtset%ne_qFD,dtset%nh_qFD,dtset%nkpt,dtset%nspinor,&
2715 : & dtset%nsppol,occ,dtset%occopt,prtvol,dtset%tphysel,dtset%tsmear,dtset%wtk,&
2716 : & prtstm=dtset%prtstm,stmbias=dtset%stmbias)
2717 :
2718 : ! Copy occupations and efermi to BigDFT variables
2719 : call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,1,wvl%wfs)
2720 :
2721 : #if defined HAVE_BIGDFT
2722 : ! Copy Fermi level to BigDFT variable:
2723 : wvl%wfs%ks%orbs%efermi=energies%e_fermie
2724 : #endif
2725 :
2726 : DBG_EXIT("COLL")
2727 :
2728 : end subroutine wvl_occ
2729 : !!***
2730 :
2731 : !!****f* ABINIT/wvl_occ_bigdft
2732 : !! NAME
2733 : !! wvl_occ_bigdft
2734 : !!
2735 : !! FUNCTION
2736 : !! Computes occupations for the wavelet case
2737 : !! Using BigDFT routines
2738 : !!
2739 : !! NOTES
2740 : !! for the wvlbigdft case, see the routine 'wvl_occ_bigdft'
2741 : !!
2742 : !! SOURCE
2743 :
2744 : subroutine wvl_occ_bigdft()
2745 :
2746 : ! *************************************************************************
2747 :
2748 : DBG_ENTER("COLL")
2749 :
2750 : ! Transfer occopt from ABINIT to BigDFT
2751 : #if defined HAVE_BIGDFT
2752 : occopt_bigdft=dtset%occopt
2753 : call wvl_occopt_abi2big(occopt_bigdft,occopt_bigdft,1)
2754 :
2755 : !Calculate occupations using BigDFT routine
2756 : call evaltoocc(mpi_enreg%me_wvl, mpi_enreg%nproc_wvl, .false., &
2757 : & dtset%tsmear, wvl%wfs%ks%orbs, occopt_bigdft)
2758 :
2759 : !Pass e_fermi from BigDFT object to ABINIT variable:
2760 : energies%e_fermie = wvl%wfs%ks%orbs%efermi
2761 :
2762 : !Copy occupations from BigDFT to ABINIT variables
2763 : call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,2,wvl%wfs)
2764 : #endif
2765 :
2766 : DBG_EXIT("COLL")
2767 :
2768 : end subroutine wvl_occ_bigdft
2769 : !!***
2770 :
2771 : !!****f* ABINIT/wvl_comm_eigen
2772 : !! NAME
2773 : !! wvl_comm_eigen
2774 : !!
2775 : !! FUNCTION
2776 : !! Computes occupations for the wavelet case
2777 : !! Using BigDFT routines
2778 : !!
2779 : !! NOTES
2780 : !! for the wvlbigdft case, see the routine 'wvl_occ_bigdft'
2781 : !!
2782 : !! SOURCE
2783 :
2784 : subroutine wvl_comm_eigen()
2785 :
2786 : !Arguments ------------------------------------
2787 :
2788 : !Local variables-------------------------------
2789 : #if defined HAVE_BIGDFT
2790 : integer:: ikpt,norb,shift
2791 : #endif
2792 : ! *************************************************************************
2793 :
2794 : DBG_ENTER("COLL")
2795 :
2796 : #if defined HAVE_BIGDFT
2797 : if(wvlbigdft) then
2798 : ! Communicates eigenvalues to all procs.
2799 : ! This will print out the eigenvalues and Fermi level.
2800 : call eigensystem_info(mpi_enreg%me_wvl, mpi_enreg%nproc_wvl,0.d0,&
2801 : & wvl%wfs%ks%Lzd%Glr%wfd%nvctr_c+7*wvl%wfs%ks%Lzd%Glr%wfd%nvctr_f,&
2802 : & wvl%wfs%ks%orbs,wvl%wfs%ks%psi)
2803 : else
2804 : ! Send all eigenvalues to all procs.
2805 : ! I simply communicate eigenvalues: I do not print them into screen, nor calculate Fermi-level.
2806 : norb=wvl%wfs%ks%orbs%norb
2807 : if (mpi_enreg%nproc_wvl > 1) then
2808 : shift=1
2809 : do ikpt = 1, wvl%wfs%ks%orbs%nkpts
2810 : call xmpi_bcast(wvl%wfs%ks%orbs%eval(shift:shift+norb-1),wvl%wfs%ks%orbs%ikptproc(ikpt),mpi_enreg%comm_wvl,ierr)
2811 : shift=shift+norb
2812 : end do
2813 : end if
2814 : end if
2815 :
2816 : !Copy eigenvalues from BigDFT object to "eigen"
2817 : call wvl_eigen_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,eigen,2,wvl%wfs)
2818 :
2819 : #else
2820 : BIGDFT_NOTENABLED_ERROR()
2821 : #endif
2822 :
2823 : DBG_EXIT("COLL")
2824 :
2825 : end subroutine wvl_comm_eigen
2826 :
2827 : end subroutine vtorho
2828 : !!***
2829 :
2830 : !!****f* ABINIT/cgq_builder
2831 : !! NAME
2832 : !! cgq_builder
2833 : !!
2834 : !! FUNCTION
2835 : !! This routine locates cgq for efield calculations, especially for parallel case
2836 : !!
2837 : !! INPUTS
2838 : !! berryflag = logical flag determining use of electric field variables
2839 : !! cg(2,mcg)=planewave coefficients of wavefunctions.
2840 : !! dtset <type(dataset_type)>=all input variables for this dataset
2841 : !! ikpt=index of current k kpt
2842 : !! ikpt_loc=index of k point on current processor (see vtorho.F90)
2843 : !! isspol=value of spin polarization currently treated
2844 : !! me_distrb=current value from spaceComm_distrb (see vtorho.F90)
2845 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
2846 : !! mcgq=size of cgq array (see vtorho.F90)
2847 : !! mkgq=size of pwnsfacq array (see vtorho.F90)
2848 : !! my_nspinor=nspinor value determined by current // set up
2849 : !! nband_k=number of bands at each k point
2850 : !! nproc_distrb=nproc from spaceComm_distrb (see vtorho.F90)
2851 : !! npwarr(nkpt)=number of planewaves in basis at this k point
2852 : !! pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
2853 : !! (see initberry.f)
2854 : !! pwind_alloc = first dimension of pwind
2855 : !! spaceComm_distrb=comm_cell from mpi_enreg
2856 : !!
2857 : !! OUTPUT
2858 : !! cgq(2,mcgq)=planewave coefficients of wavenfunctions adjacent to cg at ikpt
2859 : !! pwnsfacq(2,mkgq)=phase factors for non-symmorphic translations for cg's adjacent to cg(ikpt)
2860 : !!
2861 : !! SIDE EFFECTS
2862 : !! Input/Output
2863 : !! dtefield <type(efield_type)> = efield variables
2864 : !! mpi_enreg=information about MPI parallelization
2865 : !!
2866 : !! SOURCE
2867 :
2868 6474 : subroutine cgq_builder(berryflag,cg,cgq,dtefield,dtset,ikpt,ikpt_loc,isppol,mcg,mcgq,&
2869 : me_distrb,mkgq,mpi_enreg,my_nspinor,nband_k,nproc_distrb,&
2870 6474 : npwarr,pwnsfac,pwnsfacq,pwind_alloc,spaceComm_distrb)
2871 :
2872 : !Arguments ------------------------------------
2873 : integer,intent(in) :: ikpt,ikpt_loc,isppol,me_distrb,mcg,mcgq,mkgq,my_nspinor,nband_k
2874 : integer,intent(in) :: nproc_distrb,pwind_alloc,spaceComm_distrb
2875 : logical,intent(in) :: berryflag
2876 : type(dataset_type), intent(in) :: dtset
2877 : type(efield_type), intent(inout) :: dtefield
2878 : type(MPI_type), intent(in) :: mpi_enreg
2879 : !arrays
2880 : integer,intent(in) :: npwarr(dtset%nkpt)
2881 : real(dp),intent(in) :: cg(2,mcg),pwnsfac(2,pwind_alloc)
2882 : real(dp),intent(out) :: cgq(2,mcgq),pwnsfacq(2,mkgq)
2883 :
2884 : !Local variables -------------------------
2885 : !scalars
2886 : integer :: count,count1,icg1,icg2,dest,his_source
2887 : integer :: idir,ierr,ifor,ikg1,ikg2,ikptf,ikpt1f,ikpt1i
2888 : integer :: jkpt,jkpt1i,jkptf,jkpt1f,jsppol,my_source,npw_k1,tag
2889 : !arrays
2890 6474 : integer,allocatable :: flag_send(:,:), flag_receive(:)
2891 : real(dp) :: tsec(2)
2892 6474 : real(dp),allocatable :: buffer(:,:)
2893 : ! *************************************************************************
2894 :
2895 6474 : if (mcgq==0.or.mkgq==0) return
2896 :
2897 3056 : call timab(983,1,tsec)
2898 :
2899 : !Test compatbility of berryflag
2900 3056 : if (berryflag) then
2901 12224 : ABI_MALLOC(flag_send,(0:nproc_distrb-1,dtefield%fnkpt))
2902 : end if
2903 9168 : ABI_MALLOC(flag_receive,(dtset%nkpt))
2904 23450608 : flag_send(:,:) = 0
2905 589680 : flag_receive(:) = 0
2906 :
2907 3056 : if (berryflag) ikptf = dtefield%i2fbz(ikpt)
2908 :
2909 12224 : do idir = 1, 3
2910 :
2911 : ! skip idir values for which efield_dot(idir) = 0
2912 9168 : if (berryflag .and. abs(dtefield%efield_dot(idir)) < tol12 ) cycle
2913 :
2914 30560 : do ifor = 1, 2
2915 :
2916 18336 : if(berryflag) then
2917 127680 : dtefield%sflag(:,ikpt + dtset%nkpt*(isppol - 1),ifor,idir) = 0
2918 18336 : ikpt1f = dtefield%ikpt_dk(ikptf,ifor,idir)
2919 18336 : ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1)
2920 : end if
2921 :
2922 18336 : npw_k1 = npwarr(ikpt1i)
2923 18336 : count = npw_k1*my_nspinor*nband_k
2924 18336 : my_source = mpi_enreg%proc_distrb(ikpt1i,1,isppol)
2925 :
2926 82512 : do dest = 0, nproc_distrb-1
2927 :
2928 73344 : if ((dest==me_distrb).and.(ikpt_loc <= dtset%mkmem)) then
2929 : ! I am dest and have something to do
2930 :
2931 18336 : if ( my_source == me_distrb ) then
2932 : ! I am destination and source
2933 :
2934 11725 : if(berryflag) then
2935 11725 : ikg1 = dtefield%fkgindex(ikpt1f)
2936 11725 : ikg2 = dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*dtset%nkpt)
2937 11725 : icg1 = dtefield%cgindex(ikpt1i,isppol)
2938 11725 : icg2 = dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*dtset%nkpt)
2939 : end if
2940 :
2941 2980315 : pwnsfacq(:,ikg2 + 1:ikg2 + npw_k1) = pwnsfac(:,ikg1 + 1:ikg1 + npw_k1)
2942 21529135 : cgq(:,icg2 + 1:icg2 + count) = cg(:,icg1 + 1:icg1 + count)
2943 :
2944 : else ! I am the destination but not the source -> receive
2945 : ! receive pwnsfacq
2946 6611 : if(berryflag) then
2947 6611 : tag = ikpt1f + (isppol - 1)*dtefield%fnkpt
2948 6611 : ikg1 = dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*dtset%nkpt)
2949 : end if
2950 19833 : ABI_MALLOC(buffer,(2,npw_k1))
2951 6611 : call xmpi_recv(buffer,my_source,tag,spaceComm_distrb,ierr)
2952 1738799 : pwnsfacq(:,ikg1+1:ikg1+npw_k1) = buffer(:,1:npw_k1)
2953 6611 : ABI_FREE(buffer)
2954 :
2955 : ! receive cgq if necessary
2956 13222 : if(flag_receive(ikpt1i) == 0) then
2957 17748 : ABI_MALLOC(buffer,(2,count))
2958 5916 : tag = ikpt1i + (isppol - 1)*dtset%nkpt
2959 5916 : call xmpi_recv(buffer,my_source,tag,spaceComm_distrb,ierr)
2960 5916 : if(berryflag) icg1 = dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*dtset%nkpt)
2961 13360398 : cgq(:,icg1+1:icg1+count) = buffer(:,1:count)
2962 5916 : ABI_FREE(buffer)
2963 11832 : flag_receive(ikpt1i) = 1
2964 : end if ! end if flag_receive == 0
2965 : end if ! end tasks if I am the destination
2966 :
2967 36672 : else if (ikpt_loc <= mpi_enreg%mkmem(dest)) then ! dest != me and the dest has a k-point to treat
2968 :
2969 : ! jkpt is the kpt which is being treated by dest (in ibz)
2970 : ! jsppol is his isppol
2971 36672 : jkpt = mpi_enreg%kpt_loc2ibz_sp(dest, ikpt_loc,1)
2972 36672 : jsppol = mpi_enreg%kpt_loc2ibz_sp(dest, ikpt_loc,2)
2973 :
2974 36672 : if(jkpt > 0 .and. jsppol > 0) then
2975 :
2976 36672 : if(berryflag) then
2977 36672 : jkptf = dtefield%i2fbz(jkpt)
2978 36672 : jkpt1f = dtefield%ikpt_dk(jkptf,ifor,idir)
2979 36672 : jkpt1i = dtefield%indkk_f2ibz(jkpt1f,1)
2980 : end if
2981 36672 : his_source = mpi_enreg%proc_distrb(jkpt1i,1,jsppol)
2982 :
2983 36672 : if (his_source == me_distrb) then
2984 :
2985 : ! send
2986 : ! pwnsfacq
2987 6611 : if(berryflag) then
2988 6611 : ikg1 = dtefield%fkgindex(jkpt1f)
2989 6611 : tag = jkpt1f + (jsppol - 1)*dtefield%fnkpt
2990 : end if
2991 6611 : count1 = npwarr(jkpt1i)
2992 19833 : ABI_MALLOC(buffer,(2,count1))
2993 1738799 : buffer(:,1:count1) = pwnsfac(:,ikg1+1:ikg1+count1)
2994 6611 : call xmpi_send(buffer,dest,tag,spaceComm_distrb,ierr)
2995 6611 : ABI_FREE(buffer)
2996 :
2997 : ! send cgq if necessary
2998 13222 : if(flag_send(dest, jkpt1i)==0) then
2999 5916 : if(berryflag) icg1 = dtefield%cgindex(jkpt1i,jsppol)
3000 5916 : tag = jkpt1i + (jsppol - 1)*dtset%nkpt
3001 5916 : count1 = npwarr(jkpt1i)*nband_k*my_nspinor
3002 17748 : ABI_MALLOC(buffer,(2,count1))
3003 13360398 : buffer(:,1:count1) = cg(:,icg1+1:icg1+count1)
3004 5916 : call xmpi_send(buffer,dest,tag,spaceComm_distrb,ierr)
3005 5916 : ABI_FREE(buffer)
3006 11832 : flag_send(dest, jkpt1i)=1
3007 : end if ! if send cgq
3008 :
3009 : end if ! end check that his_source == me
3010 : end if ! end check on jkpt > 0 and jsppol > 0
3011 : end if ! end check on me = dest else if me != dest
3012 : end do ! end loop over dest = 0, nproc-1
3013 : end do !end loop over ifor
3014 : end do !end loop over idir
3015 :
3016 3056 : call timab(983,2,tsec)
3017 :
3018 3056 : ABI_FREE(flag_send)
3019 3056 : ABI_FREE(flag_receive)
3020 :
3021 3056 : end subroutine cgq_builder
3022 : !!***
3023 :
3024 : end module m_vtorho
3025 : !!***
|
