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