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 343291 : subroutine vtowfk(cg,cgq,cprj,cpus,dphase_k,dtefield,dtfil,dtset,&
180 343291 : & eig_k,ek_k,ek_k_nd,end_k,enlx_k,fixed_occ,grnl_k,gs_hamk,&
181 343291 : & ibg,icg,ikpt,iscf,isppol,kg_k,kinpw,mband_cprj,mcg,mcgq,mcprj,mkgq,mpi_enreg,&
182 343291 : & mpw,natom,nband_k,nbdbuf,nkpt,istep,nnsclo_now,npw_k,npwarr,occ_k,optforces,prtvol,&
183 343291 : & 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 343291 : 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 343291 : real(dp),allocatable,target :: cwavef(:,:)
244 343291 : real(dp),allocatable :: wfraug(:,:,:,:)
245 : #endif
246 343291 : real(dp),allocatable :: eig_save(:),enlout(:),evec(:,:),gsc(:,:),ghc_vectornd(:,:)
247 343291 : real(dp),allocatable :: subham(:),subovl(:),subvnlx(:),totvnlx(:,:)
248 343291 : real(dp),pointer :: cg_k(:,:),cg_k_block(:,:),grnl_k_block(:,:),eig_k_block(:)
249 343291 : real(dp),contiguous, pointer :: cwavef_iband(:,:), kinpw_ptr(:)
250 343291 : real(dp), allocatable :: weight_t(:) ! only allocated and used with GPU fourwf
251 343291 : 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 343291 : call timab(28,1,tsec) ! Keep track of total time spent in "vtowfk"
259 :
260 : !Structured debugging if prtvol==-level
261 343291 : 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 343291 : nkpt_max=50; if(xmpi_paral==1)nkpt_max=-1
271 :
272 343291 : wfoptalg=mod(dtset%wfoptalg,100); wfopta10=mod(wfoptalg,10)
273 343291 : xg_diago = dtset%wfoptalg == 114 .or. dtset%wfoptalg == 112 .or. dtset%wfoptalg == 111
274 343291 : istwf_k=gs_hamk%istwf_k
275 343291 : has_fock=(associated(gs_hamk%fockcommon))
276 343291 : quit=0
277 343291 : igsc=0
278 :
279 : !Parallelization over spinors management
280 343291 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
281 343291 : if (mpi_enreg%paral_spinor==0) then
282 342731 : ispinor_index=1
283 342731 : nspinor1TreatedByThisProc=.true.
284 342731 : 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 343291 : 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 343291 : blocksize=nband_k/nblockbd
299 :
300 : !Save eshift
301 343291 : 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 343291 : 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 343291 : use_rmm_diis = .False.
313 343291 : 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 343291 : mgsc=0
320 : igsc=0
321 343291 : 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 43959 : ABI_MALLOC(gsc,(0,0))
327 : end if
328 :
329 343291 : 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 343291 : 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 343291 : 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 343291 : 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 343291 : call timab(39,1,tsec) ! "vtowfk (loop)"
386 :
387 343291 : cg_k => cg(:,1+icg:npw_k*my_nspinor*nband_k+icg)
388 :
389 343291 : 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 818326 : do inonsc=1,nnsclo_now
403 : ABI_NVTX_START_RANGE(NVTX_VTOWFK_EXTRA1)
404 497291 : if (iscf < 0 .and. (inonsc <= enough .or. mod(inonsc, 10) == 0)) call cwtime(cpu, wall, gflops, "start")
405 :
406 497291 : 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 497291 : 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 4371944 : do iband=1,nband_k
427 3874653 : iwavef=(iband-1)*npw_k*my_nspinor+icg
428 3874653 : cwavef_iband => cg(:,1+iwavef:npw_k*my_nspinor+iwavef)
429 3874653 : update_cprj=.False.
430 8231898 : do ispinor=1,my_nspinor
431 4357245 : iispinor=ispinor; if (mpi_enreg%paral_spinor==1) iispinor=mpi_enreg%me_spinor+1
432 4357245 : igs=(ispinor-1)*npw_k
433 4357245 : kinpw_ptr => kinpw
434 4357245 : 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 1157736091 : do ipw=1+igs,npw_k+igs
440 1153861438 : 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 4371944 : 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 497291 : if(dtset%optdriver/=RUNL_GWLS) then
458 :
459 497283 : if(wfopta10==4.or.wfopta10==1) then
460 :
461 55217 : 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 55217 : if (wfopta10==4) then
470 :
471 30777 : 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 30349 : 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 28099 : 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 24063 : 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 24440 : else if (wfopta10 == 1) then
517 24440 : 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 24392 : 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 9024 : & 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 497291 : if (nbdbuf >= 0) then
622 4506800 : 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 497291 : 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 566785 : call wrtout(std_out,msg,'PERS')
634 : end do
635 : end if
636 :
637 : ! =========================================================================
638 : ! ========== DIAGONALIZATION OF HAMILTONIAN IN WFs SUBSPACE ===============
639 : ! =========================================================================
640 497291 : do_subdiago = .not. wfopta10 == 1 .and. .not. xg_diago
641 497291 : 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 496638 : 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 497291 : 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 566785 : 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 497291 : call timab(583,1,tsec) ! "vtowfk(pw_orthon)"
690 497291 : 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 497291 : do_ortho = (wfoptalg/=14 .and. wfoptalg /= 1 .and. wfoptalg /= 12 .and. wfoptalg /= 11 .and. dtset%ortalg>-10) .or. dtset%ortalg > 0
694 497291 : if (xg_diago) do_ortho = .false.
695 497291 : if (use_rmm_diis) do_ortho = .False.
696 :
697 496638 : 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 497291 : 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 497291 : if (xmpi_paral/=1 .or. mpi_enreg%paral_kgb/=1) then
721 : !call wrtout(std_out, "Calling cgtk_fixphase")
722 453459 : 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 10789 : 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 497291 : 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 818326 : if (max_resid < dtset%tolwfr) then
756 22256 : 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 343291 : 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 343291 : call timab(39,2,tsec)
769 343291 : 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 343291 : ndat=1;if (mpi_enreg%paral_kgb==1) ndat=mpi_enreg%bandpp
780 343291 : if(iscf>0 .and. fixed_occ) then
781 229430 : ndat_fft=ndat; if(mpi_enreg%paral_kgb==0) ndat_fft=blocksize
782 229430 : chunk = ndat_fft/gs_hamk%nfft_blocks ! Divide by 2 to construct chunk of even number of bands
783 229430 : residuchunk = ndat_fft - gs_hamk%nfft_blocks*chunk
784 229430 : 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 229430 : fftalg = gs_hamk%ngfft(7)
790 229430 : if (gs_hamk%gpu_option==ABI_GPU_DISABLED.and.fftalg/=401) then
791 1124995 : ABI_MALLOC(wfraug,(2,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6))
792 : else
793 22155 : 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 343291 : nnlout=3*natom*optforces
800 343291 : signs=1;idir=0
801 343291 : if (gs_hamk%usepaw==0) then
802 235587 : choice=1+optforces
803 235587 : paw_opt=0;cpopt=-1;tim_nonlop=2
804 : else
805 107704 : choice=2*optforces
806 107704 : paw_opt=2;cpopt=0;tim_nonlop=10-8*optforces
807 107704 : if (dtset%cprj_in_memory==2) cpopt=2 ! cprj are in memory (but not the derivatives)
808 107704 : 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 1029873 : ABI_MALLOC(enlout,(nnlout*blocksize))
817 :
818 : ! Allocation of memory space for one block of waveforms containing blocksize waveforms
819 343291 : 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 1029873 : ABI_MALLOC(cwavef, (2,npw_k*my_nspinor*blocksize))
825 : end if
826 :
827 343291 : if (dtset%cprj_in_memory/=2) then
828 343291 : if (gs_hamk%usepaw==1.and.(iscf>0.or.gs_hamk%usecprj==1)) then
829 107704 : iorder_cprj=0
830 107704 : nband_k_cprj=nband_k*(mband_cprj/dtset%mband)
831 107704 : bandpp_cprj=mpi_enreg%bandpp
832 1341864 : ABI_MALLOC(cwaveprj,(natom,my_nspinor*bandpp_cprj))
833 107704 : ncpgr=0;if (cpopt==1) ncpgr=cprj(1,1)%ncpgr
834 107704 : 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 343291 : 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 343291 : gpu_option_tmp=gs_hamk%gpu_option
857 343291 : 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 2493832 : do iblock=1,nblockbd
887 6914559 : 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 2150541 : 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 2150541 : 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 2150541 : 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 2124877 : & 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 23213501 : if(ANY(ABS(dtset%nucdipmom)>tol8)) then
922 195300 : ABI_MALLOC(ghc_vectornd,(2,npw_k*my_nspinor*blocksize))
923 : call getghc_nucdip(cwavef,ghc_vectornd,gs_hamk%gbound_k,gs_hamk%istwf_k,kg_k,gs_hamk%kpt_k,&
924 : & gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%ngfft,npw_k,gs_hamk%nvloc,&
925 : & gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,my_nspinor,gs_hamk%vectornd,gs_hamk%vlocal,&
926 65100 : & gs_hamk%zora,gs_hamk%gpu_option)
927 130200 : do iblocksize=1,blocksize
928 65100 : iband=(iblock-1)*blocksize+iblocksize
929 : end_k(iband)=DOT_PRODUCT(cwavef(1,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
930 : & ghc_vectornd(1,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor))+ &
931 : & DOT_PRODUCT(cwavef(2,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor),&
932 92500512 : & ghc_vectornd(2,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor))
933 : end do
934 65100 : ABI_FREE(ghc_vectornd)
935 : end if
936 :
937 : end if
938 :
939 4764018 : do iblocksize=1,blocksize
940 2613477 : iband=(iblock-1)*blocksize+iblocksize
941 :
942 2613477 : cwavef_iband => cg(:,1+(iband-1)*npw_k*my_nspinor+icg:iband*npw_k*my_nspinor+icg)
943 :
944 2613477 : if (gs_hamk%use_gbt /= 0) then
945 : ! Treat up and down components separately.
946 : ! Note filter 1. Also: this won't work if paral_kgb 1 and/or spinor parallelism
947 25664 : filter = 1
948 : call meanvalue_g(ar, gs_hamk%kinpw_k, filter,istwf_k,mpi_enreg,npw_k,1,&
949 : & cwavef_iband, cwavef_iband, 0,&
950 25664 : & gpu_thread_limit=dtset%gpu_thread_limit)
951 :
952 : call meanvalue_g(ar2,gs_hamk%kinpw_kp,filter,istwf_k,mpi_enreg,npw_k,1,&
953 : & cwavef_iband(:,npw_k+1:), cwavef_iband(:,npw_k+1:),0,&
954 25664 : & gpu_thread_limit=dtset%gpu_thread_limit)
955 :
956 25664 : ek_k(iband) = ar + ar2
957 : end if
958 :
959 4764018 : if(paw_dmft%use_dmft==1) then
960 529272 : do iband1=1,nband_k
961 : call meanvalue_g(ar,kinpw,0,istwf_k,mpi_enreg,npw_k,my_nspinor,&
962 : cg(:,1+(iband -1)*npw_k*my_nspinor+icg:iband *npw_k*my_nspinor+icg),&
963 : cg(:,1+(iband1-1)*npw_k*my_nspinor+icg:iband1*npw_k*my_nspinor+icg),&
964 510096 : paw_dmft%use_dmft,ar_im=ar_im,gpu_thread_limit=dtset%gpu_thread_limit)
965 510096 : ek_k_nd(1,iband,iband1)=ar
966 529272 : ek_k_nd(2,iband,iband1)=ar_im
967 : end do
968 : end if
969 : end do
970 :
971 2150541 : if (iscf>0) then
972 :
973 : ABI_NVTX_START_RANGE(NVTX_VTOWFK_FOURWF)
974 : ! In case of fixed occupation numbers, accumulates the partial density
975 2026219 : if (fixed_occ .and. mpi_enreg%paral_kgb/=1) then
976 :
977 : ! treat all bands at once on GPU
978 858866 : if (dtset%gpu_option /= ABI_GPU_DISABLED) then
979 :
980 0 : ABI_MALLOC(weight_t,(blocksize))
981 : ! compute weights
982 0 : do iblocksize=1,blocksize
983 0 : iband=(iblock-1)*blocksize+iblocksize
984 0 : weight_t(iblocksize) = occ_k(iband) * wtk / gs_hamk%ucvol
985 0 : if (abs(occ_k(iband)) < tol8) weight_t(iblocksize) = zero
986 : end do
987 :
988 0 : if(dtset%nspinor==1) then
989 : call fourwf_optmem(1,rhoaug(:,:,:,1),cwavef(:,:),dummy,wfraug,&
990 : gs_hamk%gbound_k,gs_hamk%gbound_k,istwf_k,kg_k,kg_k,&
991 : gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%nfft_blocks,gs_hamk%ngfft,&
992 : npw_k,1,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,tim_fourwf,weight,weight,&
993 : weight_array_r=weight_t,weight_array_i=weight_t,&
994 0 : gpu_option=dtset%gpu_option)
995 :
996 0 : else if(dtset%nspinor==2) then
997 0 : ABI_MALLOC(cwavefb,(2,npw_k*blocksize,2))
998 0 : ibs=(iblock-1)*npw_k*my_nspinor*blocksize
999 0 : if(gs_hamk%gpu_option==ABI_GPU_OPENMP) then
1000 : #ifdef HAVE_OPENMP_OFFLOAD
1001 : !$OMP TARGET ENTER DATA MAP(alloc:cwavefb)
1002 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cwavefb,cwavef) PRIVATE(iband)
1003 : do iband=1,blocksize
1004 : !$OMP PARALLEL DO PRIVATE(ipw)
1005 : do ipw=1,npw_k
1006 : cwavefb(1,(iband-1)*npw_k+ipw,1)=cwavef(1,ipw+(2*iband-2)*npw_k)
1007 : cwavefb(2,(iband-1)*npw_k+ipw,1)=cwavef(2,ipw+(2*iband-2)*npw_k)
1008 :
1009 : cwavefb(1,(iband-1)*npw_k+ipw,2)=cwavef(1,ipw+(2*iband-1)*npw_k)
1010 : cwavefb(2,(iband-1)*npw_k+ipw,2)=cwavef(2,ipw+(2*iband-1)*npw_k)
1011 : end do
1012 : end do
1013 : #endif
1014 : else
1015 0 : do iband=1,blocksize
1016 0 : cwavefb(:,(iband-1)*npw_k+1:iband*npw_k,1)=cwavef(:,1+(2*iband-2)*npw_k:(iband*2-1)*npw_k)
1017 0 : cwavefb(:,(iband-1)*npw_k+1:iband*npw_k,2)=cwavef(:,1+(2*iband-1)*npw_k:iband*2*npw_k)
1018 : end do
1019 : end if
1020 :
1021 : call fourwf_optmem(1,rhoaug(:,:,:,1),cwavefb(:,:,1),dummy,wfraug,&
1022 : gs_hamk%gbound_k,gs_hamk%gbound_k,istwf_k,kg_k,kg_k,&
1023 : gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%nfft_blocks,gs_hamk%ngfft,&
1024 : npw_k,1,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,tim_fourwf,weight,weight,&
1025 : weight_array_r=weight_t,weight_array_i=weight_t,&
1026 0 : gpu_option=dtset%gpu_option)
1027 :
1028 0 : if(dtset%nspden==1) then
1029 : call fourwf_optmem(1,rhoaug(:,:,:,1),cwavefb(:,:,2),dummy,wfraug,&
1030 : gs_hamk%gbound_k,gs_hamk%gbound_k,istwf_k,kg_k,kg_k,&
1031 : gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%nfft_blocks,gs_hamk%ngfft,&
1032 : npw_k,1,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,tim_fourwf,weight,weight,&
1033 : weight_array_r=weight_t,weight_array_i=weight_t,&
1034 0 : gpu_option=dtset%gpu_option)
1035 :
1036 0 : else if (dtset%nspden==4) then
1037 : ! Build the four components of rho. We use only norm quantities and, so fourwf.
1038 : ! $\sum_{n} f_n \Psi^{* \alpha}_n \Psi^{\alpha}_n =\rho^{\alpha \alpha}$
1039 : ! $\sum_{n} f_n (\Psi^{1}+\Psi^{2})^*_n (\Psi^{1}+\Psi^{2})_n=rho+m_x$
1040 : ! $\sum_{n} f_n (\Psi^{1}-i \Psi^{2})^*_n (\Psi^{1}-i \Psi^{2})_n=rho+m_y$
1041 0 : ABI_MALLOC(cwavef_x,(2,npw_k*blocksize))
1042 0 : ABI_MALLOC(cwavef_y,(2,npw_k*blocksize))
1043 :
1044 0 : if(gs_hamk%gpu_option==ABI_GPU_OPENMP) then
1045 : #ifdef HAVE_OPENMP_OFFLOAD
1046 : !$OMP TARGET ENTER DATA MAP(alloc:cwavef_x,cwavef_y)
1047 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cwavef_x,cwavef_y,cwavefb) PRIVATE(iband)
1048 : do iband=1,blocksize
1049 : !$OMP PARALLEL DO PRIVATE(ipw)
1050 : do ipw=1,npw_k
1051 : !$(\Psi^{1}+\Psi^{2})$
1052 : cwavef_x(1,(iband-1)*npw_k+ipw)=cwavefb(1,(iband-1)*npw_k+ipw,1)+cwavefb(1,(iband-1)*npw_k+ipw,2)
1053 : cwavef_x(2,(iband-1)*npw_k+ipw)=cwavefb(2,(iband-1)*npw_k+ipw,1)+cwavefb(2,(iband-1)*npw_k+ipw,2)
1054 : !$(\Psi^{1}-i \Psi^{2})$
1055 : cwavef_y(1,(iband-1)*npw_k+ipw)=cwavefb(1,(iband-1)*npw_k+ipw,1)+cwavefb(2,(iband-1)*npw_k+ipw,2)
1056 : cwavef_y(2,(iband-1)*npw_k+ipw)=cwavefb(2,(iband-1)*npw_k+ipw,1)-cwavefb(1,(iband-1)*npw_k+ipw,2)
1057 : end do
1058 : end do
1059 : #endif
1060 : else
1061 : !$(\Psi^{1}+\Psi^{2})$
1062 0 : cwavef_x(:,:)=cwavefb(:,1:npw_k*blocksize,1)+cwavefb(:,1:npw_k*blocksize,2)
1063 : !$(\Psi^{1}-i \Psi^{2})$
1064 0 : cwavef_y(1,:)=cwavefb(1,1:npw_k*blocksize,1)+cwavefb(2,1:npw_k*blocksize,2)
1065 0 : cwavef_y(2,:)=cwavefb(2,1:npw_k*blocksize,1)-cwavefb(1,1:npw_k*blocksize,2)
1066 : end if
1067 :
1068 : ! z component
1069 : call fourwf_optmem(1,rhoaug(:,:,:,4),cwavefb(:,:,2),dummy,wfraug,&
1070 : gs_hamk%gbound_k,gs_hamk%gbound_k,istwf_k,kg_k,kg_k,&
1071 : gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%nfft_blocks,gs_hamk%ngfft,&
1072 : npw_k,1,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,tim_fourwf,weight,weight,&
1073 : weight_array_r=weight_t,weight_array_i=weight_t,&
1074 0 : gpu_option=dtset%gpu_option)
1075 : ! x component
1076 : call fourwf_optmem(1,rhoaug(:,:,:,2),cwavef_x(:,:),dummy,wfraug,&
1077 : gs_hamk%gbound_k,gs_hamk%gbound_k,istwf_k,kg_k,kg_k,&
1078 : gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%nfft_blocks,gs_hamk%ngfft,&
1079 : npw_k,1,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,tim_fourwf,weight,weight,&
1080 : weight_array_r=weight_t,weight_array_i=weight_t,&
1081 0 : gpu_option=dtset%gpu_option)
1082 : ! y component
1083 : call fourwf_optmem(1,rhoaug(:,:,:,3),cwavef_y(:,:),dummy,wfraug,&
1084 : gs_hamk%gbound_k,gs_hamk%gbound_k,istwf_k,kg_k,kg_k,&
1085 : gs_hamk%mgfft,mpi_enreg,blocksize,gs_hamk%nfft_blocks,gs_hamk%ngfft,&
1086 : npw_k,1,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,tim_fourwf,weight,weight,&
1087 : weight_array_r=weight_t,weight_array_i=weight_t,&
1088 0 : gpu_option=dtset%gpu_option)
1089 : #ifdef HAVE_OPENMP_OFFLOAD
1090 : !$OMP TARGET EXIT DATA MAP(delete:cwavef_x,cwavef_y) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
1091 : #endif
1092 0 : ABI_FREE(cwavef_x)
1093 0 : ABI_FREE(cwavef_y)
1094 : end if
1095 :
1096 : #ifdef HAVE_OPENMP_OFFLOAD
1097 : !$OMP TARGET EXIT DATA MAP(delete:cwavefb) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
1098 : #endif
1099 0 : ABI_FREE(cwavefb)
1100 : end if
1101 :
1102 0 : ABI_FREE(weight_t)
1103 :
1104 : else
1105 1792812 : do iblocksize=1,blocksize
1106 933946 : iband=(iblock-1)*blocksize+iblocksize
1107 933946 : cwavef_iband => cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor)
1108 :
1109 1792812 : if (abs(occ_k(iband))>=tol8) then
1110 703285 : weight = occ_k(iband) * wtk / gs_hamk%ucvol
1111 :
1112 : ! Accumulate charge density in real space in array rhoaug
1113 : ! The same section of code is also found in mkrho.F90 : should be rationalized !
1114 : call fourwf(1,rhoaug(:,:,:,1),cwavef_iband,dummy,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
1115 : istwf_k,gs_hamk%kg_k,gs_hamk%kg_k,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%ngfft,npw_k,1,&
1116 : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,&
1117 703285 : tim_fourwf,weight,weight,gpu_option=dtset%gpu_option)
1118 :
1119 703285 : if(dtset%nspinor==2)then
1120 244866 : ABI_MALLOC(cwavef1,(2,npw_k))
1121 93039964 : cwavef1(:,:)=cwavef_iband(:,1+npw_k:2*npw_k) ! EB FR spin dn part and used for m_z component (cwavef_z)
1122 :
1123 81622 : if(dtset%nspden==1) then
1124 :
1125 : call fourwf(1,rhoaug(:,:,:,1),cwavef1,dummy,wfraug,&
1126 : gs_hamk%gbound_k,gs_hamk%gbound_k,&
1127 : istwf_k,gs_hamk%kg_k,gs_hamk%kg_k,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%ngfft,npw_k,1,&
1128 : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,&
1129 40292 : tim_fourwf,weight,weight,gpu_option=dtset%gpu_option)
1130 :
1131 41330 : else if(dtset%nspden==4) then
1132 : ! Build the four components of rho. We use only norm quantities and, so fourwf.
1133 : ! $\sum_{n} f_n \Psi^{* \alpha}_n \Psi^{\alpha}_n =\rho^{\alpha \alpha}$
1134 : ! $\sum_{n} f_n (\Psi^{1}+\Psi^{2})^*_n (\Psi^{1}+\Psi^{2})_n=rho+m_x$
1135 : ! $\sum_{n} f_n (\Psi^{1}-i \Psi^{2})^*_n (\Psi^{1}-i \Psi^{2})_n=rho+m_y$
1136 123990 : ABI_MALLOC(cwavef_x,(2,npw_k))
1137 82660 : ABI_MALLOC(cwavef_y,(2,npw_k))
1138 : !$(\Psi^{1}+\Psi^{2})$
1139 26454140 : cwavef_x(:,:)=cwavef_iband(:,1:npw_k)+cwavef1(:,1:npw_k)
1140 : !$(\Psi^{1}-i \Psi^{2})$
1141 8845600 : cwavef_y(1,:)=cwavef_iband(1,1:npw_k)+cwavef1(2,1:npw_k)
1142 8845600 : cwavef_y(2,:)=cwavef_iband(2,1:npw_k)-cwavef1(1,1:npw_k)
1143 : ! z component
1144 : call fourwf(1,rhoaug(:,:,:,4),cwavef1,dummy,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
1145 : istwf_k,gs_hamk%kg_k,gs_hamk%kg_k,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%ngfft,npw_k,1,&
1146 : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,&
1147 41330 : tim_fourwf,weight,weight,gpu_option=dtset%gpu_option)
1148 : ! x component
1149 : call fourwf(1,rhoaug(:,:,:,2),cwavef_x,dummy,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
1150 : istwf_k,gs_hamk%kg_k,gs_hamk%kg_k,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%ngfft,npw_k,1,&
1151 : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,&
1152 41330 : tim_fourwf,weight,weight,gpu_option=dtset%gpu_option)
1153 : ! y component
1154 : call fourwf(1,rhoaug(:,:,:,3),cwavef_y,dummy,wfraug,gs_hamk%gbound_k,gs_hamk%gbound_k,&
1155 : istwf_k,gs_hamk%kg_k,gs_hamk%kg_k,gs_hamk%mgfft,mpi_enreg,1,gs_hamk%ngfft,npw_k,1,&
1156 : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,1,&
1157 41330 : tim_fourwf,weight,weight,gpu_option=dtset%gpu_option)
1158 :
1159 41330 : ABI_FREE(cwavef_x)
1160 41330 : ABI_FREE(cwavef_y)
1161 :
1162 : end if ! dtset%nspden/=4
1163 81622 : ABI_FREE(cwavef1)
1164 : end if
1165 : else
1166 230661 : nskip=nskip+1
1167 : end if
1168 : end do ! Loop inside a block of bands
1169 :
1170 : end if ! dtset%gpu_option
1171 :
1172 :
1173 1167353 : else if (fixed_occ .and. mpi_enreg%paral_kgb==1) then
1174 : ! In case of fixed occupation numbers,in bandFFT mode accumulates the partial density
1175 :
1176 68990 : if (dtset%nspinor==1) then
1177 48270 : call timab(537,1,tsec) ! "prep_fourwf%vtow"
1178 : call prep_fourwf(rhoaug(:,:,:,1),blocksize,cwavef,wfraug,iblock,istwf_k,&
1179 : gs_hamk%mgfft,mpi_enreg,nband_k,ndat,gs_hamk%ngfft,npw_k,&
1180 : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,occ_k,&
1181 48270 : 1,gs_hamk%ucvol,wtk,gs_hamk%nfft_blocks,gpu_option=dtset%gpu_option)
1182 48270 : call timab(537,2,tsec)
1183 20720 : else if (dtset%nspinor==2) then
1184 82880 : ABI_MALLOC(cwavefb,(2,npw_k*blocksize,2))
1185 20720 : if(gs_hamk%gpu_option==ABI_GPU_OPENMP) then
1186 : ! --- No parallelization over spinors ---
1187 0 : if (mpi_enreg%paral_spinor==0) then
1188 : #ifdef HAVE_OPENMP_OFFLOAD
1189 : !$OMP TARGET ENTER DATA MAP(alloc:cwavefb)
1190 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cwavefb,cwavef) PRIVATE(iband)
1191 : do iband=1,blocksize
1192 : !$OMP PARALLEL DO PRIVATE(ipw)
1193 : do ipw=1,npw_k
1194 : cwavefb(1,(iband-1)*npw_k+ipw,1)=cwavef(1,ipw+(2*iband-2)*npw_k)
1195 : cwavefb(2,(iband-1)*npw_k+ipw,1)=cwavef(2,ipw+(2*iband-2)*npw_k)
1196 :
1197 : cwavefb(1,(iband-1)*npw_k+ipw,2)=cwavef(1,ipw+(2*iband-1)*npw_k)
1198 : cwavefb(2,(iband-1)*npw_k+ipw,2)=cwavef(2,ipw+(2*iband-1)*npw_k)
1199 : end do
1200 : end do
1201 : #endif
1202 : else
1203 0 : ABI_BUG("Parallelisation on spinor isn't supported with OpenMP GPU")
1204 : end if
1205 : else
1206 : ! --- No parallelization over spinors ---
1207 20720 : if (mpi_enreg%paral_spinor==0) then
1208 95680 : do iband=1,blocksize
1209 9033380 : cwavefb(:,(iband-1)*npw_k+1:iband*npw_k,1)=cwavef(:,1+(2*iband-2)*npw_k:(iband*2-1)*npw_k)
1210 9045940 : cwavefb(:,(iband-1)*npw_k+1:iband*npw_k,2)=cwavef(:,1+(2*iband-1)*npw_k:iband*2*npw_k)
1211 : end do
1212 : else
1213 : ! --- Parallelization over spinors ---
1214 : ! (split the work between 2 procs)
1215 9855216 : cwavefb(:,:,3-ispinor_index)=zero
1216 24832 : do iband=1,blocksize
1217 9871888 : cwavefb(:,(iband-1)*npw_k+1:iband*npw_k,ispinor_index) = cwavef(:,1+(iband-1)*npw_k:iband*npw_k)
1218 : end do
1219 8160 : call xmpi_sum(cwavefb,mpi_enreg%comm_spinor,ierr)
1220 : end if
1221 : end if !gpu_option
1222 :
1223 20720 : call timab(537,1,tsec) !"prep_fourwf%vtow"
1224 20720 : if (nspinor1TreatedByThisProc) then
1225 : call prep_fourwf(rhoaug(:,:,:,1),blocksize,cwavefb(:,:,1),wfraug,iblock,&
1226 : istwf_k,gs_hamk%mgfft,mpi_enreg,nband_k,ndat,gs_hamk%ngfft,npw_k,&
1227 : gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,occ_k,1,gs_hamk%ucvol,wtk,gs_hamk%nfft_blocks,&
1228 16640 : gpu_option=dtset%gpu_option)
1229 : end if
1230 20720 : if(dtset%nspden==1) then
1231 9376 : if (nspinor2TreatedByThisProc) then
1232 : call prep_fourwf(rhoaug(:,:,:,1),blocksize,cwavefb(:,:,2),wfraug,&
1233 : iblock,istwf_k,gs_hamk%mgfft,mpi_enreg,nband_k,ndat,&
1234 : gs_hamk%ngfft,npw_k,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,occ_k,1,&
1235 8768 : gs_hamk%ucvol,wtk,gs_hamk%nfft_blocks,gpu_option=dtset%gpu_option)
1236 : end if
1237 11344 : else if(dtset%nspden==4) then
1238 34032 : ABI_MALLOC(cwavef_x,(2,npw_k*blocksize))
1239 22688 : ABI_MALLOC(cwavef_y,(2,npw_k*blocksize))
1240 11344 : if(gs_hamk%gpu_option==ABI_GPU_OPENMP) then
1241 : #ifdef HAVE_OPENMP_OFFLOAD
1242 : !$OMP TARGET ENTER DATA MAP(alloc:cwavef_x,cwavef_y)
1243 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cwavef_x,cwavef_y,cwavefb) PRIVATE(iband)
1244 : do iband=1,blocksize
1245 : !$OMP PARALLEL DO PRIVATE(ipw)
1246 : do ipw=1,npw_k
1247 : !$(\Psi^{1}+\Psi^{2})$
1248 : cwavef_x(1,(iband-1)*npw_k+ipw)=cwavefb(1,(iband-1)*npw_k+ipw,1)+cwavefb(1,(iband-1)*npw_k+ipw,2)
1249 : cwavef_x(2,(iband-1)*npw_k+ipw)=cwavefb(2,(iband-1)*npw_k+ipw,1)+cwavefb(2,(iband-1)*npw_k+ipw,2)
1250 : !$(\Psi^{1}-i \Psi^{2})$
1251 : cwavef_y(1,(iband-1)*npw_k+ipw)=cwavefb(1,(iband-1)*npw_k+ipw,1)+cwavefb(2,(iband-1)*npw_k+ipw,2)
1252 : cwavef_y(2,(iband-1)*npw_k+ipw)=cwavefb(2,(iband-1)*npw_k+ipw,1)-cwavefb(1,(iband-1)*npw_k+ipw,2)
1253 : end do
1254 : end do
1255 : #endif
1256 : else
1257 5978296 : cwavef_x(:,:)=cwavefb(:,1:npw_k*blocksize,1)+cwavefb(:,:,2)
1258 2000328 : cwavef_y(1,:)=cwavefb(1,1:npw_k*blocksize,1)+cwavefb(2,:,2)
1259 2000328 : cwavef_y(2,:)=cwavefb(2,:,1)-cwavefb(1,:,2)
1260 : end if
1261 11344 : if (nspinor1TreatedByThisProc) then
1262 : call prep_fourwf(rhoaug(:,:,:,4),blocksize,cwavefb(:,:,2),wfraug,&
1263 : iblock,istwf_k,gs_hamk%mgfft,mpi_enreg,nband_k,ndat,gs_hamk%ngfft,&
1264 : npw_k,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,occ_k,1,gs_hamk%ucvol,wtk,gs_hamk%nfft_blocks,&
1265 7872 : gpu_option=dtset%gpu_option)
1266 : end if
1267 11344 : if (nspinor2TreatedByThisProc) then
1268 : call prep_fourwf(rhoaug(:,:,:,2),blocksize,cwavef_x,wfraug,&
1269 : iblock,istwf_k,gs_hamk%mgfft,mpi_enreg,nband_k,ndat,gs_hamk%ngfft,&
1270 : npw_k,gs_hamk%n4,gs_hamk%n5,gs_hamk%n6,occ_k,1,gs_hamk%ucvol,wtk,gs_hamk%nfft_blocks,&
1271 7872 : gpu_option=dtset%gpu_option)
1272 : call prep_fourwf(rhoaug(:,:,:,3),blocksize,cwavef_y,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 : end if
1277 : #ifdef HAVE_OPENMP_OFFLOAD
1278 : !$OMP TARGET EXIT DATA MAP(delete:cwavef_x,cwavef_y) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
1279 : #endif
1280 11344 : ABI_FREE(cwavef_x)
1281 11344 : ABI_FREE(cwavef_y)
1282 : end if
1283 20720 : call timab(537,2,tsec)
1284 : #ifdef HAVE_OPENMP_OFFLOAD
1285 : !$OMP TARGET EXIT DATA MAP(delete:cwavefb) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
1286 : #endif
1287 20720 : ABI_FREE(cwavefb)
1288 : end if
1289 : end if
1290 : ABI_NVTX_END_RANGE()
1291 : end if ! End of SCF calculation
1292 :
1293 : ! Call to nonlocal operator:
1294 : ! - Compute nonlocal forces from most recent wfs
1295 : ! - PAW: compute projections of WF onto NL projectors (cprj)
1296 : ABI_NVTX_START_RANGE(NVTX_VTOWFK_NONLOP)
1297 2150541 : eig_k_block => eig_k(1+(iblock-1)*blocksize:iblock*blocksize)
1298 2150541 : cg_k_block => cg_k(:,1+(iblock-1)*blocksize*my_nspinor*npw_k:iblock*blocksize*my_nspinor*npw_k)
1299 2493832 : if (dtset%cprj_in_memory==2) then
1300 0 : if (optforces>0) then
1301 0 : call timab(554,1,tsec) ! "vtowfk:rhoij"
1302 : ! Treat all wavefunctions in case of PAW
1303 0 : cwaveprj => cprj(:,1+(iblock-1)*my_nspinor*blocksize+ibg:iblock*my_nspinor*blocksize+ibg)
1304 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,eig_k_block,&
1305 0 : mpi_enreg,blocksize,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,nonlop_dum2)
1306 : ! Accumulate forces
1307 0 : iband=(iblock-1)*blocksize
1308 0 : do iblocksize=1,blocksize
1309 0 : ii=0
1310 0 : if (nnlout>3*natom) ii=6
1311 0 : iband=iband+1;ibs=ii+nnlout*(iblocksize-1)
1312 0 : grnl_k(1:nnlout,iband)=enlout(ibs+1:ibs+nnlout)
1313 : end do
1314 0 : call timab(554,2,tsec) ! "vtowfk:rhoij"
1315 : end if ! PAW or forces
1316 : else
1317 2150541 : if(iscf>0.or.gs_hamk%usecprj==1)then
1318 2033702 : if (gs_hamk%usepaw==1.or.optforces/=0) then
1319 : ! Treat all wavefunctions in case of varying occupation numbers or PAW
1320 : ! Only treat occupied bands in case of fixed occupation numbers and NCPP
1321 1201196 : if(fixed_occ.and.abs(occblock)<=tol8.and.gs_hamk%usepaw==0) then
1322 242199 : if (optforces>0) grnl_k(:,(iblock-1)*blocksize+1:iblock*blocksize)=zero
1323 : else
1324 1188251 : if (dtset%cprj_in_memory/=1) then
1325 1163499 : if(gs_hamk%usepaw==1) then
1326 968254 : call timab(554,1,tsec) ! "vtowfk:rhoij"
1327 : end if
1328 1163499 : if(cpopt==1) then
1329 916 : iband=1+(iblock-1)*bandpp_cprj
1330 916 : call pawcprj_copy(cprj(:,1+(iblock-1)*my_nspinor*blocksize+ibg:iblock*my_nspinor*blocksize+ibg),cwaveprj)
1331 : end if
1332 1163499 : if (mpi_enreg%paral_kgb==1) then
1333 85604 : call timab(572,1,tsec) ! 'prep_nonlop%vtowfk'
1334 : call prep_nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir, &
1335 : eig_k_block,blocksize,&
1336 85604 : mpi_enreg,nnlout,paw_opt,signs,nonlop_dum,tim_nonlop_prep,cwavef,nonlop_dum2,already_transposed=.false.)
1337 85604 : call timab(572,2,tsec)
1338 : else
1339 : call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,idir,eig_k_block,&
1340 : mpi_enreg,blocksize,nnlout,&
1341 1077895 : paw_opt,signs,nonlop_dum,tim_nonlop,cwavef,cwavef)
1342 : end if
1343 1163499 : if(gs_hamk%usepaw==1) then
1344 968254 : call timab(554,2,tsec)
1345 : end if
1346 : ! Acccumulate forces
1347 1163499 : if (optforces>0) then
1348 323191 : iband=(iblock-1)*blocksize
1349 766844 : do iblocksize=1,blocksize
1350 443653 : ii=0
1351 443653 : if (nnlout>3*natom) ii=6
1352 443653 : iband=iband+1;ibs=ii+nnlout*(iblocksize-1)
1353 5098886 : grnl_k(1:nnlout,iband)=enlout(ibs+1:ibs+nnlout)
1354 : end do
1355 : end if
1356 : ! Store cprj (<Pnl|Psi>)
1357 1163499 : if (gs_hamk%usepaw==1.and.gs_hamk%usecprj==1) then
1358 968254 : iband=1+(iblock-1)*bandpp_cprj
1359 : call pawcprj_put(gs_hamk%atindx,cwaveprj,cprj,natom,iband,ibg,ikpt,iorder_cprj,isppol,&
1360 : mband_cprj,dtset%mkmem,natom,bandpp_cprj,nband_k_cprj,gs_hamk%dimcprj,my_nspinor,&
1361 968254 : dtset%nsppol,dtfil%unpaw,mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1362 : end if
1363 :
1364 : else ! cprj_in_memory==1
1365 :
1366 24752 : call timab(222,1,tsec) ! 'nonlop%vtowfk'
1367 :
1368 24752 : if ( gs_hamk%istwf_k > 1 ) then ! Real only
1369 10154 : space = SPACE_CR
1370 : else ! complex
1371 14598 : space = SPACE_C
1372 : end if
1373 24752 : me_g0 = -1
1374 24752 : if (space==SPACE_CR) then
1375 10154 : me_g0 = 0
1376 10154 : if (gs_hamk%istwf_k == 2) then
1377 1544 : if (mpi_enreg%me_g0 == 1) me_g0 = 1
1378 : end if
1379 : end if
1380 : call xgBlock_map(xgx0,cg_k_block,space,npw_k*my_nspinor,blocksize,comm=mpi_enreg%comm_band,me_g0=me_g0,&
1381 24752 : gpu_option=dtset%gpu_option)
1382 24752 : call xgBlock_map_1d(xgeigen,eig_k_block,SPACE_R,blocksize)
1383 :
1384 24752 : if (optforces/=0.or.gs_hamk%usepaw==1) then
1385 24752 : call xg_nonlop_getcprj(xg_nonlop,xgx0,cprj_xgx0%self,cprj_work%self)
1386 : end if
1387 :
1388 24752 : if (optforces/=0) then
1389 17752 : grnl_k_block => grnl_k(:,1+(iblock-1)*blocksize:iblock*blocksize)
1390 17752 : call xgBlock_map(xgforces,grnl_k_block,SPACE_R,3*natom,blocksize)
1391 17752 : call xg_nonlop_forces_stress(xg_nonlop,xgx0,cprj_xgx0%self,cprj_work%self,xgeigen,forces=xgforces)
1392 : end if
1393 :
1394 24752 : call timab(222,2,tsec) ! 'nonlop%vtowfk'
1395 :
1396 24752 : if (gs_hamk%usepaw==1) then
1397 23072 : cprj_cwavef_bands => cprj(:,1+ibg+(iblock-1)*ncols_cprj:iblock*ncols_cprj+ibg)
1398 23072 : call xg_cprj_copy(cprj_cwavef_bands,cprj_xgx0%self,xg_nonlop,XG_TO_CPRJ)
1399 : end if
1400 :
1401 : end if
1402 : end if
1403 : end if ! PAW or forces
1404 : end if ! iscf>0 or iscf=-3
1405 : end if
1406 : ABI_NVTX_END_RANGE()
1407 : end do ! End of loop on blocks
1408 :
1409 : #ifdef HAVE_OPENMP_OFFLOAD
1410 : !$OMP TARGET EXIT DATA MAP(delete:cwavef) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
1411 : !$OMP TARGET EXIT DATA MAP(delete:kinpw) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP)
1412 : #endif
1413 :
1414 : ! restore safe value related to GEMM nonlop slicing and GPU in case of forces compute
1415 343291 : if(optforces==1 .and. gpu_option_tmp==ABI_GPU_OPENMP) then
1416 0 : gs_hamk%gpu_option = gpu_option_tmp
1417 0 : gemm_nonlop_block_size = blksize_gemm_nonlop_tmp
1418 0 : gemm_nonlop_is_distributed = is_distrib_tmp
1419 0 : gs_hamk%nfft_blocks = nfft_blocks_tmp
1420 : end if
1421 :
1422 343291 : if (dtset%cprj_in_memory==1) then
1423 17084 : call xg_free(cprj_xgx0)
1424 17084 : call xg_free(cprj_work)
1425 : end if
1426 :
1427 : !call cwtime_report(" Block loop", cpu, wall, gflops)
1428 :
1429 343291 : if(dtset%gpu_option==ABI_GPU_KOKKOS) then
1430 : #if defined HAVE_GPU && defined HAVE_YAKL
1431 : ABI_FREE_MANAGED(cwavef)
1432 : #endif
1433 : else
1434 343291 : ABI_FREE(cwavef)
1435 : end if
1436 :
1437 343291 : ABI_FREE(enlout)
1438 :
1439 343291 : if (dtset%cprj_in_memory/=2) then
1440 343291 : if (gs_hamk%usepaw==1.and.(iscf>0.or.gs_hamk%usecprj==1)) then
1441 107704 : call pawcprj_free(cwaveprj)
1442 : end if
1443 343291 : ABI_FREE(cwaveprj)
1444 : else
1445 0 : nullify(cwaveprj)
1446 : end if
1447 :
1448 343291 : if (fixed_occ.and.iscf>0) then
1449 229430 : if(dtset%gpu_option==ABI_GPU_KOKKOS) then
1450 : #if defined HAVE_GPU && defined HAVE_YAKL
1451 : ABI_FREE_MANAGED(wfraug)
1452 : #endif
1453 : else
1454 229430 : ABI_FREE(wfraug)
1455 : end if
1456 : end if
1457 :
1458 : !Write the number of one-way 3D ffts skipped until now (in case of fixed occupation numbers
1459 343291 : if(iscf>0 .and. fixed_occ .and. (prtvol>2 .or. ikpt<=nkpt_max) )then
1460 2090 : write(msg,'(a,i0)')' vtowfk: number of one-way 3D ffts skipped in vtowfk until now =',nskip
1461 2090 : call wrtout(std_out,msg,'PERS')
1462 : end if
1463 :
1464 : ! Norm-conserving or FockACE: Compute nonlocal+FockACE part of total energy: rotate subvnlx elements
1465 : ! Note the two calls. For (old) lobpcgwf we have a (nband_k, nband_k) matrix, whereas cgwf
1466 : ! returns results in packed form.
1467 : ! CHEBYSHEV, NEW LOBPCG and RMM-DIIS do not need this
1468 : !
1469 343291 : rotate_subvnlx = gs_hamk%usepaw == 0 .and. wfopta10 /= 1 .and. .not. xg_diago
1470 343291 : if (use_rmm_diis) rotate_subvnlx = .False.
1471 :
1472 342743 : if (rotate_subvnlx) then
1473 223951 : call timab(586,1,tsec) ! 'vtowfk(nonlocalpart)'
1474 223951 : if (wfopta10==4) then
1475 612 : call cg_hrotate_and_get_diag(istwf_k, nband_k, totvnlx, evec, enlx_k)
1476 : else
1477 223339 : call cg_hprotate_and_get_diag(nband_k, subvnlx, evec, enlx_k)
1478 : end if
1479 223951 : call timab(586,2,tsec)
1480 : end if
1481 :
1482 : !###################################################################
1483 :
1484 343291 : if (iscf<=0 .and. max_resid > dtset%tolwfr) then
1485 : write(msg,'(2(a,i0),a,es13.5)')&
1486 710 : "Wavefunctions not converged for ikpt: ", ikpt, ", nnsclo: ",nnsclo_now,', max resid: ',max_resid
1487 710 : ABI_WARNING(msg)
1488 : end if
1489 :
1490 : ! Print out eigenvalues (hartree)
1491 343291 : if (mod(dtset%wfoptalg,10)==1) then
1492 21504 : niter=dtset%mdeg_filter
1493 21504 : iter_name='as the degree of the polynomial filter'
1494 : else
1495 321787 : niter=dtset%nline
1496 321787 : iter_name='CG line minimizations'
1497 : end if
1498 343291 : if (prtvol/=5.and.(prtvol>2 .or. ikpt<=nkpt_max)) then
1499 : write(msg, '(5x,a,i5,2x,a,a,a,i4,a,i4,2a)' ) &
1500 14055 : 'eigenvalues (hartree) for',nband_k,'bands',ch10,&
1501 28110 : ' after ',inonsc,' non-SCF iterations with ',niter,' ',trim(iter_name)
1502 14055 : call wrtout(std_out,msg,'PERS')
1503 56224 : do ii=0,(nband_k-1)/6
1504 42169 : write(msg, '(1p,6e12.4)' ) (eig_k(iband),iband=1+6*ii,min(6+6*ii,nband_k))
1505 56224 : call wrtout(std_out,msg,'PERS')
1506 : end do
1507 329236 : else if(ikpt==nkpt_max+1)then
1508 0 : call wrtout(std_out,' vtowfk : prtvol=0 or 1, do not print more k-points.','PERS')
1509 : end if
1510 :
1511 : ! Print out decomposition of eigenvalues in the non-selfconsistent case or if prtvol>=10
1512 343291 : if( (iscf<0 .and. (prtvol>2 .or. ikpt<=nkpt_max)) .or. prtvol>=10)then
1513 : write(msg, '(5x,a,i5,2x,a,a,a,i4,a,i4,2a)' ) &
1514 4449 : ' mean kinetic energy (hartree) for ',nband_k,' bands',ch10,&
1515 8898 : ' after ',inonsc,' non-SCF iterations with ',niter,' ',trim(iter_name)
1516 4449 : call wrtout(std_out,msg,'PERS')
1517 :
1518 13012 : do ii=0,(nband_k-1)/6
1519 8563 : write(msg, '(1p,6e12.4)' ) (ek_k(iband),iband=1+6*ii,min(6+6*ii,nband_k))
1520 13012 : call wrtout(std_out,msg,'PERS')
1521 : end do
1522 :
1523 4449 : if (gs_hamk%usepaw==0) then
1524 : write(msg, '(5x,a,i5,2x,a,a,a,i4,a,i4,2a)' ) &
1525 3076 : ' mean NL+Fock-type energy (hartree) for ',nband_k,' bands',ch10,&
1526 6152 : ' after ',inonsc,' non-SCF iterations with ',niter,' ',trim(iter_name)
1527 3076 : call wrtout(std_out,msg,'PERS')
1528 :
1529 8029 : do ii=0,(nband_k-1)/6
1530 4953 : write(msg,'(1p,6e12.4)') (enlx_k(iband),iband=1+6*ii,min(6+6*ii,nband_k))
1531 8029 : call wrtout(std_out,msg,'PERS')
1532 : end do
1533 : end if
1534 : end if
1535 :
1536 : ! Hamiltonian constructor for gwls_sternheimer
1537 343291 : if (dtset%optdriver==RUNL_GWLS) call build_H(dtset,mpi_enreg,cpopt,cg,gs_hamk,kg_k,kinpw)
1538 :
1539 343291 : if (dtset%cprj_in_memory==2) nullify(cprj_cwavef_bands)
1540 :
1541 : #ifdef HAVE_OPENMP_OFFLOAD
1542 : if(transfer_cg) then
1543 : !$OMP TARGET EXIT DATA MAP(from:cg_k) IF(gs_hamk%gpu_option==ABI_GPU_OPENMP .and. xg_diago)
1544 : end if
1545 : #endif
1546 343291 : if(wfopta10 /= 1 .and. .not. xg_diago) then
1547 298640 : ABI_FREE(evec)
1548 298640 : ABI_FREE(subham)
1549 298640 : ABI_FREE(totvnlx)
1550 298640 : ABI_FREE(subvnlx)
1551 298640 : ABI_FREE(subovl)
1552 : end if
1553 :
1554 343291 : ABI_SFREE(gsc)
1555 :
1556 343291 : if(wfoptalg==3) then
1557 6 : ABI_FREE(eig_save)
1558 : end if
1559 :
1560 343291 : if (prtvol==-level) then
1561 : ! Structured debugging: if prtvol=-level, stop here.
1562 0 : write(msg,'(3a,i0,a)')' vtowfk: exit ',ch10,' prtvol=-',level,', debugging mode => stop '
1563 0 : ABI_ERROR(msg)
1564 : end if
1565 :
1566 343291 : call timab(30,2,tsec)
1567 343291 : call timab(28,2,tsec)
1568 :
1569 : DBG_EXIT("COLL")
1570 :
1571 686582 : end subroutine vtowfk
1572 : !!***
1573 :
1574 : end module m_vtowfk
1575 : !!***
|
