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