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