Line data Source code
1 : !!****m* ABINIT/m_getghc
2 : !! NAME
3 : !! m_getghc
4 : !!
5 : !! FUNCTION
6 : !! Compute <G|H|C> for input vector |C> expressed in reciprocal space;
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, LSI, MT, JB, JWZ)
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_getghc
26 :
27 : use, intrinsic :: iso_c_binding
28 : use defs_basis
29 : use m_errors
30 : use m_abicore
31 : use m_xmpi
32 : use m_xomp
33 : use m_gputk
34 : use m_abi_linalg
35 :
36 : use defs_abitypes, only : mpi_type
37 : use m_time, only : timab
38 : use m_fstrings, only : sjoin, itoa
39 : use m_cgtools, only : cg_copy_spin, cg_put_spin
40 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_getdim, pawcprj_copy
41 : use m_bandfft_kpt, only : bandfft_kpt, bandfft_kpt_get_ikpt
42 : use m_hamiltonian, only : gs_hamiltonian_type, KPRIME_H_K, K_H_KPRIME, K_H_K, KPRIME_H_KPRIME
43 : use m_fock, only : fock_common_type, fock_get_getghc_call
44 : use m_fock_getghc, only : fock_getghc, fock_ACE_getghc
45 : use m_nonlop, only : nonlop
46 : use m_gemm_nonlop_projectors, only : gemm_nonlop_use_gemm
47 : use m_fft, only : fourwf,fourwf_optmem
48 : use m_ompgpu_fourwf, only : ompgpu_fourwf_work_mem
49 : use m_gemm_nonlop, only : gemm_nonlop_ompgpu_work_mem
50 :
51 : #if defined(HAVE_GPU_MARKERS)
52 : use m_nvtx_data
53 : #endif
54 : #ifdef HAVE_FFTW3_THREADS
55 : use m_fftw3, only : fftw3_spawn_threads_here, fftw3_use_lib_threads
56 : #endif
57 : #if defined HAVE_GPU_CUDA
58 : use m_gpu_toolbox
59 : #endif
60 : #if defined HAVE_YAKL
61 : use gator_mod
62 : #endif
63 : #ifdef HAVE_KOKKOS
64 : use m_manage_kokkos, only : assemble_energy_contribution_kokkos
65 : #endif
66 :
67 : implicit none
68 :
69 : private
70 : !!***
71 :
72 : public :: getghc ! Compute <G|H|C> for input vector |C> expressed in reciprocal space
73 : public :: getgsc ! Compute <G|S|C> for all input vectors |Cnk> at a given k-point
74 : public :: getghc_mGGA
75 : public :: multithreaded_getghc
76 : public :: getghc_nucdip ! compute <G|H_nucdip|C> for input vector |C> expressed in recip space
77 : public :: getghc_ompgpu_work_mem ! assess GPU memory requirements for running getghc with OpenMP GPU
78 : !!***
79 :
80 : contains
81 : !!***
82 :
83 : !!****f* ABINIT/getghc_ompgpu_work_mem
84 : !! NAME
85 : !! getghc_ompgpu_work_mem
86 : !!
87 : !! FUNCTION
88 : !! Returns work memory requirement for getghc_ompgpu
89 : !!
90 : !! INPUTS
91 : !! gs_ham <type(gs_hamiltonian_type)>=contains dimensions of FFT domain
92 : !! ndat=size of batch for fourwf and nonlop processing
93 : !!
94 : !! OUTPUT
95 : !! req_mem=amount in bytes of required memory for getghc_ompgpu
96 :
97 0 : function getghc_ompgpu_work_mem(gs_ham, ndat, nfft_blocks) result(req_mem)
98 :
99 : type(gs_hamiltonian_type),intent(in),target :: gs_ham
100 : integer, intent(in) :: ndat, nfft_blocks
101 : integer(kind=c_size_t) :: req_mem, ghc_mem, nonlop_mem
102 :
103 : ! getghc use a GPU work buffer only when using fourwf
104 : ! Therefore, max GPU memory required by getghc is either:
105 : ! - the sum of getghc and fourwf work buffers memory requirements
106 : ! - the amount of memory required by gemm_nonlop_ompgpu work buffers
107 0 : ghc_mem = 0
108 0 : ghc_mem = int(2, c_size_t) * dp * gs_ham%n4 * gs_ham%n5 * gs_ham%n6 * (ndat/nfft_blocks + modulo(ndat,nfft_blocks))
109 : !ghc_mem = ghc_mem + ompgpu_fourwf_work_mem(gs_ham%ngfft, ndat)
110 :
111 : nonlop_mem = gemm_nonlop_ompgpu_work_mem(gs_ham%istwf_k, ndat, 0, gs_ham%npw_fft_k,&
112 0 : & gs_ham%indlmn, gs_ham%nattyp, gs_ham%ntypat, gs_ham%lmnmax, 2, 11)
113 :
114 0 : req_mem = MAX(ghc_mem, nonlop_mem)
115 :
116 0 : end function getghc_ompgpu_work_mem
117 : !!***
118 :
119 : !!****f* ABINIT/getghc
120 : !! NAME
121 : !! getghc
122 : !!
123 : !! FUNCTION
124 : !! Compute <G|H|C> for input vector |C> expressed in reciprocal space;
125 : !! Result is put in array ghc.
126 : !! <G|Vnonlocal + VfockACE|C> is also returned in gvnlxc if either NLoc NCPP or FockACE.
127 : !! If required, <G|S|C> is returned in gsc (S=overlap - PAW only)
128 : !! Note that left and right k points can be different, i.e. ghc=<k^prime+G|H|C_k>.
129 : !!
130 : !! INPUTS
131 : !! cpopt=flag defining the status of cwaveprj%cp(:)=<Proj_i|Cnk> scalars (PAW only)
132 : !! (same meaning as in nonlop.F90 routine)
133 : !! if cpopt=-1, <p_lmn|in> (and derivatives) are computed here (and not saved)
134 : !! if cpopt= 0, <p_lmn|in> are computed here and saved
135 : !! if cpopt= 1, <p_lmn|in> and first derivatives are computed here and saved
136 : !! if cpopt= 2 <p_lmn|in> are already in memory;
137 : !! if cpopt= 3 <p_lmn|in> are already in memory; first derivatives are computed here and saved
138 : !! if cpopt= 4 <p_lmn|in> and first derivatives are already in memory;
139 : !! cwavef(2,npw*my_nspinor*ndat)=planewave coefficients of wavefunction.
140 : !! cwavef_r(2,n4,n5,n6,nspinor) = wave function in real space
141 : !! gs_ham <type(gs_hamiltonian_type)>=all data for the Hamiltonian to be applied
142 : !! lambda=factor to be used when computing <G|H-lambda.S|C> - only for sij_opt=-1
143 : !! Typically lambda is the eigenvalue (or its guess)
144 : !! mpi_enreg=information about MPI parallelization
145 : !! ndat=number of FFT to do in parallel
146 : !! prtvol=control print volume and debugging output
147 : !! sij_opt= -PAW ONLY- if 0, only matrix elements <G|H|C> have to be computed
148 : !! (S=overlap) if 1, matrix elements <G|S|C> have to be computed in gsc in addition to ghc
149 : !! if -1, matrix elements <G|H-lambda.S|C> have to be computed in ghc (gsc not used)
150 : !! tim_getghc=timing code of the calling subroutine(can be set to 0 if not attributed)
151 : !! type_calc= option governing which part of Hamitonian is to be applied:
152 : ! 0: whole Hamiltonian
153 : !! 1: local part only
154 : !! 2: non-local+Fock+kinetic only (added to the existing Hamiltonian)
155 : !! 3: local + kinetic only (added to the existing Hamiltonian)
156 : !! ===== Optional inputs =====
157 : !! [kg_fft_k(3,:)]=optional, (k+G) vector coordinates to be used for the FFT transformation
158 : !! instead of the one contained in gs_ham datastructure.
159 : !! Typically used for real WF (in parallel) which are FFT-transformed 2 by 2.
160 : !! [kg_fft_kp(3,:)]=optional, (k^prime+G) vector coordinates to be used for the FFT transformation
161 : !! [select_k]=optional, option governing the choice of k points to be used.
162 : !! gs_ham datastructure contains quantities needed to apply Hamiltonian
163 : !! in reciprocal space between 2 kpoints, k and k^prime (equal in most cases);
164 : !! if select_k=1, <k^prime|H|k> is applied [default]
165 : !! if select_k=2, <k|H|k^prime> is applied
166 : !! if select_k=3, <k|H|k> is applied
167 : !! if select_k=4, <k^prime|H|k^prime> is applied
168 : !!
169 : !! OUTPUT
170 : !! ghc(2,npw*my_nspinor*ndat)=matrix elements <G|H|C> (if sij_opt>=0)
171 : !! or <G|H-lambda.S|C> (if sij_opt=-1)
172 : !! gvnlxc(2,npw*my_nspinor*ndat)=matrix elements <G|Vnonlocal+VFockACE|C> (if sij_opt>=0)
173 : !! or <G|Vnonlocal+VFockACE-lambda.S|C> (if sij_opt=-1)
174 : !! include Vnonlocal if NCPP and non-local Fock if associated(gs_ham%fockcommon)
175 : !! if (sij_opt=1)
176 : !! gsc(2,npw*my_nspinor*ndat)=matrix elements <G|S|C> (S=overlap).
177 : !!
178 : !! SIDE EFFECTS
179 : !! cwaveprj(natom,my_nspinor*(1+cpopt)*ndat)= wave function projected on nl projectors (PAW only)
180 : !!
181 : !! SOURCE
182 :
183 34965650 : subroutine getghc(cpopt,cwavef,cwaveprj,ghc,gsc,gs_ham,gvnlxc,lambda,mpi_enreg,ndat,&
184 : prtvol,sij_opt,tim_getghc,type_calc,&
185 34965650 : kg_fft_k,kg_fft_kp,select_k,cwavef_r,filter_dilatmx_loc) ! optional arguments
186 :
187 : !Arguments ------------------------------------
188 : !scalars
189 : logical,intent(in),optional :: filter_dilatmx_loc
190 : integer,intent(in) :: cpopt,ndat, prtvol
191 : integer,intent(in) :: sij_opt,tim_getghc,type_calc
192 : integer,intent(in),optional :: select_k
193 : real(dp),intent(in) :: lambda
194 : type(MPI_type),intent(in) :: mpi_enreg
195 : type(gs_hamiltonian_type),intent(inout),target :: gs_ham
196 : !arrays
197 : integer,intent(in),optional,target :: kg_fft_k(:,:),kg_fft_kp(:,:)
198 : real(dp),intent(out),target :: gsc(:,:)
199 : real(dp),intent(inout), target :: cwavef(:,:)
200 : real(dp),optional,intent(inout) :: cwavef_r(:,:,:,:,:)
201 : real(dp),intent(out), target :: ghc(:,:)
202 : real(dp),intent(out),target :: gvnlxc(:,:)
203 : type(pawcprj_type),intent(inout),target :: cwaveprj(:,:)
204 : !MG: Passing these arrays assumed-shape has the drawback that client code is
205 : !forced to use vec(2, npw*ndat) instead of the more user-friendly vec(2,npw,ndat)
206 :
207 : !Local variables-------------------------------
208 : !scalars
209 : integer,parameter :: level=114, tim_fourwf=1
210 : integer :: choice,cplex,cpopt_here,fftalg,i1,i2,i3,idat,idir,ierr,i0
211 : integer :: ig,igspinor,istwf_k_,ii,iispinor,ikpt_this_proc,ipw,ispinor,my_nspinor
212 : integer :: n4,n5,n6,ndat_,nnlout,npw_fft,npw_k1,npw_k2,nspinortot,option_fft
213 : integer :: paw_opt,select_k_,shift1,shift2,signs,tim_nonlop
214 : integer :: chunk,residuchunk,nfft_blocks
215 : logical(kind=c_bool) :: k1_eq_k2
216 : logical :: double_rfft_trick,have_to_reequilibrate,has_fock,local_gvnlxc
217 : logical :: nspinor1TreatedByThisProc,nspinor2TreatedByThisProc,use_cwavef_r, filter_dilatmx_loc_
218 : real(dp) :: ghcim,ghcre,weight !, kscale
219 : #ifdef HAVE_OPENMP_OFFLOAD
220 : complex(dp), parameter :: cminusone = (-1._dp,0._dp)
221 : #endif
222 : character(len=500) :: msg
223 : !arrays
224 34965650 : integer, contiguous, pointer :: gbound_k1(:,:), gbound_k2(:,:)
225 34965650 : integer, contiguous, pointer :: kg_k1(:,:), kg_k2(:,:)
226 34965650 : integer, contiguous, pointer :: indices_pw_fft(:), kg_k_fft(:,:)
227 34965650 : integer, contiguous, pointer :: recvcount_fft(:), recvdisp_fft(:)
228 34965650 : integer, contiguous, pointer :: sendcount_fft(:), senddisp_fft(:)
229 34965650 : integer, allocatable:: dimcprj(:)
230 104896950 : real(dp) :: enlout(ndat), lambda_ndat(ndat), tsec(2)
231 : real(dp), target :: nonlop_dum(1,1)
232 34965650 : real(dp), allocatable :: buff_wf(:,:)
233 34965650 : real(dp), allocatable :: cwavef1(:,:)
234 34965650 : real(dp), allocatable :: cwavef2(:,:)
235 34965650 : real(dp), allocatable :: cwavef_fft(:,:)
236 34965650 : real(dp), allocatable :: cwavef_fft_tr(:,:)
237 34965650 : real(dp), allocatable :: ghc1(:,:)
238 34965650 : real(dp), allocatable :: ghc2(:,:)
239 34965650 : real(dp), allocatable :: ghc3(:,:)
240 34965650 : real(dp), allocatable :: ghc4(:,:)
241 34965650 : real(dp), allocatable :: ghc_mGGA(:,:)
242 34965650 : real(dp), allocatable :: ghc_mGGA_fft(:,:)
243 34965650 : real(dp), allocatable :: ghc_vectornd(:,:)
244 34965650 : real(dp), allocatable :: cwavef_spin(:,:), gvnlxc_spin(:,:)
245 :
246 : #if defined HAVE_GPU && defined HAVE_YAKL
247 : real(c_double), contiguous, pointer :: gvnlc(:,:)
248 : real(c_double), contiguous, pointer :: gvnlxc_(:,:)
249 : #else
250 34965650 : real(dp), target, allocatable :: gvnlc(:,:)
251 34965650 : real(dp), contiguous, pointer :: gvnlxc_(:,:)
252 : #endif
253 34965650 : real(dp), allocatable :: vlocal_tmp(:,:,:), work(:,:,:,:)
254 34965650 : real(dp), contiguous, pointer :: kinpw_k1(:), kinpw_k2(:), kpt_k1(:), kpt_k2(:)
255 34965650 : real(dp), contiguous, pointer :: gsc_ptr(:,:)
256 : type(fock_common_type),pointer :: fock
257 34965650 : type(pawcprj_type),pointer :: cwaveprj_fock(:,:),cwaveprj_idat(:,:),cwaveprj_nonlop(:,:)
258 : logical :: transfer_ghc,transfer_gsc,transfer_cwavef,transfer_gvnlxc
259 : real(c_double), parameter :: hugevalue = huge(zero)*1.d-11
260 : ! *********************************************************************
261 :
262 : DBG_ENTER("COLL")
263 :
264 : !Keep track of total time spent in getghc:
265 34965650 : call timab(350+tim_getghc,1,tsec)
266 : ABI_NVTX_START_RANGE(NVTX_GETGHC)
267 :
268 : !Structured debugging if prtvol==-level
269 34965650 : if(prtvol==-level)then
270 0 : write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' getghc : enter, debugging '
271 0 : call wrtout(std_out,msg,'PERS')
272 : end if
273 :
274 : !Select k-dependent objects according to select_k input parameter
275 34965650 : select_k_=KPRIME_H_K;if (present(select_k)) select_k_=select_k
276 14561182 : select case (select_k_)
277 : case (KPRIME_H_K)
278 : ! <k^prime|H|k>
279 14561182 : npw_k1 = gs_ham%npw_fft_k ; npw_k2 = gs_ham%npw_fft_kp
280 14561182 : kpt_k1 => gs_ham%kpt_k ; kpt_k2 => gs_ham%kpt_kp
281 14561182 : kg_k1 => gs_ham%kg_k ; kg_k2 => gs_ham%kg_kp
282 14561182 : gbound_k1 => gs_ham%gbound_k ; gbound_k2 => gs_ham%gbound_kp
283 14561182 : kinpw_k1 => gs_ham%kinpw_k ; kinpw_k2 => gs_ham%kinpw_kp
284 : case (K_H_KPRIME)
285 : ! <k|H|k^prime>
286 0 : npw_k1 = gs_ham%npw_fft_kp; npw_k2 = gs_ham%npw_fft_k
287 0 : kpt_k1 => gs_ham%kpt_kp ; kpt_k2 => gs_ham%kpt_k
288 0 : kg_k1 => gs_ham%kg_kp ; kg_k2 => gs_ham%kg_k
289 0 : gbound_k1 => gs_ham%gbound_kp ; gbound_k2 => gs_ham%gbound_k
290 0 : kinpw_k1 => gs_ham%kinpw_kp ; kinpw_k2 => gs_ham%kinpw_k
291 : case (K_H_K)
292 : ! <k|H|k>
293 0 : npw_k1 = gs_ham%npw_fft_k ; npw_k2 = gs_ham%npw_fft_k
294 0 : kpt_k1 => gs_ham%kpt_k ; kpt_k2 => gs_ham%kpt_k
295 0 : kg_k1 => gs_ham%kg_k ; kg_k2 => gs_ham%kg_k
296 0 : gbound_k1 => gs_ham%gbound_k ; gbound_k2 => gs_ham%gbound_k
297 0 : kinpw_k1 => gs_ham%kinpw_k ; kinpw_k2 => gs_ham%kinpw_k
298 : case (KPRIME_H_KPRIME)
299 : ! <k^prime|H|k^prime>
300 20404468 : npw_k1 = gs_ham%npw_fft_kp; npw_k2 = gs_ham%npw_fft_kp
301 20404468 : kpt_k1 => gs_ham%kpt_kp ; kpt_k2 => gs_ham%kpt_kp
302 20404468 : kg_k1 => gs_ham%kg_kp ; kg_k2 => gs_ham%kg_kp
303 20404468 : gbound_k1 => gs_ham%gbound_kp ; gbound_k2 => gs_ham%gbound_kp
304 20404468 : kinpw_k1 => gs_ham%kinpw_kp ; kinpw_k2 => gs_ham%kinpw_kp
305 : case default
306 34965650 : ABI_ERROR(sjoin("Invalid select_k: ", itoa(select_k_)))
307 : end select
308 :
309 139485385 : k1_eq_k2=(all(abs(kpt_k1(:)-kpt_k2(:))<tol8))
310 :
311 : ! Check sizes
312 34965650 : my_nspinor=max(1,gs_ham%nspinor/mpi_enreg%nproc_spinor)
313 104896950 : ABI_CHECK_IGEQ(size(cwavef), 2*npw_k1*my_nspinor*ndat, 'wrong size for cwavef!')
314 104896950 : ABI_CHECK_IGEQ(size(ghc), 2*npw_k2*my_nspinor*ndat, 'wrong size for ghc!')
315 34965650 : if (sij_opt==1) then
316 24460182 : ABI_CHECK_IGEQ(size(gsc), 2*npw_k2*my_nspinor*ndat, 'wrong size for gsc!')
317 : end if
318 34965650 : if (gs_ham%usepaw==1.and.cpopt>=0) then
319 7688310 : ABI_CHECK_IGEQ(size(cwaveprj), gs_ham%natom*my_nspinor*ndat, 'wrong size for cwaveprj!')
320 : end if
321 34965650 : if (any(type_calc == [0, 2, 3])) then
322 104273178 : local_gvnlxc = size(gvnlxc)<=1
323 34757726 : if (local_gvnlxc) then
324 880108 : if(gs_ham%gpu_option==ABI_GPU_KOKKOS) then
325 : #if defined HAVE_GPU && defined HAVE_YAKL
326 : ABI_MALLOC_MANAGED(gvnlxc_, (/2,npw_k2*my_nspinor*ndat/))
327 : #endif
328 : else
329 2640324 : ABI_MALLOC(gvnlxc_,(2,npw_k2*my_nspinor*ndat))
330 : end if
331 : else
332 33877618 : gvnlxc_ => gvnlxc
333 : end if
334 104273178 : ABI_CHECK_IGEQ(size(gvnlxc_), 2*npw_k2*my_nspinor*ndat, 'wrong size for gvnlxc!')
335 : end if
336 :
337 34965650 : use_cwavef_r=present(cwavef_r)
338 34965650 : n4=gs_ham%n4 ; n5=gs_ham%n5 ; n6=gs_ham%n6
339 34965650 : nspinortot=gs_ham%nspinor
340 :
341 34965650 : if (use_cwavef_r) then
342 0 : ABI_CHECK_IEQ(size(cwavef_r,1), 2, 'wrong size for cwavef_r (dimension 1)')
343 0 : ABI_CHECK_IEQ(size(cwavef_r,2), n4, 'wrong size for cwavef_r (dimension 2)')
344 0 : ABI_CHECK_IEQ(size(cwavef_r,3), n5, 'wrong size for cwavef_r (dimension 3)')
345 0 : ABI_CHECK_IEQ(size(cwavef_r,4), n6, 'wrong size for cwavef_r (dimension 4)')
346 0 : ABI_CHECK_IEQ(size(cwavef_r,5), nspinortot, 'wrong size for cwavef_r (dimension 5)')
347 : end if
348 :
349 : !Eventually overwrite plane waves data for FFT
350 34965650 : if (present(kg_fft_k)) then
351 3274 : kg_k1 => kg_fft_k ; kg_k2 => kg_fft_k
352 3274 : npw_k1=size(kg_k1,2) ; npw_k2=size(kg_k2,2)
353 : end if
354 34965650 : if (present(kg_fft_kp)) then
355 0 : kg_k2 => kg_fft_kp ; npw_k2=size(kg_k2,2)
356 : end if
357 :
358 : !paral_kgb constraint
359 34965650 : if (mpi_enreg%paral_kgb==1.and.(.not.k1_eq_k2)) then
360 0 : ABI_BUG('paral_kgb=1 not allowed for k/=k_^prime!')
361 : end if
362 :
363 : !Do we add Fock exchange term ?
364 34965650 : has_fock = associated(gs_ham%fockcommon)
365 34965650 : if (has_fock) fock => gs_ham%fockcommon
366 :
367 : !Parallelization over spinors management
368 34965650 : if (mpi_enreg%paral_spinor==0) then
369 34914286 : shift1=npw_k1;shift2=npw_k2
370 34914286 : nspinor1TreatedByThisProc=.true.
371 34914286 : nspinor2TreatedByThisProc=(nspinortot==2)
372 : else
373 51364 : shift1=0;shift2=0
374 51364 : nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
375 51364 : nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
376 : end if
377 :
378 34965650 : filter_dilatmx_loc_ = .true.; if ( present(filter_dilatmx_loc) ) filter_dilatmx_loc_ = filter_dilatmx_loc
379 34965650 : transfer_ghc = .false.; transfer_gsc = .false.; transfer_gvnlxc = .false.; transfer_cwavef = .false.
380 : if(gs_ham%gpu_option == ABI_GPU_OPENMP) then
381 : #ifdef HAVE_OPENMP_OFFLOAD
382 : transfer_ghc = .not. xomp_target_is_present(c_loc(ghc))
383 : transfer_gsc = .not. xomp_target_is_present(c_loc(gsc)) .and. sij_opt==1
384 : transfer_gvnlxc = .not. xomp_target_is_present(c_loc(gvnlxc_))
385 : transfer_cwavef = .not. xomp_target_is_present(c_loc(cwavef))
386 :
387 : !$OMP TARGET ENTER DATA MAP(alloc:ghc) IF(transfer_ghc)
388 : !$OMP TARGET ENTER DATA MAP(alloc:gsc) IF(transfer_gsc)
389 : !$OMP TARGET ENTER DATA MAP(alloc:gvnlxc_) IF(transfer_gvnlxc)
390 : !$OMP TARGET ENTER DATA MAP(to:cwavef) IF(transfer_cwavef)
391 : if (type_calc == 2) then
392 : !$OMP TARGET UPDATE TO(ghc) IF(transfer_ghc)
393 : !$OMP TARGET UPDATE TO(gsc) IF(transfer_gsc .and. gs_ham%usepaw==1)
394 : end if
395 : #endif
396 : end if
397 :
398 : !============================================================
399 : ! Application of the local potential
400 : !============================================================
401 :
402 : ABI_NVTX_START_RANGE(NVTX_GETGHC_LOCPOT)
403 :
404 34965650 : if (any(type_calc == [0, 1, 3])) then
405 :
406 : ! Need a Vlocal
407 34962376 : ABI_CHECK(associated(gs_ham%vlocal), "We need vlocal in gs_ham!")
408 :
409 : ! fourwf can only process with one value of istwf_k
410 34962376 : if (gs_ham%use_gbt == 0) then
411 34823257 : ABI_CHECK(k1_eq_k2, 'vlocal (fourwf) cannot be computed with k/=k^prime!')
412 : end if
413 :
414 34962376 : nfft_blocks = 1;
415 34962376 : if(gs_ham%nfft_blocks > 1 .and. gs_ham%gpu_option/=ABI_GPU_DISABLED) then
416 0 : nfft_blocks = gs_ham%nfft_blocks
417 : end if
418 34962376 : chunk = ndat/nfft_blocks ! Divide by 2 to construct chunk of even number of bands
419 34962376 : residuchunk = ndat - nfft_blocks*chunk
420 :
421 : ! Eventually adjust load balancing for FFT (by changing FFT distrib)
422 : ! Not that have_to_reequilibrate can be true only if mpi_enreg%nproc_fft>1
423 34962376 : have_to_reequilibrate=.false.
424 34962376 : if (mpi_enreg%paral_kgb==1) then
425 919678 : ikpt_this_proc=bandfft_kpt_get_ikpt()
426 919678 : have_to_reequilibrate=bandfft_kpt(ikpt_this_proc)%have_to_reequilibrate
427 : end if
428 34962376 : ndat_ = ndat
429 34962376 : istwf_k_ = gs_ham%istwf_k
430 34962376 : double_rfft_trick = istwf_k_==2.and.ndat>1.and.mpi_enreg%paral_kgb==1.and.gs_ham%gpu_option==ABI_GPU_DISABLED
431 : ! LB-08-2024 : double_rfft_trick works only for paral_kgb=1, but I don't know why...
432 : ! Note that the trick can be activated only if nspinortot=1 (if =2 then istwf_k=1), so gs_ham%nvloc=1 too
433 : if (double_rfft_trick) then
434 16752 : ABI_CHECK(mpi_enreg%nproc_fft == 1, 'double_rfft_trick inside getghc is not implemented for npfft>1')
435 16752 : istwf_k_ = 1
436 16752 : ndat_ = ndat / 2
437 16752 : if (modulo(ndat,2)/=0) ndat_=ndat_+1
438 16752 : npw_fft=2*npw_k1
439 16752 : i0=1
440 16752 : if (mpi_enreg%me_g0_fft==1) then ! Do not include G=(0,0,0) twice
441 16752 : npw_fft=npw_fft-1
442 16752 : i0=2
443 : end if
444 50256 : ABI_MALLOC(kg_k_fft,(3,npw_fft))
445 50256 : ABI_MALLOC(cwavef_fft,(2,npw_fft*ndat_))
446 100676368 : kg_k_fft(:,1:npw_k1) = kg_k1(:,1:npw_k1)
447 100542352 : kg_k_fft(:,npw_k1+1:npw_fft) = -kg_k1(:,i0:npw_k1)
448 16752 : call cwavef_double_rfft_trick_pack(cwavef,cwavef_fft,mpi_enreg%me_g0_fft,ndat,npw_k1)
449 : end if
450 :
451 16752 : if (have_to_reequilibrate.and.double_rfft_trick) then
452 0 : ABI_ERROR("In getghc: have_to_reequilibrate cannot be activated with double_rfft_trick")
453 : end if
454 :
455 34962376 : if (have_to_reequilibrate) then
456 : ! Note: for this case we have ndat_=ndat
457 0 : npw_fft = bandfft_kpt(ikpt_this_proc)%npw_fft
458 0 : sendcount_fft => bandfft_kpt(ikpt_this_proc)%sendcount_fft(:)
459 0 : recvcount_fft => bandfft_kpt(ikpt_this_proc)%recvcount_fft(:)
460 0 : senddisp_fft => bandfft_kpt(ikpt_this_proc)%senddisp_fft(:)
461 0 : recvdisp_fft => bandfft_kpt(ikpt_this_proc)%recvdisp_fft(:)
462 0 : indices_pw_fft => bandfft_kpt(ikpt_this_proc)%indices_pw_fft(:)
463 0 : kg_k_fft => bandfft_kpt(ikpt_this_proc)%kg_k_fft(:,:)
464 0 : ABI_MALLOC(buff_wf,(2,npw_k1*ndat) )
465 0 : ABI_MALLOC(cwavef_fft,(2,npw_fft*ndat) )
466 0 : if(ndat>1) then
467 0 : ABI_MALLOC(cwavef_fft_tr, (2,npw_fft*ndat))
468 : end if
469 0 : do idat=1, ndat
470 0 : do ipw = 1 ,npw_k1
471 0 : buff_wf(1:2, idat + ndat*(indices_pw_fft(ipw)-1) ) = cwavef(1:2,ipw + npw_k1*(idat-1))
472 : end do
473 : end do
474 0 : if(ndat > 1) then
475 : call xmpi_alltoallv(buff_wf,2*ndat*sendcount_fft,2*ndat*senddisp_fft, &
476 0 : cwavef_fft_tr,2*ndat*recvcount_fft, 2*ndat*recvdisp_fft, mpi_enreg%comm_fft,ierr)
477 : ! We need to transpose data
478 0 : do idat=1,ndat
479 0 : do ipw = 1 ,npw_fft
480 0 : cwavef_fft(1:2, ipw + npw_fft*(idat-1)) = cwavef_fft_tr(1:2, idat + ndat*(ipw-1))
481 : end do
482 : end do
483 : else
484 : call xmpi_alltoallv(buff_wf,2*sendcount_fft,2*senddisp_fft, &
485 0 : cwavef_fft,2*recvcount_fft, 2*recvdisp_fft, mpi_enreg%comm_fft,ierr)
486 : end if
487 : end if
488 :
489 : ! Apply the local potential to the wavefunction
490 : ! Start from wavefunction in reciprocal space cwavef
491 : ! End with function ghc in reciprocal space also.
492 34962376 : fftalg = gs_ham%ngfft(7)
493 34962376 : if (gs_ham%gpu_option==ABI_GPU_DISABLED.and.fftalg/=401) then
494 171755270 : ABI_MALLOC(work,(2,gs_ham%n4,gs_ham%n5,gs_ham%n6))
495 611322 : else if(nfft_blocks==1) then
496 3056610 : ABI_MALLOC(work,(2,gs_ham%n4,gs_ham%n5,gs_ham%n6*ndat))
497 : else
498 0 : ABI_MALLOC(work,(2,gs_ham%n4,gs_ham%n5,gs_ham%n6*(chunk+residuchunk)))
499 : end if
500 : #ifdef HAVE_OPENMP_OFFLOAD
501 : !$OMP TARGET ENTER DATA MAP(alloc:work) IF(gs_ham%gpu_option==ABI_GPU_OPENMP)
502 : #endif
503 34962376 : weight=one
504 34962376 : if (.not.use_cwavef_r) then
505 34962376 : option_fft=2
506 34962376 : if (nspinortot==2) then
507 : ! Note: for this case we have ndat_=ndat
508 7488423 : ABI_MALLOC(cwavef1,(2,npw_k1*ndat))
509 4992282 : ABI_MALLOC(cwavef2,(2,npw_k1*ndat))
510 2496141 : if(gs_ham%gpu_option == ABI_GPU_OPENMP) then
511 : #ifdef HAVE_OPENMP_OFFLOAD
512 : !$OMP TARGET ENTER DATA MAP(alloc:cwavef1,cwavef2)
513 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:cwavef1,cwavef2,cwavef)
514 : do idat=1,ndat
515 : !$OMP PARALLEL DO COLLAPSE(2)
516 : do ipw=1,npw_k1
517 : do ig=1,2
518 : cwavef1(ig,ipw+(idat-1)*npw_k1)=cwavef(ig,ipw+(idat-1)*my_nspinor*npw_k1)
519 : cwavef2(ig,ipw+(idat-1)*npw_k1)=cwavef(ig,ipw+(idat-1)*my_nspinor*npw_k1+shift1)
520 : end do
521 : end do
522 : end do
523 : #endif
524 : else
525 5299641 : do idat=1,ndat
526 636002639 : do ipw=1,npw_k1
527 1892108994 : cwavef1(1:2,ipw+(idat-1)*npw_k1)=cwavef(1:2,ipw+(idat-1)*my_nspinor*npw_k1)
528 1894912494 : cwavef2(1:2,ipw+(idat-1)*npw_k1)=cwavef(1:2,ipw+(idat-1)*my_nspinor*npw_k1+shift1)
529 : end do
530 : end do
531 : end if
532 : end if
533 : else
534 0 : option_fft=3
535 0 : nfft_blocks=1
536 0 : if (nspinortot==2) then
537 0 : ABI_MALLOC(cwavef1,(0,0))
538 0 : ABI_MALLOC(cwavef2,(0,0))
539 : end if
540 : end if
541 :
542 34962376 : if (gs_ham%nvloc==1) then
543 : ! Treat scalar local potentials
544 :
545 32989271 : if (nspinortot==1) then
546 :
547 32466235 : if (use_cwavef_r) then
548 0 : do i3=1,n6
549 0 : do i2=1,n5
550 0 : do i1=1,n4
551 0 : work(1,i1,i2,i3) = gs_ham%vlocal(i1,i2,i3,1)*cwavef_r(1,i1,i2,i3,1)
552 0 : work(2,i1,i2,i3) = gs_ham%vlocal(i1,i2,i3,1)*cwavef_r(2,i1,i2,i3,1)
553 : end do
554 : end do
555 : end do
556 : end if
557 32466235 : if (have_to_reequilibrate.or.double_rfft_trick) then
558 : call fourwf(1,gs_ham%vlocal,cwavef_fft,cwavef_fft,work,gbound_k1,gbound_k2,&
559 : istwf_k_,kg_k_fft,kg_k_fft,gs_ham%mgfft,mpi_enreg,ndat_,gs_ham%ngfft,&
560 : npw_fft,npw_fft,gs_ham%n4,gs_ham%n5,gs_ham%n6,option_fft,tim_fourwf,&
561 16752 : weight,weight,gpu_option=gs_ham%gpu_option)
562 : else
563 : call fourwf_optmem(1,gs_ham%vlocal,cwavef,ghc,work,gbound_k1,gbound_k2,&
564 : istwf_k_,kg_k1,kg_k2,gs_ham%mgfft,mpi_enreg,ndat_,nfft_blocks,gs_ham%ngfft,&
565 : npw_k1,npw_k2,gs_ham%n4,gs_ham%n5,gs_ham%n6,option_fft,tim_fourwf,&
566 32449483 : weight,weight,gpu_option=gs_ham%gpu_option)
567 : end if
568 :
569 : else
570 : ! nspinortot==2
571 :
572 : ! Note: for this case we have ndat_=ndat
573 523036 : if (nspinor1TreatedByThisProc) then
574 519274 : if (use_cwavef_r) then
575 0 : do i3=1,n6
576 0 : do i2=1,n5
577 0 : do i1=1,n4
578 0 : work(1,i1,i2,i3) = gs_ham%vlocal(i1,i2,i3,1)*cwavef_r(1,i1,i2,i3,1)
579 0 : work(2,i1,i2,i3) = gs_ham%vlocal(i1,i2,i3,1)*cwavef_r(2,i1,i2,i3,1)
580 : end do
581 : end do
582 : end do
583 : #ifdef HAVE_OPENMP_OFFLOAD
584 : !$OMP TARGET UPDATE TO(work) IF(gs_ham%gpu_option==ABI_GPU_OPENMP)
585 : #endif
586 : end if
587 1557822 : ABI_MALLOC(ghc1,(2,npw_k2*ndat))
588 : #ifdef HAVE_OPENMP_OFFLOAD
589 : !$OMP TARGET ENTER DATA MAP(alloc:ghc1) IF(gs_ham%gpu_option==ABI_GPU_OPENMP)
590 : #endif
591 : call fourwf_optmem(1,gs_ham%vlocal,cwavef1,ghc1,work,gbound_k1,gbound_k2,&
592 : istwf_k_,kg_k1,kg_k2,gs_ham%mgfft,mpi_enreg,ndat,nfft_blocks,gs_ham%ngfft,&
593 : npw_k1,npw_k2,gs_ham%n4,gs_ham%n5,gs_ham%n6,option_fft,tim_fourwf,&
594 519274 : weight,weight,gpu_option=gs_ham%gpu_option)
595 519274 : if(gs_ham%gpu_option==ABI_GPU_OPENMP) then
596 : #ifdef HAVE_OPENMP_OFFLOAD
597 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:ghc,ghc1)
598 : do idat=1,ndat
599 : !$OMP PARALLEL DO COLLAPSE(2)
600 : do ipw =1, npw_k2
601 : do ig =1,2
602 : ghc(ig,ipw+(idat-1)*my_nspinor*npw_k2)=ghc1(ig,ipw+(idat-1)*npw_k2)
603 : end do
604 : end do
605 : end do
606 : #endif
607 : else
608 1344221 : do idat=1,ndat
609 258114800 : do ipw =1, npw_k2
610 771136684 : ghc(1:2,ipw+(idat-1)*my_nspinor*npw_k2)=ghc1(1:2,ipw+(idat-1)*npw_k2)
611 : end do
612 : end do
613 : end if
614 : #ifdef HAVE_OPENMP_OFFLOAD
615 : !$OMP TARGET EXIT DATA MAP(delete:ghc1) IF(gs_ham%gpu_option==ABI_GPU_OPENMP)
616 : #endif
617 519274 : ABI_FREE(ghc1)
618 : end if ! spin 1 treated by this proc
619 :
620 523036 : if (nspinor2TreatedByThisProc) then
621 519274 : if (use_cwavef_r) then
622 0 : do i3=1,n6
623 0 : do i2=1,n5
624 0 : do i1=1,n4
625 0 : work(1,i1,i2,i3) = gs_ham%vlocal(i1,i2,i3,1)*cwavef_r(1,i1,i2,i3,2)
626 0 : work(2,i1,i2,i3) = gs_ham%vlocal(i1,i2,i3,1)*cwavef_r(2,i1,i2,i3,2)
627 : end do
628 : end do
629 : end do
630 : end if
631 1557822 : ABI_MALLOC(ghc2,(2,npw_k2*ndat))
632 : #ifdef HAVE_OPENMP_OFFLOAD
633 : !$OMP TARGET ENTER DATA MAP(alloc:ghc2) IF(gs_ham%gpu_option==ABI_GPU_OPENMP)
634 : #endif
635 : call fourwf_optmem(1,gs_ham%vlocal,cwavef2,ghc2,work,gbound_k1,gbound_k2,&
636 : istwf_k_,kg_k1,kg_k2,gs_ham%mgfft,mpi_enreg,ndat,nfft_blocks,gs_ham%ngfft,&
637 : npw_k1,npw_k2,gs_ham%n4,gs_ham%n5,gs_ham%n6,option_fft,tim_fourwf,weight,weight,&
638 519274 : gpu_option=gs_ham%gpu_option)
639 :
640 519274 : if(gs_ham%gpu_option==ABI_GPU_OPENMP) then
641 : #ifdef HAVE_OPENMP_OFFLOAD
642 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:ghc,ghc2)
643 : do idat=1,ndat
644 : !$OMP PARALLEL DO COLLAPSE(2)
645 : do ipw=1,npw_k2
646 : do ig =1,2
647 : ghc(ig,ipw+(idat-1)*my_nspinor*npw_k2+shift2)=ghc2(ig,ipw+(idat-1)*npw_k2)
648 : end do
649 : end do
650 : end do
651 : #endif
652 : else
653 1344221 : do idat=1,ndat
654 258114800 : do ipw=1,npw_k2
655 771136684 : ghc(1:2,ipw+(idat-1)*my_nspinor*npw_k2+shift2)=ghc2(1:2,ipw+(idat-1)*npw_k2)
656 : end do
657 : end do
658 : end if
659 : #ifdef HAVE_OPENMP_OFFLOAD
660 : !$OMP TARGET EXIT DATA MAP(delete:ghc2) IF(gs_ham%gpu_option==ABI_GPU_OPENMP)
661 : #endif
662 519274 : ABI_FREE(ghc2)
663 : end if ! spin 2 treated by this proc
664 :
665 : end if ! nspinortot
666 :
667 1973105 : else if (gs_ham%nvloc==4) then
668 : ! Treat non-collinear local potentials
669 : ! Note: for this case we have ndat_=ndat
670 5919315 : ABI_MALLOC(ghc1,(2,npw_k2*ndat))
671 3946210 : ABI_MALLOC(ghc2,(2,npw_k2*ndat))
672 3946210 : ABI_MALLOC(ghc3,(2,npw_k2*ndat))
673 3946210 : ABI_MALLOC(ghc4,(2,npw_k2*ndat))
674 1973105 : if(gs_ham%gpu_option == ABI_GPU_OPENMP) then
675 : #ifdef HAVE_OPENMP_OFFLOAD
676 : !$OMP TARGET ENTER DATA MAP(alloc:ghc1,ghc2,ghc3,ghc4)
677 : call gpu_set_to_zero(ghc1, int(2,c_size_t)*npw_k2*ndat)
678 : call gpu_set_to_zero(ghc2, int(2,c_size_t)*npw_k2*ndat)
679 : call gpu_set_to_zero(ghc3, int(2,c_size_t)*npw_k2*ndat)
680 : call gpu_set_to_zero(ghc4, int(2,c_size_t)*npw_k2*ndat)
681 : #endif
682 : else
683 4430500181 : ghc1(:,:)=zero; ghc2(:,:)=zero; ghc3(:,:)=zero ; ghc4(:,:)=zero
684 : end if
685 1973105 : if (use_cwavef_r) then
686 0 : ABI_MALLOC(vlocal_tmp,(0,0,0))
687 : else
688 9865525 : ABI_MALLOC(vlocal_tmp,(gs_ham%n4,gs_ham%n5,gs_ham%n6))
689 : end if
690 : ! ghc1=v11*phi1
691 1973105 : if (nspinor1TreatedByThisProc) then
692 1951185 : if (use_cwavef_r) then
693 0 : do i3=1,n6
694 0 : do i2=1,n5
695 0 : do i1=1,n4
696 0 : work(1,i1,i2,i3) = gs_ham%vlocal(i1,i2,i3,1)*cwavef_r(1,i1,i2,i3,1)
697 0 : work(2,i1,i2,i3) = gs_ham%vlocal(i1,i2,i3,1)*cwavef_r(2,i1,i2,i3,1)
698 : end do
699 : end do
700 : end do
701 : #ifdef HAVE_OPENMP_OFFLOAD
702 : !$OMP TARGET UPDATE TO(work) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
703 : #endif
704 : else
705 : ! LB,07/22:
706 : ! Weird segmentation fault encountered here if called with multithreaded_getghc for big systems.
707 : ! Using an explicit loop instead of fortran syntax seems to solve the problem, I don't understand why...
708 : !vlocal_tmp(:,:,:)=gs_ham%vlocal(:,:,:,1)
709 33265230 : do i3=1,n6
710 570470291 : do i2=1,n5
711 10235254527 : do i1=1,n4
712 10203940482 : vlocal_tmp(i1,i2,i3) = gs_ham%vlocal(i1,i2,i3,1)
713 : end do
714 : end do
715 : end do
716 : end if
717 : call fourwf_optmem(1,vlocal_tmp,cwavef1,ghc1,work,gbound_k1,gbound_k2,&
718 : istwf_k_,kg_k1,kg_k2,gs_ham%mgfft,mpi_enreg,ndat,nfft_blocks,gs_ham%ngfft,&
719 : npw_k1,npw_k2,gs_ham%n4,gs_ham%n5,gs_ham%n6,option_fft,tim_fourwf,weight,weight,&
720 1951185 : gpu_option=gs_ham%gpu_option)
721 : end if
722 : ! ghc2=v22*phi2
723 1973105 : if (nspinor2TreatedByThisProc) then
724 1951185 : if (use_cwavef_r) then
725 0 : do i3=1,n6
726 0 : do i2=1,n5
727 0 : do i1=1,n4
728 0 : work(1,i1,i2,i3) = gs_ham%vlocal(i1,i2,i3,2)*cwavef_r(1,i1,i2,i3,2)
729 0 : work(2,i1,i2,i3) = gs_ham%vlocal(i1,i2,i3,2)*cwavef_r(2,i1,i2,i3,2)
730 : end do
731 : end do
732 : end do
733 : #ifdef HAVE_OPENMP_OFFLOAD
734 : !$OMP TARGET UPDATE TO(work) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
735 : #endif
736 : else
737 : ! LB,07/22:
738 : ! Weird segmentation fault encountered here if called with multithreaded_getghc for big systems.
739 : ! Using an explicit loop instead of fortran syntax seems to solve the problem, I don't understand why...
740 : !vlocal_tmp(:,:,:)=gs_ham%vlocal(:,:,:,2)
741 33265230 : do i3=1,n6
742 570470291 : do i2=1,n5
743 10235254527 : do i1=1,n4
744 10203940482 : vlocal_tmp(i1,i2,i3) = gs_ham%vlocal(i1,i2,i3,2)
745 : end do
746 : end do
747 : end do
748 : #ifdef HAVE_OPENMP_OFFLOAD
749 : !$OMP TARGET UPDATE TO(work) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
750 : #endif
751 : end if
752 : call fourwf_optmem(1,vlocal_tmp,cwavef2,ghc2,work,gbound_k1,gbound_k2,&
753 : istwf_k_,kg_k1,kg_k2,gs_ham%mgfft,mpi_enreg,ndat,nfft_blocks,gs_ham%ngfft,&
754 : npw_k1,npw_k2,gs_ham%n4,gs_ham%n5,gs_ham%n6,option_fft,tim_fourwf,weight,weight,&
755 1951185 : gpu_option=gs_ham%gpu_option)
756 : end if
757 1973105 : ABI_FREE(vlocal_tmp)
758 1973105 : cplex=2
759 1973105 : if (use_cwavef_r) then
760 0 : ABI_MALLOC(vlocal_tmp,(0,0,0))
761 : else
762 9865525 : ABI_MALLOC(vlocal_tmp,(cplex*gs_ham%n4,gs_ham%n5,gs_ham%n6))
763 : end if
764 : ! ghc3=(re(v12)-im(v12))*phi1
765 1973105 : if (nspinor1TreatedByThisProc) then
766 1951185 : if (use_cwavef_r) then
767 0 : do i3=1,n6
768 0 : do i2=1,n5
769 0 : do i1=1,n4
770 0 : work(1,i1,i2,i3)= gs_ham%vlocal(i1,i2,i3,3)*cwavef_r(1,i1,i2,i3,1)+gs_ham%vlocal(i1,i2,i3,4)*cwavef_r(2,i1,i2,i3,1)
771 0 : work(2,i1,i2,i3)=-gs_ham%vlocal(i1,i2,i3,4)*cwavef_r(1,i1,i2,i3,1)+gs_ham%vlocal(i1,i2,i3,3)*cwavef_r(2,i1,i2,i3,1)
772 : end do
773 : end do
774 : end do
775 : #ifdef HAVE_OPENMP_OFFLOAD
776 : !$OMP TARGET UPDATE TO(work) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
777 : #endif
778 : else
779 33265230 : do i3=1,gs_ham%n6
780 570470291 : do i2=1,gs_ham%n5
781 10235254527 : do i1=1,gs_ham%n4
782 9666735421 : vlocal_tmp(2*i1-1,i2,i3)= gs_ham%vlocal(i1,i2,i3,3)
783 10203940482 : vlocal_tmp(2*i1 ,i2,i3)=-gs_ham%vlocal(i1,i2,i3,4)
784 : end do
785 : end do
786 : end do
787 : end if
788 : call fourwf_optmem(cplex,vlocal_tmp,cwavef1,ghc3,work,gbound_k1,gbound_k2,&
789 : istwf_k_,kg_k1,kg_k2,gs_ham%mgfft,mpi_enreg,ndat,nfft_blocks,gs_ham%ngfft,&
790 : npw_k1,npw_k2,gs_ham%n4,gs_ham%n5,gs_ham%n6,option_fft,tim_fourwf,weight,weight,&
791 1951185 : gpu_option=gs_ham%gpu_option)
792 : end if
793 : ! ghc4=(re(v12)+im(v12))*phi2
794 1973105 : if (nspinor2TreatedByThisProc) then
795 1951185 : if (use_cwavef_r) then
796 0 : do i3=1,n6
797 0 : do i2=1,n5
798 0 : do i1=1,n4
799 0 : work(1,i1,i2,i3)= gs_ham%vlocal(i1,i2,i3,3)*cwavef_r(1,i1,i2,i3,2)-gs_ham%vlocal(i1,i2,i3,4)*cwavef_r(2,i1,i2,i3,2)
800 0 : work(2,i1,i2,i3)= gs_ham%vlocal(i1,i2,i3,4)*cwavef_r(1,i1,i2,i3,2)+gs_ham%vlocal(i1,i2,i3,3)*cwavef_r(2,i1,i2,i3,2)
801 : end do
802 : end do
803 : end do
804 : #ifdef HAVE_OPENMP_OFFLOAD
805 : !$OMP TARGET UPDATE TO(work) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
806 : #endif
807 : else
808 33265230 : do i3=1,gs_ham%n6
809 570470291 : do i2=1,gs_ham%n5
810 10235254527 : do i1=1,gs_ham%n4
811 9666735421 : vlocal_tmp(2*i1-1,i2,i3)= gs_ham%vlocal(i1,i2,i3,3)
812 10203940482 : vlocal_tmp(2*i1 ,i2,i3)= gs_ham%vlocal(i1,i2,i3,4)
813 : end do
814 : end do
815 : end do
816 : end if
817 : call fourwf_optmem(cplex,vlocal_tmp,cwavef2,ghc4,work,gbound_k1,gbound_k2,&
818 : istwf_k_,kg_k1,kg_k2,gs_ham%mgfft,mpi_enreg,ndat,nfft_blocks,gs_ham%ngfft,&
819 : npw_k1,npw_k2,gs_ham%n4,gs_ham%n5,gs_ham%n6,option_fft,tim_fourwf,weight,weight,&
820 1951185 : gpu_option=gs_ham%gpu_option)
821 : end if
822 1973105 : ABI_FREE(vlocal_tmp)
823 : ! Build ghc from pieces
824 : ! (v11,v22,Re(v12)+iIm(v12);Re(v12)-iIm(v12))(psi1;psi2): matrix product
825 1973105 : if (mpi_enreg%paral_spinor==0) then
826 1929265 : if(gs_ham%gpu_option == ABI_GPU_OPENMP) then
827 : #ifdef HAVE_OPENMP_OFFLOAD
828 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:ghc,ghc1,ghc2,ghc3,ghc4)
829 : do idat=1,ndat
830 : !$OMP PARALLEL DO
831 : do ipw=1,npw_k2
832 : ghc(1,ipw+(idat-1)*my_nspinor*npw_k2) =ghc1(1,ipw+(idat-1)*npw_k2)+ghc4(1,ipw+(idat-1)*npw_k2)
833 : ghc(2,ipw+(idat-1)*my_nspinor*npw_k2) =ghc1(2,ipw+(idat-1)*npw_k2)+ghc4(2,ipw+(idat-1)*npw_k2)
834 : ghc(1,ipw+(idat-1)*my_nspinor*npw_k2+shift2)=ghc3(1,ipw+(idat-1)*npw_k2)+ghc2(1,ipw+(idat-1)*npw_k2)
835 : ghc(2,ipw+(idat-1)*my_nspinor*npw_k2+shift2)=ghc3(2,ipw+(idat-1)*npw_k2)+ghc2(2,ipw+(idat-1)*npw_k2)
836 : end do
837 : end do
838 : #endif
839 : else
840 3858530 : do idat=1,ndat
841 363120213 : do ipw=1,npw_k2
842 1077785049 : ghc(1:2,ipw+(idat-1)*my_nspinor*npw_k2) =ghc1(1:2,ipw+(idat-1)*npw_k2)+ghc4(1:2,ipw+(idat-1)*npw_k2)
843 1079714314 : ghc(1:2,ipw+(idat-1)*my_nspinor*npw_k2+shift2)=ghc3(1:2,ipw+(idat-1)*npw_k2)+ghc2(1:2,ipw+(idat-1)*npw_k2)
844 : end do
845 : end do
846 : end if
847 : else
848 43840 : call xmpi_sum(ghc4,mpi_enreg%comm_spinor,ierr)
849 43840 : call xmpi_sum(ghc3,mpi_enreg%comm_spinor,ierr)
850 43840 : if(gs_ham%gpu_option == ABI_GPU_OPENMP) then
851 : #ifdef HAVE_OPENMP_OFFLOAD
852 : if (nspinor1TreatedByThisProc) then
853 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:ghc,ghc1,ghc4)
854 : do idat=1,ndat
855 : !$OMP PARALLEL DO
856 : do ipw=1,npw_k2
857 : ghc(1,ipw+(idat-1)*my_nspinor*npw_k2)=ghc1(1,ipw+(idat-1)*npw_k2)+ghc4(1,ipw+(idat-1)*npw_k2)
858 : ghc(2,ipw+(idat-1)*my_nspinor*npw_k2)=ghc1(2,ipw+(idat-1)*npw_k2)+ghc4(2,ipw+(idat-1)*npw_k2)
859 : end do
860 : end do
861 : else if (nspinor2TreatedByThisProc) then
862 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:ghc,ghc2,ghc3)
863 : do idat=1,ndat
864 : !$OMP PARALLEL DO
865 : do ipw=1,npw_k2
866 : ghc(1,ipw+(idat-1)*my_nspinor*npw_k2+shift2)=ghc3(1,ipw+(idat-1)*npw_k2)+ghc2(1,ipw+(idat-1)*npw_k2)
867 : ghc(2,ipw+(idat-1)*my_nspinor*npw_k2+shift2)=ghc3(2,ipw+(idat-1)*npw_k2)+ghc2(2,ipw+(idat-1)*npw_k2)
868 : end do
869 : end do
870 : end if
871 : #endif
872 : else
873 43840 : if (nspinor1TreatedByThisProc) then
874 43840 : do idat=1,ndat
875 4934960 : do ipw=1,npw_k2
876 14695280 : ghc(1:2,ipw+(idat-1)*my_nspinor*npw_k2)=ghc1(1:2,ipw+(idat-1)*npw_k2)+ghc4(1:2,ipw+(idat-1)*npw_k2)
877 : end do
878 : end do
879 21920 : else if (nspinor2TreatedByThisProc) then
880 43840 : do idat=1,ndat
881 4934960 : do ipw=1,npw_k2
882 14695280 : ghc(1:2,ipw+(idat-1)*my_nspinor*npw_k2+shift2)=ghc3(1:2,ipw+(idat-1)*npw_k2)+ghc2(1:2,ipw+(idat-1)*npw_k2)
883 : end do
884 : end do
885 : end if
886 : end if
887 : end if
888 : #ifdef HAVE_OPENMP_OFFLOAD
889 : !$OMP TARGET EXIT DATA MAP(delete:ghc1,ghc2,ghc3,ghc4) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
890 : #endif
891 1973105 : ABI_FREE(ghc1)
892 1973105 : ABI_FREE(ghc2)
893 1973105 : ABI_FREE(ghc3)
894 1973105 : ABI_FREE(ghc4)
895 : end if ! nvloc
896 :
897 34962376 : if (nspinortot==2) then
898 : #ifdef HAVE_OPENMP_OFFLOAD
899 : !$OMP TARGET EXIT DATA MAP(delete:cwavef1,cwavef2) IF (gs_ham%gpu_option == ABI_GPU_OPENMP .and. .not.use_cwavef_r)
900 : #endif
901 2496141 : ABI_FREE(cwavef1)
902 2496141 : ABI_FREE(cwavef2)
903 : end if
904 :
905 : #ifdef HAVE_OPENMP_OFFLOAD
906 : !$OMP TARGET EXIT DATA MAP(delete:work) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
907 : #endif
908 34962376 : ABI_FREE(work)
909 :
910 34962376 : if (double_rfft_trick) then
911 16752 : call cwavef_double_rfft_trick_unpack(ghc,cwavef_fft,mpi_enreg%me_g0_fft,ndat,npw_k1)
912 16752 : ABI_FREE(kg_k_fft)
913 16752 : ABI_FREE(cwavef_fft)
914 : end if
915 :
916 : ! Retrieve eventually original FFT distrib
917 34962376 : if(have_to_reequilibrate) then
918 : ! Note: for this case we have ndat_=ndat
919 0 : if(ndat > 1 ) then
920 0 : do idat=1,ndat
921 0 : do ipw = 1 ,npw_fft
922 0 : cwavef_fft_tr(1:2, idat + ndat*(ipw-1)) = cwavef_fft(1:2, ipw + npw_fft*(idat-1))
923 : end do
924 : end do
925 : call xmpi_alltoallv(cwavef_fft_tr,2*ndat*recvcount_fft, 2*ndat*recvdisp_fft, &
926 0 : buff_wf,2*ndat*sendcount_fft,2*ndat*senddisp_fft, mpi_enreg%comm_fft,ierr)
927 : else
928 : call xmpi_alltoallv(cwavef_fft,2*recvcount_fft, 2*recvdisp_fft, &
929 0 : buff_wf,2*sendcount_fft,2*senddisp_fft, mpi_enreg%comm_fft,ierr)
930 : end if
931 0 : do idat=1,ndat
932 0 : do ipw = 1 ,npw_k2
933 0 : ghc(1:2,ipw + npw_k2*(idat-1)) = buff_wf(1:2, idat + ndat*(indices_pw_fft(ipw)-1))
934 : end do
935 : end do
936 0 : ABI_FREE(buff_wf)
937 0 : ABI_FREE(cwavef_fft)
938 0 : if(ndat > 1) then
939 0 : ABI_FREE(cwavef_fft_tr)
940 : end if
941 : end if
942 :
943 : ! Add metaGGA contribution
944 34962376 : if (associated(gs_ham%vxctaulocal)) then
945 90464 : ABI_CHECK(k1_eq_k2, 'metaGGA not allowed for k/=k_^prime!')
946 542784 : ABI_CHECK_IEQ(size(gs_ham%vxctaulocal), gs_ham%n4*gs_ham%n5*gs_ham%n6*gs_ham%nvloc*4, 'wrong sizes for vxctaulocal!')
947 :
948 271392 : ABI_MALLOC(ghc_mGGA,(2,npw_k1*my_nspinor*ndat))
949 : #ifdef HAVE_OPENMP_OFFLOAD
950 : !$OMP TARGET UPDATE FROM(cwavef) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
951 : #endif
952 90464 : if (double_rfft_trick) then
953 480 : ABI_MALLOC(kg_k_fft,(3,npw_fft))
954 480 : ABI_MALLOC(cwavef_fft,(2,npw_fft*ndat_))
955 320 : ABI_MALLOC(ghc_mGGA_fft,(2,npw_fft*ndat_))
956 413600 : kg_k_fft(:,1:npw_k1) = kg_k1(:,1:npw_k1)
957 412320 : kg_k_fft(:,npw_k1+1:npw_fft) = -kg_k1(:,i0:npw_k1)
958 160 : call cwavef_double_rfft_trick_pack(cwavef,cwavef_fft,mpi_enreg%me_g0_fft,ndat,npw_k1)
959 : call getghc_mGGA(cwavef_fft,ghc_mGGA_fft,gbound_k1,gs_ham%gprimd,istwf_k_,kg_k_fft,kpt_k1,&
960 : & gs_ham%mgfft,mpi_enreg,ndat_,gs_ham%ngfft,npw_fft,gs_ham%nvloc,&
961 160 : & gs_ham%n4,gs_ham%n5,gs_ham%n6,my_nspinor,gs_ham%vxctaulocal,gs_ham%gpu_option)
962 160 : call cwavef_double_rfft_trick_unpack(ghc_mGGA,ghc_mGGA_fft,mpi_enreg%me_g0_fft,ndat,npw_k1)
963 160 : ABI_FREE(kg_k_fft)
964 160 : ABI_FREE(cwavef_fft)
965 160 : ABI_FREE(ghc_mGGA_fft)
966 : else
967 : call getghc_mGGA(cwavef,ghc_mGGA,gbound_k1,gs_ham%gprimd,istwf_k_,kg_k1,kpt_k1,&
968 : & gs_ham%mgfft,mpi_enreg,ndat,gs_ham%ngfft,npw_k1,gs_ham%nvloc,&
969 90304 : & gs_ham%n4,gs_ham%n5,gs_ham%n6,my_nspinor,gs_ham%vxctaulocal,gs_ham%gpu_option)
970 : end if
971 : #ifdef HAVE_OPENMP_OFFLOAD
972 : !$OMP TARGET UPDATE FROM(ghc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
973 : #endif
974 : !LB-2025-12-16: intel19 crashes here when npw_k1 is too big. A solution is to write the loop explicitely
975 : !ghc(1:2,1:npw_k2*my_nspinor*ndat)=ghc(1:2,1:npw_k2*my_nspinor*ndat)+ghc_mGGA(1:2,1:npw_k2*my_nspinor*ndat)
976 56805312 : do idat=1,npw_k1*my_nspinor*ndat
977 170235008 : do ig=1,2
978 170144544 : ghc(ig,idat)=ghc(ig,idat)+ghc_mGGA(ig,idat)
979 : end do
980 : end do
981 : #ifdef HAVE_OPENMP_OFFLOAD
982 : !$OMP TARGET UPDATE TO(ghc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
983 : #endif
984 90464 : ABI_FREE(ghc_mGGA)
985 : end if
986 :
987 : ! Add nuclear dipole moment contribution
988 34962376 : if (associated(gs_ham%vectornd)) then
989 673202 : ABI_CHECK(k1_eq_k2, 'nuclear dipole vector potential not allowed for k/=k_^prime!')
990 4039212 : if (size(gs_ham%vectornd)/=gs_ham%n4*gs_ham%n5*gs_ham%n6*gs_ham%nvloc*3) then
991 0 : ABI_BUG('wrong sizes for vectornd in getghc!')
992 : end if
993 2019606 : ABI_MALLOC(ghc_vectornd,(2,npw_k1*my_nspinor*ndat))
994 939408950 : ghc_vectornd=zero
995 : #ifdef HAVE_OPENMP_OFFLOAD
996 : !$OMP TARGET UPDATE FROM(cwavef) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
997 : #endif
998 : call getghc_nucdip(cwavef,ghc_vectornd,gbound_k1,istwf_k_,kg_k1,kpt_k1,&
999 : gs_ham%mgfft,mpi_enreg,ndat,gs_ham%ngfft,npw_k1,gs_ham%nvloc,&
1000 673202 : gs_ham%n4,gs_ham%n5,gs_ham%n6,my_nspinor,gs_ham%vectornd,gs_ham%vlocal,gs_ham%zora,gs_ham%gpu_option)
1001 :
1002 : #ifdef HAVE_OPENMP_OFFLOAD
1003 : !$OMP TARGET UPDATE FROM(ghc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1004 : #endif
1005 : !LB-2025-12-16: intel19 crashes here when npw_k1 is too big. A solution is to write the loop explicitely
1006 : !ghc(1:2,1:npw_k2*my_nspinor*ndat)=ghc(1:2,1:npw_k2*my_nspinor*ndat)+ghc_vectornd(1:2,1:npw_k2*my_nspinor*ndat)
1007 313585118 : do idat=1,npw_k1*my_nspinor*ndat
1008 939408950 : do ig=1,2
1009 938735748 : ghc(ig,idat)=ghc(ig,idat)+ghc_vectornd(ig,idat)
1010 : end do
1011 : end do
1012 : #ifdef HAVE_OPENMP_OFFLOAD
1013 : !$OMP TARGET UPDATE TO(ghc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1014 : #endif
1015 :
1016 673202 : ABI_FREE(ghc_vectornd)
1017 : end if
1018 :
1019 : ! If only local part is applied, still we have to filter the result because of dilatmx.
1020 : ! Otherwise it is done when adding kinetic term.
1021 34962376 : if (type_calc==1.and.filter_dilatmx_loc_) then
1022 :
1023 204650 : ABI_CHECK(gs_ham%use_gbt == 0, "use_gbt not coded!")
1024 :
1025 204650 : if(gs_ham%gpu_option == ABI_GPU_OPENMP) then
1026 : #ifdef HAVE_OPENMP_OFFLOAD
1027 : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:ghc,kinpw_k2)
1028 : do idat=1,ndat
1029 : !$OMP PARALLEL DO PRIVATE(igspinor) COLLAPSE(2)
1030 : do ispinor=1,my_nspinor
1031 : do ig=1,npw_k2
1032 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1033 : if(kinpw_k2(ig)>hugevalue) ghc(:,igspinor)=zero
1034 : end do ! ig
1035 : end do ! ispinor
1036 : end do
1037 : #endif
1038 : else
1039 1014410 : do idat=1,ndat
1040 : !$OMP PARALLEL DO PRIVATE(igspinor) COLLAPSE(2) IF(gemm_nonlop_use_gemm)
1041 2016282 : do ispinor=1,my_nspinor
1042 125631396 : do ig=1,npw_k2
1043 123819764 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1044 129219672 : if(kinpw_k2(ig)>hugevalue) ghc(:,igspinor)=zero
1045 : end do ! ig
1046 : end do ! ispinor
1047 : end do
1048 : end if
1049 : end if
1050 :
1051 : end if ! type_calc
1052 :
1053 : ABI_NVTX_END_RANGE()
1054 :
1055 34965650 : if (any(type_calc == [0, 2, 3])) then
1056 : !============================================================
1057 : ! Application of the non-local potential and the Fock potential
1058 : !============================================================
1059 :
1060 : ABI_NVTX_START_RANGE(NVTX_GETGHC_NLOCPOT)
1061 :
1062 34757726 : if (type_calc==0 .or. type_calc==2) then
1063 34757726 : signs=2 ; choice=1 ; nnlout=1 ; idir=0 ; tim_nonlop=1
1064 34757726 : cpopt_here=-1;if (gs_ham%usepaw==1) cpopt_here=cpopt
1065 34757726 : if (has_fock) then
1066 142256 : if (gs_ham%usepaw==1) then
1067 13114 : cpopt_here=max(cpopt,0)
1068 13114 : if (cpopt<2) then
1069 91798 : ABI_MALLOC(cwaveprj_fock,(gs_ham%natom,my_nspinor*ndat))
1070 39342 : ABI_MALLOC(dimcprj,(gs_ham%natom))
1071 13114 : call pawcprj_getdim(dimcprj,gs_ham%natom,gs_ham%nattyp,gs_ham%ntypat,gs_ham%typat,fock%pawtab,'O')
1072 13114 : call pawcprj_alloc(cwaveprj_fock,0,dimcprj)
1073 13114 : ABI_FREE(dimcprj)
1074 : else
1075 0 : cwaveprj_fock=>cwaveprj
1076 : end if
1077 13114 : cwaveprj_nonlop=>cwaveprj_fock
1078 : else
1079 129142 : cwaveprj_nonlop=>cwaveprj
1080 129142 : cwaveprj_fock=>cwaveprj
1081 : end if
1082 : else
1083 34615470 : cwaveprj_nonlop=>cwaveprj
1084 : end if
1085 34757726 : paw_opt=gs_ham%usepaw ; if (sij_opt/=0) paw_opt=sij_opt+3
1086 70970152 : lambda_ndat = lambda
1087 :
1088 34757726 : if (gs_ham%use_gbt == 0) then
1089 34618607 : if (gs_ham%usepaw==0) gsc_ptr => nonlop_dum
1090 34618607 : if (gs_ham%usepaw==1) gsc_ptr => gsc
1091 :
1092 : call nonlop(choice,cpopt_here,cwaveprj_nonlop,enlout,gs_ham,idir,lambda_ndat,mpi_enreg,ndat,&
1093 34618607 : nnlout,paw_opt,signs,gsc_ptr,tim_nonlop,cwavef,gvnlxc_,select_k=select_k_)
1094 : else
1095 : ! GBT case. Treat up and down components separately.
1096 139119 : gs_ham%nspinor = 1
1097 :
1098 : ! Split cwavef and gvnlxc
1099 417357 : ABI_MALLOC(cwavef_spin, (2, npw_k1*ndat))
1100 278238 : ABI_MALLOC(gvnlxc_spin, (2, npw_k1*ndat))
1101 :
1102 139119 : if (gs_ham%usepaw==0) gsc_ptr => nonlop_dum
1103 139119 : if (gs_ham%usepaw==1) gsc_ptr => gsc
1104 :
1105 : ! Apply Vnl{k-q/2} to u^up
1106 139119 : gs_ham%ispin_gbt = 1
1107 139119 : call cg_copy_spin(1, npw_k1, nspinortot, ndat, cwavef, cwavef_spin)
1108 : call nonlop(choice, cpopt_here, cwaveprj_nonlop, enlout, gs_ham, idir, lambda_ndat, mpi_enreg, ndat, &
1109 139119 : nnlout, paw_opt, signs, gsc_ptr, tim_nonlop, cwavef_spin, gvnlxc_spin, select_k=K_H_K)
1110 : ! Insert results in the right position.
1111 139119 : call cg_put_spin(1, npw_k1, nspinortot, ndat, gvnlxc_spin, gvnlxc_)
1112 :
1113 : ! Apply H_{k+q/2} to u^down
1114 : !gvnlxc_spin = zero
1115 139119 : gs_ham%ispin_gbt = 2
1116 139119 : call cg_copy_spin(2, npw_k1, nspinortot, ndat, cwavef, cwavef_spin)
1117 : call nonlop(choice, cpopt_here, cwaveprj_nonlop, enlout, gs_ham, idir, lambda_ndat, mpi_enreg, ndat, &
1118 139119 : nnlout, paw_opt, signs, gsc_ptr, tim_nonlop, cwavef_spin, gvnlxc_spin, select_k=KPRIME_H_KPRIME)
1119 : ! Insert results in the right position.
1120 139119 : call cg_put_spin(2, npw_k1, nspinortot, ndat, gvnlxc_spin, gvnlxc_)
1121 :
1122 139119 : gs_ham%nspinor = 2
1123 139119 : ABI_FREE(cwavef_spin)
1124 139119 : ABI_FREE(gvnlxc_spin)
1125 : end if ! use_gbt
1126 :
1127 34757726 : if (gs_ham%usepaw==1 .and. has_fock)then
1128 13114 : if (fock_get_getghc_call(fock)==1) then
1129 13114 : if(gs_ham%gpu_option==ABI_GPU_KOKKOS) then
1130 : #if defined HAVE_GPU && defined HAVE_YAKL
1131 : ABI_MALLOC_MANAGED(gvnlc, (/2,npw_k2*my_nspinor*ndat/))
1132 : #endif
1133 : else
1134 39342 : ABI_MALLOC(gvnlc, (2,npw_k2*my_nspinor*ndat))
1135 : #if defined HAVE_OPENMP_OFFLOAD
1136 : !$OMP TARGET ENTER DATA MAP(to:gvnlc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1137 : #endif
1138 : end if
1139 13114 : if(gs_ham%gpu_option==ABI_GPU_OPENMP) then
1140 0 : call gpu_copy(gvnlc, gvnlxc_, int(2,c_size_t)*npw_k2*my_nspinor*ndat)
1141 : else
1142 48885164 : gvnlc=gvnlxc_
1143 : end if
1144 : endif
1145 : endif
1146 :
1147 : ! Calculation of the Fock exact exchange contribution from the Fock or ACE operator
1148 26152392 : if (has_fock) then
1149 142256 : if (fock_get_getghc_call(fock)==1) then
1150 : if (gs_ham%usepaw==0) cwaveprj_idat => cwaveprj
1151 142256 : if (fock%use_ACE==0) then
1152 : #if defined HAVE_OPENMP_OFFLOAD
1153 : !$OMP TARGET UPDATE FROM(cwavef,gvnlxc_) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1154 : #endif
1155 0 : call timab(360,1,tsec)
1156 0 : call fock_getghc(cwavef,cwaveprj,gvnlxc_,gs_ham,mpi_enreg,ndat)
1157 0 : call timab(360,2,tsec)
1158 : #if defined HAVE_OPENMP_OFFLOAD
1159 : !$OMP TARGET UPDATE TO(cwavef,gvnlxc_) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1160 : #endif
1161 : else
1162 142256 : call fock_ACE_getghc(cwavef,gvnlxc_,gs_ham,mpi_enreg,ndat)
1163 : end if
1164 : end if
1165 : end if
1166 :
1167 0 : else if (type_calc == 3) then
1168 : ! for kinetic and local only, nonlocal and vfock should be zero
1169 0 : if(gs_ham%gpu_option==ABI_GPU_OPENMP) then
1170 0 : call gpu_set_to_zero(gvnlxc_, int(2,c_size_t)*npw_k2*my_nspinor*ndat)
1171 : else
1172 0 : gvnlxc_(:,:) = zero
1173 : end if
1174 : end if ! if(type_calc...
1175 :
1176 : ABI_NVTX_END_RANGE()
1177 :
1178 : !============================================================
1179 : ! Assemble kinetic, local, nonlocal and Fock contributions
1180 : !============================================================
1181 :
1182 : ABI_NVTX_START_RANGE(NVTX_GETGHC_KIN)
1183 :
1184 : #ifdef FC_NVHPC
1185 : !FIXME This Kokkos kernel seems to cause issues under NVHPC so it is disabled
1186 : if (.false.) then
1187 : #else
1188 34757726 : if (gs_ham%gpu_option == ABI_GPU_KOKKOS) then
1189 : #endif
1190 :
1191 : #if defined(HAVE_GPU_CUDA) && defined(HAVE_KOKKOS)
1192 : call assemble_energy_contribution_kokkos(c_loc(ghc), &
1193 : & c_loc(gsc), c_loc(kinpw_k2), c_loc(cwavef), c_loc(gvnlxc_), &
1194 : & ndat, my_nspinor, npw_k2, sij_opt, k1_eq_k2, hugevalue)
1195 : ! sync device so that data can be reused safely on host
1196 : ! will probably be moved elsewhere once all the scf loop runs on device
1197 : call gpu_device_synchronize()
1198 : #endif
1199 :
1200 : else
1201 :
1202 : #ifdef FC_NVHPC
1203 : #if defined(HAVE_GPU_CUDA) && defined(HAVE_KOKKOS)
1204 : !Related to FIXME above
1205 : if (gs_ham%gpu_option == ABI_GPU_KOKKOS) call gpu_device_synchronize()
1206 : #endif
1207 : #endif
1208 :
1209 : ! Assemble modified kinetic, local and nonlocal contributions
1210 : ! to <G|H|C(n,k)>. Take also into account built-in debugging.
1211 34757726 : if(prtvol/=-level)then
1212 :
1213 34757726 : if (gs_ham%gpu_option == ABI_GPU_OPENMP) then
1214 : ! OpenMP GPU
1215 : #ifdef HAVE_OPENMP_OFFLOAD
1216 : if (k1_eq_k2) then
1217 :
1218 : !MG 20260102: With nvfortran 23.11-0, this kernel is a bottleneck due to the filter on kinpw_k2 and the update of ghc.
1219 : ! Solution: branch-free mask + manual loop unrolling. loop unrolling is crucial.
1220 : ! The version with COLLAPSE(3) is faster.
1221 :
1222 : ! !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:ghc,kinpw_k2,gvnlxc_,gsc,cwavef)
1223 : ! !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(3) PRIVATE(igspinor, kscale) MAP(to:ghc,kinpw_k2,gvnlxc_,gsc,cwavef)
1224 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:ghc,kinpw_k2,gvnlxc_,cwavef)
1225 : do idat=1,ndat
1226 : do ispinor=1,my_nspinor
1227 : ! !$OMP PARALLEL DO PRIVATE(igspinor, kscale)
1228 : !$OMP PARALLEL DO PRIVATE(igspinor)
1229 : do ig=1,npw_k2
1230 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1231 :
1232 : ! ! New version with branch-free mask + manual loop unrolling.
1233 : ! kscale = merge(one, zero, kinpw_k2(ig) < hugevalue)
1234 :
1235 : ! ghc(1,igspinor) = kscale * (ghc(1,igspinor) + kinpw_k2(ig)*cwavef(1,igspinor) + gvnlxc_(1,igspinor))
1236 : ! ghc(2,igspinor) = kscale * (ghc(2,igspinor) + kinpw_k2(ig)*cwavef(2,igspinor) + gvnlxc_(2,igspinor))
1237 : ! if (sij_opt == 1) then
1238 : ! gsc(1,igspinor) = kscale * gsc(1,igspinor)
1239 : ! gsc(2,igspinor) = kscale * gsc(2,igspinor)
1240 : ! end if
1241 : if(kinpw_k2(ig)<huge(zero)*1.d-11)then
1242 : ghc(1,igspinor) = ghc(1,igspinor) + kinpw_k2(ig)*cwavef(1,igspinor) + gvnlxc_(1,igspinor)
1243 : ghc(2,igspinor) = ghc(2,igspinor) + kinpw_k2(ig)*cwavef(2,igspinor) + gvnlxc_(2,igspinor)
1244 : else
1245 : ghc(1,igspinor)=zero
1246 : ghc(2,igspinor)=zero
1247 : end if
1248 : end do ! ig
1249 : end do ! ispinor
1250 : end do ! idat
1251 : if(sij_opt==1) then
1252 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:gsc,kinpw_k2)
1253 : do idat=1,ndat
1254 : do ispinor=1,my_nspinor
1255 : !$OMP PARALLEL DO PRIVATE(igspinor)
1256 : do ig=1,npw_k2
1257 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1258 : if(kinpw_k2(ig)>=huge(zero)*1.d-11)then
1259 : gsc(1,igspinor) = zero
1260 : gsc(2,igspinor) = zero
1261 : end if
1262 : end do ! ig
1263 : end do ! ispinor
1264 : end do ! idat
1265 : end if
1266 : else
1267 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:ghc,gvnlxc_,kinpw_k2)
1268 : do idat=1,ndat
1269 : do ispinor=1,my_nspinor
1270 : !$OMP PARALLEL DO PRIVATE(igspinor)
1271 : do ig=1,npw_k2
1272 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1273 : if(kinpw_k2(ig)<hugevalue)then
1274 : ghc(1,igspinor)= ghc(1,igspinor) + gvnlxc_(1,igspinor)
1275 : ghc(2,igspinor)= ghc(2,igspinor) + gvnlxc_(2,igspinor)
1276 : else
1277 : ghc(1,igspinor)=zero
1278 : ghc(2,igspinor)=zero
1279 : end if
1280 : end do ! ig
1281 : end do ! ispinor
1282 : end do ! idat
1283 : if(sij_opt==1) then
1284 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:gsc,kinpw_k2)
1285 : do idat=1,ndat
1286 : do ispinor=1,my_nspinor
1287 : !$OMP PARALLEL DO PRIVATE(igspinor)
1288 : do ig=1,npw_k2
1289 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1290 : if(kinpw_k2(ig)>=huge(zero)*1.d-11)then
1291 : gsc(1,igspinor)=zero
1292 : gsc(2,igspinor)=zero
1293 : end if
1294 : end do ! ig
1295 : end do ! ispinor
1296 : end do ! idat
1297 : end if
1298 : end if
1299 : #endif
1300 :
1301 : else
1302 :
1303 : !CPU (+ Kokkos eventually)
1304 34757726 : if (k1_eq_k2) then
1305 : !$OMP PARALLEL DO PRIVATE(igspinor) COLLAPSE(2) IF(gemm_nonlop_use_gemm)
1306 70691914 : do idat=1,ndat
1307 109182754 : do ispinor=1,my_nspinor
1308 10812790917 : do ig=1,npw_k2
1309 10738226770 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1310 10776717610 : if(kinpw_k2(ig)<hugevalue)then
1311 53193739480 : ghc(:,igspinor) = ghc(:,igspinor) + kinpw_k2(ig)*cwavef(:,igspinor) + gvnlxc_(:,igspinor)
1312 : else
1313 298436622 : ghc(:,igspinor)=zero
1314 142274576 : if (sij_opt==1) gsc(:,igspinor)=zero
1315 : end if
1316 : end do ! ig
1317 : end do ! ispinor
1318 : end do ! idat
1319 : else
1320 139119 : if (gs_ham%use_gbt == 0) then
1321 : !$OMP PARALLEL DO PRIVATE(igspinor) COLLAPSE(2) IF(gemm_nonlop_use_gemm)
1322 0 : do idat=1,ndat
1323 0 : do ispinor=1,my_nspinor
1324 0 : do ig=1,npw_k2
1325 0 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1326 0 : if(kinpw_k2(ig)<hugevalue)then
1327 0 : ghc(:,igspinor)= ghc(:,igspinor) + gvnlxc_(:,igspinor)
1328 : else
1329 0 : ghc(:,igspinor)=zero
1330 0 : if (sij_opt==1) gsc(:,igspinor)=zero
1331 : end if
1332 : end do ! ig
1333 : end do ! ispinor
1334 : end do ! idat
1335 : else
1336 : ! GBT. use different kinetic energies for up and down components.
1337 : !$OMP PARALLEL DO PRIVATE(igspinor, iispinor) COLLAPSE(2) IF(gemm_nonlop_use_gemm)
1338 278238 : do idat=1,ndat
1339 556476 : do ispinor=1,my_nspinor
1340 278238 : iispinor=ispinor; if (mpi_enreg%paral_spinor==1) iispinor=mpi_enreg%me_spinor+1
1341 417357 : if (iispinor == 1) then
1342 31864613 : do ig=1,npw_k2
1343 31725494 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1344 31864613 : if(gs_ham%kinpw_k(ig)<hugevalue)then
1345 153976375 : ghc(:,igspinor) = ghc(:,igspinor) + gs_ham%kinpw_k(ig)*cwavef(:,igspinor) + gvnlxc_(:,igspinor)
1346 : else
1347 2790657 : ghc(:,igspinor)=zero
1348 930219 : if (sij_opt==1) gsc(:,igspinor)=zero
1349 : end if
1350 : end do ! ig
1351 : else
1352 31864613 : do ig=1,npw_k2
1353 31725494 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1354 31864613 : if(gs_ham%kinpw_kp(ig)<hugevalue)then
1355 153976375 : ghc(:,igspinor) = ghc(:,igspinor) + gs_ham%kinpw_kp(ig)*cwavef(:,igspinor) + gvnlxc_(:,igspinor)
1356 : else
1357 2790657 : ghc(:,igspinor)=zero
1358 930219 : if (sij_opt==1) gsc(:,igspinor)=zero
1359 : end if
1360 : end do ! ig
1361 : end if
1362 : end do ! ispinor
1363 : end do ! idat
1364 : end if
1365 : end if
1366 :
1367 : end if ! gs_ham%gpu_option
1368 :
1369 : else
1370 : ! Here, debugging section
1371 : #ifdef HAVE_OPENMP_OFFLOAD
1372 : !$OMP TARGET UPDATE FROM(ghc,gsc,cwavef,gvnlxc_) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1373 : #endif
1374 0 : call wrtout(std_out,' getghc : components of ghc ','PERS')
1375 : write(msg,'(a)')&
1376 0 : 'icp ig ispinor igspinor re/im ghc kinpw cwavef glocc gvnlxc gsc'
1377 0 : call wrtout(std_out,msg,'PERS')
1378 0 : do idat=1,ndat
1379 0 : do ispinor=1,my_nspinor
1380 0 : do ig=1,npw_k2
1381 0 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1382 0 : if(kinpw_k2(ig)<hugevalue)then
1383 0 : if (k1_eq_k2) then
1384 0 : ghcre=kinpw_k2(ig)*cwavef(1,igspinor)+ghc(1,igspinor)+gvnlxc_(1,igspinor)
1385 0 : ghcim=kinpw_k2(ig)*cwavef(2,igspinor)+ghc(2,igspinor)+gvnlxc_(2,igspinor)
1386 : else
1387 0 : ghcre=ghc(1,igspinor)+gvnlxc_(1,igspinor)
1388 0 : ghcim=ghc(2,igspinor)+gvnlxc_(2,igspinor)
1389 : end if
1390 : else
1391 0 : ghcre=zero
1392 0 : ghcim=zero
1393 0 : if (sij_opt==1) gsc(:,igspinor)=zero
1394 : end if
1395 0 : iispinor=ispinor; if (mpi_enreg%paral_spinor==1) iispinor=mpi_enreg%me_spinor+1
1396 0 : if (sij_opt == 1) then
1397 0 : write(msg,'(a,3(1x,i5),6(1x,es13.6))') ' 1 ', ig, iispinor, igspinor,ghcre,&
1398 0 : kinpw_k2(ig),cwavef(1,igspinor),ghc(1,igspinor),gvnlxc_(1,igspinor), gsc(1,igspinor)
1399 0 : call wrtout(std_out,msg,'PERS')
1400 0 : write(msg,'(a,3(1x,i5),6(1x,es13.6))') ' 2 ', ig, iispinor, igspinor,ghcim,&
1401 0 : kinpw_k2(ig),cwavef(2,igspinor),ghc(2,igspinor),gvnlxc_(2,igspinor), gsc(2,igspinor)
1402 0 : call wrtout(std_out,msg,'PERS')
1403 : else
1404 0 : write(msg,'(a,3(1x,i5),6(1x,es13.6))') ' 1 ', ig, iispinor, igspinor,ghcre,&
1405 0 : kinpw_k2(ig),cwavef(1,igspinor),ghc(1,igspinor),gvnlxc_(1,igspinor)
1406 0 : call wrtout(std_out,msg,'PERS')
1407 0 : write(msg,'(a,3(1x,i5),6(1x,es13.6))') ' 2 ', ig, iispinor, igspinor,ghcim,&
1408 0 : kinpw_k2(ig),cwavef(2,igspinor),ghc(2,igspinor),gvnlxc_(2,igspinor)
1409 0 : call wrtout(std_out,msg,'PERS')
1410 : end if
1411 0 : ghc(:,igspinor) = [ghcre, ghcim]
1412 : end do ! ig
1413 : end do ! ispinor
1414 : end do ! idat
1415 : #ifdef HAVE_OPENMP_OFFLOAD
1416 : !$OMP TARGET UPDATE TO(ghc,gsc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1417 : #endif
1418 : end if
1419 : end if ! gs_ham%gpu_option
1420 :
1421 : ABI_NVTX_END_RANGE()
1422 :
1423 : ! Special case of PAW + Fock : only return Fock operator contribution in gvnlxc_
1424 34757726 : if (gs_ham%usepaw==1 .and. has_fock) then
1425 13114 : if(gs_ham%gpu_option == ABI_GPU_OPENMP) then
1426 : #ifdef HAVE_OPENMP_OFFLOAD
1427 : call abi_xaxpy(2*npw_k2*my_nspinor*ndat, cminusone, gvnlc, 1, gvnlxc_, 1, x_cplx=1, gpu_option=gs_ham%gpu_option)
1428 : !$OMP TARGET EXIT DATA MAP(delete:gvnlc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1429 : #endif
1430 : else
1431 48872050 : gvnlxc_=gvnlxc_-gvnlc
1432 : end if
1433 13114 : if(gs_ham%gpu_option==ABI_GPU_KOKKOS) then
1434 : #if defined HAVE_GPU && defined HAVE_YAKL
1435 : ABI_FREE_MANAGED(gvnlc)
1436 : #endif
1437 : else
1438 13114 : ABI_FREE(gvnlc)
1439 : end if
1440 : endif
1441 :
1442 : #ifdef HAVE_OPENMP_OFFLOAD
1443 : !$OMP TARGET UPDATE FROM(ghc) IF(transfer_ghc)
1444 : !$OMP TARGET UPDATE FROM(gsc) IF(transfer_gsc)
1445 : !$OMP TARGET UPDATE FROM(cwavef) IF(transfer_cwavef)
1446 : !$OMP TARGET UPDATE FROM(gvnlxc_) IF(transfer_gvnlxc .and. .not. local_gvnlxc)
1447 : !$OMP TARGET EXIT DATA MAP(delete:gvnlxc_) IF(transfer_gvnlxc)
1448 : #endif
1449 34757726 : if (local_gvnlxc) then
1450 880108 : if(gs_ham%gpu_option==ABI_GPU_KOKKOS) then
1451 : #if defined HAVE_GPU && defined HAVE_YAKL
1452 : ABI_FREE_MANAGED(gvnlxc_)
1453 : #endif
1454 : else
1455 880108 : ABI_FREE(gvnlxc_)
1456 : end if
1457 : end if
1458 :
1459 : ! Structured debugging: if prtvol=-level, stop here.
1460 34757726 : if (prtvol == -level) then
1461 0 : ABI_ERROR(sjoin(' getghc: exit prtvol=-',itoa(level),', debugging mode => stop '))
1462 : end if
1463 :
1464 34757726 : if (type_calc==0.or.type_calc==2) then
1465 34757726 : if (has_fock.and.gs_ham%usepaw==1.and.cpopt<2) then
1466 13114 : call pawcprj_free(cwaveprj_fock)
1467 13114 : ABI_FREE(cwaveprj_fock)
1468 : end if
1469 : end if
1470 :
1471 : end if ! type_calc
1472 :
1473 : #ifdef HAVE_OPENMP_OFFLOAD
1474 : !$OMP TARGET EXIT DATA MAP(delete:ghc) IF(transfer_ghc)
1475 : !$OMP TARGET EXIT DATA MAP(delete:gsc) IF(transfer_gsc)
1476 : !$OMP TARGET EXIT DATA MAP(delete:cwavef) IF(transfer_cwavef)
1477 : #endif
1478 34965650 : call timab(350+tim_getghc,2,tsec)
1479 :
1480 : DBG_EXIT("COLL")
1481 : ABI_NVTX_END_RANGE()
1482 :
1483 69931300 : end subroutine getghc
1484 : !!***
1485 :
1486 : !----------------------------------------------------------------------
1487 :
1488 : !!****f* ABINIT/getghc_nucdip
1489 : !!
1490 : !! NAME
1491 : !! getghc_nucdip
1492 : !!
1493 : !! FUNCTION
1494 : !! Compute magnetic nuclear dipole moment contribution to <G|H|C>
1495 : !! for input vector |C> expressed in reciprocal space.
1496 : !!
1497 : !! INPUTS
1498 : !! cwavef(2,npw_k*my_nspinor*ndat)=planewave coefficients of wavefunction.
1499 : !! gbound_k(2*mgfft+4)=sphere boundary info
1500 : !! gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1)
1501 : !! istwf_k=input parameter that describes the storage of wfs
1502 : !! kg_k(3,npw_k)=G vec coordinates wrt recip lattice transl.
1503 : !! kpt(3)=current k point
1504 : !! mgfft=maximum single fft dimension
1505 : !! mpi_enreg=information about MPI parallelization
1506 : !! my_nspinor=number of spinorial components of the wavefunctions (on current proc)
1507 : !! ndat=number of FFTs to perform in parall
1508 : !! ngfft(18)=contain all needed information about 3D FFT
1509 : !! npw_k=number of planewaves in basis for given k point.
1510 : !! nvloc=number of spin components of vxctaulocal
1511 : !! n4,n5,n6=for dimensioning of vxctaulocal
1512 : !! gpu_option= GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
1513 : !! vectornd(n4,n5,n6,nvloc,3)= local potential corresponding to the vector potential of the array
1514 : !! of nuclear magnetic dipoles, in real space, on the augmented fft grid.
1515 : !!
1516 : !! OUTPUT
1517 : !! ghc_vectornd(2,npw_k*my_nspinor*ndat)=A.p contribution to <G|H|C> for array of nuclear dipoles
1518 : !!
1519 : !! SIDE EFFECTS
1520 : !!
1521 : !! NOTES
1522 : !! this code is a copied, simplified version of getghc_mGGA (see below) and should eventually be
1523 : !! integrated into that code, to simplify maintenance
1524 : !!
1525 : !! SOURCE
1526 :
1527 812570 : subroutine getghc_nucdip(cwavef,ghc_vectornd,gbound_k,istwf_k,kg_k,kpt,mgfft,mpi_enreg,&
1528 812570 : & ndat,ngfft,npw_k,nvloc,n4,n5,n6,my_nspinor,vectornd,vlocal,zora,gpu_option)
1529 :
1530 : !Arguments ------------------------------------
1531 : !scalars
1532 : integer,intent(in) :: istwf_k,mgfft,my_nspinor,ndat,npw_k,nvloc,n4,n5,n6,gpu_option,zora
1533 : type(MPI_type),intent(in) :: mpi_enreg
1534 : !arrays
1535 : integer,intent(in) :: gbound_k(2*mgfft+4),kg_k(3,npw_k),ngfft(18)
1536 : real(dp),intent(in) :: kpt(3)
1537 : real(dp),intent(inout) :: cwavef(2,npw_k*my_nspinor*ndat)
1538 : real(dp),intent(inout) :: ghc_vectornd(2,npw_k*my_nspinor*ndat)
1539 : real(dp),intent(inout) :: vectornd(n4,n5,n6,nvloc,3),vlocal(n4,n5,n6,nvloc)
1540 :
1541 : !Local variables-------------------------------
1542 : !scalars
1543 : integer,parameter :: tim_fourwf=1
1544 : integer :: idat,idir,ipw,iv1,iv2,nspinortot,shift
1545 : logical :: nspinor1TreatedByThisProc,nspinor2TreatedByThisProc,usezora
1546 : real(dp), parameter :: HalfFineStruct2=half/InvFineStruct**2
1547 : real(dp) :: weight=one
1548 : !arrays
1549 812570 : real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:)
1550 812570 : real(dp),allocatable :: gcwavef(:,:,:),gcwavef1(:,:,:),gcwavef2(:,:,:)
1551 812570 : real(dp),allocatable :: ghc1(:,:),ghc2(:,:),kgkpk(:,:),vectornd_dir(:,:,:,:)
1552 812570 : real(dp),allocatable :: work(:,:,:,:),zk(:,:,:,:)
1553 : ! *********************************************************************
1554 :
1555 1127843222 : ghc_vectornd(:,:)=zero
1556 :
1557 : !! JWZ debug initial code was only for nvloc==1 case
1558 812570 : if (nvloc/=1) return
1559 :
1560 812570 : nspinortot=min(2,(1+mpi_enreg%paral_spinor)*my_nspinor)
1561 812570 : if (mpi_enreg%paral_spinor==0) then
1562 812570 : shift=npw_k
1563 812570 : nspinor1TreatedByThisProc=.true.
1564 812570 : nspinor2TreatedByThisProc=(nspinortot==2)
1565 : else
1566 0 : shift=0
1567 0 : nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
1568 0 : nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
1569 : end if
1570 :
1571 812570 : usezora=((zora.EQ.1).OR.(zora.EQ.3))
1572 812570 : if(usezora) then
1573 0 : ABI_MALLOC(zk,(n4,n5,n6,nvloc))
1574 0 : zk(1:n4,1:1:n5,1:n6,1:nvloc)=1.0/(1.0-HalfFineStruct2*vlocal(1:n4,1:n5,1:n6,1:nvloc))
1575 : end if
1576 :
1577 4062850 : ABI_MALLOC(work,(2,n4,n5,n6*ndat))
1578 :
1579 812570 : if (nspinortot==1) then
1580 :
1581 1754592 : ABI_MALLOC(ghc1,(2,npw_k*ndat))
1582 :
1583 : ! Do it in 2 STEPs:
1584 : ! STEP1: Compute grad of cwavef
1585 2339456 : ABI_MALLOC(gcwavef,(2,npw_k*ndat,3))
1586 :
1587 1057250534 : gcwavef = zero
1588 :
1589 : ! compute k + G. Note these are in reduced coords
1590 1754592 : ABI_MALLOC(kgkpk,(npw_k,3))
1591 117797206 : do ipw = 1, npw_k
1592 469434232 : kgkpk(ipw,:) = kpt(:) + kg_k(:,ipw)
1593 : end do
1594 :
1595 : ! make 2\pi(k+G)c(G)|G> by element-wise multiplication
1596 2339456 : do idir = 1, 3
1597 4094048 : do idat = 1, ndat
1598 355146210 : do ipw=1,npw_k
1599 351637026 : gcwavef(1,ipw+(idat-1)*npw_k,idir) = cwavef(1,ipw+(idat-1)*npw_k)*kgkpk(ipw,idir)
1600 353391618 : gcwavef(2,ipw+(idat-1)*npw_k,idir) = cwavef(2,ipw+(idat-1)*npw_k)*kgkpk(ipw,idir)
1601 : end do
1602 : end do
1603 : end do
1604 584864 : ABI_FREE(kgkpk)
1605 1057250534 : gcwavef = gcwavef*two_pi
1606 :
1607 : ! STEP2: Compute sum of (grad components of vectornd)*(grad components of cwavef)
1608 3509184 : ABI_MALLOC(vectornd_dir,(n4,n5,n6,nvloc))
1609 2339456 : do idir=1,3
1610 1754592 : if (usezora) then
1611 0 : vectornd_dir(1:n4,1:n5,1:n6,1:nvloc)=zk(1:n4,1:n5,1:n6,1:nvloc)*vectornd(1:n4,1:n5,1:n6,1:nvloc,idir)
1612 : else
1613 10314799284 : vectornd_dir(1:n4,1:n5,1:n6,1:nvloc)=vectornd(1:n4,1:n5,1:n6,1:nvloc,idir)
1614 : end if
1615 : call fourwf(1,vectornd_dir,gcwavef(:,:,idir),ghc1,work,gbound_k,gbound_k,&
1616 : istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1617 1754592 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1618 : !!$OMP PARALLEL DO
1619 4094048 : do idat=1,ndat
1620 355146210 : do ipw=1,npw_k
1621 1056665670 : ghc_vectornd(:,ipw+(idat-1)*npw_k)=ghc_vectornd(:,ipw+(idat-1)*npw_k)+ghc1(:,ipw+(idat-1)*npw_k)
1622 : end do
1623 : end do
1624 : !ghc_vectornd=ghc_vectornd+ghc1
1625 : end do ! idir
1626 584864 : ABI_FREE(vectornd_dir)
1627 584864 : ABI_FREE(gcwavef)
1628 584864 : ABI_FREE(ghc1)
1629 :
1630 : else ! nspinortot==2
1631 :
1632 683118 : ABI_MALLOC(cwavef1,(2,npw_k*ndat))
1633 455412 : ABI_MALLOC(cwavef2,(2,npw_k*ndat))
1634 455412 : do idat=1,ndat
1635 227706 : iv1=1+(idat-1)*npw_k; iv2=-1+iv1+npw_k
1636 387924519 : cwavef1(1:2,iv1:iv2) = cwavef(1:2,1+(idat-1)*my_nspinor*npw_k:npw_k+(idat-1)*my_nspinor*npw_k)
1637 : cwavef2(1:2,iv1:iv2) = &
1638 388152225 : & cwavef(1:2,1+(idat-1)*my_nspinor*npw_k+shift:npw_k+(idat-1)*my_nspinor*npw_k+shift)
1639 : end do
1640 :
1641 : ! compute k + G. Note these are in reduced coords
1642 683118 : ABI_MALLOC(kgkpk,(npw_k,3))
1643 129459977 : do ipw = 1, npw_k
1644 517156790 : kgkpk(ipw,:) = kpt(:) + kg_k(:,ipw)
1645 : end do
1646 :
1647 227706 : if (nspinor1TreatedByThisProc) then
1648 :
1649 455412 : ABI_MALLOC(ghc1,(2,npw_k*ndat))
1650 :
1651 : ! Do it in 2 STEPs:
1652 : ! STEP1: Compute grad of cwavef
1653 910824 : ABI_MALLOC(gcwavef1,(2,npw_k*ndat,3))
1654 :
1655 1164001263 : gcwavef1 = zero
1656 : ! make 2\pi(k+G)c(G)|G> by element-wise multiplication
1657 910824 : do idir = 1, 3
1658 1593942 : do idat = 1, ndat
1659 683118 : iv1=1+(idat-1)*npw_k; iv2=-1+iv1+npw_k
1660 388379931 : gcwavef1(1,iv1:iv2,idir) = cwavef1(1,iv1:iv2)*kgkpk(1:npw_k,idir)
1661 389063049 : gcwavef1(2,iv1:iv2,idir) = cwavef1(2,iv1:iv2)*kgkpk(1:npw_k,idir)
1662 : end do
1663 : end do
1664 1164001263 : gcwavef1 = gcwavef1*two_pi
1665 :
1666 : ! STEP2: Compute sum of (grad components of vectornd)*(grad components of cwavef)
1667 1366236 : ABI_MALLOC(vectornd_dir,(n4,n5,n6,nvloc))
1668 910824 : do idir=1,3
1669 683118 : if (usezora) then
1670 0 : vectornd_dir(1:n4,1:n5,1:n6,1:nvloc)=zk(1:n4,1:n5,1:n6,1:nvloc)*vectornd(1:n4,1:n5,1:n6,1:nvloc,idir)
1671 : else
1672 8256805533 : vectornd_dir(1:n4,1:n5,1:n6,1:nvloc)=vectornd(1:n4,1:n5,1:n6,1:nvloc,idir)
1673 : end if
1674 : call fourwf(1,vectornd_dir,gcwavef1(:,:,idir),ghc1,work,gbound_k,gbound_k,&
1675 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1676 683118 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1677 : !!$OMP PARALLEL DO
1678 1593942 : do idat=1,ndat
1679 683118 : iv1=1+(idat-1)*npw_k; iv2=-1+iv1+npw_k
1680 : ghc_vectornd(1:2,iv1:iv2)=ghc_vectornd(1:2,iv1:iv2)+&
1681 1164456675 : & ghc1(1:2,iv1:iv2)
1682 : end do
1683 : end do ! idir
1684 227706 : ABI_FREE(vectornd_dir)
1685 227706 : ABI_FREE(gcwavef1)
1686 227706 : ABI_FREE(ghc1)
1687 :
1688 : end if ! end spinor 1
1689 :
1690 227706 : if (nspinor2TreatedByThisProc) then
1691 :
1692 455412 : ABI_MALLOC(ghc2,(2,npw_k*ndat))
1693 :
1694 : ! Do it in 2 STEPs:
1695 : ! STEP1: Compute grad of cwavef
1696 910824 : ABI_MALLOC(gcwavef2,(2,npw_k*ndat,3))
1697 1164001263 : gcwavef2 = zero
1698 : ! make 2\pi(k+G)c(G)|G> by element-wise multiplication
1699 910824 : do idir = 1, 3
1700 1593942 : do idat = 1, ndat
1701 683118 : iv1=1+(idat-1)*npw_k; iv2=-1+iv1+npw_k
1702 388379931 : gcwavef2(1,iv1:iv2,idir) = cwavef2(1,iv1:iv2)*kgkpk(1:npw_k,idir)
1703 389063049 : gcwavef2(2,iv1:iv2,idir) = cwavef2(2,iv1:iv2)*kgkpk(1:npw_k,idir)
1704 : end do
1705 : end do
1706 1164001263 : gcwavef2 = gcwavef2*two_pi
1707 :
1708 : ! STEP2: Compute sum of (grad components of vectornd)*(grad components of cwavef)
1709 1366236 : ABI_MALLOC(vectornd_dir,(n4,n5,n6,nvloc))
1710 910824 : do idir=1,3
1711 683118 : if (usezora) then
1712 0 : vectornd_dir(1:n4,1:n5,1:n6,1:nvloc)=zk(1:n4,1:n5,1:n6,1:nvloc)*vectornd(1:n4,1:n5,1:n6,1:nvloc,idir)
1713 : else
1714 8256805533 : vectornd_dir(1:n4,1:n5,1:n6,1:nvloc)=vectornd(1:n4,1:n5,1:n6,1:nvloc,idir)
1715 : end if
1716 : call fourwf(1,vectornd_dir,gcwavef2(:,:,idir),ghc2,work,gbound_k,gbound_k,&
1717 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1718 683118 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1719 : !!$OMP PARALLEL DO
1720 1593942 : do idat=1,ndat
1721 683118 : iv1=1+(idat-1)*npw_k; iv2=-1+iv1+npw_k
1722 : ghc_vectornd(1:2,iv1+shift:iv2+shift)=ghc_vectornd(1:2,iv1+shift:iv2+shift)+&
1723 1164456675 : & ghc2(1:2,iv1:iv2)
1724 : end do
1725 : end do ! idir
1726 227706 : ABI_FREE(vectornd_dir)
1727 227706 : ABI_FREE(gcwavef2)
1728 227706 : ABI_FREE(ghc2)
1729 :
1730 : end if ! end spinor 2
1731 :
1732 227706 : ABI_FREE(cwavef1)
1733 227706 : ABI_FREE(cwavef2)
1734 227706 : ABI_FREE(kgkpk)
1735 :
1736 : end if ! nspinortot
1737 :
1738 812570 : ABI_FREE(work)
1739 812570 : if(usezora) then
1740 0 : ABI_FREE(zk)
1741 : end if
1742 :
1743 812570 : end subroutine getghc_nucdip
1744 : !!***
1745 :
1746 : !!****f* ABINIT/getghc_mGGA
1747 : !!
1748 : !! NAME
1749 : !! getghc_mGGA
1750 : !!
1751 : !! FUNCTION
1752 : !! Compute metaGGA contribution to <G|H|C> for input vector |C> expressed in reciprocal space.
1753 : !!
1754 : !! INPUTS
1755 : !! cwavef(2,npw_k*my_nspinor*ndat)=planewave coefficients of wavefunction.
1756 : !! gbound_k(2*mgfft+4)=sphere boundary info
1757 : !! gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1)
1758 : !! istwf_k=input parameter that describes the storage of wfs
1759 : !! kg_k(3,npw_k)=G vec coordinates wrt recip lattice transl.
1760 : !! kpt(3)=current k point
1761 : !! mgfft=maximum single fft dimension
1762 : !! mpi_enreg=information about MPI parallelization
1763 : !! my_nspinor=number of spinorial components of the wavefunctions (on current proc)
1764 : !! ndat=number of FFTs to perform in parall
1765 : !! ngfft(18)=contain all needed information about 3D FFT
1766 : !! npw_k=number of planewaves in basis for given k point.
1767 : !! nvloc=number of spin components of vxctaulocal
1768 : !! n4,n5,n6=for dimensionning of vxctaulocal
1769 : !! gpu_option= GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
1770 : !! vxctaulocal(n4,n5,n6,nvloc,4)= local potential corresponding to the derivative of XC energy with respect to
1771 : !! kinetic energy density, in real space, on the augmented fft grid.
1772 : !! This array contains also the gradient of vxctaulocal (gvxctaulocal) in vxctaulocal(:,:,:,:,2:4).
1773 : !!
1774 : !! OUTPUT
1775 : !! ghc_mGGA(2,npw_k*my_nspinor*ndat)=metaGGA contribution to <G|H|C>
1776 : !!
1777 : !! SIDE EFFECTS
1778 : !!
1779 : !! SOURCE
1780 :
1781 98528 : subroutine getghc_mGGA(cwavef,ghc_mGGA,gbound_k,gprimd,istwf_k,kg_k,kpt,mgfft,mpi_enreg,&
1782 98528 : & ndat,ngfft,npw_k,nvloc,n4,n5,n6,my_nspinor,vxctaulocal,gpu_option)
1783 :
1784 : !Arguments ------------------------------------
1785 : !scalars
1786 : integer,intent(in) :: istwf_k,mgfft,my_nspinor,ndat,npw_k,nvloc,n4,n5,n6,gpu_option
1787 : type(MPI_type),intent(in) :: mpi_enreg
1788 : !arrays
1789 : integer,intent(in) :: gbound_k(2*mgfft+4),kg_k(3,npw_k),ngfft(18)
1790 : real(dp),intent(in) :: gprimd(3,3),kpt(3)
1791 : real(dp),intent(inout) :: cwavef(2,npw_k*my_nspinor*ndat)
1792 : real(dp),intent(inout) :: ghc_mGGA(2,npw_k*my_nspinor*ndat)
1793 : real(dp),intent(inout) :: vxctaulocal(n4,n5,n6,nvloc,4)
1794 :
1795 : !Local variables-------------------------------
1796 : !scalars
1797 : integer,parameter :: tim_fourwf=1
1798 : integer :: idat,idir,ipw,nspinortot,shift
1799 : logical :: nspinor1TreatedByThisProc,nspinor2TreatedByThisProc
1800 : real(dp) :: weight=one
1801 : !arrays
1802 : real(dp) :: kg_k_cart_vec(3)
1803 98528 : real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:)
1804 98528 : real(dp),allocatable :: gcwavef(:,:,:),gcwavef1(:,:,:),gcwavef2(:,:,:)
1805 98528 : real(dp),allocatable :: ghc1(:,:),ghc2(:,:)
1806 98528 : real(dp),allocatable :: lcwavef(:,:),lcwavef1(:,:),lcwavef2(:,:)
1807 98528 : real(dp),allocatable :: work(:,:,:,:)
1808 : ! *********************************************************************
1809 :
1810 176135072 : ghc_mGGA(:,:)=zero
1811 :
1812 98528 : if (nvloc/=1) return
1813 :
1814 98528 : nspinortot=min(2,(1+mpi_enreg%paral_spinor)*my_nspinor)
1815 98528 : if (mpi_enreg%paral_spinor==0) then
1816 98528 : shift=npw_k
1817 98528 : nspinor1TreatedByThisProc=.true.
1818 98528 : nspinor2TreatedByThisProc=(nspinortot==2)
1819 : else
1820 0 : shift=0
1821 0 : nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
1822 0 : nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
1823 : end if
1824 :
1825 492640 : ABI_MALLOC(work,(2,n4,n5,n6*ndat))
1826 :
1827 98528 : if (nspinortot==1) then
1828 :
1829 136392 : ABI_MALLOC(ghc1,(2,npw_k*ndat))
1830 :
1831 : ! Do it in 3 STEPs:
1832 : ! STEP1: Compute grad of cwavef and Laplacian of cwavef
1833 181856 : ABI_MALLOC(gcwavef,(2,npw_k*ndat,3))
1834 90928 : ABI_MALLOC(lcwavef,(2,npw_k*ndat))
1835 : !!$OMP PARALLEL DO
1836 387384872 : gcwavef = zero; lcwavef = zero
1837 95408 : do idat=1,ndat
1838 32362326 : do ipw=1,npw_k
1839 : ! convert k + G from reduced coords to Cartesian
1840 613071442 : kg_k_cart_vec = two_pi*MATMUL(gprimd,kpt(1:3)+kg_k(1:3,ipw))
1841 : ! form \grad\psi = i(k + G) \psi in Cartesian frame
1842 129067672 : gcwavef(1,ipw+(idat-1)*npw_k,1:3)= cwavef(2,ipw+(idat-1)*npw_k)*kg_k_cart_vec(1:3)
1843 129067672 : gcwavef(2,ipw+(idat-1)*npw_k,1:3)= -cwavef(1,ipw+(idat-1)*npw_k)*kg_k_cart_vec(1:3)
1844 : ! form \nabla^2\psi = -|k + G|^2 \psi
1845 : lcwavef(1:2,ipw+(idat-1)*npw_k)=&
1846 193651452 : &lcwavef(1:2,ipw+(idat-1)*npw_k)-cwavef(1:2,ipw+(idat-1)*npw_k)*DOT_PRODUCT(kg_k_cart_vec,kg_k_cart_vec)
1847 : end do
1848 : end do
1849 : ! STEP2: Compute (vxctaulocal)*(Laplacian of cwavef) and add it to ghc
1850 : call fourwf(1,vxctaulocal(:,:,:,:,1),lcwavef,ghc1,work,gbound_k,gbound_k,&
1851 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1852 45464 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1853 : !!$OMP PARALLEL DO
1854 95408 : do idat=1,ndat
1855 32362326 : do ipw=1,npw_k
1856 96850698 : ghc_mGGA(:,ipw+(idat-1)*npw_k)=ghc_mGGA(:,ipw+(idat-1)*npw_k)-half*ghc1(:,ipw+(idat-1)*npw_k)
1857 : end do
1858 : end do
1859 45464 : ABI_FREE(lcwavef)
1860 : ! STEP3: Compute sum of (grad components of vxctaulocal)*(grad components of cwavef)
1861 : ! note: since grad cwavef is in Cart frame, evidently grad vxc is also
1862 181856 : do idir=1,3
1863 : call fourwf(1,vxctaulocal(:,:,:,:,1+idir),gcwavef(:,:,idir),ghc1,work,gbound_k,gbound_k,&
1864 : istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1865 136392 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1866 : !!$OMP PARALLEL DO
1867 331688 : do idat=1,ndat
1868 97086978 : do ipw=1,npw_k
1869 290552094 : ghc_mGGA(:,ipw+(idat-1)*npw_k)=ghc_mGGA(:,ipw+(idat-1)*npw_k)-half*ghc1(:,ipw+(idat-1)*npw_k)
1870 : end do
1871 : end do
1872 : end do ! idir
1873 45464 : ABI_FREE(gcwavef)
1874 45464 : ABI_FREE(ghc1)
1875 :
1876 : else ! nspinortot==2
1877 :
1878 159192 : ABI_MALLOC(cwavef1,(2,npw_k*ndat))
1879 106128 : ABI_MALLOC(cwavef2,(2,npw_k*ndat))
1880 109248 : do idat=1,ndat
1881 13315213 : do ipw=1,npw_k
1882 39617895 : cwavef1(1:2,ipw+(idat-1)*npw_k)=cwavef(1:2,ipw+(idat-1)*my_nspinor*npw_k)
1883 39674079 : cwavef2(1:2,ipw+(idat-1)*npw_k)=cwavef(1:2,ipw+(idat-1)*my_nspinor*npw_k+shift)
1884 : end do
1885 : end do
1886 : ! call cg_zcopy(npw*ndat,cwavef(1,1),cwavef1)
1887 : ! call cg_zcopy(npw*ndat,cwavef(1,1+shift),cwavef2)
1888 :
1889 :
1890 53064 : if (nspinor1TreatedByThisProc) then
1891 :
1892 106128 : ABI_MALLOC(ghc1,(2,npw_k*ndat))
1893 :
1894 : ! Do it in 3 STEPs:
1895 : ! STEP1: Compute grad of cwavef and Laplacian of cwavef
1896 212256 : ABI_MALLOC(gcwavef1,(2,npw_k*ndat,3))
1897 106128 : ABI_MALLOC(lcwavef1,(2,npw_k*ndat))
1898 158683836 : gcwavef1 = zero; lcwavef1 = zero
1899 : !!$OMP PARALLEL DO
1900 109248 : do idat=1,ndat
1901 13315213 : do ipw=1,npw_k
1902 : ! convert k + G from reduced coords to Cartesian
1903 250913335 : kg_k_cart_vec = two_pi*MATMUL(gprimd,kpt(1:3)+kg_k(1:3,ipw))
1904 : ! form \grad\psi = i(k + G) \psi in Cartesian frame
1905 52823860 : gcwavef1(1,ipw+(idat-1)*npw_k,1:3)= cwavef1(2,ipw+(idat-1)*npw_k)*kg_k_cart_vec(1:3)
1906 52823860 : gcwavef1(2,ipw+(idat-1)*npw_k,1:3)= -cwavef1(1,ipw+(idat-1)*npw_k)*kg_k_cart_vec(1:3)
1907 : ! form \nabla^2\psi = -|k + G|^2 \psi
1908 : lcwavef1(1:2,ipw+(idat-1)*npw_k)=&
1909 79291974 : &lcwavef1(1:2,ipw+(idat-1)*npw_k)-cwavef1(1:2,ipw+(idat-1)*npw_k)*DOT_PRODUCT(kg_k_cart_vec,kg_k_cart_vec)
1910 : end do
1911 : end do
1912 : ! STEP2: Compute (vxctaulocal)*(Laplacian of cwavef) and add it to ghc
1913 : call fourwf(1,vxctaulocal(:,:,:,:,1),lcwavef1,ghc1,work,gbound_k,gbound_k,&
1914 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1915 53064 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1916 : !!$OMP PARALLEL DO
1917 109248 : do idat=1,ndat
1918 13315213 : do ipw=1,npw_k
1919 39674079 : ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)=ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)-half*ghc1(:,ipw+(idat-1)*npw_k)
1920 : end do
1921 : end do
1922 53064 : ABI_FREE(lcwavef1)
1923 : ! STEP3: Compute (grad components of vxctaulocal)*(grad components of cwavef)
1924 212256 : do idir=1,3
1925 : call fourwf(1,vxctaulocal(:,:,:,:,1+idir),gcwavef1(:,:,idir),ghc1,work,gbound_k,gbound_k,&
1926 : istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1927 159192 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1928 : !!$OMP PARALLEL DO
1929 380808 : do idat=1,ndat
1930 39945639 : do ipw=1,npw_k
1931 119022237 : ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k) = ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)-half*ghc1(:,ipw+(idat-1)*npw_k)
1932 : end do
1933 : end do
1934 : end do ! idir
1935 53064 : ABI_FREE(gcwavef1)
1936 53064 : ABI_FREE(ghc1)
1937 :
1938 : end if ! spin 1 treated by this proc
1939 :
1940 53064 : if (nspinor2TreatedByThisProc) then
1941 :
1942 106128 : ABI_MALLOC(ghc2,(2,npw_k*ndat))
1943 :
1944 : ! Do it in 3 STEPs:
1945 : ! STEP1: Compute grad of cwavef and Laplacian of cwavef
1946 212256 : ABI_MALLOC(gcwavef2,(2,npw_k*ndat,3))
1947 106128 : ABI_MALLOC(lcwavef2,(2,npw_k*ndat))
1948 : !!$OMP PARALLEL DO
1949 158683836 : gcwavef2 = zero; lcwavef2 = zero
1950 109248 : do idat=1,ndat
1951 13315213 : do ipw=1,npw_k
1952 : ! convert k + G from reduced coords to Cartesian
1953 250913335 : kg_k_cart_vec = two_pi*MATMUL(gprimd,kpt(1:3)+kg_k(1:3,ipw))
1954 : ! form \grad\psi = i(k + G) \psi in Cartesian frame
1955 52823860 : gcwavef2(1,ipw+(idat-1)*npw_k,1:3)= cwavef2(2,ipw+(idat-1)*npw_k)*kg_k_cart_vec(1:3)
1956 52823860 : gcwavef2(2,ipw+(idat-1)*npw_k,1:3)= -cwavef2(1,ipw+(idat-1)*npw_k)*kg_k_cart_vec(1:3)
1957 : ! form \nabla^2\psi = -|k + G|^2 \psi
1958 : lcwavef2(1:2,ipw+(idat-1)*npw_k)=&
1959 79291974 : &lcwavef2(1:2,ipw+(idat-1)*npw_k)-cwavef2(1:2,ipw+(idat-1)*npw_k)*DOT_PRODUCT(kg_k_cart_vec,kg_k_cart_vec)
1960 : end do
1961 : end do
1962 : ! STEP2: Compute (vxctaulocal)*(Laplacian of cwavef) and add it to ghc
1963 : call fourwf(1,vxctaulocal(:,:,:,:,1),lcwavef2,ghc2,work,gbound_k,gbound_k,&
1964 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1965 53064 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1966 : !!$OMP PARALLEL DO
1967 109248 : do idat=1,ndat
1968 13315213 : do ipw=1,npw_k
1969 : ! original code
1970 : ! ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)=ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)-half*ghc2(:,ipw+(idat-1)*npw_k)
1971 : ! but this stores the spinor2 result in the spinor1 location. Should be stored with shift
1972 : ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k+shift)=ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k+shift)&
1973 39674079 : & -half*ghc2(:,ipw+(idat-1)*npw_k)
1974 : end do
1975 : end do
1976 53064 : ABI_FREE(lcwavef2)
1977 : ! STEP3: Compute sum of (grad components of vxctaulocal)*(grad components of cwavef)
1978 212256 : do idir=1,3
1979 : call fourwf(1,vxctaulocal(:,:,:,:,1+idir),gcwavef2(:,:,idir),ghc2,work,gbound_k,gbound_k,&
1980 : istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1981 159192 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1982 : !!$OMP PARALLEL DO
1983 380808 : do idat=1,ndat
1984 39945639 : do ipw=1,npw_k
1985 : ! original code
1986 : ! ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)=ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)-half*ghc2(:,ipw+(idat-1)*npw_k)
1987 : ! but this stores the spinor2 result in the spinor1 location. Should be stored with shift
1988 : ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k+shift)=ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k+shift)&
1989 119022237 : & -half*ghc2(:,ipw+(idat-1)*npw_k)
1990 : end do
1991 : end do
1992 : end do ! idir
1993 :
1994 53064 : ABI_FREE(gcwavef2)
1995 53064 : ABI_FREE(ghc2)
1996 :
1997 : end if ! spin 2 treated by this proc
1998 :
1999 53064 : ABI_FREE(cwavef1)
2000 53064 : ABI_FREE(cwavef2)
2001 :
2002 : end if ! nspinortot
2003 :
2004 98528 : ABI_FREE(work)
2005 :
2006 : end subroutine getghc_mGGA
2007 : !!***
2008 :
2009 : !!****f* ABINIT/cwavef_double_rfft_trick_pack
2010 : !!
2011 : !! NAME
2012 : !! cwavef_double_rfft_trick_pack
2013 : !!
2014 : !! FUNCTION
2015 : !!
2016 : !! We have C(G)=C(-G)^* and D(G)=D(-G)^* and only G components are in memory (istwfk=2)
2017 : !! The Fourier transform is:
2018 : !! C(r) = sum_G e^(iGr) C(G) = sum_(G_z>=0,G/=0) 2 Re[e^(iGr) C(G)] + C(0)
2019 : !! so C(r) is a real function (same for D)
2020 : !! Here we construct:
2021 : !! E( G) = C(G) + i D(G)
2022 : !! E(-G) = C(G)^* + i D(G)^* (G/=0)
2023 : !! so:
2024 : !! E(r) = C(r) + i D(r)
2025 : !! In short, one can do only one FFT on E (with istwfk=1) and obtains the FFT of C and D (istwfk=2)
2026 : !!
2027 : !! INPUTS
2028 : !! cwavef(2,npw_k*my_nspinor*(2*ndat))=planewave coefficients of wavefunctioni (istwfk=2)
2029 : !! mpi_enreg=information about MPI parallelization
2030 : !! ndat=number of FFTs to perform in parall
2031 : !! npw_k=number of planewaves in basis for given k point (istwfk=2).
2032 : !!
2033 : !! OUTPUT
2034 : !! cwavef_fft(2,npw_fft*my_nspinor*ndat)=planewave coefficients of wavefunction (istwfk=1)
2035 : !!
2036 : !! SOURCE
2037 : !!
2038 16912 : subroutine cwavef_double_rfft_trick_pack(cwavef,cwavef_fft,me_g0,ndat,npw_k)
2039 :
2040 : !Arguments ------------------------------------
2041 : !scalars
2042 : integer,intent(in) :: ndat,npw_k,me_g0
2043 : !arrays
2044 : real(dp),intent(in) :: cwavef(:,:)
2045 : real(dp),intent(out) :: cwavef_fft(:,:)
2046 : !Local variables-------------------------------
2047 : !scalars
2048 : integer :: idat,i0,ib1,ib2,npw_fft,ndat_
2049 : logical :: ndat_is_odd
2050 :
2051 16912 : ndat_ = ndat / 2
2052 16912 : ndat_is_odd=.False.
2053 16912 : if (modulo(ndat,2)/=0) then
2054 0 : ndat_=ndat_+1
2055 0 : ndat_is_odd=.True.
2056 : end if
2057 :
2058 16912 : npw_fft=2*npw_k
2059 16912 : i0=1
2060 16912 : if (me_g0==1) then! Do not include G=(0,0,0) twice
2061 16912 : npw_fft=npw_fft-1
2062 16912 : i0=2
2063 : end if
2064 :
2065 16912 : if (size(cwavef,1)/=2) then
2066 0 : ABI_BUG('wrong size for cwavef (dim 1)')
2067 : end if
2068 16912 : if (ndat_is_odd) then
2069 0 : if (size(cwavef,2)/=npw_k*(2*ndat_-1)) then
2070 0 : ABI_BUG('wrong size for cwavef (dim 2) (odd)')
2071 : end if
2072 : else
2073 16912 : if (size(cwavef,2)/=npw_k*2*ndat_) then
2074 0 : ABI_BUG('wrong size for cwavef (dim 2)')
2075 : end if
2076 : end if
2077 16912 : if (size(cwavef_fft,1)/=2.or.size(cwavef_fft,2)/=npw_fft*ndat_) then
2078 0 : ABI_BUG('wrong size for cwavef_fft')
2079 : end if
2080 :
2081 106164 : do idat=1,ndat_
2082 89252 : ib1=(idat-1)*npw_fft ! band shift for cwavef_fft
2083 89252 : ib2=(idat-1)*2*npw_k ! band shift for cwavef
2084 106164 : if (.not.ndat_is_odd.or.idat<ndat_) then
2085 : ! E(G) = C(G) + i D(G)
2086 : cwavef_fft(1,1+ib1:npw_k+ib1) = cwavef(1,1+ib2 : npw_k+ib2) &
2087 68960932 : -cwavef(2,1+npw_k+ib2:2*npw_k+ib2)
2088 : cwavef_fft(2,1+ib1:npw_k+ib1) = cwavef(2,1+ib2 : npw_k+ib2) &
2089 68960932 : +cwavef(1,1+npw_k+ib2:2*npw_k+ib2)
2090 : ! E(-G) = C(G)^* + i D(G)^* (G/=0)
2091 : cwavef_fft(1,1+npw_k+ib1:npw_fft+ib1) = cwavef(1,i0 +ib2: npw_k+ib2) &
2092 68871680 : +cwavef(2,i0+npw_k+ib2:2*npw_k+ib2)
2093 : cwavef_fft(2,1+npw_k+ib1:npw_fft+ib1) =-cwavef(2,i0 +ib2: npw_k+ib2) &
2094 68871680 : +cwavef(1,i0+npw_k+ib2:2*npw_k+ib2)
2095 : else ! idat=ndat_ and ndat_is_odd : the vector D does not exist
2096 : ! E(G) = C(G)
2097 0 : cwavef_fft(1,1+ib1:npw_k+ib1) = cwavef(1,1+ib2:npw_k+ib2)
2098 0 : cwavef_fft(2,1+ib1:npw_k+ib1) = cwavef(2,1+ib2:npw_k+ib2)
2099 : ! E(-G) = C(G)^* (G/=0)
2100 0 : cwavef_fft(1,1+npw_k+ib1:npw_fft+ib1) = cwavef(1,i0+ib2:npw_k+ib2)
2101 0 : cwavef_fft(2,1+npw_k+ib1:npw_fft+ib1) =-cwavef(2,i0+ib2:npw_k+ib2)
2102 : end if
2103 : end do
2104 :
2105 16912 : end subroutine cwavef_double_rfft_trick_pack
2106 : !!***
2107 :
2108 : !!****f* ABINIT/cwavef_double_rfft_trick_unpack
2109 : !!
2110 : !! NAME
2111 : !! cwavef_double_rfft_trick_unpack
2112 : !!
2113 : !! FUNCTION
2114 : !!
2115 : !! From the "cwavef_double_rfft_trick_pack" routine we have:
2116 : !! E( G) = C(G) + i D(G)
2117 : !! E(-G) = C(G)^* + i D(G)^* (G/=0)
2118 : !! Here we compute:
2119 : !! C(G) = ( E(G) + E(-G)^* ) / 2 (G/=0)
2120 : !! D(G) = ( iE(-G)^* - iE(G)) / 2 (G/=0)
2121 : !! and:
2122 : !! C(0) = Re(E(0))
2123 : !! D(0) = Im(E(0))
2124 : !!
2125 : !! INPUTS
2126 : !! cwavef_fft(2,npw_fft*my_nspinor*ndat_)=planewave coefficients of wavefunction (istwfk=1)
2127 : !! mpi_enreg=information about MPI parallelization
2128 : !! ndat=number of FFTs to perform in parall
2129 : !! npw_k=number of planewaves in basis for given k point (istwfk=2).
2130 : !!
2131 : !! OUTPUT
2132 : !! cwavef(2,npw_fft*my_nspinor*(2*ndat))=planewave coefficients of wavefunction (istwfk=2)
2133 : !!
2134 : !! SOURCE
2135 : !!
2136 16912 : subroutine cwavef_double_rfft_trick_unpack(cwavef,cwavef_fft,me_g0,ndat,npw_k)
2137 :
2138 : !Arguments ------------------------------------
2139 : !scalars
2140 : integer,intent(in) :: ndat,npw_k,me_g0
2141 : !arrays
2142 : real(dp),intent(out) :: cwavef(:,:)
2143 : real(dp),intent(in) :: cwavef_fft(:,:)
2144 : !Local variables-------------------------------
2145 : !scalars
2146 : integer :: idat,i0,ib1,ib2,npw_fft,ndat_
2147 : logical :: ndat_is_odd
2148 :
2149 16912 : ndat_ = ndat / 2
2150 16912 : ndat_is_odd=.False.
2151 16912 : if (modulo(ndat,2)/=0) then
2152 0 : ndat_=ndat_+1
2153 0 : ndat_is_odd=.True.
2154 : end if
2155 :
2156 16912 : npw_fft=2*npw_k
2157 16912 : i0=1
2158 16912 : if (me_g0==1) then! Do not include G=(0,0,0) twice
2159 16912 : npw_fft=npw_fft-1
2160 16912 : i0=2
2161 : end if
2162 :
2163 16912 : if (size(cwavef,1)/=2) then
2164 0 : ABI_BUG('wrong size for cwavef (dim 1)')
2165 : end if
2166 16912 : if (ndat_is_odd) then
2167 0 : if (size(cwavef,2)/=npw_k*(2*ndat_-1)) then
2168 0 : ABI_BUG('wrong size for cwavef (dim 2) (odd)')
2169 : end if
2170 : else
2171 16912 : if (size(cwavef,2)/=npw_k*2*ndat_) then
2172 0 : ABI_BUG('wrong size for cwavef (dim 2)')
2173 : end if
2174 : end if
2175 16912 : if (size(cwavef_fft,1)/=2.or.size(cwavef_fft,2)/=npw_fft*ndat_) then
2176 0 : ABI_BUG('wrong size for cwavef_fft')
2177 : end if
2178 :
2179 106164 : do idat=1,ndat_
2180 89252 : ib1=(idat-1)*npw_fft ! band shift for cwavef_fft
2181 89252 : ib2=(idat-1)*2*npw_k ! band shift for cwavef
2182 : ! C(G) = ( E(G) + E(-G)^* ) / 2 (factor 1/2 will be applied later)
2183 : cwavef(1,i0+ib2:npw_k+ib2) = cwavef_fft(1,i0 +ib1:npw_k +ib1) & !+Re(E( G))
2184 68871680 : +cwavef_fft(1,1 +npw_k+ib1:npw_fft+ib1) !+Re(E(-G))
2185 : cwavef(2,i0+ib2:npw_k+ib2) = cwavef_fft(2,i0 +ib1:npw_k +ib1) & !+Im(E( G))
2186 68871680 : -cwavef_fft(2,1 +npw_k+ib1:npw_fft+ib1) !-Im(E(-G))
2187 89252 : if (.not.ndat_is_odd.or.idat<ndat_) then
2188 : ! D(G) = ( iE(-G)^* - iE(G) ) / 2 (factor 1/2 will be applied later)
2189 : cwavef(1,i0+npw_k+ib2:2*npw_k+ib2) = cwavef_fft(2,i0 +ib1:npw_k +ib1) & !+Im(E( G))
2190 68871680 : +cwavef_fft(2,1 +npw_k+ib1:npw_fft+ib1) !+Im(E(-G))
2191 : cwavef(2,i0+npw_k+ib2:2*npw_k+ib2) =-cwavef_fft(1,i0 +ib1:npw_k +ib1) & !-Re(E( G))
2192 68871680 : +cwavef_fft(1,1 +npw_k+ib1:npw_fft+ib1) !+Re(E(-G))
2193 : end if
2194 106164 : if (me_g0==1) then
2195 : ! Compute C(G=0) and D(G=0) and multiply by 2 as we apply 1/2 to the whole array shortly afterwards
2196 : ! C(G=0) = Re(E(G=0))
2197 89252 : cwavef(1,1+ib2) = two*cwavef_fft(1,1+ib1)
2198 89252 : cwavef(2,1+ib2) = zero
2199 89252 : if (.not.ndat_is_odd.or.idat<ndat_) then
2200 : ! D(G=0) = Im(E(G=0))
2201 89252 : cwavef(1,1+npw_k+ib2) = two*cwavef_fft(2,1+ib1)
2202 89252 : cwavef(2,1+npw_k+ib2) = zero
2203 : end if
2204 : end if
2205 :
2206 : end do
2207 :
2208 413246992 : cwavef(:,:) = half*cwavef(:,:)
2209 :
2210 16912 : end subroutine cwavef_double_rfft_trick_unpack
2211 : !!***
2212 :
2213 : !!****f* ABINIT/getgsc
2214 : !! NAME
2215 : !! getgsc
2216 : !!
2217 : !! FUNCTION
2218 : !! Compute <G|S|C> for all input vectors |Cnk> at a given k-point,
2219 : !! OR for one input vector |Cnk>.
2220 : !! |Cnk> are expressed in reciprocal space.
2221 : !! S is the overlap operator between |Cnk> (used for PAW).
2222 : !!
2223 : !! INPUTS
2224 : !! cg(2,mcg)=planewave coefficients of wavefunctions
2225 : !! cprj(natom,mcprj)= wave functions projected with non-local projectors: cprj=<p_i|Cnk>
2226 : !! gs_ham <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
2227 : !! ibg=shift to be applied on the location of data in the array cprj (beginning of current k-point)
2228 : !! icg=shift to be applied on the location of data in the array cg (beginning of current k-point)
2229 : !! igsc=shift to be applied on the location of data in the array gsc (beginning of current k-point)
2230 : !! ikpt,isppol=indexes of current (spin.kpoint)
2231 : !! mcg=second dimension of the cg array
2232 : !! mcprj=second dimension of the cprj array
2233 : !! mgsc=second dimension of the gsc array
2234 : !! mpi_enreg=information about MPI parallelization
2235 : !! ndat=number of bands to compute in parallel
2236 : !! natom=number of atoms in unit cell.
2237 : !! nband= if positive: number of bands at this k point for that spin polarization
2238 : !! if negative: abs(nband) is the index of the only band to be computed
2239 : !! npw_k=number of planewaves in basis for given k point.
2240 : !! nspinor=number of spinorial components of the wavefunctions
2241 : !! [select_k]=optional, option governing the choice of k points to be used.
2242 : !! gs_ham datastructure contains quantities needed to apply overlap operator
2243 : !! in reciprocal space between 2 kpoints, k and k^prime (equal in most cases);
2244 : !! if select_k=1, <k^prime|S|k> is applied [default]
2245 : !! if select_k=2, <k|S|k^prime> is applied
2246 : !! if select_k=3, <k|S|k> is applied
2247 : !! if select_k=4, <k^prime|S|k^prime> is applied
2248 : !!
2249 : !! OUTPUT
2250 : !! gsc(2,mgsc)= <g|S|Cnk> or <g|S^(1)|Cnk> (S=overlap)
2251 : !!
2252 : !! SOURCE
2253 :
2254 179581 : subroutine getgsc(cg,cprj,gs_ham,gsc,ibg,icg,igsc,ikpt,isppol,&
2255 : & mcg,mcprj,mgsc,mpi_enreg,ndat,natom,nband,npw_k,nspinor,select_k)
2256 :
2257 : !Arguments ------------------------------------
2258 : !scalars
2259 : integer,intent(in) :: ibg,icg,igsc,ikpt,isppol,mcg,mcprj
2260 : integer,intent(in) :: mgsc,natom,nband,npw_k,nspinor,ndat
2261 : !TODO : may be needed to distribute cprj over band procs
2262 : ! integer,intent(in) :: mband_mem
2263 : integer,intent(in),optional :: select_k
2264 : type(MPI_type),intent(in) :: mpi_enreg
2265 : type(gs_hamiltonian_type),intent(inout),target :: gs_ham
2266 : !arrays
2267 : real(dp),intent(in), target :: cg(2,mcg)
2268 : real(dp),intent(out), target :: gsc(2,mgsc)
2269 : type(pawcprj_type),intent(in) :: cprj(natom,mcprj)
2270 :
2271 : !Local variables-------------------------------
2272 : !scalars
2273 : integer :: choice,cpopt,dimenl1,dimenl2,iband,iband1,iband2,index_cg,index_cprj
2274 : integer :: index_gsc,me,my_nspinor,my_ndat,paw_opt,select_k_,signs,tim_nonlop,useylm
2275 : !character(len=500) :: msg
2276 : !arrays
2277 359162 : real(dp) :: enlout_dum(ndat),tsec(2)
2278 179581 : real(dp), contiguous, pointer :: cwavef(:,:),scwavef(:,:)
2279 179581 : type(pawcprj_type),allocatable :: cwaveprj(:,:)
2280 : ! *********************************************************************
2281 :
2282 : DBG_ENTER("COLL")
2283 :
2284 : !Compatibility tests
2285 179581 : my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
2286 66816 : my_ndat=ndat; if(nband<1) my_ndat=1
2287 179581 : if(gs_ham%usepaw==0) then
2288 0 : ABI_BUG('Only compatible with PAW (usepaw=1) !')
2289 : end if
2290 179581 : if(nband<0.and.(mcg<npw_k*my_nspinor.or.mgsc<npw_k*my_nspinor.or.mcprj<my_nspinor)) then
2291 0 : ABI_BUG('Invalid value for mcg, mgsc or mcprj !')
2292 : end if
2293 :
2294 : !Keep track of total time spent in getgsc:
2295 179581 : call timab(565,1,tsec)
2296 :
2297 179581 : if(gs_ham%gpu_option==ABI_GPU_DISABLED) then
2298 273119845 : gsc = zero
2299 : else if(gs_ham%gpu_option==ABI_GPU_DISABLED) then
2300 : #ifdef HAVE_OPENMP_OFFLOAD
2301 : call gpu_set_to_zero(gsc,int(2,c_size_t)*mgsc)
2302 : #endif
2303 : end if
2304 :
2305 : !Prepare some data
2306 179581 : if (gs_ham%usecprj==1) then
2307 804010 : ABI_MALLOC(cwaveprj,(natom,my_nspinor*my_ndat))
2308 112765 : call pawcprj_alloc(cwaveprj,0,gs_ham%dimcprj)
2309 : else
2310 66816 : ABI_MALLOC(cwaveprj,(0,0))
2311 : end if
2312 179581 : dimenl1=gs_ham%dimekb1;dimenl2=natom;tim_nonlop=0
2313 179581 : choice=1;signs=2;cpopt=-1+3*gs_ham%usecprj;paw_opt=3;useylm=1
2314 179581 : select_k_=1;if (present(select_k)) select_k_=select_k
2315 179581 : me=mpi_enreg%me_kpt
2316 :
2317 : !Loop over bands
2318 179581 : index_cprj=ibg;index_cg=icg;index_gsc=igsc
2319 179581 : if (nband>0) then
2320 : iband1=1;iband2=nband
2321 : !only do 1 band in case nband < 0 (the |nband|th one)
2322 66816 : else if (nband<0) then
2323 66816 : iband1=abs(nband);iband2=iband1
2324 66816 : index_cprj=index_cprj+(iband1-1)*my_nspinor
2325 66816 : index_cg =index_cg +(iband1-1)*npw_k*my_nspinor
2326 66816 : index_gsc =index_gsc +(iband1-1)*npw_k*my_nspinor
2327 : end if
2328 :
2329 848129 : do iband=iband1,iband2,my_ndat
2330 :
2331 668548 : if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me.and.nband>0) then
2332 : ! No longer needed 28/03/2020 to parallelize memory
2333 : ! gsc(:,1+index_gsc:npw_k*my_nspinor+index_gsc)=zero
2334 : ! index_gsc=index_gsc+npw_k*my_nspinor
2335 : !index_cprj=index_cprj+my_nspinor
2336 : !index_cg=index_cg+npw_k*my_nspinor
2337 :
2338 : cycle
2339 : end if
2340 :
2341 : ! Retrieve WF at (n,k)
2342 665380 : cwavef(1:2,1:npw_k*my_nspinor*my_ndat) => cg(:,1+index_cg:npw_k*my_nspinor*my_ndat+index_cg)
2343 665380 : scwavef(1:2,1:npw_k*my_nspinor*my_ndat) => gsc(:,1+index_gsc:npw_k*my_nspinor*my_ndat+index_gsc)
2344 665380 : if (gs_ham%usecprj==1) then
2345 598564 : call pawcprj_copy(cprj(:,1+index_cprj:my_nspinor*my_ndat+index_cprj),cwaveprj)
2346 : end if
2347 :
2348 : ! Compute <g|S|Cnk>
2349 : call nonlop(choice,cpopt,cwaveprj,enlout_dum,gs_ham,0,(/zero/),mpi_enreg,my_ndat,1,paw_opt,&
2350 665380 : & signs,scwavef,tim_nonlop,cwavef,cwavef,select_k=select_k_)
2351 :
2352 :
2353 : ! End of loop over bands
2354 665380 : index_cprj=index_cprj+my_nspinor*my_ndat
2355 665380 : index_cg=index_cg+npw_k*my_nspinor*my_ndat
2356 848129 : index_gsc=index_gsc+npw_k*my_nspinor*my_ndat
2357 : end do
2358 :
2359 : !Memory deallocation
2360 179581 : if (gs_ham%usecprj==1) then
2361 112765 : call pawcprj_free(cwaveprj)
2362 : end if
2363 416652 : ABI_FREE(cwaveprj)
2364 :
2365 179581 : call timab(565,2,tsec)
2366 :
2367 : DBG_EXIT("COLL")
2368 :
2369 359162 : end subroutine getgsc
2370 : !!***
2371 :
2372 : !!****f* ABINIT/multithreaded_getghc
2373 : !!
2374 : !! NAME
2375 : !! multithreaded_getghc
2376 : !!
2377 : !! FUNCTION
2378 : !!
2379 : !! INPUTS
2380 : !! cpopt=flag defining the status of cwaveprj%cp(:)=<Proj_i|Cnk> scalars (PAW only)
2381 : !! (same meaning as in nonlop.F90 routine)
2382 : !! if cpopt=-1, <p_lmn|in> (and derivatives) are computed here (and not saved)
2383 : !! if cpopt= 0, <p_lmn|in> are computed here and saved
2384 : !! if cpopt= 1, <p_lmn|in> and first derivatives are computed here and saved
2385 : !! if cpopt= 2 <p_lmn|in> are already in memory;
2386 : !! if cpopt= 3 <p_lmn|in> are already in memory; first derivatives are computed here and saved
2387 : !! if cpopt= 4 <p_lmn|in> and first derivatives are already in memory;
2388 : !! cwavef(2,npw*my_nspinor*ndat)=planewave coefficients of wavefunction.
2389 : !! gs_ham <type(gs_hamiltonian_type)>=all data for the Hamiltonian to be applied
2390 : !! lambda=factor to be used when computing <G|H-lambda.S|C> - only for sij_opt=-1
2391 : !! Typically lambda is the eigenvalue (or its guess)
2392 : !! mpi_enreg=information about MPI parallelization
2393 : !! ndat=number of FFT to do in parallel
2394 : !! prtvol=control print volume and debugging output
2395 : !! sij_opt= -PAW ONLY- if 0, only matrix elements <G|H|C> have to be computed
2396 : !! (S=overlap) if 1, matrix elements <G|S|C> have to be computed in gsc in addition to ghc
2397 : !! if -1, matrix elements <G|H-lambda.S|C> have to be computed in ghc (gsc not used)
2398 : !! tim_getghc=timing code of the calling subroutine(can be set to 0 if not attributed)
2399 : !! type_calc= option governing which part of Hamitonian is to be applied:
2400 : !! 0: whole Hamiltonian
2401 : !! 1: local part only
2402 : !! 2: non-local+kinetic only (added to the existing Hamiltonian)
2403 : !! 3: local + kinetic only (added to the existing Hamiltonian)
2404 : !! ===== Optional inputs =====
2405 : !! [kg_fft_k(3,:)]=optional, (k+G) vector coordinates to be used for the FFT transformation
2406 : !! instead of the one contained in gs_ham datastructure.
2407 : !! Typically used for real WF (in parallel) which are FFT-transformed 2 by 2.
2408 : !! [kg_fft_kp(3,:)]=optional, (k^prime+G) vector coordinates to be used for the FFT transformation
2409 : !! [select_k]=optional, option governing the choice of k points to be used.
2410 : !! gs_ham datastructure contains quantities needed to apply Hamiltonian
2411 : !! in reciprocal space between 2 kpoints, k and k^prime (equal in most cases);
2412 : !! if select_k=1, <k^prime|H|k> is applied [default]
2413 : !! if select_k=2, <k|H|k^prime> is applied
2414 : !! if select_k=3, <k|H|k> is applied
2415 : !! if select_k=4, <k^prime|H|k^prime> is applied
2416 : !!
2417 : !! OUTPUT
2418 : !! ghc(2,npw*my_nspinor*ndat)=matrix elements <G|H|C> (if sij_opt>=0)
2419 : !! or <G|H-lambda.S|C> (if sij_opt=-1)
2420 : !! gvnlxc(2,npw*my_nspinor*ndat)=matrix elements <G|Vnonlocal|C> (if sij_opt>=0)
2421 : !! or <G|Vnonlocal-lambda.S|C> (if sij_opt=-1)
2422 : !! if (sij_opt=1)
2423 : !! gsc(2,npw*my_nspinor*ndat)=matrix elements <G|S|C> (S=overlap).
2424 : !!
2425 : !! SIDE EFFECTS
2426 : !! cwaveprj(natom,my_nspinor*(1+cpopt)*ndat)= wave function projected on nl projectors (PAW only)
2427 : !!
2428 : !! SOURCE
2429 :
2430 2238536 : subroutine multithreaded_getghc(cpopt,cwavef,cwaveprj,ghc,gsc,gs_ham,gvnlxc,lambda,mpi_enreg,ndat,&
2431 : & prtvol,sij_opt,tim_getghc,type_calc,&
2432 1119268 : & kg_fft_k,kg_fft_kp,select_k,filter_dilatmx_loc) ! optional arguments
2433 :
2434 : #ifdef HAVE_OPENMP
2435 : use omp_lib
2436 : #endif
2437 :
2438 : !Arguments ------------------------------------
2439 : !scalars
2440 : logical,intent(in),optional :: filter_dilatmx_loc
2441 : integer,intent(in) :: cpopt,ndat, prtvol
2442 : integer,intent(in) :: sij_opt,tim_getghc,type_calc
2443 : integer,intent(in),optional :: select_k
2444 : real(dp),intent(in) :: lambda
2445 : type(MPI_type),intent(in) :: mpi_enreg
2446 : type(gs_hamiltonian_type),intent(inout),target :: gs_ham
2447 : !arrays
2448 : integer,intent(in),optional,target :: kg_fft_k(:,:),kg_fft_kp(:,:)
2449 : real(dp),intent(out),target :: gsc(:,:)
2450 : real(dp),intent(inout) :: cwavef(:,:)
2451 : real(dp),intent(out) :: ghc(:,:),gvnlxc(:,:)
2452 : type(pawcprj_type),intent(inout),target :: cwaveprj(:,:)
2453 :
2454 : !Local variables-------------------------------
2455 : !scalars
2456 : logical :: filter_dilatmx_loc_
2457 : integer :: firstelt, firstprj, lastelt, lastprj,usegvnlxc,usegsc
2458 : integer :: nthreads,fftalga
2459 : integer :: ithread
2460 : integer :: chunk
2461 : integer :: residuchunk
2462 : integer :: firstband
2463 : integer :: lastband
2464 : integer :: spacedim, spacedim_prj
2465 : logical :: fftw3_use_lib_threads_sav
2466 : integer :: select_k_default
2467 : ! *************************************************************************
2468 :
2469 1119268 : select_k_default = 1; if (present(select_k)) select_k_default = select_k
2470 1119268 : filter_dilatmx_loc_ = .true.; if ( present(filter_dilatmx_loc) ) filter_dilatmx_loc_ = filter_dilatmx_loc
2471 :
2472 1119268 : spacedim = size(cwavef ,dim=2)/ndat
2473 1119268 : spacedim_prj = size(cwaveprj,dim=2)/ndat
2474 :
2475 1119268 : nthreads = xomp_get_num_threads(open_parallel=.True.)
2476 1119268 : fftalga = gs_ham%ngfft(7)/100
2477 1119268 : if (fftalga==FFT_SG.and.nthreads>1.and.ndat>1) then
2478 0 : ABI_ERROR("fftalg=1XX is not thread-safe, so it cannot be used in multi-threaded hamiltonian with nthreads>1 and ndat>1.")
2479 : end if
2480 :
2481 : ! Disabling multithreading for GPU variants (getghc_ompgpu is not thread-safe for now)
2482 : !$omp parallel default (none) &
2483 : !$omp& private(ithread,nthreads,chunk,firstband,lastband,residuchunk,firstelt,lastelt), &
2484 : !$omp& private(firstprj,lastprj,usegvnlxc,usegsc,fftw3_use_lib_threads_sav), &
2485 : !$omp& shared(cwavef,ghc,gsc, gvnlxc,spacedim,spacedim_prj,ndat,kg_fft_k,kg_fft_kp,gs_ham,cwaveprj,mpi_enreg), &
2486 : !$omp& shared(gemm_nonlop_use_gemm), &
2487 : !$omp& firstprivate(cpopt,lambda,prtvol,sij_opt,tim_getghc,type_calc,select_k_default,filter_dilatmx_loc_) &
2488 : !$omp& IF(gs_ham%gpu_option==ABI_GPU_DISABLED .and. .not. gemm_nonlop_use_gemm)
2489 1119268 : ithread = 0
2490 1119268 : nthreads = 1
2491 1119268 : fftw3_use_lib_threads_sav = .false.
2492 : if(gs_ham%gpu_option==ABI_GPU_DISABLED .and. .not. gemm_nonlop_use_gemm) then
2493 : #ifdef HAVE_OPENMP
2494 : ithread = omp_get_thread_num()
2495 : nthreads = omp_get_num_threads()
2496 : !Ensure that libs are used without threads (mkl, openblas, fftw3, ...)
2497 : #ifdef HAVE_LINALG_MKL_THREADS
2498 : call mkl_set_num_threads(1)
2499 : #endif
2500 : !LB-23/07/24: OpenBLAS detects parallel sections automatically. To comment this line improves performances for some cases.
2501 : !#ifdef HAVE_LINALG_OPENBLAS_THREADS
2502 : ! call openblas_set_num_threads(1)
2503 : !#endif
2504 : #ifdef HAVE_LINALG_NVPL_THREADS
2505 : call nvpl_blas_set_num_threads(1)
2506 : #endif
2507 : #ifdef HAVE_FFTW3_THREADS
2508 : fftw3_use_lib_threads_sav=(.not.fftw3_spawn_threads_here(nthreads,nthreads))
2509 : call fftw3_use_lib_threads(.false.)
2510 : #endif
2511 : #endif
2512 : end if
2513 1119268 : chunk = ndat/nthreads ! Divide by 2 to construct chunk of even number of bands
2514 1119268 : residuchunk = ndat - nthreads*chunk
2515 : if ( ithread < nthreads-residuchunk ) then
2516 : firstband = ithread*chunk+1
2517 : lastband = (ithread+1)*chunk
2518 : else
2519 : firstband = (nthreads-residuchunk)*chunk + ( ithread -(nthreads-residuchunk) )*(chunk+1) +1
2520 : lastband = firstband+chunk
2521 : end if
2522 1119268 : usegvnlxc=1
2523 3357804 : if (size(gvnlxc)<=1) usegvnlxc=0
2524 1119268 : usegsc=0
2525 1119268 : if (gs_ham%usepaw==1.and.sij_opt==1) usegsc=1
2526 :
2527 1119268 : if ( lastband /= 0 ) then
2528 1119268 : firstelt = (firstband-1)*spacedim+1
2529 1119268 : firstprj = (firstband-1)*spacedim_prj+1
2530 1119268 : lastelt = lastband*spacedim
2531 1119268 : lastprj = lastband*spacedim_prj
2532 : ! Don't know how to manage optional arguments .... :(
2533 1119268 : if ( present(kg_fft_k) ) then
2534 3274 : if (present(kg_fft_kp)) then
2535 : call getghc(cpopt,cwavef(:,firstelt:lastelt),cwaveprj(:,firstprj:lastprj),&
2536 : & ghc(:,firstelt:lastelt),gsc(:,firstelt:lastelt*usegsc),&
2537 : & gs_ham,gvnlxc(:,firstelt:lastelt*usegvnlxc),lambda, mpi_enreg,lastband-firstband+1,&
2538 : & prtvol,sij_opt,tim_getghc,type_calc,&
2539 0 : & select_k=select_k_default,kg_fft_k=kg_fft_k,kg_fft_kp=kg_fft_kp,filter_dilatmx_loc=filter_dilatmx_loc_)
2540 : else
2541 : call getghc(cpopt,cwavef(:,firstelt:lastelt),cwaveprj(:,firstprj:lastprj),&
2542 : & ghc(:,firstelt:lastelt),gsc(:,firstelt:lastelt*usegsc),&
2543 : & gs_ham,gvnlxc(:,firstelt:lastelt*usegvnlxc),lambda, mpi_enreg,lastband-firstband+1,&
2544 : & prtvol,sij_opt,tim_getghc,type_calc,&
2545 3274 : & select_k=select_k_default,kg_fft_k=kg_fft_k,filter_dilatmx_loc=filter_dilatmx_loc_)
2546 : end if
2547 : else
2548 1115994 : if (present(kg_fft_kp)) then
2549 : call getghc(cpopt,cwavef(:,firstelt:lastelt),cwaveprj(:,firstprj:lastprj),&
2550 : & ghc(:,firstelt:lastelt),gsc(:,firstelt:lastelt*usegsc),&
2551 : & gs_ham,gvnlxc(:,firstelt:lastelt*usegvnlxc),lambda, mpi_enreg,lastband-firstband+1,&
2552 : & prtvol,sij_opt,tim_getghc,type_calc,&
2553 0 : & select_k=select_k_default,kg_fft_kp=kg_fft_kp,filter_dilatmx_loc=filter_dilatmx_loc_)
2554 : else
2555 : call getghc(cpopt,cwavef(:,firstelt:lastelt),cwaveprj(:,firstprj:lastprj),&
2556 : & ghc(:,firstelt:lastelt),gsc(:,firstelt:lastelt*usegsc),&
2557 : & gs_ham,gvnlxc(:,firstelt:lastelt*usegvnlxc),lambda, mpi_enreg,lastband-firstband+1,&
2558 : & prtvol,sij_opt,tim_getghc,type_calc,&
2559 1115994 : & select_K=select_k_default,filter_dilatmx_loc=filter_dilatmx_loc_)
2560 : end if
2561 : end if
2562 : end if
2563 : if(gs_ham%gpu_option==ABI_GPU_DISABLED .and. .not. gemm_nonlop_use_gemm) then
2564 : #ifdef HAVE_OPENMP
2565 : !Restore libs behavior (mkl, openblas, fftw3, ...)
2566 : #ifdef HAVE_LINALG_MKL_THREADS
2567 : call mkl_set_num_threads(nthreads)
2568 : #endif
2569 : !LB-23/07/24: OpenBLAS detects parallel sections automatically. To comment this line improves performances for some cases.
2570 : !#ifdef HAVE_LINALG_OPENBLAS_THREADS
2571 : ! call openblas_set_num_threads(nthreads)
2572 : !#endif
2573 : #ifdef HAVE_LINALG_NVPL_THREADS
2574 : call nvpl_blas_set_num_threads(nthreads)
2575 : #endif
2576 : #ifdef HAVE_FFTW3_THREADS
2577 : call fftw3_use_lib_threads(fftw3_use_lib_threads_sav)
2578 : #endif
2579 : #endif
2580 : end if
2581 : !$omp end parallel
2582 :
2583 1119268 : end subroutine multithreaded_getghc
2584 : !!***
2585 :
2586 : end module m_getghc
2587 : !!***
|
