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