Branch data Line data Source code
1 : : !!****m* ABINIT/m_vtowfk
2 : : !! NAME
3 : : !! m_vtowfk
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, MT)
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_vtowfk
26 : :
27 : : use, intrinsic :: iso_fortran_env, only: int32, int64, real32, real64
28 : : use, intrinsic :: iso_c_binding, only: c_size_t, c_loc
29 : :
30 : : use defs_basis
31 : : use m_abicore
32 : : use m_errors
33 : : use m_xmpi
34 : : use m_efield
35 : : use m_linalg_interfaces
36 : : use m_cgtools
37 : : use m_dtset
38 : : use m_dtfil
39 : : use m_xomp
40 : : use m_xg
41 : : use m_xg_nonlop
42 : :
43 : : use defs_abitypes, only : MPI_type
44 : : use m_time, only : timab, cwtime, cwtime_report, sec2str
45 : : use m_fstrings, only : sjoin, itoa, ftoa
46 : : use m_hamiltonian, only : gs_hamiltonian_type, K_H_K, KPRIME_H_KPRIME
47 : : use m_getghc, only : getghc_nucdip
48 : : use m_paw_dmft, only : paw_dmft_type
49 : : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_put,pawcprj_copy
50 : : use m_paw_dmft, only : paw_dmft_type
51 : : use m_gwls_hamiltonian, only : build_H
52 : : use m_fftcore, only : fftcore_set_mixprec, fftcore_mixprec
53 : : use m_cgwf, only : cgwf
54 : : use m_cgwf_cprj, only : cgwf_cprj,mksubovl,cprj_update,cprj_update_oneband
55 : : use m_lobpcgwf_old,only : lobpcgwf
56 : : use m_lobpcgwf, only : lobpcgwf2
57 : : use m_chebfiwf, only : chebfiwf2
58 : : use m_chebfiwf_cprj,only : chebfiwf2_cprj
59 : : use m_lobpcgwf_cprj,only : lobpcgwf2_cprj
60 : : use m_slicewf, only : slicewf
61 : : use m_slicewf_cprj, only : slicewf_cprj
62 : : use m_spacepar, only : meanvalue_g, meanvalue_g_batch
63 : : use m_chebfi, only : chebfi
64 : : use m_rmm_diis, only : rmm_diis
65 : : use m_nonlop, only : nonlop !, nonlop_counter
66 : : use m_prep_kgb, only : prep_nonlop, prep_fourwf
67 : : use m_cgprj, only : cprj_rotate,xg_cprj_copy,XG_TO_CPRJ
68 : : use m_fft, only : fourwf, fourwf_optmem
69 : : use m_cgtk, only : cgtk_fixphase
70 : : use m_common, only : get_gemm_nonlop_ompgpu_blocksize
71 : : use m_gemm_nonlop_projectors, only : gemm_nonlop_block_size, gemm_nonlop_is_distributed
72 : : use m_gputk, only : gpu_copy
73 : : #if defined HAVE_YAKL
74 : : use gator_mod
75 : : #endif
76 : : #if defined(HAVE_GPU_MARKERS)
77 : : use m_nvtx_data
78 : : #endif
79 : :
80 : : implicit none
81 : :
82 : : private
83 : : !!***
84 : :
85 : : public :: vtowfk
86 : : !!***
87 : :
88 : : contains
89 : : !!***
90 : :
91 : : !!****f* ABINIT/vtowfk
92 : : !! NAME
93 : : !! vtowfk
94 : : !!
95 : : !! FUNCTION
96 : : !! This routine compute the partial density at a given k-point,
97 : : !! for a given spin-polarization, from a fixed Hamiltonian
98 : : !! but might also simply compute eigenvectors and eigenvalues at this k point
99 : : !!
100 : : !! INPUTS
101 : : !! cgq = array that holds the WF of the nearest neighbours of
102 : : !! the current k-point (electric field, MPI //)
103 : : !! cpus= cpu time limit in seconds
104 : : !! dtefield <type(efield_type)> = variables related to Berry phase
105 : : !! calculations (see initberry.f)
106 : : !! dtfil <type(datafiles_type)>=variables related to files
107 : : !! dtset <type(dataset_type)>=all input variables for this dataset
108 : : !! fixed_occ=true if electronic occupations are fixed (occopt<3)
109 : : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
110 : : !! ibg=shift to be applied on the location of data in the array cprj
111 : : !! icg=shift to be applied on the location of data in the array cg
112 : : !! ikpt=number of the k-point
113 : : !! iscf=(<= 0 =>non-SCF), >0 => SCF
114 : : !! isppol= 1 for unpolarized, 2 for spin-polarized
115 : : !! kg_k(3,npw_k)=reduced planewave coordinates.
116 : : !! kinpw(npw_k)=(modified) kinetic energy for each plane wave (Hartree)
117 : : !! mcg=second dimension of the cg array
118 : : !! mcgq=second dimension of the cgq array (electric field, MPI //)
119 : : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
120 : : !! mkgq = second dimension of pwnsfacq
121 : : !! mpi_enreg=information about MPI parallelization
122 : : !! mpw=maximum dimensioned size of npw
123 : : !! natom=number of atoms in cell.
124 : : !! nband_k=number of bands at this k point for that spin polarization
125 : : !! nkpt=number of k points.
126 : : !! istep=index of the number of steps in the routine scfcv
127 : : !! nnsclo_now=number of non-self-consistent loops for the current vtrial
128 : : !! (often 1 for SCF calculation, =nstep for non-SCF calculations)
129 : : !! npw_k=number of plane waves at this k point
130 : : !! npwarr(nkpt)=number of planewaves in basis at this k point
131 : : !! occ_k(nband_k)=occupation number for each band (usually 2) for each k.
132 : : !! optforces=option for the computation of forces
133 : : !! prtvol=control print volume and debugging output
134 : : !! pwind(pwind_alloc,2,3)= array used to compute
135 : : !! the overlap matrix smat between k-points (see initberry.f)
136 : : !! pwind_alloc= first dimension of pwind
137 : : !! pwnsfac(2,pwind_alloc)= phase factors for non-symmorphic translations
138 : : !! (see initberry.f)
139 : : !! pwnsfacq(2,mkgq)= phase factors for the nearest neighbours of the
140 : : !! current k-point (electric field, MPI //)
141 : : !! usebanfft=flag for band-fft parallelism
142 : : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
143 : : !! wtk=weight assigned to the k point.
144 : : !! zshift(nband_k)=energy shifts for the squared shifted hamiltonian algorithm
145 : : !!
146 : : !! OUTPUT
147 : : !! dphase_k(3)=change in Zak phase for the current k-point
148 : : !! eig_k(nband_k)=array for holding eigenvalues (hartree)
149 : : !! ek_k(nband_k)=contribution from each band to kinetic energy, at this k-point
150 : : !! ek_k_nd(2,nband_k,nband_k*use_dmft)=contribution to kinetic energy,
151 : : !! including non-diagonal terms, at this k-point (useful if use_dmft)
152 : : !! end_k(nband_k)=contribution from each band to nuclear dipole energy, at this k-point
153 : : !! resid_k(nband_k)=residuals for each band over all k points, BEFORE the band rotation.
154 : : !! In input: previous residuals.
155 : : !! ==== if optforces>0 ====
156 : : !! grnl_k(3*natom,nband_k)=nonlocal gradients, at this k-point
157 : : !! ==== if gs_hamk%usepaw==0 ====
158 : : !! enlx_k(nband_k)=contribution from each band to
159 : : !! nonlocal pseudopotential + Fock-type part of total energy, at this k-point
160 : : !! ==== if (gs_hamk%usepaw==1) ====
161 : : !! cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors:
162 : : !! cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector.
163 : : !!
164 : : !! SIDE EFFECTS
165 : : !! cg(2,mcg)=updated wavefunctions
166 : : !! rhoaug(n4,n5,n6,nvloc)= density in electrons/bohr**3, on the augmented fft grid.
167 : : !! (cumulative, so input as well as output). Update only
168 : : !! for occopt<3 (fixed occupation numbers)
169 : : !! rmm_diis_status: Status of the RMM-DIIS eigensolver. See m_rmm_diis.
170 : : !!
171 : : !! NOTES
172 : : !! The cprj are distributed over band and spinors processors.
173 : : !! One processor doesn't know all the cprj.
174 : : !! Only the mod((iband-1)/mpi_enreg%bandpp,mpi_enreg%nproc_band) projectors
175 : : !! are stored on each proc.
176 : : !!
177 : : !! SOURCE
178 : :
179 : 693002 : subroutine vtowfk(cg,cgq,cprj,cpus,dphase_k,dtefield,dtfil,dtset,&
180 : 346501 : & eig_k,ek_k,ek_k_nd,end_k,enlx_k,fixed_occ,grnl_k,gs_hamk,&
181 : 693002 : & ibg,icg,ikpt,iscf,isppol,kg_k,kinpw,kinpw_kp,mband_cprj,mcg,mcgq,mcprj,mkgq,mpi_enreg,&
182 : 346501 : & mpw,natom,nband_k,nbdbuf,nkpt,istep,nnsclo_now,npw_k,npwarr,occ_k,optforces,prtvol,&
183 [ - + ]: 346501 : & pwind,pwind_alloc,pwnsfac,pwnsfacq,resid_k,rhoaug,paw_dmft,wtk,xg_nonlop,zshift,rmm_diis_status)
184 : :
185 : : !Arguments ------------------------------------
186 : : integer, intent(in) :: ibg,icg,ikpt,iscf,isppol,mband_cprj,mcg,mcgq,mcprj,mkgq,mpw
187 : : integer, intent(in) :: natom,nband_k,nbdbuf,nkpt,nnsclo_now,npw_k,optforces
188 : : integer, intent(in) :: prtvol,pwind_alloc,istep
189 : : logical,intent(in) :: fixed_occ
190 : : real(dp), intent(in) :: cpus,wtk
191 : : type(datafiles_type), intent(in) :: dtfil
192 : : type(efield_type), intent(inout) :: dtefield
193 : : type(dataset_type), intent(in) :: dtset
194 : : type(gs_hamiltonian_type), intent(inout) :: gs_hamk
195 : : type(MPI_type), intent(inout) :: mpi_enreg
196 : : type(paw_dmft_type), intent(in) :: paw_dmft
197 : : integer, intent(in) :: kg_k(3,npw_k)
198 : : integer, intent(in) :: npwarr(nkpt),pwind(pwind_alloc,2,3)
199 : : integer, intent(inout) :: rmm_diis_status(2)
200 : : real(dp), intent(in) :: cgq(2,mcgq),occ_k(nband_k)
201 : : real(dp),target,intent(in) :: kinpw(npw_k)
202 : : real(dp),target,intent(in) :: kinpw_kp(:)
203 : : real(dp), intent(in) :: pwnsfac(2,pwind_alloc),pwnsfacq(2,mkgq)
204 : : real(dp), intent(in) :: zshift(nband_k)
205 : : real(dp), target, intent(out) :: eig_k(nband_k)
206 : : real(dp), intent(out) :: ek_k(nband_k),dphase_k(3),ek_k_nd(2,nband_k,nband_k*paw_dmft%use_dmft)
207 : : real(dp), intent(out) :: end_k(nband_k),enlx_k(nband_k)
208 : : real(dp), intent(out),target :: grnl_k(3*natom,nband_k*optforces)
209 : : real(dp), intent(inout) :: resid_k(nband_k)
210 : : real(dp), intent(inout),target :: cg(2,mcg)
211 : : real(dp), intent(inout) :: rhoaug(gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,gs_hamk%nvloc)
212 : : type(pawcprj_type),intent(inout),target :: cprj(natom,mcprj*gs_hamk%usecprj)
213 : : type(xg_nonlop_t),intent(in) :: xg_nonlop
214 : :
215 : : !Local variables-------------------------------
216 : : logical :: has_fock,xg_diago,update_cprj
217 : : logical :: do_subdiago,do_ortho,rotate_subvnlx,use_rmm_diis,is_distrib_tmp
218 : : integer,parameter :: level=112,tim_fourwf=2,tim_nonlop_prep=11,enough=3,tim_getcprj=5
219 : : integer,save :: nskip=0
220 : : ! Flag use_subovl: 1 if "subovl" array is computed (see below)
221 : : ! subovl should be Identity (in that case we should use use_subovl=0)
222 : : ! But this is true only if conjugate gradient algo. converges
223 : : integer :: use_subovl=0, use_subvnlx=0, use_totvnlx=0
224 : : integer :: bandpp_cprj,blocksize,choice,cpopt,fftalg,iband,iband1,filter
225 : : integer :: nstep_mixed
226 : : integer :: iblock,iblocksize,ibs,idir,ierr,igs,igsc,ii,inonsc
227 : : integer :: iorder_cprj,ipw,ispinor,iispinor,ispinor_index,istwf_k,iwavef,me_g0,mgsc,my_nspinor,n1,n2,n3 !kk
228 : : integer :: nband_k_cprj,ncols_cprj,nblockbd,ncpgr,ndat,niter,nkpt_max,nnlout,ortalgo,ndat_fft
229 : : integer :: paw_opt,quit,signs,space,spaceComm,tim_nonlop,wfoptalg,wfopta10
230 : : integer :: gpu_option_tmp,nblk_gemm_nonlop,blksize_gemm_nonlop_tmp,nfft_blocks_tmp
231 : : integer :: chunk,residuchunk
232 : : logical :: nspinor1TreatedByThisProc,nspinor2TreatedByThisProc
233 : : logical :: transfer_cg
234 : : real(dp) :: ar,ar2,ar_im,eshift,occblock,norm
235 : : real(dp) :: max_resid,weight,cpu,wall,gflops
236 : : character(len=50) :: iter_name
237 : : character(len=500) :: msg
238 : : real(dp) :: dummy(2,1),nonlop_dum(1,1),nonlop_dum2(1,1),tsec(2)
239 : 346501 : real(dp),allocatable :: cwavef1(:,:),cwavef_spin(:,:),cwavef_x(:,:),cwavef_y(:,:),cwavefb(:,:,:)
240 : : #if defined HAVE_GPU && defined HAVE_YAKL
241 : : real(real64), contiguous, pointer :: cwavef(:,:) => null()
242 : : real(real64), contiguous, pointer :: wfraug(:,:,:,:)
243 : : #else
244 : 346501 : real(dp),allocatable,target :: cwavef(:,:)
245 : 346501 : real(dp),allocatable :: wfraug(:,:,:,:)
246 : : #endif
247 : 346501 : real(dp),allocatable :: eig_save(:),enlout(:),enlout_spin(:),evec(:,:),gsc(:,:),ghc_vectornd(:,:)
248 : 346501 : real(dp),allocatable :: subham(:),subovl(:),subvnlx(:),totvnlx(:,:)
249 : 346501 : real(dp),pointer :: cg_k(:,:),cg_k_block(:,:),grnl_k_block(:,:),eig_k_block(:)
250 : 346501 : real(dp),contiguous, pointer :: cwavef_iband(:,:), kinpw_ptr(:)
251 : 346501 : real(dp), allocatable :: weight_t(:) ! only allocated and used with GPU fourwf
252 : 346501 : type(pawcprj_type),pointer :: cwaveprj(:,:), cprj_cwavef_bands(:,:),cprj_cwavef(:,:)
253 : : type(xgBlock_t) :: xgx0,xgeigen,xgforces
254 : : type(xg_t) :: cprj_xgx0,cprj_work
255 : : ! **********************************************************************
256 : :
257 : : DBG_ENTER("COLL")
258 : :
259 : 346501 : call timab(28,1,tsec) ! Keep track of total time spent in "vtowfk"
260 : :
261 : : !Structured debugging if prtvol==-level
262 [ - + ]: 346501 : if(prtvol==-level)then
263 [ # # # # ]: 0 : write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,'vtowfk: enter'
264 : 0 : call wrtout(std_out,msg,'PERS')
265 : : end if
266 : :
267 : : !=========================================================================
268 : : !============= INITIALIZATIONS AND ALLOCATIONS ===========================
269 : : !=========================================================================
270 : :
271 : 346501 : nkpt_max=50; if(xmpi_paral==1)nkpt_max=-1
272 : :
273 : 346501 : wfoptalg=mod(dtset%wfoptalg,100); wfopta10=mod(wfoptalg,10)
274 [ + + + + ]: 346501 : xg_diago = dtset%wfoptalg == 114 .or. dtset%wfoptalg == 112 .or. dtset%wfoptalg == 111
275 : 346501 : istwf_k=gs_hamk%istwf_k
276 : 346501 : has_fock=(associated(gs_hamk%fockcommon))
277 : 346501 : quit=0
278 : 346501 : igsc=0
279 : :
280 : : !Parallelization over spinors management
281 : 346501 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
282 [ + + ]: 346501 : if (mpi_enreg%paral_spinor==0) then
283 : 345941 : ispinor_index=1
284 : 345941 : nspinor1TreatedByThisProc=.true.
285 : 345941 : nspinor2TreatedByThisProc=(dtset%nspinor==2)
286 : : else
287 : 560 : ispinor_index=mpi_enreg%me_spinor+1
288 : 560 : nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
289 : 560 : nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
290 : : end if
291 : :
292 : : !Parallelism over FFT and/or bands: define sizes and tabs
293 : : !if (mpi_enreg%paral_kgb==1) then
294 : 346501 : nblockbd=nband_k/(mpi_enreg%nproc_band*mpi_enreg%bandpp)
295 : : !else
296 : : ! nblockbd=nband_k/mpi_enreg%nproc_fft
297 : : ! if (nband_k/=nblockbd*mpi_enreg%nproc_fft) nblockbd=nblockbd+1
298 : : !end if
299 : 346501 : blocksize=nband_k/nblockbd
300 : :
301 : : !Save eshift
302 [ + + ]: 346501 : if(wfoptalg==3)then
303 : 6 : eshift=zshift(1)
304 [ - + + - : 18 : ABI_MALLOC(eig_save,(nband_k))
+ - + - -
+ - + ]
305 [ + + ]: 30 : eig_save(:)=eshift
306 : : end if
307 : :
308 : 346501 : n1=gs_hamk%ngfft(1); n2=gs_hamk%ngfft(2); n3=gs_hamk%ngfft(3)
309 : :
310 : : ! Decide whether RMM-DIIS eigensolver should be activated.
311 : : ! rmm_diis > 0 --> Activate it after (3 + rmm_diis) iterations with wfoptalg algorithm.
312 : : ! rmm_diis < 0 --> Start with RMM-DIIS directly (risky)
313 : 346501 : use_rmm_diis = .False.
314 [ + + + + ]: 346501 : if (dtset%rmm_diis /= 0 .and. iscf > 0) then
315 : 880 : use_rmm_diis = istep > 3 + dtset%rmm_diis
316 : : !if (use_rmm_diis) call wrtout(std_out, " Activating RMM-DIIS eigensolver in SCF mode.")
317 : : end if
318 : : !nonlop_counter = 0
319 : :
320 : 346501 : mgsc=0
321 : : igsc=0
322 [ + + - + : 346501 : if ((.not. xg_diago .and. dtset%cprj_in_memory==0) .or. dtset%rmm_diis /= 0) then
+ + ]
323 : 300482 : mgsc=nband_k*npw_k*my_nspinor*gs_hamk%usepaw
324 [ + + + - : 901446 : ABI_MALLOC_OR_DIE(gsc,(2,mgsc), ierr)
+ - + + +
- - + -
+ ]
325 [ + + + + ]: 986143307 : gsc=zero
326 : : else
327 [ - + ]: 46019 : ABI_MALLOC(gsc,(0,0))
328 : : end if
329 : :
330 [ + + ]: 346501 : if(wfopta10 /= 1 .and. .not. xg_diago) then
331 : : !chebfi already does this stuff inside
332 [ + - - + : 1498950 : ABI_MALLOC(evec,(2*nband_k,nband_k))
+ - - + +
- + - + -
- + - + ]
333 [ - + + - : 899370 : ABI_MALLOC(subham,(nband_k*(nband_k+1)))
+ - + - -
+ - + ]
334 : :
335 [ - + ]: 299790 : ABI_MALLOC(subvnlx,(0))
336 [ - + ]: 299790 : ABI_MALLOC(totvnlx,(0,0))
337 [ + + ]: 299790 : if (wfopta10==4) then
338 : : ! Later, will have to generalize to Fock case, like when wfopta10/=4
339 [ + + ]: 1966 : if (gs_hamk%usepaw==0) then
340 [ - + ]: 612 : ABI_FREE(totvnlx)
341 [ + + ]: 612 : if (istwf_k==1) then
342 [ - + + - : 1608 : ABI_MALLOC(totvnlx,(2*nband_k,nband_k))
- + + - -
+ - + ]
343 [ + - ]: 76 : else if (istwf_k==2) then
344 [ - + + - : 304 : ABI_MALLOC(totvnlx,(nband_k,nband_k))
- + + - +
- + - - +
- + ]
345 : : end if
346 : 612 : use_totvnlx=1
347 : : endif
348 : : else
349 [ + + + + ]: 297824 : if (gs_hamk%usepaw==0 .or. has_fock) then
350 [ - + ]: 224818 : ABI_FREE(subvnlx)
351 [ - + + - : 449636 : ABI_MALLOC(subvnlx,(nband_k*(nband_k+1)))
- + - + ]
352 : 224818 : use_subvnlx=1
353 : : end if
354 : : end if
355 : :
356 [ - + ]: 299790 : if (use_subovl==1) then
357 [ # # # # : 0 : ABI_MALLOC(subovl,(nband_k*(nband_k+1)))
# # # # ]
358 : : else
359 [ - + ]: 299790 : ABI_MALLOC(subovl,(0))
360 : : end if
361 : : end if
362 : :
363 : : ! Carry out UP TO dtset%nline (or dtset%mdeg_filter) steps, or until resid for every band is < dtset%tolwfr
364 [ + + + + : 346501 : if (prtvol/=5 .and. (prtvol>2 .or. ikpt <= nkpt_max)) then
- + ]
365 : 14055 : write(msg,'(a,i5,2x,a,3f9.5,2x,a)')' non-scf iterations; kpt # ',ikpt,', k= (',gs_hamk%kpt_k,'), band residuals:'
366 : 14055 : call wrtout(std_out,msg,'PERS')
367 : : end if
368 : :
369 [ - + ]: 346501 : if (dtset%cprj_in_memory==2) then
370 [ # # ]: 0 : if (ikpt==1) then
371 : 0 : write(msg,'(a,i3)') ' In vtowfk : use of cprj in memory with cprj_update_lvl=',dtset%cprj_update_lvl
372 : 0 : call wrtout(std_out,msg,'COLL')
373 : : end if
374 : 0 : cprj_cwavef_bands => cprj(:,1+ibg:nband_k/mpi_enreg%nproc_band*my_nspinor+ibg)
375 : : end if
376 : :
377 : : !Electric field: initialize dphase_k
378 : 346501 : dphase_k(:) = zero
379 : :
380 : : !=========================================================================
381 : : !==================== NON-SELF-CONSISTENT LOOP ===========================
382 : : !=========================================================================
383 : :
384 : : !nnsclo_now=number of non-self-consistent loops for the current vtrial
385 : : !(often 1 for SCF calculation, =nstep for non-SCF calculations)
386 : 346501 : call timab(39,1,tsec) ! "vtowfk (loop)"
387 : :
388 : 346501 : cg_k => cg(:,1+icg:npw_k*my_nspinor*nband_k+icg)
389 : :
390 : 346501 : transfer_cg = .false.
391 : : #ifdef HAVE_OPENMP_OFFLOAD
392 : : transfer_cg = .not. xomp_target_is_present(c_loc(cg))
393 : : if(transfer_cg) then
394 : : if(xg_diago) then
395 : : !$OMP TARGET ENTER DATA MAP(alloc:cg_k) IF(dtset%gpu_option==ABI_GPU_OPENMP)
396 : : !$OMP TARGET UPDATE TO(cg_k) IF(dtset%gpu_option==ABI_GPU_OPENMP .and. .not. use_rmm_diis)
397 : : end if
398 : : else if(istep == 1) then
399 : : !$OMP TARGET UPDATE FROM(cg_k) IF(dtset%gpu_option==ABI_GPU_OPENMP .and. use_rmm_diis)
400 : : end if
401 : : #endif
402 : :
403 [ + + ]: 825734 : do inonsc=1,nnsclo_now
404 : : ABI_NVTX_START_RANGE(NVTX_VTOWFK_EXTRA1)
405 [ + + + + : 501686 : if (iscf < 0 .and. (inonsc <= enough .or. mod(inonsc, 10) == 0)) call cwtime(cpu, wall, gflops, "start")
+ + ]
406 : :
407 [ + + + + ]: 501686 : if (dtset%rmm_diis /= 0 .and. iscf < 0) then
408 : 65 : use_rmm_diis = inonsc > 3 + dtset%rmm_diis
409 : : !if (use_rmm_diis) call wrtout(std_out, " Activating RMM-DIIS eigensolver in NSCF mode.")
410 : : end if
411 : :
412 : : ! This initialisation is needed for the MPI-parallelisation (gathering using sum)
413 [ + + ]: 501686 : if(wfopta10 /= 1 .and. .not. xg_diago) then
414 [ + + ]: 62176601 : subham(:)=zero
415 [ + + ]: 446357 : if (gs_hamk%usepaw==0) then
416 [ + + ]: 346763 : if (wfopta10==4) then
417 [ + + + + ]: 87328 : totvnlx(:,:)=zero
418 : : else
419 [ + + ]: 34674373 : subvnlx(:)=zero
420 : : end if
421 : : end if
422 [ - + - - ]: 446357 : if (use_subovl==1)subovl(:)=zero
423 : : end if
424 : :
425 : : ! Filter the WFs when modified kinetic energy is too large (see routine mkkin.f)
426 : : ! !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(igs,iwavef)
427 [ + + ]: 4429200 : do iband=1,nband_k
428 : 3927514 : iwavef=(iband-1)*npw_k*my_nspinor+icg
429 : 3927514 : cwavef_iband => cg(:,1+iwavef:npw_k*my_nspinor+iwavef)
430 : 3927514 : update_cprj=.False.
431 [ + + ]: 8353236 : do ispinor=1,my_nspinor
432 [ + + ]: 4425722 : iispinor=ispinor; if (mpi_enreg%paral_spinor==1) iispinor=mpi_enreg%me_spinor+1
433 : 4425722 : igs=(ispinor-1)*npw_k
434 : 4425722 : kinpw_ptr => kinpw
435 [ + + ]: 4425722 : if (gs_hamk%use_gbt /= 0) then
436 : : ! Use different kinetic energies for up and down component (this won't work if paral_kgb /= 0)
437 : : ! cg is still in the distributed band layout here, so its G ordering
438 : : ! matches the local kinetic table passed by vtorho, not the gathered
439 : : ! table installed in gs_hamk for prep_getghc.
440 : : if (iispinor == 1) kinpw_ptr => kinpw
441 [ + + ]: 86912 : if (iispinor == 2) kinpw_ptr => kinpw_kp
442 : : end if
443 [ + + ]: 1184941552 : do ipw=1+igs,npw_k+igs
444 [ + + ]: 1181014038 : if(kinpw_ptr(ipw-igs)>huge(zero)*1.d-11)then
445 : 22215259 : norm=cwavef_iband(1,ipw)**2+cwavef_iband(2,ipw)**2
446 [ + + ]: 22215259 : if (norm>tol15*tol15) update_cprj=.True.
447 [ + + ]: 66645777 : cwavef_iband(:,ipw)=zero
448 : : end if
449 : : end do
450 : : end do
451 [ - + - - ]: 4429200 : if (dtset%cprj_in_memory==2.and.update_cprj) then
452 : 0 : cprj_cwavef => cprj_cwavef_bands(:,my_nspinor*(iband-1)+1:my_nspinor*iband)
453 : 0 : call cprj_update_oneband(cwavef_iband,cprj_cwavef,gs_hamk,mpi_enreg,tim_getcprj)
454 : : end if
455 : : end do
456 : : ABI_NVTX_END_RANGE()
457 : :
458 : : ! JLJ 17/10/2014: If it is a GWLS calculation, construct the hamiltonian
459 : : ! as in a usual GS calc., but skip any minimisation procedure.
460 : : ! This would be equivalent to nstep=0, if the latter did work.
461 [ + + ]: 501686 : if(dtset%optdriver/=RUNL_GWLS) then
462 : :
463 [ + + ]: 501678 : if(wfopta10==4.or.wfopta10==1) then
464 : :
465 [ - + ]: 57557 : if (dtset%gpu_option==ABI_GPU_KOKKOS) then
466 : : ! Kokkos GPU branch is not OpenMP thread-safe, setting OpenMP num threads to 1
467 : 0 : call xomp_set_num_threads(1)
468 : : end if
469 : :
470 : : ! =========================================================================
471 : : ! ============ MINIMIZATION OF BANDS: LOBPCG ==============================
472 : : ! =========================================================================
473 [ + + ]: 57557 : if (wfopta10==4) then
474 : :
475 [ + + ]: 30885 : if (use_rmm_diis) then
476 : : call rmm_diis(istep, ikpt, isppol, cg_k, dtset, eig_k, occ_k, enlx_k, gs_hamk, kinpw, gsc, &
477 [ - + ]: 428 : mpi_enreg, nband_k, npw_k, my_nspinor, resid_k, rmm_diis_status)
478 : : #ifdef HAVE_OPENMP_OFFLOAD
479 : : !$OMP TARGET UPDATE TO(cg_k) IF(dtset%gpu_option==ABI_GPU_OPENMP .and. xg_diago)
480 : : #endif
481 : : else
482 : :
483 [ + + ]: 30457 : if ( .not. xg_diago ) then
484 : :
485 : : ABI_NVTX_START_RANGE(NVTX_LOBPCG1)
486 : : call lobpcgwf(cg,dtset,gs_hamk,gsc,icg,igsc,kinpw,kinpw_kp,mcg,mgsc,mpi_enreg,&
487 : 2250 : nband_k,nblockbd,npw_k,prtvol,resid_k,subham,totvnlx,use_totvnlx)
488 : : ! In case of FFT parallelism, exchange subspace arrays
489 : 2250 : spaceComm=mpi_enreg%comm_bandspinorfft
490 : 2250 : call xmpi_sum(subham,spaceComm,ierr)
491 [ + + ]: 2250 : if (gs_hamk%usepaw==0) then
492 : 704 : if (wfopta10==4) then
493 : 704 : call xmpi_sum(totvnlx,spaceComm,ierr)
494 : : else
495 : : call xmpi_sum(subvnlx,spaceComm,ierr)
496 : : end if
497 : : end if
498 [ - + ]: 2250 : if (use_subovl==1) call xmpi_sum(subovl,spaceComm,ierr)
499 : : ABI_NVTX_END_RANGE()
500 : :
501 : : else
502 : :
503 : : ABI_NVTX_START_RANGE(NVTX_LOBPCG2)
504 [ + + ]: 28207 : if (dtset%cprj_in_memory==1) then
505 : : call lobpcgwf2_cprj(cg_k,dtset,eig_k,occ_k,enlx_k,gs_hamk,isppol,ikpt,inonsc,istep,&
506 [ - + ]: 4036 : kinpw,mpi_enreg,nband_k,npw_k,my_nspinor,prtvol,resid_k,nbdbuf,xg_nonlop)
507 : : else
508 : : call lobpcgwf2(cg_k,dtset,eig_k,occ_k,enlx_k,gs_hamk,isppol,ikpt,inonsc,istep,kinpw,kinpw_kp,mpi_enreg,&
509 [ - + ]: 24171 : nband_k,npw_k,my_nspinor,prtvol,resid_k,nbdbuf)
510 : : end if
511 : : ABI_NVTX_END_RANGE()
512 : :
513 : : end if
514 : :
515 : : end if
516 : :
517 : : ! =========================================================================
518 : : ! ============ MINIMIZATION OF BANDS: CHEBYSHEV FILTERING =================
519 : : ! =========================================================================
520 [ + - ]: 26672 : else if (wfopta10 == 1) then
521 [ + + ]: 26672 : if ( .not. xg_diago) then
522 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI1)
523 : : call chebfi(cg_k,dtset,eig_k,enlx_k,gs_hamk,gsc,kinpw,&
524 [ - + ]: 48 : & mpi_enreg,nband_k,npw_k,my_nspinor,prtvol,resid_k)
525 : : ABI_NVTX_END_RANGE()
526 [ + + ]: 26624 : else if (dtset%cprj_in_memory==1) then
527 : : call chebfiwf2_cprj(cg_k,dtset,eig_k,occ_k,enlx_k,gs_hamk,&
528 [ - + ]: 15368 : mpi_enreg,nband_k,npw_k,my_nspinor,prtvol,resid_k,xg_nonlop)
529 : : else
530 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2)
531 : : call chebfiwf2(cg_k,dtset,eig_k,occ_k,enlx_k,gs_hamk,&
532 [ - + ]: 11256 : & mpi_enreg,nband_k,npw_k,my_nspinor,prtvol,resid_k)
533 : : ABI_NVTX_END_RANGE()
534 : : end if
535 : : end if
536 : :
537 : : ! =========================================================================
538 : : ! ============ MINIMIZATION OF BANDS: SPECTRUM SLICING == =================
539 : : ! =========================================================================
540 [ + + ]: 444121 : else if (wfopta10 == 2) then
541 : 18 : nstep_mixed = dtset%nstep_mixed ! below which perform chebfi
542 : 18 : write(std_out,'(a,i0)') 'running vtowfk for nstep_mixed=', nstep_mixed
543 [ + - + - ]: 18 : if ( xg_diago .and. dtset%cprj_in_memory == 0 ) then
544 [ + - ]: 18 : if (istep > nstep_mixed) then
545 : 18 : write(std_out,'(a,i0)') 'entering slicewf'
546 : : !ABI_NVTX_START_RANGE(NVTX_SPESLI)
547 : : call slicewf(cg_k,dtset,eig_k,enlx_k,gs_hamk,mpi_enreg,&
548 [ - + ]: 18 : & nband_k,npw_k,my_nspinor,prtvol,resid_k)
549 : : !ABI_NVTX_END_RANGE()
550 : : else
551 : 0 : write(std_out,'(a,i0)') 'entering chebfiwf2'
552 : : ABI_NVTX_START_RANGE(NVTX_CHEBFI2)
553 : : call chebfiwf2(cg_k,dtset,eig_k,occ_k,enlx_k,gs_hamk,&
554 [ # # ]: 0 : & mpi_enreg,nband_k,npw_k,my_nspinor,prtvol,resid_k)
555 : : ABI_NVTX_END_RANGE()
556 : : end if
557 : : else
558 [ # # ]: 0 : if (istep > nstep_mixed) then
559 : 0 : write(std_out,'(a,i0)') 'entering slicewf_cprj'
560 : : ! ITEST
561 : 0 : write(901,*)
562 : 0 : write(901,*) '**'
563 : 0 : write(901,*) 'SCF iteration=', istep
564 : 0 : write(901,*) '**'
565 : 0 : write(901,*)
566 : 0 : flush(901)
567 : : ! ITEST
568 : : !ABI_NVTX_START_RANGE(NVTX_SPESLI)
569 : : call slicewf_cprj(cg_k,dtset,eig_k,occ_k,enlx_k,gs_hamk,mpi_enreg,&
570 [ # # ]: 0 : & nband_k,npw_k,my_nspinor,prtvol,resid_k,xg_nonlop)
571 : : !ABI_NVTX_END_RANGE()
572 : : else
573 : 0 : write(std_out,'(a,i0)') 'entering chebfiwf2_cprj'
574 : : call chebfiwf2_cprj(cg_k,dtset,eig_k,occ_k,enlx_k,gs_hamk,&
575 [ # # ]: 0 : mpi_enreg,nband_k,npw_k,my_nspinor,prtvol,resid_k,xg_nonlop)
576 : : end if
577 : : end if
578 : :
579 : : ! =========================================================================
580 : : ! ======== MINIMIZATION OF BANDS: CONJUGATE GRADIENT (Teter et al.) =======
581 : : ! =========================================================================
582 : : else
583 : : ! use_subvnlx=0; if (gs_hamk%usepaw==0 .or. associated(gs_hamk%fockcommon)) use_subvnlx=1
584 : : ! use_subvnlx=0; if (gs_hamk%usepaw==0) use_subvnlx=1
585 : :
586 [ + + ]: 444103 : if (.not. use_rmm_diis) then
587 : :
588 [ + + + + : 443878 : if (isppol==1.and.ikpt==1.and.inonsc==1.and.istep==1) then
+ + + + ]
589 [ + + ]: 4991 : if (dtset%tolwfr_diago/=zero) then
590 : 1304 : write(msg, '(a,es16.6)' ) ' cgwf: tolwfr_diago=',dtset%tolwfr_diago
591 : 1304 : call wrtout(std_out,msg,'COLL')
592 : : end if
593 : : end if
594 : :
595 [ - + ]: 443878 : if (dtset%cprj_in_memory==2) then
596 : : call cgwf_cprj(cg,cprj_cwavef_bands,dtset%cprj_update_lvl,eig_k,&
597 : : gs_hamk,icg,mcg,mpi_enreg,nband_k,dtset%nline,&
598 : 0 : dtset%ortalg,prtvol,quit,resid_k,subham,dtset%tolrde,dtset%tolwfr_diago,wfoptalg)
599 : : else
600 : : call cgwf(dtset%berryopt,cg,cgq,dtset%chkexit,cpus,dphase_k,dtefield,dtfil%filnam_ds(1),&
601 : : gsc,gs_hamk,icg,igsc,ikpt,inonsc,isppol,dtset%mband,mcg,mcgq,mgsc,mkgq,&
602 : : mpi_enreg,mpw,nband_k,dtset%nbdblock,nkpt,dtset%nline,npw_k,npwarr,my_nspinor,&
603 : : dtset%nsppol,dtset%ortalg,prtvol,pwind,pwind_alloc,pwnsfac,pwnsfacq,quit,resid_k,&
604 : 443878 : subham,subovl,subvnlx,dtset%tolrde,dtset%tolwfr_diago,use_subovl,use_subvnlx,wfoptalg,zshift)
605 : : end if
606 : : else
607 : : call rmm_diis(istep, ikpt, isppol, cg(:,icg+1:), dtset, eig_k, occ_k, enlx_k, gs_hamk, kinpw, gsc, &
608 : 225 : mpi_enreg, nband_k, npw_k, my_nspinor, resid_k, rmm_diis_status)
609 : : end if
610 : :
611 [ - + ]: 444103 : if (dtset%gpu_option==ABI_GPU_KOKKOS) then
612 : : ! Kokkos GPU branch is not OpenMp thread-safe, restoring OpenMP threads num
613 : 0 : call xomp_set_num_threads(dtset%gpu_kokkos_nthrd)
614 : : end if
615 : : end if
616 : : end if
617 : :
618 : : ! =========================================================================
619 : : ! ===================== FIND LARGEST RESIDUAL =============================
620 : : ! =========================================================================
621 : :
622 : : ! Find largest resid over bands at this k point
623 : : ! Note that this operation is done BEFORE rotation of bands:
624 : : ! it would be time-consuming to recompute the residuals after.
625 [ + + ]: 501686 : if (nbdbuf >= 0) then
626 [ - + + - : 4557084 : max_resid = maxval(resid_k(1:max(1,nband_k-nbdbuf)))
+ + ]
627 [ + - ]: 1948 : else if (nbdbuf==-101) then
628 [ + - - + : 22712 : max_resid = maxval(occ_k(1:nband_k)*resid_k(1:nband_k))
- - + + ]
629 : : else
630 : 0 : ABI_ERROR(sjoin('Bad value of nbdbuf:', itoa(nbdbuf)))
631 : : end if
632 : :
633 : : ! Print residuals
634 [ + + + + : 501686 : if(prtvol/=5.and.(prtvol>2 .or. ikpt<=nkpt_max))then
+ - ]
635 [ + + ]: 95908 : do ii=0,(nband_k-1)/8
636 : 69494 : write(msg,'(a,8es10.2)')' res:',(resid_k(iband),iband=1+ii*8,min(nband_k,8+ii*8))
637 : 571180 : call wrtout(std_out,msg,'PERS')
638 : : end do
639 : : end if
640 : :
641 : : ! =========================================================================
642 : : ! ========== DIAGONALIZATION OF HAMILTONIAN IN WFs SUBSPACE ===============
643 : : ! =========================================================================
644 : 501686 : do_subdiago = .not. wfopta10 == 1 .and. .not. xg_diago
645 [ + + ]: 501686 : if (use_rmm_diis) do_subdiago = .False. ! subdiago is already performed before RMM-DIIS.
646 : :
647 : : ABI_NVTX_START_RANGE(NVTX_SUB_SPC_DIAGO)
648 [ + + ]: 501033 : if (do_subdiago) then
649 [ + + ]: 446132 : if (prtvol > 1) call wrtout(std_out, " Performing subspace diagonalization.")
650 : 446132 : call timab(585,1,tsec) !"vtowfk(subdiago)"
651 [ - + ]: 446132 : if (dtset%cprj_in_memory==2) then
652 : : call subdiago_low_memory(cg,eig_k,evec,icg,istwf_k,&
653 : 0 : mcg,nband_k,npw_k,my_nspinor,dtset%paral_kgb,subham)
654 : 0 : call timab(585,2,tsec)
655 : 0 : call timab(578,1,tsec)
656 [ # # ]: 0 : call cprj_rotate(cprj_cwavef_bands,evec,gs_hamk%dimcprj,natom,nband_k,gs_hamk%nspinor)
657 : 0 : call timab(578,2,tsec)
658 : : else
659 : : call subdiago(cg, eig_k, evec, gsc, icg, igsc, istwf_k, &
660 : : mcg, mgsc, nband_k, npw_k, my_nspinor, dtset%paral_kgb, &
661 : 446132 : subham, subovl, use_subovl, gs_hamk%usepaw, mpi_enreg%me_g0)
662 : 446132 : call timab(585,2,tsec)
663 : : end if
664 : : end if
665 : : ABI_NVTX_END_RANGE()
666 : :
667 : : ! Print energies
668 [ + + + + : 501686 : if(prtvol/=5.and.(prtvol>2 .or. ikpt<=nkpt_max))then
+ - ]
669 [ + + ]: 95908 : do ii=0,(nband_k-1)/8
670 : 69494 : write(msg, '(a,8es10.2)' )' ene:',(eig_k(iband),iband=1+ii*8,min(nband_k,8+ii*8))
671 : 571180 : call wrtout(std_out,msg,'PERS')
672 : : end do
673 : : end if
674 : :
675 : : ! THIS CHANGE OF SHIFT DOES NOT WORK WELL
676 : : ! Update zshift in the case of wfoptalg==3
677 : : ! if(wfoptalg==3 .and. inonsc/=1)then
678 : : ! do iband=1,nband_k
679 : : ! if(eig_k(iband)<eshift .and. eig_save(iband)<eshift) zshift(iband)=max(eig_k(iband),eig_save(iband))
680 : : ! if(eig_k(iband)>eshift .and. eig_save(iband)>eshift) zshift(iband)=min(eig_k(iband),eig_save(iband))
681 : : ! end do
682 : : ! eig_save(:)=eig_k(:)
683 : : ! end if
684 : :
685 : : ! =========================================================================
686 : : ! =============== ORTHOGONALIZATION OF WFs (if needed) ====================
687 : : ! =========================================================================
688 : :
689 : : ! Re-orthonormalize the wavefunctions at this k point.
690 : : ! this step is redundant but is performed to combat rounding error in wavefunction orthogonality.
691 : : ! This step is performed inside rmm_diis if RMM-DIIS is activated.
692 : :
693 : 501686 : call timab(583,1,tsec) ! "vtowfk(pw_orthon)"
694 : 501686 : ortalgo = mpi_enreg%paral_kgb
695 : : ! The orthogonalization is completely disabled with ortalg<=-10.
696 : : ! This option is useful for testing only and is not documented.
697 [ + + + + : 501686 : do_ortho = (wfoptalg/=14 .and. wfoptalg /= 1 .and. wfoptalg /= 12 .and. wfoptalg /= 11 .and. dtset%ortalg>-10) .or. dtset%ortalg > 0
+ + - + -
+ ]
698 [ + + ]: 501686 : if (xg_diago) do_ortho = .false.
699 [ + + ]: 501686 : if (use_rmm_diis) do_ortho = .False.
700 : :
701 [ + + ]: 501033 : if (do_ortho) then
702 : :
703 : : ABI_NVTX_START_RANGE(NVTX_ORTHO_WF)
704 : :
705 [ + + ]: 443996 : if (prtvol > 0) call wrtout(std_out, " Calling pw_orthon to orthonormalize bands.")
706 [ - + ]: 443996 : if (dtset%cprj_in_memory==2) then
707 [ # # ]: 0 : ABI_FREE(subovl)
708 [ # # # # : 0 : ABI_MALLOC(subovl,(nband_k*(nband_k+1)))
# # # # #
# # # ]
709 : 0 : call mksubovl(cg,cprj_cwavef_bands,gs_hamk,icg,nband_k,subovl,mpi_enreg)
710 : 0 : call pw_orthon_cprj(icg,mcg,npw_k*my_nspinor,my_nspinor,nband_k,ortalgo,subovl,cg,cprj=cprj_cwavef_bands)
711 : : else
712 : : call pw_orthon(icg,igsc,istwf_k,mcg,mgsc,npw_k*my_nspinor,nband_k,ortalgo,gsc,gs_hamk%usepaw,cg,&
713 : 443996 : mpi_enreg%me_g0,mpi_enreg%comm_bandspinorfft)
714 : : end if
715 : :
716 : : ABI_NVTX_END_RANGE()
717 : : end if
718 : 501686 : call timab(583,2,tsec)
719 : :
720 : : ABI_NVTX_START_RANGE(NVTX_VTOWFK_EXTRA2)
721 : :
722 : : ! DEBUG seq==par comment next block
723 : : ! Fix phases of all bands
724 [ + + ]: 501686 : if (xmpi_paral/=1 .or. mpi_enreg%paral_kgb/=1) then
725 : : !call wrtout(std_out, "Calling cgtk_fixphase")
726 [ + + + - ]: 457746 : if ( (.not.xg_diago) .and. dtset%cprj_in_memory==0 ) then
727 : 444725 : call cgtk_fixphase(cg,gsc,icg,igsc,istwf_k,mcg,mgsc,mpi_enreg,nband_k,npw_k*my_nspinor,gs_hamk%usepaw)
728 : : else if ( xg_diago ) then
729 : : ! GSC is local to vtowfk and is completely useless since everything
730 : : ! is calculated in my lobpcg, we don't care about the phase of gsc !
731 : 13021 : call cgtk_fixphase(cg,gsc,icg,igsc,istwf_k,mcg,mgsc,mpi_enreg,nband_k,npw_k*my_nspinor,0)
732 : : else ! dtset%cprj_in_memory/=0 .and. .not.xg_diago
733 : : call cgtk_fixphase(cg,gsc,icg,igsc,istwf_k,mcg,mgsc,mpi_enreg,nband_k,npw_k*my_nspinor,0,&
734 : 0 : & cprj=cprj_cwavef_bands,nspinor=dtset%nspinor)
735 : : end if
736 : : end if
737 : :
738 [ + + ]: 501686 : if (iscf < 0) then
739 [ + + + + ]: 71566 : if (max_resid > dtset%tolwfr .and. max_resid < tol7) then
740 [ - + ]: 32000 : if (fftcore_mixprec == 1) call wrtout(std_out, " Approaching NSCF convergence. Activating FFT in double-precision")
741 : 32000 : ii = fftcore_set_mixprec(0)
742 : : end if
743 : :
744 : : ! Print residual and wall-time required by NSCF iteration.
745 [ + + + + ]: 71566 : if (inonsc <= enough .or. mod(inonsc, 20) == 0) then
746 : 36862 : call cwtime(cpu, wall, gflops, "stop")
747 [ + + + + ]: 36862 : if (ikpt == 1 .or. mod(ikpt, 100) == 0) then
748 [ + + ]: 1283 : if (inonsc == 1) call wrtout(std_out, sjoin(" k-point: [", itoa(ikpt), "/", itoa(nkpt), "], spin:", itoa(isppol)))
749 : : call wrtout(std_out, sjoin(" Max resid =", ftoa(max_resid, fmt="es13.5"), &
750 : : " (exclude nbdbuf bands). One NSCF iteration cpu-time:", &
751 : 1283 : sec2str(cpu), ", wall-time:", sec2str(wall)), do_flush=.True.)
752 [ + + ]: 1283 : if (inonsc == enough) call wrtout(std_out, " Printing residuals every mod(20) iterations...")
753 : : end if
754 : : end if
755 : : end if
756 : : ABI_NVTX_END_RANGE()
757 : :
758 : : ! Exit loop over inonsc if converged
759 [ + + ]: 825734 : if (max_resid < dtset%tolwfr) then
760 [ + + + + : 22453 : if (iscf < 0 .and. (ikpt == 1 .or. mod(ikpt, 100) == 0)) then
+ + ]
761 : 519 : call wrtout(std_out, sjoin(" NSCF loop completed after", itoa(inonsc), "iterations"))
762 : : end if
763 : : exit
764 : : end if
765 : : end do ! inonsc (NON SELF-CONSISTENT LOOP)
766 : :
767 [ - + ]: 346501 : if (dtset%cprj_in_memory==2) then
768 : 0 : update_cprj=dtset%cprj_update_lvl<=3.and.dtset%cprj_update_lvl/=2
769 [ # # ]: 0 : if (update_cprj) call cprj_update(cg,cprj_cwavef_bands,gs_hamk,icg,nband_k,mpi_enreg,tim_getcprj)
770 : : end if
771 : :
772 : 346501 : call timab(39,2,tsec)
773 : 346501 : call timab(30,1,tsec) ! "vtowfk (afterloop)"
774 : :
775 : : !if (dtset%prtvol > 0)
776 : : !call wrtout(std_out, sjoin(" Number of Vnl|Psi> applications:", itoa(nonlop_counter)))
777 : :
778 : : !###################################################################
779 : :
780 : : !Compute kinetic energy and non-local energy for each band, and in the SCF
781 : : !case, contribution to forces, and eventually accumulate rhoaug
782 : :
783 [ + + ]: 346501 : ndat=1;if (mpi_enreg%paral_kgb==1) ndat=mpi_enreg%bandpp
784 [ + + + + ]: 346501 : if(iscf>0 .and. fixed_occ) then
785 [ + + ]: 232077 : ndat_fft=ndat; if(mpi_enreg%paral_kgb==0) ndat_fft=blocksize
786 : 232077 : chunk = ndat_fft/gs_hamk%nfft_blocks ! Divide by 2 to construct chunk of even number of bands
787 : 232077 : residuchunk = ndat_fft - gs_hamk%nfft_blocks*chunk
788 [ + - ]: 232077 : if(dtset%gpu_option==ABI_GPU_KOKKOS) then
789 : : #if defined HAVE_GPU && defined HAVE_YAKL
790 : : ABI_MALLOC_MANAGED(wfraug,(/2,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6*ndat_fft/))
791 : : #endif
792 : : else
793 : 232077 : fftalg = gs_hamk%ngfft(7)
794 [ + - + + ]: 232077 : if (gs_hamk%gpu_option==ABI_GPU_DISABLED.and.fftalg/=401) then
795 [ - + + - : 1137730 : ABI_MALLOC(wfraug,(2,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6))
- + + - -
+ + - + -
+ - - + -
+ - + ]
796 : : else
797 [ - + + - : 22655 : ABI_MALLOC(wfraug,(2,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6*(chunk+residuchunk)))
- + + - -
+ + - + -
+ - - + -
+ - + ]
798 : : end if
799 : : end if
800 : : end if
801 : :
802 : : !"nonlop" routine input parameters
803 : 346501 : nnlout=3*natom*optforces
804 : 346501 : signs=1;idir=0
805 [ + + ]: 346501 : if (gs_hamk%usepaw==0) then
806 : 236360 : choice=1+optforces
807 : 236360 : paw_opt=0;cpopt=-1;tim_nonlop=2
808 : : else
809 : 110141 : choice=2*optforces
810 : 110141 : paw_opt=2;cpopt=0;tim_nonlop=10-8*optforces
811 [ - + ]: 110141 : if (dtset%cprj_in_memory==2) cpopt=2 ! cprj are in memory (but not the derivatives)
812 [ + + ]: 110141 : if (dtset%usefock==1) then
813 : : !if (dtset%optforces/= 0) then
814 [ + + ]: 602 : if (optforces/= 0) then
815 : 217 : choice=2;cpopt=1; nnlout=3*natom
816 : : end if
817 : : end if
818 : : end if
819 : :
820 [ + + + - : 1039503 : ABI_MALLOC(enlout,(nnlout*blocksize))
+ - + + -
+ - + -
+ ]
821 : :
822 : : ! Allocation of memory space for one block of waveforms containing blocksize waveforms
823 [ + - ]: 346501 : if(dtset%gpu_option==ABI_GPU_KOKKOS) then
824 : : #if defined HAVE_GPU && defined HAVE_YAKL
825 : : ABI_MALLOC_MANAGED(cwavef, (/2,npw_k*my_nspinor*blocksize/))
826 : : #endif
827 : : else
828 [ - + + - : 1039503 : ABI_MALLOC(cwavef, (2,npw_k*my_nspinor*blocksize))
+ - + - -
+ - + -
+ ]
829 : : end if
830 : :
831 [ + - ]: 346501 : if (dtset%cprj_in_memory/=2) then
832 [ + + + + : 346501 : if (gs_hamk%usepaw==1.and.(iscf>0.or.gs_hamk%usecprj==1)) then
+ - ]
833 : 110141 : iorder_cprj=0
834 : 110141 : nband_k_cprj=nband_k*(mband_cprj/dtset%mband)
835 : 110141 : bandpp_cprj=mpi_enreg%bandpp
836 [ - + + - : 1409463 : ABI_MALLOC(cwaveprj,(natom,my_nspinor*bandpp_cprj))
- + + - +
- + - - +
- + + + +
+ ]
837 [ + + ]: 110141 : ncpgr=0;if (cpopt==1) ncpgr=cprj(1,1)%ncpgr
838 : 110141 : call pawcprj_alloc(cwaveprj,ncpgr,gs_hamk%dimcprj)
839 : : else
840 [ - + ]: 236360 : ABI_MALLOC(cwaveprj,(0,0))
841 : : end if
842 : : end if
843 : :
844 : : !The code below is more efficient if paral_kgb==1 (less MPI communications)
845 : : !however OMP is not compatible with paral_kgb since we should define
846 : : !which threads performs the call to MPI_ALL_REDUCE.
847 : : !This problem can be easily solved by removing MPI_enreg from meanvalue_g so that
848 : : !the MPI call is done only once outside the OMP parallel region.
849 : :
850 : : !call cwtime(cpu, wall, gflops, "start")
851 : :
852 [ + + ]: 346501 : if (dtset%cprj_in_memory==1) then
853 : 17084 : ncols_cprj = blocksize*my_nspinor/mpi_enreg%nproc_band
854 : 17084 : call xg_init(cprj_xgx0,xg_nonlop%space_cprj,xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
855 : 17084 : call xg_init(cprj_work,xg_nonlop%space_cprj,xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
856 : : end if
857 : :
858 : : ! In case of GEMM nonlop distribution + force computation,
859 : : ! recompute distribution as projectors arrays are bigger in this case
860 : 346501 : gpu_option_tmp=gs_hamk%gpu_option
861 [ + + - + ]: 346501 : if(optforces==1 .and. gs_hamk%gpu_option==ABI_GPU_OPENMP) then
862 : 0 : blksize_gemm_nonlop_tmp = gemm_nonlop_block_size; is_distrib_tmp = gemm_nonlop_is_distributed
863 : 0 : nfft_blocks_tmp = gs_hamk%nfft_blocks
864 : 0 : gemm_nonlop_block_size = dtset%gpu_nl_splitsize
865 : : call get_gemm_nonlop_ompgpu_blocksize(ikpt,gs_hamk,mpi_enreg%bandpp,nband_k,&
866 : : & dtset%nspinor,dtset%nspden,mpi_enreg%paral_kgb,mpi_enreg%nproc_band,&
867 : : & optforces,0,-1,gs_hamk%gpu_option,(dtset%gpu_nl_distrib/=0),&
868 : : & gemm_nonlop_block_size,nblk_gemm_nonlop,gs_hamk%nfft_blocks,&
869 [ # # # # ]: 0 : & warn_on_fail=.true.,disable_output=(istep>1.and.ikpt>1))
870 [ # # # # ]: 0 : gemm_nonlop_is_distributed = (dtset%gpu_nl_distrib/=0 .and. nblk_gemm_nonlop > 0)
871 [ # # ]: 0 : if(nblk_gemm_nonlop==-1) then
872 : 0 : gs_hamk%gpu_option=ABI_GPU_DISABLED
873 : 0 : ABI_WARNING("GPU has been disabled for forces computation during SCF step due to memory constraints.")
874 : : end if
875 : : end if
876 : :
877 : : #ifdef HAVE_OPENMP_OFFLOAD
878 : : !$OMP TARGET ENTER DATA MAP(alloc:cwavef) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
879 : : !$OMP TARGET ENTER DATA MAP(to:kinpw) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
880 : : #endif
881 : :
882 : : ! Transferring cg back is only needed in case of DMFT case or if GBT is on
883 : : if (gs_hamk%use_gbt /= 0 .or. paw_dmft%use_dmft==1) then
884 : : #ifdef HAVE_OPENMP_OFFLOAD
885 : : !$OMP TARGET UPDATE FROM(cg_k) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP .and. xg_diago)
886 : : #endif
887 : : end if
888 : : ! Loop over bands or blocks of bands.
889 : : ! Note that in sequential mode iblock=iband, nblockbd=nband_k and blocksize=1
890 [ + + ]: 2523815 : do iblock=1,nblockbd
891 [ + - - + : 7009858 : occblock=maxval(occ_k(1+(iblock-1)*blocksize:iblock*blocksize))
- - + + ]
892 : : !cwavef(:,:)=cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg)
893 [ + - ]: 2177314 : if(gs_hamk%gpu_option==ABI_GPU_OPENMP) then
894 : : #ifdef HAVE_OPENMP_OFFLOAD
895 : : ! cg is already on GPU, simply copy it in cwavef
896 : : if(xg_diago) then
897 : : call gpu_copy(cwavef, &
898 : : & cg_k(:,1+(iblock-1)*npw_k*my_nspinor*blocksize:iblock*npw_k*my_nspinor*blocksize),&
899 : : & int(2,c_size_t)*npw_k*my_nspinor*blocksize)
900 : : else
901 : : ! cg isn't on GPU, single transfer and copy it in cwavef
902 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cg_k,cwavef) PRIVATE(iblocksize)
903 : : do iblocksize=1,blocksize*my_nspinor
904 : : !$OMP PARALLEL DO PRIVATE(ipw)
905 : : do ipw=1,npw_k
906 : : cwavef(1,ipw+(iblocksize-1)*npw_k)=cg_k(1,ipw+(iblocksize-1)*npw_k+(iblock-1)*npw_k*blocksize*my_nspinor)
907 : : cwavef(2,ipw+(iblocksize-1)*npw_k)=cg_k(2,ipw+(iblocksize-1)*npw_k+(iblock-1)*npw_k*blocksize*my_nspinor)
908 : : end do
909 : : end do
910 : : end if
911 : : #endif
912 : : else
913 : : call DCOPY(2*npw_k*my_nspinor*blocksize, &
914 [ - + ]: 2177314 : cg_k(:,1+(iblock-1)*npw_k*my_nspinor*blocksize:iblock*npw_k*my_nspinor*blocksize), 1, cwavef, 1)
915 : : end if
916 : :
917 : : ! Compute kinetic energies for all bands in this block (use_gbt==0).
918 : : ! meanvalue_g_batch handles both istwf_k==1 (phase 1) and istwf_k>=2 (phase 3).
919 [ + + ]: 2177314 : if (gs_hamk%use_gbt == 0) then
920 : : call meanvalue_g_batch(ek_k(1+(iblock-1)*blocksize:iblock*blocksize), kinpw, &
921 : : & 0, istwf_k, mpi_enreg, npw_k, my_nspinor, blocksize, &
922 : 2139746 : & cwavef, cwavef, 0, gpu_option=gs_hamk%gpu_option)
923 : :
924 : : ! Compute nuclear dipole energies for all bands in this block (use_gbt==0).
925 : : ! NOTE: this call is problematic in the npband /= 1 case. It does not crash
926 : : ! but doesn't return the correct band energies. However, the resulting end_k is
927 : : ! used only for printing.
928 : : ! To be fixed eventually--jzwanzig 20 August 2026
929 [ + + + + : 23423322 : if(ANY(ABS(dtset%nucdipmom)>tol8)) then
+ + + + ]
930 [ - + + - : 195600 : ABI_MALLOC(ghc_vectornd,(2,npw_k*my_nspinor*blocksize))
+ - + - -
+ - + ]
931 : : call getghc_nucdip(cwavef,ghc_vectornd,gs_hamk%gbound_k,gs_hamk%istwf_k,kg_k,gs_hamk%kpt_k,&
932 : : & gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%ngfft,npw_k,gs_hamk%nvloc,&
933 : : & gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,my_nspinor,gs_hamk%vectornd,gs_hamk%vlocal,&
934 [ + - - - : 65200 : & gs_hamk%zora,gs_hamk%gpu_option)
- - - - -
- - - - -
+ - - - -
- - - - -
- - - + -
- - - - -
- - - - -
+ - - - -
- - - - ]
935 [ + + ]: 130700 : do iblocksize=1,blocksize
936 : 65500 : iband=(iblock-1)*blocksize+iblocksize
937 : : end_k(iband)=DOT_PRODUCT(cwavef(1,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
938 : : & ghc_vectornd(1,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor))+ &
939 : : & DOT_PRODUCT(cwavef(2,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
940 [ + + + + ]: 96924812 : & ghc_vectornd(2,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor))
941 : : end do
942 [ - + ]: 65200 : ABI_FREE(ghc_vectornd)
943 : : end if
944 : :
945 : : end if
946 : :
947 [ + + ]: 4832544 : do iblocksize=1,blocksize
948 : 2655230 : iband=(iblock-1)*blocksize+iblocksize
949 : :
950 : 2655230 : cwavef_iband => cg(:,1+(iband-1)*npw_k*my_nspinor+icg:iband*npw_k*my_nspinor+icg)
951 : :
952 [ + + ]: 2655230 : if (gs_hamk%use_gbt /= 0) then
953 : : ! Treat up and down components separately. cwavef_iband is still in
954 : : ! the local distributed G layout here, so use the local k/k' kinetic
955 : : ! tables passed by vtorho rather than the transposed tables in gs_hamk.
956 : 37568 : filter = 1
957 : : call meanvalue_g(ar,kinpw,filter,istwf_k,mpi_enreg,npw_k,1,&
958 : : & cwavef_iband, cwavef_iband, 0,&
959 : 37568 : & gpu_thread_limit=dtset%gpu_thread_limit)
960 : :
961 : : call meanvalue_g(ar2,kinpw_kp,filter,istwf_k,mpi_enreg,npw_k,1,&
962 : : & cwavef_iband(:,npw_k+1:), cwavef_iband(:,npw_k+1:),0,&
963 [ - + - - : 37568 : & gpu_thread_limit=dtset%gpu_thread_limit)
- - - + ]
964 : :
965 : 37568 : ek_k(iband) = ar + ar2
966 : : end if
967 : :
968 [ + + ]: 4832544 : if(paw_dmft%use_dmft==1) then
969 [ + + ]: 529272 : do iband1=1,nband_k
970 : : call meanvalue_g(ar,kinpw,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
971 : : cg(:,1+(iband -1)*npw_k*my_nspinor+icg:iband *npw_k*my_nspinor+icg),&
972 : : cg(:,1+(iband1-1)*npw_k*my_nspinor+icg:iband1*npw_k*my_nspinor+icg),&
973 : 510096 : paw_dmft%use_dmft,ar_im=ar_im,gpu_thread_limit=dtset%gpu_thread_limit)
974 : 510096 : ek_k_nd(1,iband,iband1)=ar
975 : 529272 : ek_k_nd(2,iband,iband1)=ar_im
976 : : end do
977 : : end if
978 : : end do
979 : :
980 [ + + ]: 2177314 : if (iscf>0) then
981 : :
982 : : ABI_NVTX_START_RANGE(NVTX_VTOWFK_FOURWF)
983 : : ! In case of fixed occupation numbers, accumulates the partial density
984 [ + + + + ]: 2052435 : if (fixed_occ .and. mpi_enreg%paral_kgb/=1) then
985 : :
986 : : ! treat all bands at once on GPU
987 [ - + ]: 873078 : if (dtset%gpu_option /= ABI_GPU_DISABLED) then
988 : :
989 [ # # # # : 0 : ABI_MALLOC(weight_t,(blocksize))
# # # # #
# # # ]
990 : : ! compute weights
991 [ # # ]: 0 : do iblocksize=1,blocksize
992 : 0 : iband=(iblock-1)*blocksize+iblocksize
993 : 0 : weight_t(iblocksize) = occ_k(iband) * wtk / gs_hamk%ucvol
994 [ # # ]: 0 : if (abs(occ_k(iband)) < tol8) weight_t(iblocksize) = zero
995 : : end do
996 : :
997 [ # # ]: 0 : if(dtset%nspinor==1) then
998 : : call fourwf_optmem(1,rhoaug(:,:,:,1),cwavef(:,:),dummy,wfraug,&
999 : : gs_hamk%gbound_k,gs_hamk%gbound_k,istwf_k,kg_k,kg_k,&
1000 : : gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%nfft_blocks,gs_hamk%ngfft,&
1001 : : npw_k,1,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,tim_fourwf,weight,weight,&
1002 : : weight_array_r=weight_t,weight_array_i=weight_t,&
1003 : 0 : gpu_option=dtset%gpu_option)
1004 : :
1005 [ # # ]: 0 : else if(dtset%nspinor==2) then
1006 [ # # # # : 0 : ABI_MALLOC(cwavefb,(2,npw_k*blocksize,2))
# # # # #
# # # # #
# # ]
1007 : 0 : ibs=(iblock-1)*npw_k*my_nspinor*blocksize
1008 [ # # ]: 0 : if(gs_hamk%gpu_option==ABI_GPU_OPENMP) then
1009 : : #ifdef HAVE_OPENMP_OFFLOAD
1010 : : !$OMP TARGET ENTER DATA MAP(alloc:cwavefb)
1011 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cwavefb,cwavef) PRIVATE(iband)
1012 : : do iband=1,blocksize
1013 : : !$OMP PARALLEL DO PRIVATE(ipw)
1014 : : do ipw=1,npw_k
1015 : : cwavefb(1,(iband-1)*npw_k+ipw,1)=cwavef(1,ipw+(2*iband-2)*npw_k)
1016 : : cwavefb(2,(iband-1)*npw_k+ipw,1)=cwavef(2,ipw+(2*iband-2)*npw_k)
1017 : :
1018 : : cwavefb(1,(iband-1)*npw_k+ipw,2)=cwavef(1,ipw+(2*iband-1)*npw_k)
1019 : : cwavefb(2,(iband-1)*npw_k+ipw,2)=cwavef(2,ipw+(2*iband-1)*npw_k)
1020 : : end do
1021 : : end do
1022 : : #endif
1023 : : else
1024 [ # # ]: 0 : do iband=1,blocksize
1025 [ # # # # ]: 0 : cwavefb(:,(iband-1)*npw_k+1:iband*npw_k,1)=cwavef(:,1+(2*iband-2)*npw_k:(iband*2-1)*npw_k)
1026 [ # # # # ]: 0 : cwavefb(:,(iband-1)*npw_k+1:iband*npw_k,2)=cwavef(:,1+(2*iband-1)*npw_k:iband*2*npw_k)
1027 : : end do
1028 : : end if
1029 : :
1030 : : call fourwf_optmem(1,rhoaug(:,:,:,1),cwavefb(:,:,1),dummy,wfraug,&
1031 : : gs_hamk%gbound_k,gs_hamk%gbound_k,istwf_k,kg_k,kg_k,&
1032 : : gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%nfft_blocks,gs_hamk%ngfft,&
1033 : : npw_k,1,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,tim_fourwf,weight,weight,&
1034 : : weight_array_r=weight_t,weight_array_i=weight_t,&
1035 : 0 : gpu_option=dtset%gpu_option)
1036 : :
1037 [ # # ]: 0 : if(dtset%nspden==1) then
1038 : : call fourwf_optmem(1,rhoaug(:,:,:,1),cwavefb(:,:,2),dummy,wfraug,&
1039 : : gs_hamk%gbound_k,gs_hamk%gbound_k,istwf_k,kg_k,kg_k,&
1040 : : gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%nfft_blocks,gs_hamk%ngfft,&
1041 : : npw_k,1,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,tim_fourwf,weight,weight,&
1042 : : weight_array_r=weight_t,weight_array_i=weight_t,&
1043 : 0 : gpu_option=dtset%gpu_option)
1044 : :
1045 [ # # ]: 0 : else if (dtset%nspden==4) then
1046 : : ! Build the four components of rho. We use only norm quantities and, so fourwf.
1047 : : ! $\sum_{n} f_n \Psi^{* \alpha}_n \Psi^{\alpha}_n =\rho^{\alpha \alpha}$
1048 : : ! $\sum_{n} f_n (\Psi^{1}+\Psi^{2})^*_n (\Psi^{1}+\Psi^{2})_n=rho+m_x$
1049 : : ! $\sum_{n} f_n (\Psi^{1}-i \Psi^{2})^*_n (\Psi^{1}-i \Psi^{2})_n=rho+m_y$
1050 [ # # # # : 0 : ABI_MALLOC(cwavef_x,(2,npw_k*blocksize))
# # # # #
# # # ]
1051 [ # # # # : 0 : ABI_MALLOC(cwavef_y,(2,npw_k*blocksize))
# # # # ]
1052 : :
1053 [ # # ]: 0 : if(gs_hamk%gpu_option==ABI_GPU_OPENMP) then
1054 : : #ifdef HAVE_OPENMP_OFFLOAD
1055 : : !$OMP TARGET ENTER DATA MAP(alloc:cwavef_x,cwavef_y)
1056 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cwavef_x,cwavef_y,cwavefb) PRIVATE(iband)
1057 : : do iband=1,blocksize
1058 : : !$OMP PARALLEL DO PRIVATE(ipw)
1059 : : do ipw=1,npw_k
1060 : : !$(\Psi^{1}+\Psi^{2})$
1061 : : cwavef_x(1,(iband-1)*npw_k+ipw)=cwavefb(1,(iband-1)*npw_k+ipw,1)+cwavefb(1,(iband-1)*npw_k+ipw,2)
1062 : : cwavef_x(2,(iband-1)*npw_k+ipw)=cwavefb(2,(iband-1)*npw_k+ipw,1)+cwavefb(2,(iband-1)*npw_k+ipw,2)
1063 : : !$(\Psi^{1}-i \Psi^{2})$
1064 : : cwavef_y(1,(iband-1)*npw_k+ipw)=cwavefb(1,(iband-1)*npw_k+ipw,1)+cwavefb(2,(iband-1)*npw_k+ipw,2)
1065 : : cwavef_y(2,(iband-1)*npw_k+ipw)=cwavefb(2,(iband-1)*npw_k+ipw,1)-cwavefb(1,(iband-1)*npw_k+ipw,2)
1066 : : end do
1067 : : end do
1068 : : #endif
1069 : : else
1070 : : !$(\Psi^{1}+\Psi^{2})$
1071 [ # # # # ]: 0 : cwavef_x(:,:)=cwavefb(:,1:npw_k*blocksize,1)+cwavefb(:,1:npw_k*blocksize,2)
1072 : : !$(\Psi^{1}-i \Psi^{2})$
1073 [ # # ]: 0 : cwavef_y(1,:)=cwavefb(1,1:npw_k*blocksize,1)+cwavefb(2,1:npw_k*blocksize,2)
1074 [ # # ]: 0 : cwavef_y(2,:)=cwavefb(2,1:npw_k*blocksize,1)-cwavefb(1,1:npw_k*blocksize,2)
1075 : : end if
1076 : :
1077 : : ! z component
1078 : : call fourwf_optmem(1,rhoaug(:,:,:,4),cwavefb(:,:,2),dummy,wfraug,&
1079 : : gs_hamk%gbound_k,gs_hamk%gbound_k,istwf_k,kg_k,kg_k,&
1080 : : gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%nfft_blocks,gs_hamk%ngfft,&
1081 : : npw_k,1,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,tim_fourwf,weight,weight,&
1082 : : weight_array_r=weight_t,weight_array_i=weight_t,&
1083 : 0 : gpu_option=dtset%gpu_option)
1084 : : ! x component
1085 : : call fourwf_optmem(1,rhoaug(:,:,:,2),cwavef_x(:,:),dummy,wfraug,&
1086 : : gs_hamk%gbound_k,gs_hamk%gbound_k,istwf_k,kg_k,kg_k,&
1087 : : gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%nfft_blocks,gs_hamk%ngfft,&
1088 : : npw_k,1,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,tim_fourwf,weight,weight,&
1089 : : weight_array_r=weight_t,weight_array_i=weight_t,&
1090 : 0 : gpu_option=dtset%gpu_option)
1091 : : ! y component
1092 : : call fourwf_optmem(1,rhoaug(:,:,:,3),cwavef_y(:,:),dummy,wfraug,&
1093 : : gs_hamk%gbound_k,gs_hamk%gbound_k,istwf_k,kg_k,kg_k,&
1094 : : gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%nfft_blocks,gs_hamk%ngfft,&
1095 : : npw_k,1,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,tim_fourwf,weight,weight,&
1096 : : weight_array_r=weight_t,weight_array_i=weight_t,&
1097 : 0 : gpu_option=dtset%gpu_option)
1098 : : #ifdef HAVE_OPENMP_OFFLOAD
1099 : : !$OMP TARGET EXIT DATA MAP(delete:cwavef_x,cwavef_y) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
1100 : : #endif
1101 [ # # ]: 0 : ABI_FREE(cwavef_x)
1102 [ # # ]: 0 : ABI_FREE(cwavef_y)
1103 : : end if
1104 : :
1105 : : #ifdef HAVE_OPENMP_OFFLOAD
1106 : : !$OMP TARGET EXIT DATA MAP(delete:cwavefb) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
1107 : : #endif
1108 [ # # ]: 0 : ABI_FREE(cwavefb)
1109 : : end if
1110 : :
1111 [ # # ]: 0 : ABI_FREE(weight_t)
1112 : :
1113 : : else
1114 [ + + ]: 1835916 : do iblocksize=1,blocksize
1115 : 962838 : iband=(iblock-1)*blocksize+iblocksize
1116 : 962838 : cwavef_iband => cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor)
1117 : :
1118 [ + + ]: 1835916 : if (abs(occ_k(iband))>=tol8) then
1119 : 720608 : weight = occ_k(iband) * wtk / gs_hamk%ucvol
1120 : :
1121 : : ! Accumulate charge density in real space in array rhoaug
1122 : : ! The same section of code is also found in mkrho.F90 : should be rationalized !
1123 : : call fourwf(1,rhoaug(:,:,:,1),cwavef_iband,dummy,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
1124 : : istwf_k,gs_hamk%kg_k,gs_hamk%kg_k,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%ngfft,npw_k,1,&
1125 : : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,&
1126 [ + - - - : 720608 : tim_fourwf,weight,weight,gpu_option=dtset%gpu_option)
- - - - +
- - - - -
- - - + -
+ ]
1127 : :
1128 [ + + ]: 720608 : if(dtset%nspinor==2)then
1129 [ - + + - : 247746 : ABI_MALLOC(cwavef1,(2,npw_k))
+ - + - -
+ - + ]
1130 [ + + + + ]: 93425044 : cwavef1(:,:)=cwavef_iband(:,1+npw_k:2*npw_k) ! EB FR spin dn part and used for m_z component (cwavef_z)
1131 : :
1132 [ + + ]: 82582 : if(dtset%nspden==1) then
1133 : :
1134 : : call fourwf(1,rhoaug(:,:,:,1),cwavef1,dummy,wfraug,&
1135 : : gs_hamk%gbound_k,gs_hamk%gbound_k,&
1136 : : istwf_k,gs_hamk%kg_k,gs_hamk%kg_k,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%ngfft,npw_k,1,&
1137 : : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,&
1138 [ + - - - : 41252 : tim_fourwf,weight,weight,gpu_option=dtset%gpu_option)
- - - - +
- - - - -
- - - + -
+ ]
1139 : :
1140 [ + - ]: 41330 : else if(dtset%nspden==4) then
1141 : : ! Build the four components of rho. We use only norm quantities and, so fourwf.
1142 : : ! $\sum_{n} f_n \Psi^{* \alpha}_n \Psi^{\alpha}_n =\rho^{\alpha \alpha}$
1143 : : ! $\sum_{n} f_n (\Psi^{1}+\Psi^{2})^*_n (\Psi^{1}+\Psi^{2})_n=rho+m_x$
1144 : : ! $\sum_{n} f_n (\Psi^{1}-i \Psi^{2})^*_n (\Psi^{1}-i \Psi^{2})_n=rho+m_y$
1145 [ - + + - : 82660 : ABI_MALLOC(cwavef_x,(2,npw_k))
- + - + ]
1146 [ - + + - : 82660 : ABI_MALLOC(cwavef_y,(2,npw_k))
- + - + ]
1147 : : !$(\Psi^{1}+\Psi^{2})$
1148 [ + + + + ]: 26454140 : cwavef_x(:,:)=cwavef_iband(:,1:npw_k)+cwavef1(:,1:npw_k)
1149 : : !$(\Psi^{1}-i \Psi^{2})$
1150 [ + + ]: 8845600 : cwavef_y(1,:)=cwavef_iband(1,1:npw_k)+cwavef1(2,1:npw_k)
1151 [ + + ]: 8845600 : cwavef_y(2,:)=cwavef_iband(2,1:npw_k)-cwavef1(1,1:npw_k)
1152 : : ! z component
1153 : : call fourwf(1,rhoaug(:,:,:,4),cwavef1,dummy,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
1154 : : istwf_k,gs_hamk%kg_k,gs_hamk%kg_k,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%ngfft,npw_k,1,&
1155 : : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,&
1156 [ + - - - : 41330 : tim_fourwf,weight,weight,gpu_option=dtset%gpu_option)
- - - - +
- - - - -
- - - + -
+ ]
1157 : : ! x component
1158 : : call fourwf(1,rhoaug(:,:,:,2),cwavef_x,dummy,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
1159 : : istwf_k,gs_hamk%kg_k,gs_hamk%kg_k,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%ngfft,npw_k,1,&
1160 : : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,&
1161 [ + - - - : 41330 : tim_fourwf,weight,weight,gpu_option=dtset%gpu_option)
- - - - +
- - - - -
- - - + -
+ ]
1162 : : ! y component
1163 : : call fourwf(1,rhoaug(:,:,:,3),cwavef_y,dummy,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
1164 : : istwf_k,gs_hamk%kg_k,gs_hamk%kg_k,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%ngfft,npw_k,1,&
1165 : : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,&
1166 [ + - - - : 41330 : tim_fourwf,weight,weight,gpu_option=dtset%gpu_option)
- - - - +
- - - - -
- - - + -
+ ]
1167 : :
1168 [ - + ]: 41330 : ABI_FREE(cwavef_x)
1169 [ - + ]: 41330 : ABI_FREE(cwavef_y)
1170 : :
1171 : : end if ! dtset%nspden/=4
1172 [ - + ]: 82582 : ABI_FREE(cwavef1)
1173 : : end if
1174 : : else
1175 : 242230 : nskip=nskip+1
1176 : : end if
1177 : : end do ! Loop inside a block of bands
1178 : :
1179 : : end if ! dtset%gpu_option
1180 : :
1181 : :
1182 [ + + - + ]: 1179357 : else if (fixed_occ .and. mpi_enreg%paral_kgb==1) then
1183 : : ! In case of fixed occupation numbers,in bandFFT mode accumulates the partial density
1184 : :
1185 [ + + ]: 69090 : if (dtset%nspinor==1) then
1186 : 48370 : call timab(537,1,tsec) ! "prep_fourwf%vtow"
1187 : : call prep_fourwf(rhoaug(:,:,:,1),blocksize,cwavef,wfraug,iblock,istwf_k,&
1188 : : gs_hamk%mgfft,mpi_enreg,nband_k,ndat,gs_hamk%ngfft,npw_k,&
1189 : : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,occ_k,&
1190 : 48370 : 1,gs_hamk%ucvol,wtk,gs_hamk%nfft_blocks,gpu_option=dtset%gpu_option)
1191 : 48370 : call timab(537,2,tsec)
1192 [ - + ]: 20720 : else if (dtset%nspinor==2) then
1193 [ - + + - : 82880 : ABI_MALLOC(cwavefb,(2,npw_k*blocksize,2))
+ - + - +
- - + - +
- + ]
1194 [ - + ]: 20720 : if(gs_hamk%gpu_option==ABI_GPU_OPENMP) then
1195 : : ! --- No parallelization over spinors ---
1196 [ # # ]: 0 : if (mpi_enreg%paral_spinor==0) then
1197 : : #ifdef HAVE_OPENMP_OFFLOAD
1198 : : !$OMP TARGET ENTER DATA MAP(alloc:cwavefb)
1199 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cwavefb,cwavef) PRIVATE(iband)
1200 : : do iband=1,blocksize
1201 : : !$OMP PARALLEL DO PRIVATE(ipw)
1202 : : do ipw=1,npw_k
1203 : : cwavefb(1,(iband-1)*npw_k+ipw,1)=cwavef(1,ipw+(2*iband-2)*npw_k)
1204 : : cwavefb(2,(iband-1)*npw_k+ipw,1)=cwavef(2,ipw+(2*iband-2)*npw_k)
1205 : :
1206 : : cwavefb(1,(iband-1)*npw_k+ipw,2)=cwavef(1,ipw+(2*iband-1)*npw_k)
1207 : : cwavefb(2,(iband-1)*npw_k+ipw,2)=cwavef(2,ipw+(2*iband-1)*npw_k)
1208 : : end do
1209 : : end do
1210 : : #endif
1211 : : else
1212 : 0 : ABI_BUG("Parallelisation on spinor isn't supported with OpenMP GPU")
1213 : : end if
1214 : : else
1215 : : ! --- No parallelization over spinors ---
1216 [ + + ]: 20720 : if (mpi_enreg%paral_spinor==0) then
1217 [ + + ]: 95680 : do iband=1,blocksize
1218 [ + + + + ]: 9033380 : cwavefb(:,(iband-1)*npw_k+1:iband*npw_k,1)=cwavef(:,1+(2*iband-2)*npw_k:(iband*2-1)*npw_k)
1219 [ + + + + ]: 9045940 : cwavefb(:,(iband-1)*npw_k+1:iband*npw_k,2)=cwavef(:,1+(2*iband-1)*npw_k:iband*2*npw_k)
1220 : : end do
1221 : : else
1222 : : ! --- Parallelization over spinors ---
1223 : : ! (split the work between 2 procs)
1224 [ + + + + ]: 9855216 : cwavefb(:,:,3-ispinor_index)=zero
1225 [ + + ]: 24832 : do iband=1,blocksize
1226 [ + + + + ]: 9871888 : cwavefb(:,(iband-1)*npw_k+1:iband*npw_k,ispinor_index) = cwavef(:,1+(iband-1)*npw_k:iband*npw_k)
1227 : : end do
1228 : 8160 : call xmpi_sum(cwavefb,mpi_enreg%comm_spinor,ierr)
1229 : : end if
1230 : : end if !gpu_option
1231 : :
1232 : 20720 : call timab(537,1,tsec) !"prep_fourwf%vtow"
1233 [ + + ]: 20720 : if (nspinor1TreatedByThisProc) then
1234 : : call prep_fourwf(rhoaug(:,:,:,1),blocksize,cwavefb(:,:,1),wfraug,iblock,&
1235 : : istwf_k,gs_hamk%mgfft,mpi_enreg,nband_k,ndat,gs_hamk%ngfft,npw_k,&
1236 : : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,occ_k,1,gs_hamk%ucvol,wtk,gs_hamk%nfft_blocks,&
1237 : 16640 : gpu_option=dtset%gpu_option)
1238 : : end if
1239 [ + + ]: 20720 : if(dtset%nspden==1) then
1240 [ + + ]: 9376 : if (nspinor2TreatedByThisProc) then
1241 : : call prep_fourwf(rhoaug(:,:,:,1),blocksize,cwavefb(:,:,2),wfraug,&
1242 : : iblock,istwf_k,gs_hamk%mgfft,mpi_enreg,nband_k,ndat,&
1243 : : gs_hamk%ngfft,npw_k,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,occ_k,1,&
1244 : 8768 : gs_hamk%ucvol,wtk,gs_hamk%nfft_blocks,gpu_option=dtset%gpu_option)
1245 : : end if
1246 [ + - ]: 11344 : else if(dtset%nspden==4) then
1247 [ - + + - : 34032 : ABI_MALLOC(cwavef_x,(2,npw_k*blocksize))
+ - + - -
+ - + ]
1248 [ - + + - : 22688 : ABI_MALLOC(cwavef_y,(2,npw_k*blocksize))
- + - + ]
1249 [ + - ]: 11344 : if(gs_hamk%gpu_option==ABI_GPU_OPENMP) then
1250 : : #ifdef HAVE_OPENMP_OFFLOAD
1251 : : !$OMP TARGET ENTER DATA MAP(alloc:cwavef_x,cwavef_y)
1252 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cwavef_x,cwavef_y,cwavefb) PRIVATE(iband)
1253 : : do iband=1,blocksize
1254 : : !$OMP PARALLEL DO PRIVATE(ipw)
1255 : : do ipw=1,npw_k
1256 : : !$(\Psi^{1}+\Psi^{2})$
1257 : : cwavef_x(1,(iband-1)*npw_k+ipw)=cwavefb(1,(iband-1)*npw_k+ipw,1)+cwavefb(1,(iband-1)*npw_k+ipw,2)
1258 : : cwavef_x(2,(iband-1)*npw_k+ipw)=cwavefb(2,(iband-1)*npw_k+ipw,1)+cwavefb(2,(iband-1)*npw_k+ipw,2)
1259 : : !$(\Psi^{1}-i \Psi^{2})$
1260 : : cwavef_y(1,(iband-1)*npw_k+ipw)=cwavefb(1,(iband-1)*npw_k+ipw,1)+cwavefb(2,(iband-1)*npw_k+ipw,2)
1261 : : cwavef_y(2,(iband-1)*npw_k+ipw)=cwavefb(2,(iband-1)*npw_k+ipw,1)-cwavefb(1,(iband-1)*npw_k+ipw,2)
1262 : : end do
1263 : : end do
1264 : : #endif
1265 : : else
1266 [ + + + + ]: 5978296 : cwavef_x(:,:)=cwavefb(:,1:npw_k*blocksize,1)+cwavefb(:,:,2)
1267 [ + + ]: 2000328 : cwavef_y(1,:)=cwavefb(1,1:npw_k*blocksize,1)+cwavefb(2,:,2)
1268 [ + + ]: 2000328 : cwavef_y(2,:)=cwavefb(2,:,1)-cwavefb(1,:,2)
1269 : : end if
1270 [ + + ]: 11344 : if (nspinor1TreatedByThisProc) then
1271 : : call prep_fourwf(rhoaug(:,:,:,4),blocksize,cwavefb(:,:,2),wfraug,&
1272 : : iblock,istwf_k,gs_hamk%mgfft,mpi_enreg,nband_k,ndat,gs_hamk%ngfft,&
1273 : : npw_k,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,occ_k,1,gs_hamk%ucvol,wtk,gs_hamk%nfft_blocks,&
1274 : 7872 : gpu_option=dtset%gpu_option)
1275 : : end if
1276 [ + + ]: 11344 : if (nspinor2TreatedByThisProc) then
1277 : : call prep_fourwf(rhoaug(:,:,:,2),blocksize,cwavef_x,wfraug,&
1278 : : iblock,istwf_k,gs_hamk%mgfft,mpi_enreg,nband_k,ndat,gs_hamk%ngfft,&
1279 : : npw_k,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,occ_k,1,gs_hamk%ucvol,wtk,gs_hamk%nfft_blocks,&
1280 : 7872 : gpu_option=dtset%gpu_option)
1281 : : call prep_fourwf(rhoaug(:,:,:,3),blocksize,cwavef_y,wfraug,&
1282 : : iblock,istwf_k,gs_hamk%mgfft,mpi_enreg,nband_k,ndat,gs_hamk%ngfft,&
1283 : : npw_k,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,occ_k,1,gs_hamk%ucvol,wtk,gs_hamk%nfft_blocks,&
1284 : 7872 : gpu_option=dtset%gpu_option)
1285 : : end if
1286 : : #ifdef HAVE_OPENMP_OFFLOAD
1287 : : !$OMP TARGET EXIT DATA MAP(delete:cwavef_x,cwavef_y) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
1288 : : #endif
1289 [ - + ]: 11344 : ABI_FREE(cwavef_x)
1290 [ - + ]: 11344 : ABI_FREE(cwavef_y)
1291 : : end if
1292 : 20720 : call timab(537,2,tsec)
1293 : : #ifdef HAVE_OPENMP_OFFLOAD
1294 : : !$OMP TARGET EXIT DATA MAP(delete:cwavefb) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
1295 : : #endif
1296 [ - + ]: 20720 : ABI_FREE(cwavefb)
1297 : : end if
1298 : : end if
1299 : : ABI_NVTX_END_RANGE()
1300 : : end if ! End of SCF calculation
1301 : :
1302 : : ! Call to nonlocal operator:
1303 : : ! - Compute nonlocal forces from most recent wfs
1304 : : ! - PAW: compute projections of WF onto NL projectors (cprj)
1305 : : ABI_NVTX_START_RANGE(NVTX_VTOWFK_NONLOP)
1306 : 2177314 : eig_k_block => eig_k(1+(iblock-1)*blocksize:iblock*blocksize)
1307 : 2177314 : cg_k_block => cg_k(:,1+(iblock-1)*blocksize*my_nspinor*npw_k:iblock*blocksize*my_nspinor*npw_k)
1308 [ - + ]: 2523815 : if (dtset%cprj_in_memory==2) then
1309 [ # # ]: 0 : if (optforces>0) then
1310 : 0 : call timab(554,1,tsec) ! "vtowfk:rhoij"
1311 : : ! Treat all wavefunctions in case of PAW
1312 : 0 : cwaveprj => cprj(:,1+(iblock-1)*my_nspinor*blocksize+ibg:iblock*my_nspinor*blocksize+ibg)
1313 : : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,eig_k_block,&
1314 [ # # ]: 0 : mpi_enreg,blocksize,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,nonlop_dum2)
1315 : : ! Accumulate forces
1316 : 0 : iband=(iblock-1)*blocksize
1317 [ # # ]: 0 : do iblocksize=1,blocksize
1318 : 0 : ii=0
1319 [ # # ]: 0 : if (nnlout>3*natom) ii=6
1320 : 0 : iband=iband+1;ibs=ii+nnlout*(iblocksize-1)
1321 [ # # ]: 0 : grnl_k(1:nnlout,iband)=enlout(ibs+1:ibs+nnlout)
1322 : : end do
1323 : 0 : call timab(554,2,tsec) ! "vtowfk:rhoij"
1324 : : end if ! PAW or forces
1325 : : else
1326 [ + + + + ]: 2177314 : if(iscf>0.or.gs_hamk%usecprj==1)then
1327 [ + + + + ]: 2059918 : if (gs_hamk%usepaw==1.or.optforces/=0) then
1328 : : ! Treat all wavefunctions in case of varying occupation numbers or PAW
1329 : : ! Only treat occupied bands in case of fixed occupation numbers and NCPP
1330 [ + + + + : 1213702 : if(fixed_occ.and.abs(occblock)<=tol8.and.gs_hamk%usepaw==0) then
+ + ]
1331 [ + - + + : 242199 : if (optforces>0) grnl_k(:,(iblock-1)*blocksize+1:iblock*blocksize)=zero
+ + ]
1332 : : else
1333 [ + + ]: 1200757 : if (dtset%cprj_in_memory/=1) then
1334 [ + + ]: 1176005 : if(gs_hamk%usepaw==1) then
1335 : 980760 : call timab(554,1,tsec) ! "vtowfk:rhoij"
1336 : : end if
1337 [ + + ]: 1176005 : if(cpopt==1) then
1338 : 916 : iband=1+(iblock-1)*bandpp_cprj
1339 : 916 : call pawcprj_copy(cprj(:,1+(iblock-1)*my_nspinor*blocksize+ibg:iblock*my_nspinor*blocksize+ibg),cwaveprj)
1340 : : end if
1341 [ + + ]: 1176005 : if (mpi_enreg%paral_kgb==1) then
1342 : 85704 : call timab(572,1,tsec) ! 'prep_nonlop%vtowfk'
1343 [ + - ]: 85704 : if (gs_hamk%use_gbt == 0) then
1344 : : call prep_nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,eig_k_block,blocksize,&
1345 : : mpi_enreg,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop_prep,cwavef,nonlop_dum2,&
1346 [ - + ]: 85704 : already_transposed=.false.)
1347 : : else
1348 : 0 : gs_hamk%nspinor=1
1349 [ # # # # : 0 : ABI_MALLOC(cwavef_spin,(2,npw_k*blocksize))
# # # # #
# # # #
# ]
1350 [ # # # # : 0 : ABI_MALLOC(enlout_spin,(nnlout*blocksize))
# # # # #
# # # #
# ]
1351 [ # # ]: 0 : enlout_spin=zero
1352 : :
1353 : 0 : gs_hamk%ispin_gbt=1
1354 : 0 : call cg_copy_spin(1,npw_k,my_nspinor,blocksize,cwavef,cwavef_spin)
1355 : : call prep_nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,eig_k_block,blocksize,&
1356 : : mpi_enreg,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop_prep,cwavef_spin,nonlop_dum2,&
1357 [ # # ]: 0 : already_transposed=.false.,select_k=K_H_K)
1358 : :
1359 : 0 : gs_hamk%ispin_gbt=2
1360 : 0 : call cg_copy_spin(2,npw_k,my_nspinor,blocksize,cwavef,cwavef_spin)
1361 : : call prep_nonlop(choice,cpopt,cwaveprj,enlout_spin,gs_hamk,idir,eig_k_block,blocksize,&
1362 : : mpi_enreg,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop_prep,cwavef_spin,nonlop_dum2,&
1363 [ # # ]: 0 : already_transposed=.false.,select_k=KPRIME_H_KPRIME)
1364 : :
1365 [ # # # # ]: 0 : enlout=enlout+enlout_spin
1366 : 0 : gs_hamk%nspinor=2
1367 [ # # ]: 0 : ABI_FREE(cwavef_spin)
1368 [ # # ]: 0 : ABI_FREE(enlout_spin)
1369 : : end if
1370 : 85704 : call timab(572,2,tsec)
1371 : : else
1372 [ + - ]: 1090301 : if (gs_hamk%use_gbt == 0) then
1373 : : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,eig_k_block,&
1374 : : mpi_enreg,blocksize,nnlout,&
1375 [ - + ]: 1090301 : paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
1376 : : else
1377 : 0 : gs_hamk%nspinor=1
1378 [ # # # # : 0 : ABI_MALLOC(cwavef_spin,(2,npw_k*blocksize))
# # # # #
# # # #
# ]
1379 [ # # # # : 0 : ABI_MALLOC(enlout_spin,(nnlout*blocksize))
# # # # #
# # # #
# ]
1380 [ # # ]: 0 : enlout_spin(:)=zero
1381 : 0 : gs_hamk%ispin_gbt=1
1382 : 0 : call cg_copy_spin(1,npw_k,my_nspinor,blocksize,cwavef,cwavef_spin)
1383 : : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,eig_k_block,&
1384 : : mpi_enreg,blocksize,nnlout,&
1385 [ # # ]: 0 : paw_opt,signs,nonlop_dum,tim_nonlop,cwavef_spin,cwavef_spin,select_k=K_H_K)
1386 : 0 : gs_hamk%ispin_gbt=2
1387 : 0 : call cg_copy_spin(2,npw_k,my_nspinor,blocksize,cwavef,cwavef_spin)
1388 : : call nonlop(choice,cpopt,cwaveprj,enlout_spin,gs_hamk,idir,eig_k_block,&
1389 : : mpi_enreg,blocksize,nnlout,&
1390 [ # # ]: 0 : paw_opt,signs,nonlop_dum,tim_nonlop,cwavef_spin,cwavef_spin,select_k=KPRIME_H_KPRIME)
1391 [ # # ]: 0 : enlout(1:nnlout*blocksize)=enlout(1:nnlout*blocksize)+enlout_spin(1:nnlout*blocksize)
1392 : 0 : gs_hamk%nspinor=2
1393 [ # # ]: 0 : ABI_FREE(cwavef_spin)
1394 [ # # ]: 0 : ABI_FREE(enlout_spin)
1395 : : end if
1396 : : end if
1397 [ + + ]: 1176005 : if(gs_hamk%usepaw==1) then
1398 : 980760 : call timab(554,2,tsec)
1399 : : end if
1400 : : ! Acccumulate forces
1401 [ + + ]: 1176005 : if (optforces>0) then
1402 : 325151 : iband=(iblock-1)*blocksize
1403 [ + + ]: 785444 : do iblocksize=1,blocksize
1404 : 460293 : ii=0
1405 [ - + ]: 460293 : if (nnlout>3*natom) ii=6
1406 : 460293 : iband=iband+1;ibs=ii+nnlout*(iblocksize-1)
1407 [ + + ]: 5217326 : grnl_k(1:nnlout,iband)=enlout(ibs+1:ibs+nnlout)
1408 : : end do
1409 : : end if
1410 : : ! Store cprj (<Pnl|Psi>)
1411 [ + + ]: 1176005 : if (gs_hamk%usepaw==1.and.gs_hamk%usecprj==1) then
1412 : 980760 : iband=1+(iblock-1)*bandpp_cprj
1413 : : call pawcprj_put(gs_hamk%atindx,cwaveprj,cprj,natom,iband,ibg,ikpt,iorder_cprj,isppol,&
1414 : : mband_cprj,dtset%mkmem,natom,bandpp_cprj,nband_k_cprj,gs_hamk%dimcprj,my_nspinor,&
1415 [ - + ]: 980760 : dtset%nsppol,dtfil%unpaw,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1416 : : end if
1417 : :
1418 : : else ! cprj_in_memory==1
1419 : :
1420 : 24752 : call timab(222,1,tsec) ! 'nonlop%vtowfk'
1421 : :
1422 [ + + ]: 24752 : if ( gs_hamk%istwf_k > 1 ) then ! Real only
1423 : 10154 : space = SPACE_CR
1424 : : else ! complex
1425 : 14598 : space = SPACE_C
1426 : : end if
1427 : 24752 : me_g0 = -1
1428 [ + + ]: 24752 : if (space==SPACE_CR) then
1429 : 10154 : me_g0 = 0
1430 [ + + ]: 10154 : if (gs_hamk%istwf_k == 2) then
1431 [ + + ]: 1544 : if (mpi_enreg%me_g0 == 1) me_g0 = 1
1432 : : end if
1433 : : end if
1434 : : call xgBlock_map(xgx0,cg_k_block,space,npw_k*my_nspinor,blocksize,comm=mpi_enreg%comm_band,me_g0=me_g0,&
1435 : 24752 : gpu_option=dtset%gpu_option)
1436 : 24752 : call xgBlock_map_1d(xgeigen,eig_k_block,SPACE_R,blocksize)
1437 : :
1438 [ + + + - ]: 24752 : if (optforces/=0.or.gs_hamk%usepaw==1) then
1439 : 24752 : call xg_nonlop_getcprj(xg_nonlop,xgx0,cprj_xgx0%self,cprj_work%self)
1440 : : end if
1441 : :
1442 [ + + ]: 24752 : if (optforces/=0) then
1443 : 17752 : grnl_k_block => grnl_k(:,1+(iblock-1)*blocksize:iblock*blocksize)
1444 : 17752 : call xgBlock_map(xgforces,grnl_k_block,SPACE_R,3*natom,blocksize)
1445 : 17752 : call xg_nonlop_forces_stress(xg_nonlop,xgx0,cprj_xgx0%self,cprj_work%self,xgeigen,forces=xgforces)
1446 : : end if
1447 : :
1448 : 24752 : call timab(222,2,tsec) ! 'nonlop%vtowfk'
1449 : :
1450 [ + + ]: 24752 : if (gs_hamk%usepaw==1) then
1451 : 23072 : cprj_cwavef_bands => cprj(:,1+ibg+(iblock-1)*ncols_cprj:iblock*ncols_cprj+ibg)
1452 : 23072 : call xg_cprj_copy(cprj_cwavef_bands,cprj_xgx0%self,xg_nonlop,XG_TO_CPRJ)
1453 : : end if
1454 : :
1455 : : end if
1456 : : end if
1457 : : end if ! PAW or forces
1458 : : end if ! iscf>0 or iscf=-3
1459 : : end if
1460 : : ABI_NVTX_END_RANGE()
1461 : : end do ! End of loop on blocks
1462 : :
1463 : : #ifdef HAVE_OPENMP_OFFLOAD
1464 : : !$OMP TARGET EXIT DATA MAP(delete:cwavef) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
1465 : : !$OMP TARGET EXIT DATA MAP(delete:kinpw) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
1466 : : #endif
1467 : :
1468 : : ! restore safe value related to GEMM nonlop slicing and GPU in case of forces compute
1469 [ + + - + ]: 346501 : if(optforces==1 .and. gpu_option_tmp==ABI_GPU_OPENMP) then
1470 : 0 : gs_hamk%gpu_option = gpu_option_tmp
1471 : 0 : gemm_nonlop_block_size = blksize_gemm_nonlop_tmp
1472 : 0 : gemm_nonlop_is_distributed = is_distrib_tmp
1473 : 0 : gs_hamk%nfft_blocks = nfft_blocks_tmp
1474 : : end if
1475 : :
1476 [ + + ]: 346501 : if (dtset%cprj_in_memory==1) then
1477 : 17084 : call xg_free(cprj_xgx0)
1478 : 17084 : call xg_free(cprj_work)
1479 : : end if
1480 : :
1481 : : !call cwtime_report(" Block loop", cpu, wall, gflops)
1482 : :
1483 [ + - ]: 346501 : if(dtset%gpu_option==ABI_GPU_KOKKOS) then
1484 : : #if defined HAVE_GPU && defined HAVE_YAKL
1485 : : ABI_FREE_MANAGED(cwavef)
1486 : : #endif
1487 : : else
1488 [ - + ]: 346501 : ABI_FREE(cwavef)
1489 : : end if
1490 : :
1491 [ - + ]: 346501 : ABI_FREE(enlout)
1492 : :
1493 [ + - ]: 346501 : if (dtset%cprj_in_memory/=2) then
1494 [ + + + + : 346501 : if (gs_hamk%usepaw==1.and.(iscf>0.or.gs_hamk%usecprj==1)) then
+ - ]
1495 : 110141 : call pawcprj_free(cwaveprj)
1496 : : end if
1497 [ - + ]: 346501 : ABI_FREE(cwaveprj)
1498 : : else
1499 : 0 : nullify(cwaveprj)
1500 : : end if
1501 : :
1502 [ + + + + ]: 346501 : if (fixed_occ.and.iscf>0) then
1503 [ + - ]: 232077 : if(dtset%gpu_option==ABI_GPU_KOKKOS) then
1504 : : #if defined HAVE_GPU && defined HAVE_YAKL
1505 : : ABI_FREE_MANAGED(wfraug)
1506 : : #endif
1507 : : else
1508 [ - + ]: 232077 : ABI_FREE(wfraug)
1509 : : end if
1510 : : end if
1511 : :
1512 : : !Write the number of one-way 3D ffts skipped until now (in case of fixed occupation numbers
1513 [ + + + + : 346501 : if(iscf>0 .and. fixed_occ .and. (prtvol>2 .or. ikpt<=nkpt_max) )then
+ + - + ]
1514 : 2090 : write(msg,'(a,i0)')' vtowfk: number of one-way 3D ffts skipped in vtowfk until now =',nskip
1515 : 2090 : call wrtout(std_out,msg,'PERS')
1516 : : end if
1517 : :
1518 : : ! Norm-conserving or FockACE: Compute nonlocal+FockACE part of total energy: rotate subvnlx elements
1519 : : ! Note the two calls. For (old) lobpcgwf we have a (nband_k, nband_k) matrix, whereas cgwf
1520 : : ! returns results in packed form.
1521 : : ! CHEBYSHEV, NEW LOBPCG and RMM-DIIS do not need this
1522 : : !
1523 [ + + + + ]: 346501 : rotate_subvnlx = gs_hamk%usepaw == 0 .and. wfopta10 /= 1 .and. .not. xg_diago
1524 [ + + ]: 346501 : if (use_rmm_diis) rotate_subvnlx = .False.
1525 : :
1526 [ + + ]: 345953 : if (rotate_subvnlx) then
1527 : 224724 : call timab(586,1,tsec) ! 'vtowfk(nonlocalpart)'
1528 [ + + ]: 224724 : if (wfopta10==4) then
1529 : 612 : call cg_hrotate_and_get_diag(istwf_k, nband_k, totvnlx, evec, enlx_k)
1530 : : else
1531 : 224112 : call cg_hprotate_and_get_diag(nband_k, subvnlx, evec, enlx_k)
1532 : : end if
1533 : 224724 : call timab(586,2,tsec)
1534 : : end if
1535 : :
1536 : : !###################################################################
1537 : :
1538 [ + + + + ]: 346501 : if (iscf<=0 .and. max_resid > dtset%tolwfr) then
1539 : : write(msg,'(2(a,i0),a,es13.5)')&
1540 : 710 : "Wavefunctions not converged for ikpt: ", ikpt, ", nnsclo: ",nnsclo_now,', max resid: ',max_resid
1541 : 710 : ABI_WARNING(msg)
1542 : : end if
1543 : :
1544 : : ! Print out eigenvalues (hartree)
1545 [ + + ]: 346501 : if (mod(dtset%wfoptalg,10)==1) then
1546 : 23464 : niter=dtset%mdeg_filter
1547 : 23464 : iter_name='as the degree of the polynomial filter'
1548 : : else
1549 : 323037 : niter=dtset%nline
1550 : 323037 : iter_name='CG line minimizations'
1551 : : end if
1552 [ + + + + : 346501 : if (prtvol/=5.and.(prtvol>2 .or. ikpt<=nkpt_max)) then
- + ]
1553 : : write(msg, '(5x,a,i5,2x,a,a,a,i4,a,i4,2a)' ) &
1554 : 14055 : 'eigenvalues (hartree) for',nband_k,'bands',ch10,&
1555 : 28110 : ' after ',inonsc,' non-SCF iterations with ',niter,' ',trim(iter_name)
1556 : 14055 : call wrtout(std_out,msg,'PERS')
1557 [ + + ]: 56224 : do ii=0,(nband_k-1)/6
1558 : 42169 : write(msg, '(1p,6e12.4)' ) (eig_k(iband),iband=1+6*ii,min(6+6*ii,nband_k))
1559 : 56224 : call wrtout(std_out,msg,'PERS')
1560 : : end do
1561 [ - + ]: 332446 : else if(ikpt==nkpt_max+1)then
1562 : 0 : call wrtout(std_out,' vtowfk : prtvol=0 or 1, do not print more k-points.','PERS')
1563 : : end if
1564 : :
1565 : : ! Print out decomposition of eigenvalues in the non-selfconsistent case or if prtvol>=10
1566 [ + + + + : 346501 : if( (iscf<0 .and. (prtvol>2 .or. ikpt<=nkpt_max)) .or. prtvol>=10)then
+ - + + ]
1567 : : write(msg, '(5x,a,i5,2x,a,a,a,i4,a,i4,2a)' ) &
1568 : 4449 : ' mean kinetic energy (hartree) for ',nband_k,' bands',ch10,&
1569 : 8898 : ' after ',inonsc,' non-SCF iterations with ',niter,' ',trim(iter_name)
1570 : 4449 : call wrtout(std_out,msg,'PERS')
1571 : :
1572 [ + + ]: 13012 : do ii=0,(nband_k-1)/6
1573 : 8563 : write(msg, '(1p,6e12.4)' ) (ek_k(iband),iband=1+6*ii,min(6+6*ii,nband_k))
1574 : 13012 : call wrtout(std_out,msg,'PERS')
1575 : : end do
1576 : :
1577 [ + + ]: 4449 : if (gs_hamk%usepaw==0) then
1578 : : write(msg, '(5x,a,i5,2x,a,a,a,i4,a,i4,2a)' ) &
1579 : 3076 : ' mean NL+Fock-type energy (hartree) for ',nband_k,' bands',ch10,&
1580 : 6152 : ' after ',inonsc,' non-SCF iterations with ',niter,' ',trim(iter_name)
1581 : 3076 : call wrtout(std_out,msg,'PERS')
1582 : :
1583 [ + + ]: 8029 : do ii=0,(nband_k-1)/6
1584 : 4953 : write(msg,'(1p,6e12.4)') (enlx_k(iband),iband=1+6*ii,min(6+6*ii,nband_k))
1585 : 8029 : call wrtout(std_out,msg,'PERS')
1586 : : end do
1587 : : end if
1588 : : end if
1589 : :
1590 : : ! Hamiltonian constructor for gwls_sternheimer
1591 [ + + ]: 346501 : if (dtset%optdriver==RUNL_GWLS) call build_H(dtset,mpi_enreg,cpopt,cg,gs_hamk,kg_k,kinpw)
1592 : :
1593 [ - + ]: 346501 : if (dtset%cprj_in_memory==2) nullify(cprj_cwavef_bands)
1594 : :
1595 : : #ifdef HAVE_OPENMP_OFFLOAD
1596 : : if(transfer_cg) then
1597 : : !$OMP TARGET EXIT DATA MAP(from:cg_k) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP .and. xg_diago)
1598 : : end if
1599 : : #endif
1600 [ + + ]: 346501 : if(wfopta10 /= 1 .and. .not. xg_diago) then
1601 [ - + ]: 299790 : ABI_FREE(evec)
1602 [ - + ]: 299790 : ABI_FREE(subham)
1603 [ - + ]: 299790 : ABI_FREE(totvnlx)
1604 [ - + ]: 299790 : ABI_FREE(subvnlx)
1605 [ - + ]: 299790 : ABI_FREE(subovl)
1606 : : end if
1607 : :
1608 [ + - ]: 346501 : ABI_SFREE(gsc)
1609 : :
1610 [ + + ]: 346501 : if(wfoptalg==3) then
1611 [ - + ]: 6 : ABI_FREE(eig_save)
1612 : : end if
1613 : :
1614 [ - + ]: 346501 : if (prtvol==-level) then
1615 : : ! Structured debugging: if prtvol=-level, stop here.
1616 : 0 : write(msg,'(3a,i0,a)')' vtowfk: exit ',ch10,' prtvol=-',level,', debugging mode => stop '
1617 : 0 : ABI_ERROR(msg)
1618 : : end if
1619 : :
1620 : 346501 : call timab(30,2,tsec)
1621 : 346501 : call timab(28,2,tsec)
1622 : :
1623 : : DBG_EXIT("COLL")
1624 : :
1625 [ - + - + : 693002 : end subroutine vtowfk
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
1626 : : !!***
1627 : :
1628 : : end module m_vtowfk
1629 : : !!***
|
