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 [ - + - + : 35264769 : subroutine getghc(cpopt,cwavef,cwaveprj,ghc,gsc,gs_ham,gvnlxc,lambda,mpi_enreg,ndat,&
- + - + -
+ ]
184 : : prtvol,sij_opt,tim_getghc,type_calc,&
185 [ - + - - : 35264769 : 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 : 35264769 : integer, contiguous, pointer :: gbound_k1(:,:), gbound_k2(:,:)
225 : 35264769 : integer, contiguous, pointer :: kg_k1(:,:), kg_k2(:,:)
226 : 35264769 : integer, contiguous, pointer :: indices_pw_fft(:), kg_k_fft(:,:)
227 : 35264769 : integer, contiguous, pointer :: recvcount_fft(:), recvdisp_fft(:)
228 : 35264769 : integer, contiguous, pointer :: sendcount_fft(:), senddisp_fft(:)
229 : 35264769 : integer, allocatable:: dimcprj(:)
230 : 105794307 : real(dp) :: enlout(ndat), lambda_ndat(ndat), tsec(2)
231 : : real(dp), target :: nonlop_dum(1,1)
232 : 35264769 : real(dp), allocatable :: buff_wf(:,:)
233 : 35264769 : real(dp), allocatable :: cwavef1(:,:)
234 : 35264769 : real(dp), allocatable :: cwavef2(:,:)
235 : 35264769 : real(dp), allocatable :: cwavef_fft(:,:)
236 : 35264769 : real(dp), allocatable :: cwavef_fft_tr(:,:)
237 : 35264769 : real(dp), allocatable :: ghc1(:,:)
238 : 35264769 : real(dp), allocatable :: ghc2(:,:)
239 : 35264769 : real(dp), allocatable :: ghc3(:,:)
240 : 35264769 : real(dp), allocatable :: ghc4(:,:)
241 : 35264769 : real(dp), allocatable :: ghc_mGGA(:,:)
242 : 35264769 : real(dp), allocatable :: ghc_mGGA_fft(:,:)
243 : 35264769 : real(dp), allocatable :: ghc_vectornd(:,:)
244 : 35264769 : 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 : 35264769 : real(dp), target, allocatable :: gvnlc(:,:)
251 : 35264769 : real(dp), contiguous, pointer :: gvnlxc_(:,:)
252 : : #endif
253 : 35264769 : real(dp), allocatable :: vlocal_tmp(:,:,:), work(:,:,:,:)
254 : 35264769 : real(dp), contiguous, pointer :: kinpw_k1(:), kinpw_k2(:), kpt_k1(:), kpt_k2(:)
255 : 35264769 : real(dp), contiguous, pointer :: gsc_ptr(:,:)
256 : : type(fock_common_type),pointer :: fock
257 : 35264769 : 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 : 35264769 : call timab(350+tim_getghc,1,tsec)
266 : : ABI_NVTX_START_RANGE(NVTX_GETGHC)
267 : :
268 : : !Structured debugging if prtvol==-level
269 [ - + ]: 35264769 : 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 [ + + ]: 35264769 : select_k_=KPRIME_H_K;if (present(select_k)) select_k_=select_k
276 : 14803455 : select case (select_k_)
277 : : case (KPRIME_H_K)
278 : : ! <k^prime|H|k>
279 : 14803455 : npw_k1 = gs_ham%npw_fft_k ; npw_k2 = gs_ham%npw_fft_kp
280 : 14803455 : kpt_k1 => gs_ham%kpt_k ; kpt_k2 => gs_ham%kpt_kp
281 : 14803455 : kg_k1 => gs_ham%kg_k ; kg_k2 => gs_ham%kg_kp
282 : 14803455 : gbound_k1 => gs_ham%gbound_k ; gbound_k2 => gs_ham%gbound_kp
283 : 14803455 : 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 [ + - - + : 35264769 : ABI_ERROR(sjoin("Invalid select_k: ", itoa(select_k_)))
- ]
307 : : end select
308 : :
309 [ + + + + ]: 140533972 : k1_eq_k2=(all(abs(kpt_k1(:)-kpt_k2(:))<tol8))
310 : :
311 : : ! Check sizes
312 : 35264769 : my_nspinor=max(1,gs_ham%nspinor/mpi_enreg%nproc_spinor)
313 [ + + - + : 105794307 : ABI_CHECK_IGEQ(size(cwavef), 2*npw_k1*my_nspinor*ndat, 'wrong size for cwavef!')
- - ]
314 [ + + - + : 105794307 : ABI_CHECK_IGEQ(size(ghc), 2*npw_k2*my_nspinor*ndat, 'wrong size for ghc!')
- - ]
315 [ + + ]: 35264769 : if (sij_opt==1) then
316 [ + + - + : 25028241 : ABI_CHECK_IGEQ(size(gsc), 2*npw_k2*my_nspinor*ndat, 'wrong size for gsc!')
- - ]
317 : : end if
318 [ + + + + ]: 35264769 : 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 [ + + - + ]: 35264769 : if (any(type_calc == [0, 2, 3])) then
322 [ + + ]: 105170535 : local_gvnlxc = size(gvnlxc)<=1
323 [ + + ]: 35056845 : if (local_gvnlxc) then
324 [ + - ]: 912524 : 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 [ - + + - : 2737572 : ABI_MALLOC(gvnlxc_,(2,npw_k2*my_nspinor*ndat))
+ - + - -
+ - + ]
330 : : end if
331 : : else
332 : 34144321 : gvnlxc_ => gvnlxc
333 : : end if
334 [ + + + - : 105170535 : ABI_CHECK_IGEQ(size(gvnlxc_), 2*npw_k2*my_nspinor*ndat, 'wrong size for gvnlxc!')
- - ]
335 : : end if
336 : :
337 : 35264769 : use_cwavef_r=present(cwavef_r)
338 : 35264769 : n4=gs_ham%n4 ; n5=gs_ham%n5 ; n6=gs_ham%n6
339 : 35264769 : nspinortot=gs_ham%nspinor
340 : :
341 [ - + ]: 35264769 : 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 [ + + ]: 35264769 : 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 [ - + ]: 35264769 : 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 [ + + - + : 35264769 : 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 : 35264769 : has_fock = associated(gs_ham%fockcommon)
369 [ + + ]: 35264769 : if (has_fock) fock => gs_ham%fockcommon
370 : :
371 : : !Parallelization over spinors management
372 [ + + ]: 35264769 : if (mpi_enreg%paral_spinor==0) then
373 : 35201405 : shift1=npw_k1;shift2=npw_k2
374 : 35201405 : nspinor1TreatedByThisProc=.true.
375 : 35201405 : nspinor2TreatedByThisProc=(nspinortot==2)
376 : : else
377 : 63364 : shift1=0;shift2=0
378 : 63364 : nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
379 : 63364 : nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
380 : : end if
381 : :
382 [ + + ]: 35264769 : filter_dilatmx_loc_ = .true.; if ( present(filter_dilatmx_loc) ) filter_dilatmx_loc_ = filter_dilatmx_loc
383 : 35264769 : 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 [ + + ]: 35264769 : if (any(type_calc == [0, 1, 3])) then
409 : :
410 : : ! Need a Vlocal
411 [ - + ]: 35261495 : 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 [ + + ]: 35261495 : if (gs_ham%use_gbt == 0) then
415 [ - + ]: 35056849 : ABI_CHECK(k1_eq_k2, 'vlocal (fourwf) cannot be computed with k/=k^prime!')
416 : : end if
417 : :
418 : 35261495 : nfft_blocks = 1;
419 [ - + - - ]: 35261495 : 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 : 35261495 : chunk = ndat/nfft_blocks ! Divide by 2 to construct chunk of even number of bands
423 : 35261495 : 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 : 35261495 : have_to_reequilibrate=.false.
428 [ + + ]: 35261495 : if (mpi_enreg%paral_kgb==1) then
429 : 932206 : ikpt_this_proc=bandfft_kpt_get_ikpt()
430 : 932206 : have_to_reequilibrate=bandfft_kpt(ikpt_this_proc)%have_to_reequilibrate
431 : : end if
432 : 35261495 : ndat_ = ndat
433 : 35261495 : istwf_k_ = gs_ham%istwf_k
434 [ + + + + : 35261495 : 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 [ - + ]: 35261495 : 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 : 35261495 : fftalg = gs_ham%ngfft(7)
497 [ + - + + ]: 35261495 : if (gs_ham%gpu_option==ABI_GPU_DISABLED.and.fftalg/=401) then
498 [ - + + - : 173188225 : ABI_MALLOC(work,(2,gs_ham%n4,gs_ham%n5,gs_ham%n6))
- + + - -
+ + - + -
+ - - + -
+ ]
499 [ + - ]: 623850 : else if(nfft_blocks==1) then
500 [ - + + - : 3119250 : 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 : 35261495 : weight=one
508 [ + - ]: 35261495 : if (.not.use_cwavef_r) then
509 : 35261495 : option_fft=2
510 [ + + ]: 35261495 : if (nspinortot==2) then
511 : : ! Note: for this case we have ndat_=ndat
512 [ - + + - : 7933521 : ABI_MALLOC(cwavef1,(2,npw_k1*ndat))
+ - + - -
+ - + ]
513 [ - + + - : 5289014 : ABI_MALLOC(cwavef2,(2,npw_k1*ndat))
- + - + ]
514 [ + - ]: 2644507 : 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 [ + + ]: 5604933 : do idat=1,ndat
530 [ + + ]: 1019738177 : do ipw=1,npw_k1
531 [ + + ]: 3042399732 : cwavef1(1:2,ipw+(idat-1)*npw_k1)=cwavef(1:2,ipw+(idat-1)*my_nspinor*npw_k1)
532 [ + + ]: 3045360158 : 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 [ + + ]: 35261495 : if (gs_ham%nvloc==1) then
547 : : ! Treat scalar local potentials
548 : :
549 [ + + ]: 33134727 : if (nspinortot==1) then
550 : :
551 [ - + ]: 32616988 : 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 [ + + ]: 32616988 : 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 [ + - - - : 32600236 : 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 [ + + ]: 517739 : if (nspinor1TreatedByThisProc) then
578 [ - + ]: 513977 : 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 [ - + + - : 1541931 : 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 [ + - - - : 513977 : weight,weight,gpu_option=gs_ham%gpu_option)
- - - - -
- - - - +
- - - - -
- - - ]
599 [ + - ]: 513977 : 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 [ + + ]: 1342187 : do idat=1,ndat
613 [ + + ]: 256554290 : do ipw =1, npw_k2
614 [ + + ]: 766464519 : 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 [ - + ]: 513977 : ABI_FREE(ghc1)
622 : : end if ! spin 1 treated by this proc
623 : :
624 [ + + ]: 517739 : if (nspinor2TreatedByThisProc) then
625 [ - + ]: 513977 : 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 [ - + + - : 1541931 : 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 [ + - - - : 513977 : gpu_option=gs_ham%gpu_option)
- - - - -
- - - - +
- - - - -
- - - ]
643 : :
644 [ + - ]: 513977 : 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 [ + + ]: 1342187 : do idat=1,ndat
658 [ + + ]: 256554290 : do ipw=1,npw_k2
659 [ + + ]: 766464519 : 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 [ - + ]: 513977 : ABI_FREE(ghc2)
667 : : end if ! spin 2 treated by this proc
668 : :
669 : : end if ! nspinortot
670 : :
671 [ + - ]: 2126768 : else if (gs_ham%nvloc==4) then
672 : : ! Treat non-collinear local potentials
673 : : ! Note: for this case we have ndat_=ndat
674 [ - + + - : 6380304 : ABI_MALLOC(ghc1,(2,npw_k2*ndat))
+ - + - -
+ - + ]
675 [ - + + - : 4253536 : ABI_MALLOC(ghc2,(2,npw_k2*ndat))
- + - + ]
676 [ - + + - : 4253536 : ABI_MALLOC(ghc3,(2,npw_k2*ndat))
- + - + -
+ ]
677 [ - + + - : 4253536 : ABI_MALLOC(ghc4,(2,npw_k2*ndat))
- + - + -
+ ]
678 [ + - ]: 2126768 : 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 [ + + + + : 9050518508 : ghc1(:,:)=zero; ghc2(:,:)=zero; ghc3(:,:)=zero ; ghc4(:,:)=zero
+ + + + +
+ + + + +
+ + ]
688 : : end if
689 [ - + ]: 2126768 : if (use_cwavef_r) then
690 [ # # ]: 0 : ABI_MALLOC(vlocal_tmp,(0,0,0))
691 : : else
692 [ - + + - : 10633840 : ABI_MALLOC(vlocal_tmp,(gs_ham%n4,gs_ham%n5,gs_ham%n6))
- + + - -
+ + - + -
+ - - + -
+ ]
693 : : end if
694 : : ! ghc1=v11*phi1
695 [ + + ]: 2126768 : if (nspinor1TreatedByThisProc) then
696 [ - + ]: 2098848 : 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 [ + + ]: 39109011 : do i3=1,n6
714 [ + + ]: 790638374 : do i2=1,n5
715 [ + + ]: 18237739185 : do i1=1,n4
716 : 18200729022 : 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 : 2098848 : gpu_option=gs_ham%gpu_option)
725 : : end if
726 : : ! ghc2=v22*phi2
727 [ + + ]: 2126768 : if (nspinor2TreatedByThisProc) then
728 [ - + ]: 2098848 : 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 [ + + ]: 39109011 : do i3=1,n6
746 [ + + ]: 790638374 : do i2=1,n5
747 [ + + ]: 18237739185 : do i1=1,n4
748 : 18200729022 : 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 : 2098848 : gpu_option=gs_ham%gpu_option)
760 : : end if
761 [ - + ]: 2126768 : ABI_FREE(vlocal_tmp)
762 : 2126768 : cplex=2
763 [ - + ]: 2126768 : if (use_cwavef_r) then
764 [ # # ]: 0 : ABI_MALLOC(vlocal_tmp,(0,0,0))
765 : : else
766 [ - + + - : 10633840 : 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 [ + + ]: 2126768 : if (nspinor1TreatedByThisProc) then
770 [ - + ]: 2098848 : 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 [ + + ]: 39109011 : do i3=1,gs_ham%n6
784 [ + + ]: 790638374 : do i2=1,gs_ham%n5
785 [ + + ]: 18237739185 : do i1=1,gs_ham%n4
786 : 17449199659 : vlocal_tmp(2*i1-1,i2,i3)= gs_ham%vlocal(i1,i2,i3,3)
787 : 18200729022 : 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 : 2098848 : gpu_option=gs_ham%gpu_option)
796 : : end if
797 : : ! ghc4=(re(v12)+im(v12))*phi2
798 [ + + ]: 2126768 : if (nspinor2TreatedByThisProc) then
799 [ - + ]: 2098848 : 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 [ + + ]: 39109011 : do i3=1,gs_ham%n6
813 [ + + ]: 790638374 : do i2=1,gs_ham%n5
814 [ + + ]: 18237739185 : do i1=1,gs_ham%n4
815 : 17449199659 : vlocal_tmp(2*i1-1,i2,i3)= gs_ham%vlocal(i1,i2,i3,3)
816 : 18200729022 : 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 : 2098848 : gpu_option=gs_ham%gpu_option)
825 : : end if
826 [ - + ]: 2126768 : 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 [ + + ]: 2126768 : if (mpi_enreg%paral_spinor==0) then
830 [ + - ]: 2070928 : 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 [ + + ]: 4141856 : do idat=1,ndat
845 [ + + ]: 747414261 : do ipw=1,npw_k2
846 [ + + ]: 2229817215 : 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 [ + + ]: 2231888143 : 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 : 55840 : call xmpi_sum(ghc4,mpi_enreg%comm_spinor,ierr)
853 : 55840 : call xmpi_sum(ghc3,mpi_enreg%comm_spinor,ierr)
854 [ + - ]: 55840 : 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 [ + + ]: 55840 : if (nspinor1TreatedByThisProc) then
878 [ + + ]: 55840 : do idat=1,ndat
879 [ + + ]: 5435960 : do ipw=1,npw_k2
880 [ + + ]: 16168280 : 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 [ + - ]: 27920 : else if (nspinor2TreatedByThisProc) then
884 [ + + ]: 55840 : do idat=1,ndat
885 [ + + ]: 5435960 : do ipw=1,npw_k2
886 [ + + ]: 16168280 : 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 [ - + ]: 2126768 : ABI_FREE(ghc1)
896 [ - + ]: 2126768 : ABI_FREE(ghc2)
897 [ - + ]: 2126768 : ABI_FREE(ghc3)
898 [ - + ]: 2126768 : ABI_FREE(ghc4)
899 : : end if ! nvloc
900 : :
901 [ + + ]: 35261495 : 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 [ - + ]: 2644507 : ABI_FREE(cwavef1)
906 [ - + ]: 2644507 : 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 [ - + ]: 35261495 : ABI_FREE(work)
913 : :
914 [ + + ]: 35261495 : 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 [ - + ]: 35261495 : 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 [ + + ]: 35261495 : 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 : : #ifdef HAVE_OPENMP_OFFLOAD
954 : : !$OMP TARGET UPDATE FROM(cwavef) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
955 : : #endif
956 [ + + ]: 90464 : if (double_rfft_trick) then
957 [ - + + - : 480 : ABI_MALLOC(kg_k_fft,(3,npw_fft))
+ - + - -
+ - + ]
958 [ - + + - : 480 : ABI_MALLOC(cwavef_fft,(2,npw_fft*ndat_))
+ - + - -
+ - + -
+ ]
959 [ - + + - : 320 : ABI_MALLOC(ghc_mGGA_fft,(2,npw_fft*ndat_))
- + - + -
+ ]
960 [ + - + + : 413600 : kg_k_fft(:,1:npw_k1) = kg_k1(:,1:npw_k1)
+ + + + +
+ ]
961 [ + - + + : 412320 : kg_k_fft(:,npw_k1+1:npw_fft) = -kg_k1(:,i0:npw_k1)
+ + + + +
+ ]
962 : 160 : call cwavef_double_rfft_trick_pack(cwavef,cwavef_fft,mpi_enreg%me_g0_fft,ndat,npw_k1)
963 : : call getghc_mGGA(cwavef_fft,ghc_mGGA_fft,gbound_k1,gs_ham%gprimd,istwf_k_,kg_k_fft,kpt_k1,&
964 : : & gs_ham%mgfft,mpi_enreg,ndat_,gs_ham%ngfft,npw_fft,gs_ham%nvloc,&
965 [ + - - - : 160 : & gs_ham%n4,gs_ham%n5,gs_ham%n6,my_nspinor,gs_ham%vxctaulocal,gs_ham%gpu_option)
- - - - -
- - - - -
- + - - -
- - - - -
- - ]
966 : 160 : call cwavef_double_rfft_trick_unpack(ghc_mGGA,ghc_mGGA_fft,mpi_enreg%me_g0_fft,ndat,npw_k1)
967 [ - + ]: 160 : ABI_FREE(kg_k_fft)
968 [ - + ]: 160 : ABI_FREE(cwavef_fft)
969 [ - + ]: 160 : ABI_FREE(ghc_mGGA_fft)
970 : : else
971 : : call getghc_mGGA(cwavef,ghc_mGGA,gbound_k1,gs_ham%gprimd,istwf_k_,kg_k1,kpt_k1,&
972 : : & gs_ham%mgfft,mpi_enreg,ndat,gs_ham%ngfft,npw_k1,gs_ham%nvloc,&
973 [ - + - - : 90304 : & gs_ham%n4,gs_ham%n5,gs_ham%n6,my_nspinor,gs_ham%vxctaulocal,gs_ham%gpu_option)
- - - - +
- - - - -
- - - - -
- - - - +
- - - - -
+ - - - -
- - - - -
- ]
974 : : end if
975 : : #ifdef HAVE_OPENMP_OFFLOAD
976 : : !$OMP TARGET UPDATE FROM(ghc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
977 : : #endif
978 : : !LB-2025-12-16: intel19 crashes here when npw_k1 is too big. A solution is to write the loop explicitely
979 : : !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)
980 [ + + ]: 56805312 : do idat=1,npw_k1*my_nspinor*ndat
981 [ + + ]: 170235008 : do ig=1,2
982 : 170144544 : ghc(ig,idat)=ghc(ig,idat)+ghc_mGGA(ig,idat)
983 : : end do
984 : : end do
985 : : #ifdef HAVE_OPENMP_OFFLOAD
986 : : !$OMP TARGET UPDATE TO(ghc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
987 : : #endif
988 [ - + ]: 90464 : ABI_FREE(ghc_mGGA)
989 : : end if
990 : :
991 : : ! Add nuclear dipole moment contribution
992 [ + + ]: 35261495 : if (associated(gs_ham%vectornd)) then
993 [ - + ]: 674090 : ABI_CHECK(k1_eq_k2, 'nuclear dipole vector potential not allowed for k/=k_^prime!')
994 [ + + - + ]: 4044540 : if (size(gs_ham%vectornd)/=gs_ham%n4*gs_ham%n5*gs_ham%n6*gs_ham%nvloc*3) then
995 : 0 : ABI_BUG('wrong sizes for vectornd in getghc!')
996 : : end if
997 [ - + + - : 2022270 : ABI_MALLOC(ghc_vectornd,(2,npw_k1*my_nspinor*ndat))
+ - + - -
+ - + ]
998 [ + + + + ]: 998334854 : ghc_vectornd=zero
999 : : #ifdef HAVE_OPENMP_OFFLOAD
1000 : : !$OMP TARGET UPDATE FROM(cwavef) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1001 : : #endif
1002 : : call getghc_nucdip(cwavef,ghc_vectornd,gbound_k1,istwf_k_,kg_k1,kpt_k1,&
1003 : : gs_ham%mgfft,mpi_enreg,ndat,gs_ham%ngfft,npw_k1,gs_ham%nvloc,&
1004 [ - + - - : 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)
- - - - +
- - - - -
- - - - -
- - - + -
- - - - -
- - - - -
- + - - -
- - + - -
- - - - -
- - - - +
- - - - -
- - - ]
1005 : :
1006 : : #ifdef HAVE_OPENMP_OFFLOAD
1007 : : !$OMP TARGET UPDATE FROM(ghc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1008 : : #endif
1009 : : !LB-2025-12-16: intel19 crashes here when npw_k1 is too big. A solution is to write the loop explicitely
1010 : : !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)
1011 [ + + ]: 333227678 : do idat=1,npw_k1*my_nspinor*ndat
1012 [ + + ]: 998334854 : do ig=1,2
1013 : 997660764 : ghc(ig,idat)=ghc(ig,idat)+ghc_vectornd(ig,idat)
1014 : : end do
1015 : : end do
1016 : : #ifdef HAVE_OPENMP_OFFLOAD
1017 : : !$OMP TARGET UPDATE TO(ghc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1018 : : #endif
1019 : :
1020 [ - + ]: 674090 : ABI_FREE(ghc_vectornd)
1021 : : end if
1022 : :
1023 : : ! If only local part is applied, still we have to filter the result because of dilatmx.
1024 : : ! Otherwise it is done when adding kinetic term.
1025 [ + + + + ]: 35261495 : if (type_calc==1.and.filter_dilatmx_loc_) then
1026 : :
1027 [ - + ]: 204650 : ABI_CHECK(gs_ham%use_gbt == 0, "use_gbt not coded!")
1028 : :
1029 [ + - ]: 204650 : if(gs_ham%gpu_option == ABI_GPU_OPENMP) then
1030 : : #ifdef HAVE_OPENMP_OFFLOAD
1031 : : !$OMP TARGET TEAMS DISTRIBUTE MAP(to:ghc,kinpw_k2)
1032 : : do idat=1,ndat
1033 : : !$OMP PARALLEL DO PRIVATE(igspinor) COLLAPSE(2)
1034 : : do ispinor=1,my_nspinor
1035 : : do ig=1,npw_k2
1036 : : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1037 : : if(kinpw_k2(ig)>hugevalue) ghc(:,igspinor)=zero
1038 : : end do ! ig
1039 : : end do ! ispinor
1040 : : end do
1041 : : #endif
1042 : : else
1043 [ + + ]: 1014410 : do idat=1,ndat
1044 : : !$OMP PARALLEL DO PRIVATE(igspinor) COLLAPSE(2) IF(gemm_nonlop_use_gemm)
1045 [ + + ]: 2016282 : do ispinor=1,my_nspinor
1046 [ + + ]: 125631396 : do ig=1,npw_k2
1047 : 123819764 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1048 [ + + + + ]: 129219672 : if(kinpw_k2(ig)>hugevalue) ghc(:,igspinor)=zero
1049 : : end do ! ig
1050 : : end do ! ispinor
1051 : : end do
1052 : : end if
1053 : : end if
1054 : :
1055 : : end if ! type_calc
1056 : :
1057 : : ABI_NVTX_END_RANGE()
1058 : :
1059 [ + + - + ]: 35264769 : if (any(type_calc == [0, 2, 3])) then
1060 : : !============================================================
1061 : : ! Application of the non-local potential and the Fock potential
1062 : : !============================================================
1063 : :
1064 : : ABI_NVTX_START_RANGE(NVTX_GETGHC_NLOCPOT)
1065 : :
1066 [ + - ]: 35056845 : if (type_calc==0 .or. type_calc==2) then
1067 : 35056845 : signs=2 ; choice=1 ; nnlout=1 ; idir=0 ; tim_nonlop=1
1068 [ + + ]: 35056845 : cpopt_here=-1;if (gs_ham%usepaw==1) cpopt_here=cpopt
1069 [ + + ]: 35056845 : if (has_fock) then
1070 [ + + ]: 142256 : if (gs_ham%usepaw==1) then
1071 : 13114 : cpopt_here=max(cpopt,0)
1072 [ + - ]: 13114 : if (cpopt<2) then
1073 [ - + + - : 91798 : ABI_MALLOC(cwaveprj_fock,(gs_ham%natom,my_nspinor*ndat))
- + + - +
- + - - +
- + + + +
+ ]
1074 [ - + + - : 39342 : ABI_MALLOC(dimcprj,(gs_ham%natom))
+ - + - -
+ - + -
+ ]
1075 [ - + ]: 13114 : call pawcprj_getdim(dimcprj,gs_ham%natom,gs_ham%nattyp,gs_ham%ntypat,gs_ham%typat,fock%pawtab,'O')
1076 : 13114 : call pawcprj_alloc(cwaveprj_fock,0,dimcprj)
1077 [ - + ]: 13114 : ABI_FREE(dimcprj)
1078 : : else
1079 : 0 : cwaveprj_fock=>cwaveprj
1080 : : end if
1081 : 13114 : cwaveprj_nonlop=>cwaveprj_fock
1082 : : else
1083 : 129142 : cwaveprj_nonlop=>cwaveprj
1084 : 129142 : cwaveprj_fock=>cwaveprj
1085 : : end if
1086 : : else
1087 : 34914589 : cwaveprj_nonlop=>cwaveprj
1088 : : end if
1089 [ + + ]: 35056845 : paw_opt=gs_ham%usepaw ; if (sij_opt/=0) paw_opt=sij_opt+3
1090 [ + + ]: 71709526 : lambda_ndat = lambda
1091 : :
1092 [ + + ]: 35056845 : if (gs_ham%use_gbt == 0) then
1093 [ + + ]: 34852199 : if (gs_ham%usepaw==0) gsc_ptr => nonlop_dum
1094 [ + + ]: 34852199 : if (gs_ham%usepaw==1) gsc_ptr => gsc
1095 : :
1096 : : call nonlop(choice,cpopt_here,cwaveprj_nonlop,enlout,gs_ham,idir,lambda_ndat,mpi_enreg,ndat,&
1097 : 34852199 : nnlout,paw_opt,signs,gsc_ptr,tim_nonlop,cwavef,gvnlxc_,select_k=select_k_)
1098 : : else
1099 : : ! GBT case. Treat up and down components separately.
1100 : 204646 : gs_ham%nspinor = 1
1101 : :
1102 : : ! Split cwavef and gvnlxc
1103 [ - + + - : 613938 : ABI_MALLOC(cwavef_spin, (2, npw_k1*ndat))
+ - + - -
+ - + -
+ ]
1104 [ - + + - : 409292 : ABI_MALLOC(gvnlxc_spin, (2, npw_k1*ndat))
- + - + -
+ ]
1105 : :
1106 [ + - ]: 204646 : if (gs_ham%usepaw==0) gsc_ptr => nonlop_dum
1107 [ - + ]: 204646 : if (gs_ham%usepaw==1) gsc_ptr => gsc
1108 : :
1109 : : ! Apply Vnl{k-q/2} to u^up
1110 : 204646 : gs_ham%ispin_gbt = 1
1111 [ - + - - : 204646 : call cg_copy_spin(1, npw_k1, nspinortot, ndat, cwavef, cwavef_spin)
- - - - -
+ - - -
- ]
1112 : : call nonlop(choice, cpopt_here, cwaveprj_nonlop, enlout, gs_ham, idir, lambda_ndat, mpi_enreg, ndat, &
1113 : 204646 : nnlout, paw_opt, signs, gsc_ptr, tim_nonlop, cwavef_spin, gvnlxc_spin, select_k=K_H_K)
1114 : : ! Insert results in the right position.
1115 : 204646 : call cg_put_spin(1, npw_k1, nspinortot, ndat, gvnlxc_spin, gvnlxc_)
1116 : :
1117 : : ! Apply H_{k+q/2} to u^down
1118 : : !gvnlxc_spin = zero
1119 : 204646 : gs_ham%ispin_gbt = 2
1120 [ - + - - : 204646 : call cg_copy_spin(2, npw_k1, nspinortot, ndat, cwavef, cwavef_spin)
- - - - -
+ - - -
- ]
1121 : : call nonlop(choice, cpopt_here, cwaveprj_nonlop, enlout, gs_ham, idir, lambda_ndat, mpi_enreg, ndat, &
1122 : 204646 : nnlout, paw_opt, signs, gsc_ptr, tim_nonlop, cwavef_spin, gvnlxc_spin, select_k=KPRIME_H_KPRIME)
1123 : : ! Insert results in the right position.
1124 : 204646 : call cg_put_spin(2, npw_k1, nspinortot, ndat, gvnlxc_spin, gvnlxc_)
1125 : :
1126 : 204646 : gs_ham%nspinor = 2
1127 [ - + ]: 204646 : ABI_FREE(cwavef_spin)
1128 [ - + ]: 204646 : ABI_FREE(gvnlxc_spin)
1129 : : end if ! use_gbt
1130 : :
1131 [ + + + + ]: 35056845 : if (gs_ham%usepaw==1 .and. has_fock)then
1132 [ + - ]: 13114 : if (fock_get_getghc_call(fock)==1) then
1133 [ + - ]: 13114 : if(gs_ham%gpu_option==ABI_GPU_KOKKOS) then
1134 : : #if defined HAVE_GPU && defined HAVE_YAKL
1135 : : ABI_MALLOC_MANAGED(gvnlc, (/2,npw_k2*my_nspinor*ndat/))
1136 : : #endif
1137 : : else
1138 [ - + + - : 39342 : ABI_MALLOC(gvnlc, (2,npw_k2*my_nspinor*ndat))
+ - + - -
+ - + ]
1139 : : #if defined HAVE_OPENMP_OFFLOAD
1140 : : !$OMP TARGET ENTER DATA MAP(to:gvnlc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1141 : : #endif
1142 : : end if
1143 [ - + ]: 13114 : if(gs_ham%gpu_option==ABI_GPU_OPENMP) then
1144 : 0 : call gpu_copy(gvnlc, gvnlxc_, int(2,c_size_t)*npw_k2*my_nspinor*ndat)
1145 : : else
1146 [ + - + + : 48885164 : gvnlc=gvnlxc_
+ + + - +
- - + - -
- - - - +
+ + + ]
1147 : : end if
1148 : : endif
1149 : : endif
1150 : :
1151 : : ! Calculation of the Fock exact exchange contribution from the Fock or ACE operator
1152 [ + + ]: 26244398 : if (has_fock) then
1153 [ + - ]: 142256 : if (fock_get_getghc_call(fock)==1) then
1154 : : if (gs_ham%usepaw==0) cwaveprj_idat => cwaveprj
1155 [ - + ]: 142256 : if (fock%use_ACE==0) then
1156 : : #if defined HAVE_OPENMP_OFFLOAD
1157 : : !$OMP TARGET UPDATE FROM(cwavef,gvnlxc_) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1158 : : #endif
1159 : 0 : call timab(360,1,tsec)
1160 : 0 : call fock_getghc(cwavef,cwaveprj,gvnlxc_,gs_ham,mpi_enreg,ndat)
1161 : 0 : call timab(360,2,tsec)
1162 : : #if defined HAVE_OPENMP_OFFLOAD
1163 : : !$OMP TARGET UPDATE TO(cwavef,gvnlxc_) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1164 : : #endif
1165 : : else
1166 : 142256 : call fock_ACE_getghc(cwavef,gvnlxc_,gs_ham,mpi_enreg,ndat)
1167 : : end if
1168 : : end if
1169 : : end if
1170 : :
1171 [ # # ]: 0 : else if (type_calc == 3) then
1172 : : ! for kinetic and local only, nonlocal and vfock should be zero
1173 [ # # ]: 0 : if(gs_ham%gpu_option==ABI_GPU_OPENMP) then
1174 : 0 : call gpu_set_to_zero(gvnlxc_, int(2,c_size_t)*npw_k2*my_nspinor*ndat)
1175 : : else
1176 [ # # # # ]: 0 : gvnlxc_(:,:) = zero
1177 : : end if
1178 : : end if ! if(type_calc...
1179 : :
1180 : : ABI_NVTX_END_RANGE()
1181 : :
1182 : : !============================================================
1183 : : ! Assemble kinetic, local, nonlocal and Fock contributions
1184 : : !============================================================
1185 : :
1186 : : ABI_NVTX_START_RANGE(NVTX_GETGHC_KIN)
1187 : :
1188 : : #ifdef FC_NVHPC
1189 : : !FIXME This Kokkos kernel seems to cause issues under NVHPC so it is disabled
1190 : : if (.false.) then
1191 : : #else
1192 [ + - ]: 35056845 : if (gs_ham%gpu_option == ABI_GPU_KOKKOS) then
1193 : : #endif
1194 : :
1195 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_KOKKOS)
1196 : : call assemble_energy_contribution_kokkos(c_loc(ghc), &
1197 : : & c_loc(gsc), c_loc(kinpw_k2), c_loc(cwavef), c_loc(gvnlxc_), &
1198 : : & ndat, my_nspinor, npw_k2, sij_opt, k1_eq_k2, hugevalue)
1199 : : ! sync device so that data can be reused safely on host
1200 : : ! will probably be moved elsewhere once all the scf loop runs on device
1201 : : call gpu_device_synchronize()
1202 : : #endif
1203 : :
1204 : : else
1205 : :
1206 : : #ifdef FC_NVHPC
1207 : : #if defined(HAVE_GPU_CUDA) && defined(HAVE_KOKKOS)
1208 : : !Related to FIXME above
1209 : : if (gs_ham%gpu_option == ABI_GPU_KOKKOS) call gpu_device_synchronize()
1210 : : #endif
1211 : : #endif
1212 : :
1213 : : ! Assemble modified kinetic, local and nonlocal contributions
1214 : : ! to <G|H|C(n,k)>. Take also into account built-in debugging.
1215 [ + - ]: 35056845 : if(prtvol/=-level)then
1216 : :
1217 [ + - ]: 35056845 : if (gs_ham%gpu_option == ABI_GPU_OPENMP) then
1218 : : ! OpenMP GPU
1219 : : #ifdef HAVE_OPENMP_OFFLOAD
1220 : : if (k1_eq_k2) then
1221 : :
1222 : : !MG 20260102: With nvfortran 23.11-0, this kernel is a bottleneck due to the filter on kinpw_k2 and the update of ghc.
1223 : : ! Solution: branch-free mask + manual loop unrolling. loop unrolling is crucial.
1224 : : ! The version with COLLAPSE(3) is faster.
1225 : :
1226 : : ! !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:ghc,kinpw_k2,gvnlxc_,gsc,cwavef)
1227 : : ! !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(3) PRIVATE(igspinor, kscale) MAP(to:ghc,kinpw_k2,gvnlxc_,gsc,cwavef)
1228 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:ghc,kinpw_k2,gvnlxc_,cwavef)
1229 : : do idat=1,ndat
1230 : : do ispinor=1,my_nspinor
1231 : : ! !$OMP PARALLEL DO PRIVATE(igspinor, kscale)
1232 : : !$OMP PARALLEL DO PRIVATE(igspinor)
1233 : : do ig=1,npw_k2
1234 : : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1235 : :
1236 : : ! ! New version with branch-free mask + manual loop unrolling.
1237 : : ! kscale = merge(one, zero, kinpw_k2(ig) < hugevalue)
1238 : :
1239 : : ! ghc(1,igspinor) = kscale * (ghc(1,igspinor) + kinpw_k2(ig)*cwavef(1,igspinor) + gvnlxc_(1,igspinor))
1240 : : ! ghc(2,igspinor) = kscale * (ghc(2,igspinor) + kinpw_k2(ig)*cwavef(2,igspinor) + gvnlxc_(2,igspinor))
1241 : : ! if (sij_opt == 1) then
1242 : : ! gsc(1,igspinor) = kscale * gsc(1,igspinor)
1243 : : ! gsc(2,igspinor) = kscale * gsc(2,igspinor)
1244 : : ! end if
1245 : : if(kinpw_k2(ig)<huge(zero)*1.d-11)then
1246 : : ghc(1,igspinor) = ghc(1,igspinor) + kinpw_k2(ig)*cwavef(1,igspinor) + gvnlxc_(1,igspinor)
1247 : : ghc(2,igspinor) = ghc(2,igspinor) + kinpw_k2(ig)*cwavef(2,igspinor) + gvnlxc_(2,igspinor)
1248 : : else
1249 : : ghc(1,igspinor)=zero
1250 : : ghc(2,igspinor)=zero
1251 : : end if
1252 : : end do ! ig
1253 : : end do ! ispinor
1254 : : end do ! idat
1255 : : if(sij_opt==1) then
1256 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:gsc,kinpw_k2)
1257 : : do idat=1,ndat
1258 : : do ispinor=1,my_nspinor
1259 : : !$OMP PARALLEL DO PRIVATE(igspinor)
1260 : : do ig=1,npw_k2
1261 : : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1262 : : if(kinpw_k2(ig)>=huge(zero)*1.d-11)then
1263 : : gsc(1,igspinor) = zero
1264 : : gsc(2,igspinor) = zero
1265 : : end if
1266 : : end do ! ig
1267 : : end do ! ispinor
1268 : : end do ! idat
1269 : : end if
1270 : : else
1271 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:ghc,gvnlxc_,kinpw_k2)
1272 : : do idat=1,ndat
1273 : : do ispinor=1,my_nspinor
1274 : : !$OMP PARALLEL DO PRIVATE(igspinor)
1275 : : do ig=1,npw_k2
1276 : : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1277 : : if(kinpw_k2(ig)<hugevalue)then
1278 : : ghc(1,igspinor)= ghc(1,igspinor) + gvnlxc_(1,igspinor)
1279 : : ghc(2,igspinor)= ghc(2,igspinor) + gvnlxc_(2,igspinor)
1280 : : else
1281 : : ghc(1,igspinor)=zero
1282 : : ghc(2,igspinor)=zero
1283 : : end if
1284 : : end do ! ig
1285 : : end do ! ispinor
1286 : : end do ! idat
1287 : : if(sij_opt==1) then
1288 : : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:gsc,kinpw_k2)
1289 : : do idat=1,ndat
1290 : : do ispinor=1,my_nspinor
1291 : : !$OMP PARALLEL DO PRIVATE(igspinor)
1292 : : do ig=1,npw_k2
1293 : : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1294 : : if(kinpw_k2(ig)>=huge(zero)*1.d-11)then
1295 : : gsc(1,igspinor)=zero
1296 : : gsc(2,igspinor)=zero
1297 : : end if
1298 : : end do ! ig
1299 : : end do ! ispinor
1300 : : end do ! idat
1301 : : end if
1302 : : end if
1303 : : #endif
1304 : :
1305 : : else
1306 : :
1307 : : !CPU (+ Kokkos eventually)
1308 [ + + ]: 35056845 : if (k1_eq_k2) then
1309 : : !$OMP PARALLEL DO PRIVATE(igspinor) COLLAPSE(2) IF(gemm_nonlop_use_gemm)
1310 [ + + ]: 71300234 : do idat=1,ndat
1311 [ + + ]: 110245201 : do ispinor=1,my_nspinor
1312 [ + + ]: 11581183080 : do ig=1,npw_k2
1313 : 11505790078 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1314 [ + + ]: 11544735045 : if(kinpw_k2(ig)<hugevalue)then
1315 [ + - + + : 57031556020 : ghc(:,igspinor) = ghc(:,igspinor) + kinpw_k2(ig)*cwavef(:,igspinor) + gvnlxc_(:,igspinor)
+ + ]
1316 : : else
1317 [ + + ]: 298436622 : ghc(:,igspinor)=zero
1318 [ + + + + ]: 142274576 : if (sij_opt==1) gsc(:,igspinor)=zero
1319 : : end if
1320 : : end do ! ig
1321 : : end do ! ispinor
1322 : : end do ! idat
1323 : : else
1324 [ - + ]: 204646 : if (gs_ham%use_gbt == 0) then
1325 : : !$OMP PARALLEL DO PRIVATE(igspinor) COLLAPSE(2) IF(gemm_nonlop_use_gemm)
1326 [ # # ]: 0 : do idat=1,ndat
1327 [ # # ]: 0 : do ispinor=1,my_nspinor
1328 [ # # ]: 0 : do ig=1,npw_k2
1329 : 0 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1330 [ # # ]: 0 : if(kinpw_k2(ig)<hugevalue)then
1331 [ # # # # : 0 : ghc(:,igspinor)= ghc(:,igspinor) + gvnlxc_(:,igspinor)
# # ]
1332 : : else
1333 [ # # ]: 0 : ghc(:,igspinor)=zero
1334 [ # # # # ]: 0 : if (sij_opt==1) gsc(:,igspinor)=zero
1335 : : end if
1336 : : end do ! ig
1337 : : end do ! ispinor
1338 : : end do ! idat
1339 : : else
1340 : : ! GBT. use different kinetic energies for up and down components.
1341 : : !$OMP PARALLEL DO PRIVATE(igspinor, iispinor) COLLAPSE(2) IF(gemm_nonlop_use_gemm)
1342 [ + + ]: 409292 : do idat=1,ndat
1343 [ + + ]: 818584 : do ispinor=1,my_nspinor
1344 [ - + ]: 409292 : iispinor=ispinor; if (mpi_enreg%paral_spinor==1) iispinor=mpi_enreg%me_spinor+1
1345 [ + + ]: 613938 : if (iispinor == 1) then
1346 [ + + ]: 73416399 : do ig=1,npw_k2
1347 : 73211753 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1348 [ + + ]: 73416399 : if(gs_ham%kinpw_k(ig)<hugevalue)then
1349 [ + - + + : 357781890 : ghc(:,igspinor) = ghc(:,igspinor) + gs_ham%kinpw_k(ig)*cwavef(:,igspinor) + gvnlxc_(:,igspinor)
+ + ]
1350 : : else
1351 [ + + ]: 4966125 : ghc(:,igspinor)=zero
1352 [ - + - - ]: 1655375 : if (sij_opt==1) gsc(:,igspinor)=zero
1353 : : end if
1354 : : end do ! ig
1355 : : else
1356 [ + + ]: 73416399 : do ig=1,npw_k2
1357 : 73211753 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1358 [ + + ]: 73416399 : if(gs_ham%kinpw_kp(ig)<hugevalue)then
1359 [ + - + + : 357781890 : ghc(:,igspinor) = ghc(:,igspinor) + gs_ham%kinpw_kp(ig)*cwavef(:,igspinor) + gvnlxc_(:,igspinor)
+ + ]
1360 : : else
1361 [ + + ]: 4966125 : ghc(:,igspinor)=zero
1362 [ - + - - ]: 1655375 : if (sij_opt==1) gsc(:,igspinor)=zero
1363 : : end if
1364 : : end do ! ig
1365 : : end if
1366 : : end do ! ispinor
1367 : : end do ! idat
1368 : : end if
1369 : : end if
1370 : :
1371 : : end if ! gs_ham%gpu_option
1372 : :
1373 : : else
1374 : : ! Here, debugging section
1375 : : #ifdef HAVE_OPENMP_OFFLOAD
1376 : : !$OMP TARGET UPDATE FROM(ghc,gsc,cwavef,gvnlxc_) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1377 : : #endif
1378 : 0 : call wrtout(std_out,' getghc : components of ghc ','PERS')
1379 : : write(msg,'(a)')&
1380 : 0 : 'icp ig ispinor igspinor re/im ghc kinpw cwavef glocc gvnlxc gsc'
1381 : 0 : call wrtout(std_out,msg,'PERS')
1382 [ # # ]: 0 : do idat=1,ndat
1383 [ # # ]: 0 : do ispinor=1,my_nspinor
1384 [ # # ]: 0 : do ig=1,npw_k2
1385 : 0 : igspinor=ig+npw_k2*(ispinor-1)+npw_k2*my_nspinor*(idat-1)
1386 [ # # ]: 0 : if(kinpw_k2(ig)<hugevalue)then
1387 [ # # ]: 0 : if (k1_eq_k2) then
1388 : 0 : ghcre=kinpw_k2(ig)*cwavef(1,igspinor)+ghc(1,igspinor)+gvnlxc_(1,igspinor)
1389 : 0 : ghcim=kinpw_k2(ig)*cwavef(2,igspinor)+ghc(2,igspinor)+gvnlxc_(2,igspinor)
1390 : : else
1391 : 0 : ghcre=ghc(1,igspinor)+gvnlxc_(1,igspinor)
1392 : 0 : ghcim=ghc(2,igspinor)+gvnlxc_(2,igspinor)
1393 : : end if
1394 : : else
1395 : 0 : ghcre=zero
1396 : 0 : ghcim=zero
1397 [ # # # # ]: 0 : if (sij_opt==1) gsc(:,igspinor)=zero
1398 : : end if
1399 [ # # ]: 0 : iispinor=ispinor; if (mpi_enreg%paral_spinor==1) iispinor=mpi_enreg%me_spinor+1
1400 [ # # ]: 0 : if (sij_opt == 1) then
1401 : 0 : write(msg,'(a,3(1x,i5),6(1x,es13.6))') ' 1 ', ig, iispinor, igspinor,ghcre,&
1402 : 0 : kinpw_k2(ig),cwavef(1,igspinor),ghc(1,igspinor),gvnlxc_(1,igspinor), gsc(1,igspinor)
1403 : 0 : call wrtout(std_out,msg,'PERS')
1404 : 0 : write(msg,'(a,3(1x,i5),6(1x,es13.6))') ' 2 ', ig, iispinor, igspinor,ghcim,&
1405 : 0 : kinpw_k2(ig),cwavef(2,igspinor),ghc(2,igspinor),gvnlxc_(2,igspinor), gsc(2,igspinor)
1406 : 0 : call wrtout(std_out,msg,'PERS')
1407 : : else
1408 : 0 : write(msg,'(a,3(1x,i5),6(1x,es13.6))') ' 1 ', ig, iispinor, igspinor,ghcre,&
1409 : 0 : kinpw_k2(ig),cwavef(1,igspinor),ghc(1,igspinor),gvnlxc_(1,igspinor)
1410 : 0 : call wrtout(std_out,msg,'PERS')
1411 : 0 : write(msg,'(a,3(1x,i5),6(1x,es13.6))') ' 2 ', ig, iispinor, igspinor,ghcim,&
1412 : 0 : kinpw_k2(ig),cwavef(2,igspinor),ghc(2,igspinor),gvnlxc_(2,igspinor)
1413 : 0 : call wrtout(std_out,msg,'PERS')
1414 : : end if
1415 [ # # ]: 0 : ghc(:,igspinor) = [ghcre, ghcim]
1416 : : end do ! ig
1417 : : end do ! ispinor
1418 : : end do ! idat
1419 : : #ifdef HAVE_OPENMP_OFFLOAD
1420 : : !$OMP TARGET UPDATE TO(ghc,gsc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1421 : : #endif
1422 : : end if
1423 : : end if ! gs_ham%gpu_option
1424 : :
1425 : : ABI_NVTX_END_RANGE()
1426 : :
1427 : : ! Special case of PAW + Fock : only return Fock operator contribution in gvnlxc_
1428 [ + + + + ]: 35056845 : if (gs_ham%usepaw==1 .and. has_fock) then
1429 [ + - ]: 13114 : if(gs_ham%gpu_option == ABI_GPU_OPENMP) then
1430 : : #ifdef HAVE_OPENMP_OFFLOAD
1431 : : call abi_xaxpy(2*npw_k2*my_nspinor*ndat, cminusone, gvnlc, 1, gvnlxc_, 1, x_cplx=1, gpu_option=gs_ham%gpu_option)
1432 : : !$OMP TARGET EXIT DATA MAP(delete:gvnlc) IF(gs_ham%gpu_option == ABI_GPU_OPENMP)
1433 : : #endif
1434 : : else
1435 [ + - + + : 48872050 : gvnlxc_=gvnlxc_-gvnlc
+ + + + +
+ ]
1436 : : end if
1437 [ - + ]: 13114 : if(gs_ham%gpu_option==ABI_GPU_KOKKOS) then
1438 : : #if defined HAVE_GPU && defined HAVE_YAKL
1439 : : ABI_FREE_MANAGED(gvnlc)
1440 : : #endif
1441 : : else
1442 [ - + ]: 13114 : ABI_FREE(gvnlc)
1443 : : end if
1444 : : endif
1445 : :
1446 : : #ifdef HAVE_OPENMP_OFFLOAD
1447 : : !$OMP TARGET UPDATE FROM(ghc) IF(transfer_ghc)
1448 : : !$OMP TARGET UPDATE FROM(gsc) IF(transfer_gsc)
1449 : : !$OMP TARGET UPDATE FROM(cwavef) IF(transfer_cwavef)
1450 : : !$OMP TARGET UPDATE FROM(gvnlxc_) IF(transfer_gvnlxc .and. .not. local_gvnlxc)
1451 : : !$OMP TARGET EXIT DATA MAP(delete:gvnlxc_) IF(transfer_gvnlxc)
1452 : : #endif
1453 [ + + ]: 35056845 : if (local_gvnlxc) then
1454 [ + - ]: 912524 : if(gs_ham%gpu_option==ABI_GPU_KOKKOS) then
1455 : : #if defined HAVE_GPU && defined HAVE_YAKL
1456 : : ABI_FREE_MANAGED(gvnlxc_)
1457 : : #endif
1458 : : else
1459 [ - + ]: 912524 : ABI_FREE(gvnlxc_)
1460 : : end if
1461 : : end if
1462 : :
1463 : : ! Structured debugging: if prtvol=-level, stop here.
1464 [ - + ]: 35056845 : if (prtvol == -level) then
1465 : 0 : ABI_ERROR(sjoin(' getghc: exit prtvol=-',itoa(level),', debugging mode => stop '))
1466 : : end if
1467 : :
1468 [ - + ]: 35056845 : if (type_calc==0.or.type_calc==2) then
1469 [ + + + + : 35056845 : if (has_fock.and.gs_ham%usepaw==1.and.cpopt<2) then
- + ]
1470 : 13114 : call pawcprj_free(cwaveprj_fock)
1471 [ - + ]: 13114 : ABI_FREE(cwaveprj_fock)
1472 : : end if
1473 : : end if
1474 : :
1475 : : end if ! type_calc
1476 : :
1477 : : #ifdef HAVE_OPENMP_OFFLOAD
1478 : : !$OMP TARGET EXIT DATA MAP(delete:ghc) IF(transfer_ghc)
1479 : : !$OMP TARGET EXIT DATA MAP(delete:gsc) IF(transfer_gsc)
1480 : : !$OMP TARGET EXIT DATA MAP(delete:cwavef) IF(transfer_cwavef)
1481 : : #endif
1482 : 35264769 : call timab(350+tim_getghc,2,tsec)
1483 : :
1484 : : DBG_EXIT("COLL")
1485 : : ABI_NVTX_END_RANGE()
1486 : :
1487 [ - + - + : 70529538 : end subroutine getghc
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + ]
1488 : : !!***
1489 : :
1490 : : !----------------------------------------------------------------------
1491 : :
1492 : : !!****f* ABINIT/getghc_nucdip
1493 : : !!
1494 : : !! NAME
1495 : : !! getghc_nucdip
1496 : : !!
1497 : : !! FUNCTION
1498 : : !! Compute magnetic nuclear dipole moment contribution to <G|H|C>
1499 : : !! for input vector |C> expressed in reciprocal space.
1500 : : !!
1501 : : !! INPUTS
1502 : : !! cwavef(2,npw_k*my_nspinor*ndat)=planewave coefficients of wavefunction.
1503 : : !! gbound_k(2*mgfft+4)=sphere boundary info
1504 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1)
1505 : : !! istwf_k=input parameter that describes the storage of wfs
1506 : : !! kg_k(3,npw_k)=G vec coordinates wrt recip lattice transl.
1507 : : !! kpt(3)=current k point
1508 : : !! mgfft=maximum single fft dimension
1509 : : !! mpi_enreg=information about MPI parallelization
1510 : : !! my_nspinor=number of spinorial components of the wavefunctions (on current proc)
1511 : : !! ndat=number of FFTs to perform in parall
1512 : : !! ngfft(18)=contain all needed information about 3D FFT
1513 : : !! npw_k=number of planewaves in basis for given k point.
1514 : : !! nvloc=number of spin components of vxctaulocal
1515 : : !! n4,n5,n6=for dimensioning of vxctaulocal
1516 : : !! gpu_option= GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
1517 : : !! vectornd(n4,n5,n6,nvloc,3)= local potential corresponding to the vector potential of the array
1518 : : !! of nuclear magnetic dipoles, in real space, on the augmented fft grid.
1519 : : !!
1520 : : !! OUTPUT
1521 : : !! ghc_vectornd(2,npw_k*my_nspinor*ndat)=A.p contribution to <G|H|C> for array of nuclear dipoles
1522 : : !!
1523 : : !! SIDE EFFECTS
1524 : : !!
1525 : : !! NOTES
1526 : : !! this code is a copied, simplified version of getghc_mGGA (see below) and should eventually be
1527 : : !! integrated into that code, to simplify maintenance
1528 : : !!
1529 : : !! SOURCE
1530 : :
1531 : 813630 : subroutine getghc_nucdip(cwavef,ghc_vectornd,gbound_k,istwf_k,kg_k,kpt,mgfft,mpi_enreg,&
1532 : 813630 : & ndat,ngfft,npw_k,nvloc,n4,n5,n6,my_nspinor,vectornd,vlocal,zora,gpu_option)
1533 : :
1534 : : !Arguments ------------------------------------
1535 : : !scalars
1536 : : integer,intent(in) :: istwf_k,mgfft,my_nspinor,ndat,npw_k,nvloc,n4,n5,n6,gpu_option,zora
1537 : : type(MPI_type),intent(in) :: mpi_enreg
1538 : : !arrays
1539 : : integer,intent(in) :: gbound_k(2*mgfft+4),kg_k(3,npw_k),ngfft(18)
1540 : : real(dp),intent(in) :: kpt(3)
1541 : : real(dp),intent(inout) :: cwavef(2,npw_k*my_nspinor*ndat)
1542 : : real(dp),intent(inout) :: ghc_vectornd(2,npw_k*my_nspinor*ndat)
1543 : : real(dp),intent(inout) :: vectornd(n4,n5,n6,nvloc,3),vlocal(n4,n5,n6,nvloc)
1544 : :
1545 : : !Local variables-------------------------------
1546 : : !scalars
1547 : : integer,parameter :: tim_fourwf=1
1548 : : integer :: idat,idir,ipw,iv1,iv2,nspinortot,shift
1549 : : logical :: nspinor1TreatedByThisProc,nspinor2TreatedByThisProc,usezora
1550 : : real(dp), parameter :: HalfFineStruct2=half/InvFineStruct**2
1551 : : real(dp) :: weight=one
1552 : : !arrays
1553 : 813630 : real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:)
1554 : 813630 : real(dp),allocatable :: gcwavef(:,:,:),gcwavef1(:,:,:),gcwavef2(:,:,:)
1555 : 813630 : real(dp),allocatable :: ghc1(:,:),ghc2(:,:),kgkpk(:,:),vectornd_dir(:,:,:,:)
1556 : 813630 : real(dp),allocatable :: work(:,:,:,:),zk(:,:,:,:)
1557 : : ! *********************************************************************
1558 : :
1559 [ + + + + ]: 1198182702 : ghc_vectornd(:,:)=zero
1560 : :
1561 : : !! JWZ debug initial code was only for nvloc==1 case
1562 [ + - ]: 813630 : if (nvloc/=1) return
1563 : :
1564 : 813630 : nspinortot=min(2,(1+mpi_enreg%paral_spinor)*my_nspinor)
1565 [ + - ]: 813630 : if (mpi_enreg%paral_spinor==0) then
1566 : 813630 : shift=npw_k
1567 : 813630 : nspinor1TreatedByThisProc=.true.
1568 : 813630 : nspinor2TreatedByThisProc=(nspinortot==2)
1569 : : else
1570 : 0 : shift=0
1571 : 0 : nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
1572 : 0 : nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
1573 : : end if
1574 : :
1575 : 813630 : usezora=((zora.EQ.1).OR.(zora.EQ.3))
1576 [ - + ]: 813630 : if(usezora) then
1577 [ # # # # : 0 : ABI_MALLOC(zk,(n4,n5,n6,nvloc))
# # # # #
# # # # #
# # # # #
# # # #
# ]
1578 [ # # # # : 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))
# # # # ]
1579 : : end if
1580 : :
1581 [ - + + - : 4068150 : ABI_MALLOC(work,(2,n4,n5,n6*ndat))
- + + - -
+ + - + -
+ - - + -
+ ]
1582 : :
1583 [ + + ]: 813630 : if (nspinortot==1) then
1584 : :
1585 [ - + + - : 1757772 : ABI_MALLOC(ghc1,(2,npw_k*ndat))
+ - + - -
+ - + ]
1586 : :
1587 : : ! Do it in 2 STEPs:
1588 : : ! STEP1: Compute grad of cwavef
1589 [ - + + - : 2343696 : ABI_MALLOC(gcwavef,(2,npw_k*ndat,3))
+ - + - +
- - + -
+ ]
1590 : :
1591 [ + + + + : 1268270034 : gcwavef = zero
+ + ]
1592 : :
1593 : : ! compute k + G. Note these are in reduced coords
1594 [ - + + - : 1757772 : ABI_MALLOC(kgkpk,(npw_k,3))
+ - + - -
+ - + ]
1595 [ + + ]: 133746065 : do ipw = 1, npw_k
1596 [ + + ]: 533226488 : kgkpk(ipw,:) = kpt(:) + kg_k(:,ipw)
1597 : : end do
1598 : :
1599 : : ! make 2\pi(k+G)c(G)|G> by element-wise multiplication
1600 [ + + ]: 2343696 : do idir = 1, 3
1601 [ + + ]: 4103952 : do idat = 1, ndat
1602 [ + + ]: 425493474 : do ipw=1,npw_k
1603 : 421975446 : gcwavef(1,ipw+(idat-1)*npw_k,idir) = cwavef(1,ipw+(idat-1)*npw_k)*kgkpk(ipw,idir)
1604 : 423735702 : gcwavef(2,ipw+(idat-1)*npw_k,idir) = cwavef(2,ipw+(idat-1)*npw_k)*kgkpk(ipw,idir)
1605 : : end do
1606 : : end do
1607 : : end do
1608 [ - + ]: 585924 : ABI_FREE(kgkpk)
1609 [ + + + + : 1268270034 : gcwavef = gcwavef*two_pi
+ + ]
1610 : :
1611 : : ! STEP2: Compute sum of (grad components of vectornd)*(grad components of cwavef)
1612 [ - + + - : 3515544 : ABI_MALLOC(vectornd_dir,(n4,n5,n6,nvloc))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
1613 [ + + ]: 2343696 : do idir=1,3
1614 [ - + ]: 1757772 : if (usezora) then
1615 [ # # # # : 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)
# # # # ]
1616 : : else
1617 [ + + + + : 11551876524 : vectornd_dir(1:n4,1:n5,1:n6,1:nvloc)=vectornd(1:n4,1:n5,1:n6,1:nvloc,idir)
+ + + + ]
1618 : : end if
1619 : : call fourwf(1,vectornd_dir,gcwavef(:,:,idir),ghc1,work,gbound_k,gbound_k,&
1620 : : istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1621 : 1757772 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1622 : : !!$OMP PARALLEL DO
1623 [ + + ]: 4103952 : do idat=1,ndat
1624 [ + + ]: 425493474 : do ipw=1,npw_k
1625 [ + + ]: 1267686594 : ghc_vectornd(:,ipw+(idat-1)*npw_k)=ghc_vectornd(:,ipw+(idat-1)*npw_k)+ghc1(:,ipw+(idat-1)*npw_k)
1626 : : end do
1627 : : end do
1628 : : !ghc_vectornd=ghc_vectornd+ghc1
1629 : : end do ! idir
1630 [ - + ]: 585924 : ABI_FREE(vectornd_dir)
1631 [ - + ]: 585924 : ABI_FREE(gcwavef)
1632 [ - + ]: 585924 : ABI_FREE(ghc1)
1633 : :
1634 : : else ! nspinortot==2
1635 : :
1636 [ - + + - : 683118 : ABI_MALLOC(cwavef1,(2,npw_k*ndat))
+ - + - -
+ - + ]
1637 [ - + + - : 455412 : ABI_MALLOC(cwavef2,(2,npw_k*ndat))
- + - + ]
1638 [ + + ]: 455412 : do idat=1,ndat
1639 : 227706 : iv1=1+(idat-1)*npw_k; iv2=-1+iv1+npw_k
1640 [ + + + + ]: 387924519 : cwavef1(1:2,iv1:iv2) = cwavef(1:2,1+(idat-1)*my_nspinor*npw_k:npw_k+(idat-1)*my_nspinor*npw_k)
1641 : : cwavef2(1:2,iv1:iv2) = &
1642 [ + + + + ]: 388152225 : & cwavef(1:2,1+(idat-1)*my_nspinor*npw_k+shift:npw_k+(idat-1)*my_nspinor*npw_k+shift)
1643 : : end do
1644 : :
1645 : : ! compute k + G. Note these are in reduced coords
1646 [ - + + - : 683118 : ABI_MALLOC(kgkpk,(npw_k,3))
+ - + - -
+ - + ]
1647 [ + + ]: 129459977 : do ipw = 1, npw_k
1648 [ + + ]: 517156790 : kgkpk(ipw,:) = kpt(:) + kg_k(:,ipw)
1649 : : end do
1650 : :
1651 [ + - ]: 227706 : if (nspinor1TreatedByThisProc) then
1652 : :
1653 [ - + + - : 455412 : ABI_MALLOC(ghc1,(2,npw_k*ndat))
- + - + ]
1654 : :
1655 : : ! Do it in 2 STEPs:
1656 : : ! STEP1: Compute grad of cwavef
1657 [ - + + - : 910824 : ABI_MALLOC(gcwavef1,(2,npw_k*ndat,3))
+ - + - +
- - + -
+ ]
1658 : :
1659 [ + + + + : 1164001263 : gcwavef1 = zero
+ + ]
1660 : : ! make 2\pi(k+G)c(G)|G> by element-wise multiplication
1661 [ + + ]: 910824 : do idir = 1, 3
1662 [ + + ]: 1593942 : do idat = 1, ndat
1663 : 683118 : iv1=1+(idat-1)*npw_k; iv2=-1+iv1+npw_k
1664 [ + + ]: 388379931 : gcwavef1(1,iv1:iv2,idir) = cwavef1(1,iv1:iv2)*kgkpk(1:npw_k,idir)
1665 [ + + ]: 389063049 : gcwavef1(2,iv1:iv2,idir) = cwavef1(2,iv1:iv2)*kgkpk(1:npw_k,idir)
1666 : : end do
1667 : : end do
1668 [ + + + + : 1164001263 : gcwavef1 = gcwavef1*two_pi
+ + ]
1669 : :
1670 : : ! STEP2: Compute sum of (grad components of vectornd)*(grad components of cwavef)
1671 [ - + + - : 1366236 : ABI_MALLOC(vectornd_dir,(n4,n5,n6,nvloc))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
1672 [ + + ]: 910824 : do idir=1,3
1673 [ - + ]: 683118 : if (usezora) then
1674 [ # # # # : 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)
# # # # ]
1675 : : else
1676 [ + + + + : 8256805533 : vectornd_dir(1:n4,1:n5,1:n6,1:nvloc)=vectornd(1:n4,1:n5,1:n6,1:nvloc,idir)
+ + + + ]
1677 : : end if
1678 : : call fourwf(1,vectornd_dir,gcwavef1(:,:,idir),ghc1,work,gbound_k,gbound_k,&
1679 : : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1680 : 683118 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1681 : : !!$OMP PARALLEL DO
1682 [ + + ]: 1593942 : do idat=1,ndat
1683 : 683118 : iv1=1+(idat-1)*npw_k; iv2=-1+iv1+npw_k
1684 : : ghc_vectornd(1:2,iv1:iv2)=ghc_vectornd(1:2,iv1:iv2)+&
1685 [ + + + + ]: 1164456675 : & ghc1(1:2,iv1:iv2)
1686 : : end do
1687 : : end do ! idir
1688 [ - + ]: 227706 : ABI_FREE(vectornd_dir)
1689 [ - + ]: 227706 : ABI_FREE(gcwavef1)
1690 [ - + ]: 227706 : ABI_FREE(ghc1)
1691 : :
1692 : : end if ! end spinor 1
1693 : :
1694 [ + - ]: 227706 : if (nspinor2TreatedByThisProc) then
1695 : :
1696 [ - + + - : 455412 : ABI_MALLOC(ghc2,(2,npw_k*ndat))
- + - + ]
1697 : :
1698 : : ! Do it in 2 STEPs:
1699 : : ! STEP1: Compute grad of cwavef
1700 [ - + + - : 910824 : ABI_MALLOC(gcwavef2,(2,npw_k*ndat,3))
+ - + - +
- - + -
+ ]
1701 [ + + + + : 1164001263 : gcwavef2 = zero
+ + ]
1702 : : ! make 2\pi(k+G)c(G)|G> by element-wise multiplication
1703 [ + + ]: 910824 : do idir = 1, 3
1704 [ + + ]: 1593942 : do idat = 1, ndat
1705 : 683118 : iv1=1+(idat-1)*npw_k; iv2=-1+iv1+npw_k
1706 [ + + ]: 388379931 : gcwavef2(1,iv1:iv2,idir) = cwavef2(1,iv1:iv2)*kgkpk(1:npw_k,idir)
1707 [ + + ]: 389063049 : gcwavef2(2,iv1:iv2,idir) = cwavef2(2,iv1:iv2)*kgkpk(1:npw_k,idir)
1708 : : end do
1709 : : end do
1710 [ + + + + : 1164001263 : gcwavef2 = gcwavef2*two_pi
+ + ]
1711 : :
1712 : : ! STEP2: Compute sum of (grad components of vectornd)*(grad components of cwavef)
1713 [ - + + - : 1366236 : ABI_MALLOC(vectornd_dir,(n4,n5,n6,nvloc))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
1714 [ + + ]: 910824 : do idir=1,3
1715 [ - + ]: 683118 : if (usezora) then
1716 [ # # # # : 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)
# # # # ]
1717 : : else
1718 [ + + + + : 8256805533 : vectornd_dir(1:n4,1:n5,1:n6,1:nvloc)=vectornd(1:n4,1:n5,1:n6,1:nvloc,idir)
+ + + + ]
1719 : : end if
1720 : : call fourwf(1,vectornd_dir,gcwavef2(:,:,idir),ghc2,work,gbound_k,gbound_k,&
1721 : : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1722 : 683118 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1723 : : !!$OMP PARALLEL DO
1724 [ + + ]: 1593942 : do idat=1,ndat
1725 : 683118 : iv1=1+(idat-1)*npw_k; iv2=-1+iv1+npw_k
1726 : : ghc_vectornd(1:2,iv1+shift:iv2+shift)=ghc_vectornd(1:2,iv1+shift:iv2+shift)+&
1727 [ + + + + ]: 1164456675 : & ghc2(1:2,iv1:iv2)
1728 : : end do
1729 : : end do ! idir
1730 [ - + ]: 227706 : ABI_FREE(vectornd_dir)
1731 [ - + ]: 227706 : ABI_FREE(gcwavef2)
1732 [ - + ]: 227706 : ABI_FREE(ghc2)
1733 : :
1734 : : end if ! end spinor 2
1735 : :
1736 [ - + ]: 227706 : ABI_FREE(cwavef1)
1737 [ - + ]: 227706 : ABI_FREE(cwavef2)
1738 [ - + ]: 227706 : ABI_FREE(kgkpk)
1739 : :
1740 : : end if ! nspinortot
1741 : :
1742 [ - + ]: 813630 : ABI_FREE(work)
1743 [ - + ]: 813630 : if(usezora) then
1744 [ # # ]: 0 : ABI_FREE(zk)
1745 : : end if
1746 : :
1747 [ - + ]: 813630 : end subroutine getghc_nucdip
1748 : : !!***
1749 : :
1750 : : !!****f* ABINIT/getghc_mGGA
1751 : : !!
1752 : : !! NAME
1753 : : !! getghc_mGGA
1754 : : !!
1755 : : !! FUNCTION
1756 : : !! Compute metaGGA contribution to <G|H|C> for input vector |C> expressed in reciprocal space.
1757 : : !!
1758 : : !! INPUTS
1759 : : !! cwavef(2,npw_k*my_nspinor*ndat)=planewave coefficients of wavefunction.
1760 : : !! gbound_k(2*mgfft+4)=sphere boundary info
1761 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1)
1762 : : !! istwf_k=input parameter that describes the storage of wfs
1763 : : !! kg_k(3,npw_k)=G vec coordinates wrt recip lattice transl.
1764 : : !! kpt(3)=current k point
1765 : : !! mgfft=maximum single fft dimension
1766 : : !! mpi_enreg=information about MPI parallelization
1767 : : !! my_nspinor=number of spinorial components of the wavefunctions (on current proc)
1768 : : !! ndat=number of FFTs to perform in parall
1769 : : !! ngfft(18)=contain all needed information about 3D FFT
1770 : : !! npw_k=number of planewaves in basis for given k point.
1771 : : !! nvloc=number of spin components of vxctaulocal
1772 : : !! n4,n5,n6=for dimensionning of vxctaulocal
1773 : : !! gpu_option= GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
1774 : : !! vxctaulocal(n4,n5,n6,nvloc,4)= local potential corresponding to the derivative of XC energy with respect to
1775 : : !! kinetic energy density, in real space, on the augmented fft grid.
1776 : : !! This array contains also the gradient of vxctaulocal (gvxctaulocal) in vxctaulocal(:,:,:,:,2:4).
1777 : : !!
1778 : : !! OUTPUT
1779 : : !! ghc_mGGA(2,npw_k*my_nspinor*ndat)=metaGGA contribution to <G|H|C>
1780 : : !!
1781 : : !! SIDE EFFECTS
1782 : : !!
1783 : : !! SOURCE
1784 : :
1785 : 98528 : subroutine getghc_mGGA(cwavef,ghc_mGGA,gbound_k,gprimd,istwf_k,kg_k,kpt,mgfft,mpi_enreg,&
1786 : 98528 : & ndat,ngfft,npw_k,nvloc,n4,n5,n6,my_nspinor,vxctaulocal,gpu_option)
1787 : :
1788 : : !Arguments ------------------------------------
1789 : : !scalars
1790 : : integer,intent(in) :: istwf_k,mgfft,my_nspinor,ndat,npw_k,nvloc,n4,n5,n6,gpu_option
1791 : : type(MPI_type),intent(in) :: mpi_enreg
1792 : : !arrays
1793 : : integer,intent(in) :: gbound_k(2*mgfft+4),kg_k(3,npw_k),ngfft(18)
1794 : : real(dp),intent(in) :: gprimd(3,3),kpt(3)
1795 : : real(dp),intent(inout) :: cwavef(2,npw_k*my_nspinor*ndat)
1796 : : real(dp),intent(inout) :: ghc_mGGA(2,npw_k*my_nspinor*ndat)
1797 : : real(dp),intent(inout) :: vxctaulocal(n4,n5,n6,nvloc,4)
1798 : :
1799 : : !Local variables-------------------------------
1800 : : !scalars
1801 : : integer,parameter :: tim_fourwf=1
1802 : : integer :: idat,idir,ipw,nspinortot,shift
1803 : : logical :: nspinor1TreatedByThisProc,nspinor2TreatedByThisProc
1804 : : real(dp) :: weight=one
1805 : : !arrays
1806 : : real(dp) :: kg_k_cart_vec(3)
1807 : 98528 : real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:)
1808 : 98528 : real(dp),allocatable :: gcwavef(:,:,:),gcwavef1(:,:,:),gcwavef2(:,:,:)
1809 : 98528 : real(dp),allocatable :: ghc1(:,:),ghc2(:,:)
1810 : 98528 : real(dp),allocatable :: lcwavef(:,:),lcwavef1(:,:),lcwavef2(:,:)
1811 : 98528 : real(dp),allocatable :: work(:,:,:,:)
1812 : : ! *********************************************************************
1813 : :
1814 [ + + + + ]: 176135072 : ghc_mGGA(:,:)=zero
1815 : :
1816 [ + - ]: 98528 : if (nvloc/=1) return
1817 : :
1818 : 98528 : nspinortot=min(2,(1+mpi_enreg%paral_spinor)*my_nspinor)
1819 [ + - ]: 98528 : if (mpi_enreg%paral_spinor==0) then
1820 : 98528 : shift=npw_k
1821 : 98528 : nspinor1TreatedByThisProc=.true.
1822 : 98528 : nspinor2TreatedByThisProc=(nspinortot==2)
1823 : : else
1824 : 0 : shift=0
1825 : 0 : nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
1826 : 0 : nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
1827 : : end if
1828 : :
1829 [ - + + - : 492640 : ABI_MALLOC(work,(2,n4,n5,n6*ndat))
- + + - -
+ + - + -
+ - - + -
+ ]
1830 : :
1831 [ + + ]: 98528 : if (nspinortot==1) then
1832 : :
1833 [ - + + - : 136392 : ABI_MALLOC(ghc1,(2,npw_k*ndat))
+ - + - -
+ - + ]
1834 : :
1835 : : ! Do it in 3 STEPs:
1836 : : ! STEP1: Compute grad of cwavef and Laplacian of cwavef
1837 [ - + + - : 181856 : ABI_MALLOC(gcwavef,(2,npw_k*ndat,3))
+ - + - +
- - + -
+ ]
1838 [ - + + - : 90928 : ABI_MALLOC(lcwavef,(2,npw_k*ndat))
- + - + ]
1839 : : !!$OMP PARALLEL DO
1840 [ + + + + : 387384872 : gcwavef = zero; lcwavef = zero
+ + + + +
+ ]
1841 [ + + ]: 95408 : do idat=1,ndat
1842 [ + + ]: 32362326 : do ipw=1,npw_k
1843 : : ! convert k + G from reduced coords to Cartesian
1844 [ + + + + : 613071442 : kg_k_cart_vec = two_pi*MATMUL(gprimd,kpt(1:3)+kg_k(1:3,ipw))
+ + + + ]
1845 : : ! form \grad\psi = i(k + G) \psi in Cartesian frame
1846 [ + + ]: 129067672 : gcwavef(1,ipw+(idat-1)*npw_k,1:3)= cwavef(2,ipw+(idat-1)*npw_k)*kg_k_cart_vec(1:3)
1847 [ + + ]: 129067672 : gcwavef(2,ipw+(idat-1)*npw_k,1:3)= -cwavef(1,ipw+(idat-1)*npw_k)*kg_k_cart_vec(1:3)
1848 : : ! form \nabla^2\psi = -|k + G|^2 \psi
1849 : : lcwavef(1:2,ipw+(idat-1)*npw_k)=&
1850 [ + + + + ]: 193651452 : &lcwavef(1:2,ipw+(idat-1)*npw_k)-cwavef(1:2,ipw+(idat-1)*npw_k)*DOT_PRODUCT(kg_k_cart_vec,kg_k_cart_vec)
1851 : : end do
1852 : : end do
1853 : : ! STEP2: Compute (vxctaulocal)*(Laplacian of cwavef) and add it to ghc
1854 : : call fourwf(1,vxctaulocal(:,:,:,:,1),lcwavef,ghc1,work,gbound_k,gbound_k,&
1855 : : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1856 : 45464 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1857 : : !!$OMP PARALLEL DO
1858 [ + + ]: 95408 : do idat=1,ndat
1859 [ + + ]: 32362326 : do ipw=1,npw_k
1860 [ + + ]: 96850698 : ghc_mGGA(:,ipw+(idat-1)*npw_k)=ghc_mGGA(:,ipw+(idat-1)*npw_k)-half*ghc1(:,ipw+(idat-1)*npw_k)
1861 : : end do
1862 : : end do
1863 [ - + ]: 45464 : ABI_FREE(lcwavef)
1864 : : ! STEP3: Compute sum of (grad components of vxctaulocal)*(grad components of cwavef)
1865 : : ! note: since grad cwavef is in Cart frame, evidently grad vxc is also
1866 [ + + ]: 181856 : do idir=1,3
1867 : : call fourwf(1,vxctaulocal(:,:,:,:,1+idir),gcwavef(:,:,idir),ghc1,work,gbound_k,gbound_k,&
1868 : : istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1869 : 136392 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1870 : : !!$OMP PARALLEL DO
1871 [ + + ]: 331688 : do idat=1,ndat
1872 [ + + ]: 97086978 : do ipw=1,npw_k
1873 [ + + ]: 290552094 : ghc_mGGA(:,ipw+(idat-1)*npw_k)=ghc_mGGA(:,ipw+(idat-1)*npw_k)-half*ghc1(:,ipw+(idat-1)*npw_k)
1874 : : end do
1875 : : end do
1876 : : end do ! idir
1877 [ - + ]: 45464 : ABI_FREE(gcwavef)
1878 [ - + ]: 45464 : ABI_FREE(ghc1)
1879 : :
1880 : : else ! nspinortot==2
1881 : :
1882 [ - + + - : 159192 : ABI_MALLOC(cwavef1,(2,npw_k*ndat))
+ - + - -
+ - + ]
1883 [ - + + - : 106128 : ABI_MALLOC(cwavef2,(2,npw_k*ndat))
- + - + ]
1884 [ + + ]: 109248 : do idat=1,ndat
1885 [ + + ]: 13315213 : do ipw=1,npw_k
1886 [ + + ]: 39617895 : cwavef1(1:2,ipw+(idat-1)*npw_k)=cwavef(1:2,ipw+(idat-1)*my_nspinor*npw_k)
1887 [ + + ]: 39674079 : cwavef2(1:2,ipw+(idat-1)*npw_k)=cwavef(1:2,ipw+(idat-1)*my_nspinor*npw_k+shift)
1888 : : end do
1889 : : end do
1890 : : ! call cg_zcopy(npw*ndat,cwavef(1,1),cwavef1)
1891 : : ! call cg_zcopy(npw*ndat,cwavef(1,1+shift),cwavef2)
1892 : :
1893 : :
1894 [ + - ]: 53064 : if (nspinor1TreatedByThisProc) then
1895 : :
1896 [ - + + - : 106128 : ABI_MALLOC(ghc1,(2,npw_k*ndat))
- + - + ]
1897 : :
1898 : : ! Do it in 3 STEPs:
1899 : : ! STEP1: Compute grad of cwavef and Laplacian of cwavef
1900 [ - + + - : 212256 : ABI_MALLOC(gcwavef1,(2,npw_k*ndat,3))
+ - + - +
- - + -
+ ]
1901 [ - + + - : 106128 : ABI_MALLOC(lcwavef1,(2,npw_k*ndat))
- + - + ]
1902 [ + + + + : 158683836 : gcwavef1 = zero; lcwavef1 = zero
+ + + + +
+ ]
1903 : : !!$OMP PARALLEL DO
1904 [ + + ]: 109248 : do idat=1,ndat
1905 [ + + ]: 13315213 : do ipw=1,npw_k
1906 : : ! convert k + G from reduced coords to Cartesian
1907 [ + + + + : 250913335 : kg_k_cart_vec = two_pi*MATMUL(gprimd,kpt(1:3)+kg_k(1:3,ipw))
+ + + + ]
1908 : : ! form \grad\psi = i(k + G) \psi in Cartesian frame
1909 [ + + ]: 52823860 : gcwavef1(1,ipw+(idat-1)*npw_k,1:3)= cwavef1(2,ipw+(idat-1)*npw_k)*kg_k_cart_vec(1:3)
1910 [ + + ]: 52823860 : gcwavef1(2,ipw+(idat-1)*npw_k,1:3)= -cwavef1(1,ipw+(idat-1)*npw_k)*kg_k_cart_vec(1:3)
1911 : : ! form \nabla^2\psi = -|k + G|^2 \psi
1912 : : lcwavef1(1:2,ipw+(idat-1)*npw_k)=&
1913 [ + + + + ]: 79291974 : &lcwavef1(1:2,ipw+(idat-1)*npw_k)-cwavef1(1:2,ipw+(idat-1)*npw_k)*DOT_PRODUCT(kg_k_cart_vec,kg_k_cart_vec)
1914 : : end do
1915 : : end do
1916 : : ! STEP2: Compute (vxctaulocal)*(Laplacian of cwavef) and add it to ghc
1917 : : call fourwf(1,vxctaulocal(:,:,:,:,1),lcwavef1,ghc1,work,gbound_k,gbound_k,&
1918 : : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1919 : 53064 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1920 : : !!$OMP PARALLEL DO
1921 [ + + ]: 109248 : do idat=1,ndat
1922 [ + + ]: 13315213 : do ipw=1,npw_k
1923 [ + + ]: 39674079 : ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)=ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)-half*ghc1(:,ipw+(idat-1)*npw_k)
1924 : : end do
1925 : : end do
1926 [ - + ]: 53064 : ABI_FREE(lcwavef1)
1927 : : ! STEP3: Compute (grad components of vxctaulocal)*(grad components of cwavef)
1928 [ + + ]: 212256 : do idir=1,3
1929 : : call fourwf(1,vxctaulocal(:,:,:,:,1+idir),gcwavef1(:,:,idir),ghc1,work,gbound_k,gbound_k,&
1930 : : istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1931 : 159192 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1932 : : !!$OMP PARALLEL DO
1933 [ + + ]: 380808 : do idat=1,ndat
1934 [ + + ]: 39945639 : do ipw=1,npw_k
1935 [ + + ]: 119022237 : ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k) = ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)-half*ghc1(:,ipw+(idat-1)*npw_k)
1936 : : end do
1937 : : end do
1938 : : end do ! idir
1939 [ - + ]: 53064 : ABI_FREE(gcwavef1)
1940 [ - + ]: 53064 : ABI_FREE(ghc1)
1941 : :
1942 : : end if ! spin 1 treated by this proc
1943 : :
1944 [ + - ]: 53064 : if (nspinor2TreatedByThisProc) then
1945 : :
1946 [ - + + - : 106128 : ABI_MALLOC(ghc2,(2,npw_k*ndat))
- + - + ]
1947 : :
1948 : : ! Do it in 3 STEPs:
1949 : : ! STEP1: Compute grad of cwavef and Laplacian of cwavef
1950 [ - + + - : 212256 : ABI_MALLOC(gcwavef2,(2,npw_k*ndat,3))
+ - + - +
- - + -
+ ]
1951 [ - + + - : 106128 : ABI_MALLOC(lcwavef2,(2,npw_k*ndat))
- + - + ]
1952 : : !!$OMP PARALLEL DO
1953 [ + + + + : 158683836 : gcwavef2 = zero; lcwavef2 = zero
+ + + + +
+ ]
1954 [ + + ]: 109248 : do idat=1,ndat
1955 [ + + ]: 13315213 : do ipw=1,npw_k
1956 : : ! convert k + G from reduced coords to Cartesian
1957 [ + + + + : 250913335 : kg_k_cart_vec = two_pi*MATMUL(gprimd,kpt(1:3)+kg_k(1:3,ipw))
+ + + + ]
1958 : : ! form \grad\psi = i(k + G) \psi in Cartesian frame
1959 [ + + ]: 52823860 : gcwavef2(1,ipw+(idat-1)*npw_k,1:3)= cwavef2(2,ipw+(idat-1)*npw_k)*kg_k_cart_vec(1:3)
1960 [ + + ]: 52823860 : gcwavef2(2,ipw+(idat-1)*npw_k,1:3)= -cwavef2(1,ipw+(idat-1)*npw_k)*kg_k_cart_vec(1:3)
1961 : : ! form \nabla^2\psi = -|k + G|^2 \psi
1962 : : lcwavef2(1:2,ipw+(idat-1)*npw_k)=&
1963 [ + + + + ]: 79291974 : &lcwavef2(1:2,ipw+(idat-1)*npw_k)-cwavef2(1:2,ipw+(idat-1)*npw_k)*DOT_PRODUCT(kg_k_cart_vec,kg_k_cart_vec)
1964 : : end do
1965 : : end do
1966 : : ! STEP2: Compute (vxctaulocal)*(Laplacian of cwavef) and add it to ghc
1967 : : call fourwf(1,vxctaulocal(:,:,:,:,1),lcwavef2,ghc2,work,gbound_k,gbound_k,&
1968 : : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1969 : 53064 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1970 : : !!$OMP PARALLEL DO
1971 [ + + ]: 109248 : do idat=1,ndat
1972 [ + + ]: 13315213 : do ipw=1,npw_k
1973 : : ! original code
1974 : : ! ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)=ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)-half*ghc2(:,ipw+(idat-1)*npw_k)
1975 : : ! but this stores the spinor2 result in the spinor1 location. Should be stored with shift
1976 : : ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k+shift)=ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k+shift)&
1977 [ + + ]: 39674079 : & -half*ghc2(:,ipw+(idat-1)*npw_k)
1978 : : end do
1979 : : end do
1980 [ - + ]: 53064 : ABI_FREE(lcwavef2)
1981 : : ! STEP3: Compute sum of (grad components of vxctaulocal)*(grad components of cwavef)
1982 [ + + ]: 212256 : do idir=1,3
1983 : : call fourwf(1,vxctaulocal(:,:,:,:,1+idir),gcwavef2(:,:,idir),ghc2,work,gbound_k,gbound_k,&
1984 : : istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,&
1985 : 159192 : & tim_fourwf,weight,weight,gpu_option=gpu_option)
1986 : : !!$OMP PARALLEL DO
1987 [ + + ]: 380808 : do idat=1,ndat
1988 [ + + ]: 39945639 : do ipw=1,npw_k
1989 : : ! original code
1990 : : ! ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)=ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k)-half*ghc2(:,ipw+(idat-1)*npw_k)
1991 : : ! but this stores the spinor2 result in the spinor1 location. Should be stored with shift
1992 : : ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k+shift)=ghc_mGGA(:,ipw+(idat-1)*my_nspinor*npw_k+shift)&
1993 [ + + ]: 119022237 : & -half*ghc2(:,ipw+(idat-1)*npw_k)
1994 : : end do
1995 : : end do
1996 : : end do ! idir
1997 : :
1998 [ - + ]: 53064 : ABI_FREE(gcwavef2)
1999 [ - + ]: 53064 : ABI_FREE(ghc2)
2000 : :
2001 : : end if ! spin 2 treated by this proc
2002 : :
2003 [ - + ]: 53064 : ABI_FREE(cwavef1)
2004 [ - + ]: 53064 : ABI_FREE(cwavef2)
2005 : :
2006 : : end if ! nspinortot
2007 : :
2008 [ - + ]: 98528 : ABI_FREE(work)
2009 : :
2010 : : end subroutine getghc_mGGA
2011 : : !!***
2012 : :
2013 : : !!****f* ABINIT/cwavef_double_rfft_trick_pack
2014 : : !!
2015 : : !! NAME
2016 : : !! cwavef_double_rfft_trick_pack
2017 : : !!
2018 : : !! FUNCTION
2019 : : !!
2020 : : !! We have C(G)=C(-G)^* and D(G)=D(-G)^* and only G components are in memory (istwfk=2)
2021 : : !! The Fourier transform is:
2022 : : !! C(r) = sum_G e^(iGr) C(G) = sum_(G_z>=0,G/=0) 2 Re[e^(iGr) C(G)] + C(0)
2023 : : !! so C(r) is a real function (same for D)
2024 : : !! Here we construct:
2025 : : !! E( G) = C(G) + i D(G)
2026 : : !! E(-G) = C(G)^* + i D(G)^* (G/=0)
2027 : : !! so:
2028 : : !! E(r) = C(r) + i D(r)
2029 : : !! In short, one can do only one FFT on E (with istwfk=1) and obtains the FFT of C and D (istwfk=2)
2030 : : !!
2031 : : !! INPUTS
2032 : : !! cwavef(2,npw_k*my_nspinor*(2*ndat))=planewave coefficients of wavefunctioni (istwfk=2)
2033 : : !! mpi_enreg=information about MPI parallelization
2034 : : !! ndat=number of FFTs to perform in parall
2035 : : !! npw_k=number of planewaves in basis for given k point (istwfk=2).
2036 : : !!
2037 : : !! OUTPUT
2038 : : !! cwavef_fft(2,npw_fft*my_nspinor*ndat)=planewave coefficients of wavefunction (istwfk=1)
2039 : : !!
2040 : : !! SOURCE
2041 : : !!
2042 [ - + - + ]: 16912 : subroutine cwavef_double_rfft_trick_pack(cwavef,cwavef_fft,me_g0,ndat,npw_k)
2043 : :
2044 : : !Arguments ------------------------------------
2045 : : !scalars
2046 : : integer,intent(in) :: ndat,npw_k,me_g0
2047 : : !arrays
2048 : : real(dp),intent(in) :: cwavef(:,:)
2049 : : real(dp),intent(out) :: cwavef_fft(:,:)
2050 : : !Local variables-------------------------------
2051 : : !scalars
2052 : : integer :: idat,i0,ib1,ib2,npw_fft,ndat_
2053 : : logical :: ndat_is_odd
2054 : :
2055 : 16912 : ndat_ = ndat / 2
2056 : 16912 : ndat_is_odd=.False.
2057 [ - + ]: 16912 : if (modulo(ndat,2)/=0) then
2058 : 0 : ndat_=ndat_+1
2059 : 0 : ndat_is_odd=.True.
2060 : : end if
2061 : :
2062 : 16912 : npw_fft=2*npw_k
2063 : 16912 : i0=1
2064 [ + - ]: 16912 : if (me_g0==1) then! Do not include G=(0,0,0) twice
2065 : 16912 : npw_fft=npw_fft-1
2066 : 16912 : i0=2
2067 : : end if
2068 : :
2069 [ - + ]: 16912 : if (size(cwavef,1)/=2) then
2070 : 0 : ABI_BUG('wrong size for cwavef (dim 1)')
2071 : : end if
2072 [ - + ]: 16912 : if (ndat_is_odd) then
2073 [ # # ]: 0 : if (size(cwavef,2)/=npw_k*(2*ndat_-1)) then
2074 : 0 : ABI_BUG('wrong size for cwavef (dim 2) (odd)')
2075 : : end if
2076 : : else
2077 [ - + ]: 16912 : if (size(cwavef,2)/=npw_k*2*ndat_) then
2078 : 0 : ABI_BUG('wrong size for cwavef (dim 2)')
2079 : : end if
2080 : : end if
2081 [ + - - + ]: 16912 : if (size(cwavef_fft,1)/=2.or.size(cwavef_fft,2)/=npw_fft*ndat_) then
2082 : 0 : ABI_BUG('wrong size for cwavef_fft')
2083 : : end if
2084 : :
2085 [ + + ]: 106164 : do idat=1,ndat_
2086 : 89252 : ib1=(idat-1)*npw_fft ! band shift for cwavef_fft
2087 : 89252 : ib2=(idat-1)*2*npw_k ! band shift for cwavef
2088 [ + - ]: 106164 : if (.not.ndat_is_odd.or.idat<ndat_) then
2089 : : ! E(G) = C(G) + i D(G)
2090 : : cwavef_fft(1,1+ib1:npw_k+ib1) = cwavef(1,1+ib2 : npw_k+ib2) &
2091 [ + + ]: 68960932 : -cwavef(2,1+npw_k+ib2:2*npw_k+ib2)
2092 : : cwavef_fft(2,1+ib1:npw_k+ib1) = cwavef(2,1+ib2 : npw_k+ib2) &
2093 [ + + ]: 68960932 : +cwavef(1,1+npw_k+ib2:2*npw_k+ib2)
2094 : : ! E(-G) = C(G)^* + i D(G)^* (G/=0)
2095 : : cwavef_fft(1,1+npw_k+ib1:npw_fft+ib1) = cwavef(1,i0 +ib2: npw_k+ib2) &
2096 [ + + ]: 68871680 : +cwavef(2,i0+npw_k+ib2:2*npw_k+ib2)
2097 : : cwavef_fft(2,1+npw_k+ib1:npw_fft+ib1) =-cwavef(2,i0 +ib2: npw_k+ib2) &
2098 [ + + ]: 68871680 : +cwavef(1,i0+npw_k+ib2:2*npw_k+ib2)
2099 : : else ! idat=ndat_ and ndat_is_odd : the vector D does not exist
2100 : : ! E(G) = C(G)
2101 [ # # ]: 0 : cwavef_fft(1,1+ib1:npw_k+ib1) = cwavef(1,1+ib2:npw_k+ib2)
2102 [ # # ]: 0 : cwavef_fft(2,1+ib1:npw_k+ib1) = cwavef(2,1+ib2:npw_k+ib2)
2103 : : ! E(-G) = C(G)^* (G/=0)
2104 [ # # ]: 0 : cwavef_fft(1,1+npw_k+ib1:npw_fft+ib1) = cwavef(1,i0+ib2:npw_k+ib2)
2105 [ # # ]: 0 : cwavef_fft(2,1+npw_k+ib1:npw_fft+ib1) =-cwavef(2,i0+ib2:npw_k+ib2)
2106 : : end if
2107 : : end do
2108 : :
2109 : 16912 : end subroutine cwavef_double_rfft_trick_pack
2110 : : !!***
2111 : :
2112 : : !!****f* ABINIT/cwavef_double_rfft_trick_unpack
2113 : : !!
2114 : : !! NAME
2115 : : !! cwavef_double_rfft_trick_unpack
2116 : : !!
2117 : : !! FUNCTION
2118 : : !!
2119 : : !! From the "cwavef_double_rfft_trick_pack" routine we have:
2120 : : !! E( G) = C(G) + i D(G)
2121 : : !! E(-G) = C(G)^* + i D(G)^* (G/=0)
2122 : : !! Here we compute:
2123 : : !! C(G) = ( E(G) + E(-G)^* ) / 2 (G/=0)
2124 : : !! D(G) = ( iE(-G)^* - iE(G)) / 2 (G/=0)
2125 : : !! and:
2126 : : !! C(0) = Re(E(0))
2127 : : !! D(0) = Im(E(0))
2128 : : !!
2129 : : !! INPUTS
2130 : : !! cwavef_fft(2,npw_fft*my_nspinor*ndat_)=planewave coefficients of wavefunction (istwfk=1)
2131 : : !! mpi_enreg=information about MPI parallelization
2132 : : !! ndat=number of FFTs to perform in parall
2133 : : !! npw_k=number of planewaves in basis for given k point (istwfk=2).
2134 : : !!
2135 : : !! OUTPUT
2136 : : !! cwavef(2,npw_fft*my_nspinor*(2*ndat))=planewave coefficients of wavefunction (istwfk=2)
2137 : : !!
2138 : : !! SOURCE
2139 : : !!
2140 [ - + - + ]: 16912 : subroutine cwavef_double_rfft_trick_unpack(cwavef,cwavef_fft,me_g0,ndat,npw_k)
2141 : :
2142 : : !Arguments ------------------------------------
2143 : : !scalars
2144 : : integer,intent(in) :: ndat,npw_k,me_g0
2145 : : !arrays
2146 : : real(dp),intent(out) :: cwavef(:,:)
2147 : : real(dp),intent(in) :: cwavef_fft(:,:)
2148 : : !Local variables-------------------------------
2149 : : !scalars
2150 : : integer :: idat,i0,ib1,ib2,npw_fft,ndat_
2151 : : logical :: ndat_is_odd
2152 : :
2153 : 16912 : ndat_ = ndat / 2
2154 : 16912 : ndat_is_odd=.False.
2155 [ - + ]: 16912 : if (modulo(ndat,2)/=0) then
2156 : 0 : ndat_=ndat_+1
2157 : 0 : ndat_is_odd=.True.
2158 : : end if
2159 : :
2160 : 16912 : npw_fft=2*npw_k
2161 : 16912 : i0=1
2162 [ + - ]: 16912 : if (me_g0==1) then! Do not include G=(0,0,0) twice
2163 : 16912 : npw_fft=npw_fft-1
2164 : 16912 : i0=2
2165 : : end if
2166 : :
2167 [ - + ]: 16912 : if (size(cwavef,1)/=2) then
2168 : 0 : ABI_BUG('wrong size for cwavef (dim 1)')
2169 : : end if
2170 [ - + ]: 16912 : if (ndat_is_odd) then
2171 [ # # ]: 0 : if (size(cwavef,2)/=npw_k*(2*ndat_-1)) then
2172 : 0 : ABI_BUG('wrong size for cwavef (dim 2) (odd)')
2173 : : end if
2174 : : else
2175 [ - + ]: 16912 : if (size(cwavef,2)/=npw_k*2*ndat_) then
2176 : 0 : ABI_BUG('wrong size for cwavef (dim 2)')
2177 : : end if
2178 : : end if
2179 [ + - - + ]: 16912 : if (size(cwavef_fft,1)/=2.or.size(cwavef_fft,2)/=npw_fft*ndat_) then
2180 : 0 : ABI_BUG('wrong size for cwavef_fft')
2181 : : end if
2182 : :
2183 [ + + ]: 106164 : do idat=1,ndat_
2184 : 89252 : ib1=(idat-1)*npw_fft ! band shift for cwavef_fft
2185 : 89252 : ib2=(idat-1)*2*npw_k ! band shift for cwavef
2186 : : ! C(G) = ( E(G) + E(-G)^* ) / 2 (factor 1/2 will be applied later)
2187 : : cwavef(1,i0+ib2:npw_k+ib2) = cwavef_fft(1,i0 +ib1:npw_k +ib1) & !+Re(E( G))
2188 [ + + ]: 68871680 : +cwavef_fft(1,1 +npw_k+ib1:npw_fft+ib1) !+Re(E(-G))
2189 : : cwavef(2,i0+ib2:npw_k+ib2) = cwavef_fft(2,i0 +ib1:npw_k +ib1) & !+Im(E( G))
2190 [ + + ]: 68871680 : -cwavef_fft(2,1 +npw_k+ib1:npw_fft+ib1) !-Im(E(-G))
2191 [ + - ]: 89252 : if (.not.ndat_is_odd.or.idat<ndat_) then
2192 : : ! D(G) = ( iE(-G)^* - iE(G) ) / 2 (factor 1/2 will be applied later)
2193 : : cwavef(1,i0+npw_k+ib2:2*npw_k+ib2) = cwavef_fft(2,i0 +ib1:npw_k +ib1) & !+Im(E( G))
2194 [ + + ]: 68871680 : +cwavef_fft(2,1 +npw_k+ib1:npw_fft+ib1) !+Im(E(-G))
2195 : : cwavef(2,i0+npw_k+ib2:2*npw_k+ib2) =-cwavef_fft(1,i0 +ib1:npw_k +ib1) & !-Re(E( G))
2196 [ + + ]: 68871680 : +cwavef_fft(1,1 +npw_k+ib1:npw_fft+ib1) !+Re(E(-G))
2197 : : end if
2198 [ + - ]: 106164 : if (me_g0==1) then
2199 : : ! Compute C(G=0) and D(G=0) and multiply by 2 as we apply 1/2 to the whole array shortly afterwards
2200 : : ! C(G=0) = Re(E(G=0))
2201 : 89252 : cwavef(1,1+ib2) = two*cwavef_fft(1,1+ib1)
2202 : 89252 : cwavef(2,1+ib2) = zero
2203 [ + - ]: 89252 : if (.not.ndat_is_odd.or.idat<ndat_) then
2204 : : ! D(G=0) = Im(E(G=0))
2205 : 89252 : cwavef(1,1+npw_k+ib2) = two*cwavef_fft(2,1+ib1)
2206 : 89252 : cwavef(2,1+npw_k+ib2) = zero
2207 : : end if
2208 : : end if
2209 : :
2210 : : end do
2211 : :
2212 [ + + + + ]: 413246992 : cwavef(:,:) = half*cwavef(:,:)
2213 : :
2214 : 16912 : end subroutine cwavef_double_rfft_trick_unpack
2215 : : !!***
2216 : :
2217 : : !!****f* ABINIT/getgsc
2218 : : !! NAME
2219 : : !! getgsc
2220 : : !!
2221 : : !! FUNCTION
2222 : : !! Compute <G|S|C> for all input vectors |Cnk> at a given k-point,
2223 : : !! OR for one input vector |Cnk>.
2224 : : !! |Cnk> are expressed in reciprocal space.
2225 : : !! S is the overlap operator between |Cnk> (used for PAW).
2226 : : !!
2227 : : !! INPUTS
2228 : : !! cg(2,mcg)=planewave coefficients of wavefunctions
2229 : : !! cprj(natom,mcprj)= wave functions projected with non-local projectors: cprj=<p_i|Cnk>
2230 : : !! gs_ham <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
2231 : : !! ibg=shift to be applied on the location of data in the array cprj (beginning of current k-point)
2232 : : !! icg=shift to be applied on the location of data in the array cg (beginning of current k-point)
2233 : : !! igsc=shift to be applied on the location of data in the array gsc (beginning of current k-point)
2234 : : !! ikpt,isppol=indexes of current (spin.kpoint)
2235 : : !! mcg=second dimension of the cg array
2236 : : !! mcprj=second dimension of the cprj array
2237 : : !! mgsc=second dimension of the gsc array
2238 : : !! mpi_enreg=information about MPI parallelization
2239 : : !! ndat=number of bands to compute in parallel
2240 : : !! natom=number of atoms in unit cell.
2241 : : !! nband= if positive: number of bands at this k point for that spin polarization
2242 : : !! if negative: abs(nband) is the index of the only band to be computed
2243 : : !! npw_k=number of planewaves in basis for given k point.
2244 : : !! nspinor=number of spinorial components of the wavefunctions
2245 : : !! [select_k]=optional, option governing the choice of k points to be used.
2246 : : !! gs_ham datastructure contains quantities needed to apply overlap operator
2247 : : !! in reciprocal space between 2 kpoints, k and k^prime (equal in most cases);
2248 : : !! if select_k=1, <k^prime|S|k> is applied [default]
2249 : : !! if select_k=2, <k|S|k^prime> is applied
2250 : : !! if select_k=3, <k|S|k> is applied
2251 : : !! if select_k=4, <k^prime|S|k^prime> is applied
2252 : : !!
2253 : : !! OUTPUT
2254 : : !! gsc(2,mgsc)= <g|S|Cnk> or <g|S^(1)|Cnk> (S=overlap)
2255 : : !!
2256 : : !! SOURCE
2257 : :
2258 : 180133 : subroutine getgsc(cg,cprj,gs_ham,gsc,ibg,icg,igsc,ikpt,isppol,&
2259 : : & mcg,mcprj,mgsc,mpi_enreg,ndat,natom,nband,npw_k,nspinor,select_k)
2260 : :
2261 : : !Arguments ------------------------------------
2262 : : !scalars
2263 : : integer,intent(in) :: ibg,icg,igsc,ikpt,isppol,mcg,mcprj
2264 : : integer,intent(in) :: mgsc,natom,nband,npw_k,nspinor,ndat
2265 : : !TODO : may be needed to distribute cprj over band procs
2266 : : ! integer,intent(in) :: mband_mem
2267 : : integer,intent(in),optional :: select_k
2268 : : type(MPI_type),intent(in) :: mpi_enreg
2269 : : type(gs_hamiltonian_type),intent(inout),target :: gs_ham
2270 : : !arrays
2271 : : real(dp),intent(in), target :: cg(2,mcg)
2272 : : real(dp),intent(out), target :: gsc(2,mgsc)
2273 : : type(pawcprj_type),intent(in) :: cprj(natom,mcprj)
2274 : :
2275 : : !Local variables-------------------------------
2276 : : !scalars
2277 : : integer :: choice,cpopt,dimenl1,dimenl2,iband,iband1,iband2,index_cg,index_cprj
2278 : : integer :: index_gsc,me,my_nspinor,my_ndat,paw_opt,select_k_,signs,tim_nonlop,useylm
2279 : : !character(len=500) :: msg
2280 : : !arrays
2281 : 360266 : real(dp) :: enlout_dum(ndat),tsec(2)
2282 : 180133 : real(dp), contiguous, pointer :: cwavef(:,:),scwavef(:,:)
2283 : 180133 : type(pawcprj_type),allocatable :: cwaveprj(:,:)
2284 : : ! *********************************************************************
2285 : :
2286 : : DBG_ENTER("COLL")
2287 : :
2288 : : !Compatibility tests
2289 : 180133 : my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
2290 : 66816 : my_ndat=ndat; if(nband<1) my_ndat=1
2291 [ - + ]: 180133 : if(gs_ham%usepaw==0) then
2292 : 0 : ABI_BUG('Only compatible with PAW (usepaw=1) !')
2293 : : end if
2294 [ + + + - : 180133 : if(nband<0.and.(mcg<npw_k*my_nspinor.or.mgsc<npw_k*my_nspinor.or.mcprj<my_nspinor)) then
+ - - + ]
2295 : 0 : ABI_BUG('Invalid value for mcg, mgsc or mcprj !')
2296 : : end if
2297 : :
2298 : : !Keep track of total time spent in getgsc:
2299 : 180133 : call timab(565,1,tsec)
2300 : :
2301 [ + - ]: 180133 : if(gs_ham%gpu_option==ABI_GPU_DISABLED) then
2302 [ + + + + ]: 299541301 : gsc = zero
2303 : : else if(gs_ham%gpu_option==ABI_GPU_DISABLED) then
2304 : : #ifdef HAVE_OPENMP_OFFLOAD
2305 : : call gpu_set_to_zero(gsc,int(2,c_size_t)*mgsc)
2306 : : #endif
2307 : : end if
2308 : :
2309 : : !Prepare some data
2310 [ + + ]: 180133 : if (gs_ham%usecprj==1) then
2311 [ - + + - : 808762 : ABI_MALLOC(cwaveprj,(natom,my_nspinor*my_ndat))
- + + - +
- + - - +
- + - + +
+ + + ]
2312 : 113317 : call pawcprj_alloc(cwaveprj,0,gs_ham%dimcprj)
2313 : : else
2314 [ - + - + ]: 66816 : ABI_MALLOC(cwaveprj,(0,0))
2315 : : end if
2316 : 180133 : dimenl1=gs_ham%dimekb1;dimenl2=natom;tim_nonlop=0
2317 : 180133 : choice=1;signs=2;cpopt=-1+3*gs_ham%usecprj;paw_opt=3;useylm=1
2318 [ + - ]: 180133 : select_k_=1;if (present(select_k)) select_k_=select_k
2319 : 180133 : me=mpi_enreg%me_kpt
2320 : :
2321 : : !Loop over bands
2322 : 180133 : index_cprj=ibg;index_cg=icg;index_gsc=igsc
2323 [ + + ]: 180133 : if (nband>0) then
2324 : : iband1=1;iband2=nband
2325 : : !only do 1 band in case nband < 0 (the |nband|th one)
2326 [ + - ]: 66816 : else if (nband<0) then
2327 : 66816 : iband1=abs(nband);iband2=iband1
2328 : 66816 : index_cprj=index_cprj+(iband1-1)*my_nspinor
2329 : 66816 : index_cg =index_cg +(iband1-1)*npw_k*my_nspinor
2330 : 66816 : index_gsc =index_gsc +(iband1-1)*npw_k*my_nspinor
2331 : : end if
2332 : :
2333 [ - + - - : 863369 : do iband=iband1,iband2,my_ndat
- + + + ]
2334 : :
2335 [ + + + - ]: 683236 : if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me.and.nband>0) then
2336 : : ! No longer needed 28/03/2020 to parallelize memory
2337 : : ! gsc(:,1+index_gsc:npw_k*my_nspinor+index_gsc)=zero
2338 : : ! index_gsc=index_gsc+npw_k*my_nspinor
2339 : : !index_cprj=index_cprj+my_nspinor
2340 : : !index_cg=index_cg+npw_k*my_nspinor
2341 : :
2342 : : cycle
2343 : : end if
2344 : :
2345 : : ! Retrieve WF at (n,k)
2346 : 679852 : cwavef(1:2,1:npw_k*my_nspinor*my_ndat) => cg(:,1+index_cg:npw_k*my_nspinor*my_ndat+index_cg)
2347 : 679852 : scwavef(1:2,1:npw_k*my_nspinor*my_ndat) => gsc(:,1+index_gsc:npw_k*my_nspinor*my_ndat+index_gsc)
2348 [ + + ]: 679852 : if (gs_ham%usecprj==1) then
2349 : 613036 : call pawcprj_copy(cprj(:,1+index_cprj:my_nspinor*my_ndat+index_cprj),cwaveprj)
2350 : : end if
2351 : :
2352 : : ! Compute <g|S|Cnk>
2353 : : call nonlop(choice,cpopt,cwaveprj,enlout_dum,gs_ham,0,(/zero/),mpi_enreg,my_ndat,1,paw_opt,&
2354 : 679852 : & signs,scwavef,tim_nonlop,cwavef,cwavef,select_k=select_k_)
2355 : :
2356 : :
2357 : : ! End of loop over bands
2358 : 679852 : index_cprj=index_cprj+my_nspinor*my_ndat
2359 : 679852 : index_cg=index_cg+npw_k*my_nspinor*my_ndat
2360 : 863369 : index_gsc=index_gsc+npw_k*my_nspinor*my_ndat
2361 : : end do
2362 : :
2363 : : !Memory deallocation
2364 [ + + ]: 180133 : if (gs_ham%usecprj==1) then
2365 : 113317 : call pawcprj_free(cwaveprj)
2366 : : end if
2367 [ + - + + : 419196 : ABI_FREE(cwaveprj)
- + - + -
+ ]
2368 : :
2369 : 180133 : call timab(565,2,tsec)
2370 : :
2371 : : DBG_EXIT("COLL")
2372 : :
2373 [ + + - + : 360266 : end subroutine getgsc
- - - - -
- ]
2374 : : !!***
2375 : :
2376 : : !!****f* ABINIT/multithreaded_getghc
2377 : : !!
2378 : : !! NAME
2379 : : !! multithreaded_getghc
2380 : : !!
2381 : : !! FUNCTION
2382 : : !!
2383 : : !! INPUTS
2384 : : !! cpopt=flag defining the status of cwaveprj%cp(:)=<Proj_i|Cnk> scalars (PAW only)
2385 : : !! (same meaning as in nonlop.F90 routine)
2386 : : !! if cpopt=-1, <p_lmn|in> (and derivatives) are computed here (and not saved)
2387 : : !! if cpopt= 0, <p_lmn|in> are computed here and saved
2388 : : !! if cpopt= 1, <p_lmn|in> and first derivatives are computed here and saved
2389 : : !! if cpopt= 2 <p_lmn|in> are already in memory;
2390 : : !! if cpopt= 3 <p_lmn|in> are already in memory; first derivatives are computed here and saved
2391 : : !! if cpopt= 4 <p_lmn|in> and first derivatives are already in memory;
2392 : : !! cwavef(2,npw*my_nspinor*ndat)=planewave coefficients of wavefunction.
2393 : : !! gs_ham <type(gs_hamiltonian_type)>=all data for the Hamiltonian to be applied
2394 : : !! lambda=factor to be used when computing <G|H-lambda.S|C> - only for sij_opt=-1
2395 : : !! Typically lambda is the eigenvalue (or its guess)
2396 : : !! mpi_enreg=information about MPI parallelization
2397 : : !! ndat=number of FFT to do in parallel
2398 : : !! prtvol=control print volume and debugging output
2399 : : !! sij_opt= -PAW ONLY- if 0, only matrix elements <G|H|C> have to be computed
2400 : : !! (S=overlap) if 1, matrix elements <G|S|C> have to be computed in gsc in addition to ghc
2401 : : !! if -1, matrix elements <G|H-lambda.S|C> have to be computed in ghc (gsc not used)
2402 : : !! tim_getghc=timing code of the calling subroutine(can be set to 0 if not attributed)
2403 : : !! type_calc= option governing which part of Hamitonian is to be applied:
2404 : : !! 0: whole Hamiltonian
2405 : : !! 1: local part only
2406 : : !! 2: non-local+kinetic only (added to the existing Hamiltonian)
2407 : : !! 3: local + kinetic only (added to the existing Hamiltonian)
2408 : : !! ===== Optional inputs =====
2409 : : !! [kg_fft_k(3,:)]=optional, (k+G) vector coordinates to be used for the FFT transformation
2410 : : !! instead of the one contained in gs_ham datastructure.
2411 : : !! Typically used for real WF (in parallel) which are FFT-transformed 2 by 2.
2412 : : !! [kg_fft_kp(3,:)]=optional, (k^prime+G) vector coordinates to be used for the FFT transformation
2413 : : !! [select_k]=optional, option governing the choice of k points to be used.
2414 : : !! gs_ham datastructure contains quantities needed to apply Hamiltonian
2415 : : !! in reciprocal space between 2 kpoints, k and k^prime (equal in most cases);
2416 : : !! if select_k=1, <k^prime|H|k> is applied [default]
2417 : : !! if select_k=2, <k|H|k^prime> is applied
2418 : : !! if select_k=3, <k|H|k> is applied
2419 : : !! if select_k=4, <k^prime|H|k^prime> is applied
2420 : : !!
2421 : : !! OUTPUT
2422 : : !! ghc(2,npw*my_nspinor*ndat)=matrix elements <G|H|C> (if sij_opt>=0)
2423 : : !! or <G|H-lambda.S|C> (if sij_opt=-1)
2424 : : !! gvnlxc(2,npw*my_nspinor*ndat)=matrix elements <G|Vnonlocal|C> (if sij_opt>=0)
2425 : : !! or <G|Vnonlocal-lambda.S|C> (if sij_opt=-1)
2426 : : !! if (sij_opt=1)
2427 : : !! gsc(2,npw*my_nspinor*ndat)=matrix elements <G|S|C> (S=overlap).
2428 : : !!
2429 : : !! SIDE EFFECTS
2430 : : !! cwaveprj(natom,my_nspinor*(1+cpopt)*ndat)= wave function projected on nl projectors (PAW only)
2431 : : !!
2432 : : !! SOURCE
2433 : :
2434 [ - + - + : 2303368 : subroutine multithreaded_getghc(cpopt,cwavef,cwaveprj,ghc,gsc,gs_ham,gvnlxc,lambda,mpi_enreg,ndat,&
- + - + -
+ ]
2435 : : & prtvol,sij_opt,tim_getghc,type_calc,&
2436 [ + + + - : 1151684 : & kg_fft_k,kg_fft_kp,select_k,filter_dilatmx_loc) ! optional arguments
- + - + -
- - - ]
2437 : :
2438 : : #ifdef HAVE_OPENMP
2439 : : use omp_lib
2440 : : #endif
2441 : :
2442 : : !Arguments ------------------------------------
2443 : : !scalars
2444 : : logical,intent(in),optional :: filter_dilatmx_loc
2445 : : integer,intent(in) :: cpopt,ndat, prtvol
2446 : : integer,intent(in) :: sij_opt,tim_getghc,type_calc
2447 : : integer,intent(in),optional :: select_k
2448 : : real(dp),intent(in) :: lambda
2449 : : type(MPI_type),intent(in) :: mpi_enreg
2450 : : type(gs_hamiltonian_type),intent(inout),target :: gs_ham
2451 : : !arrays
2452 : : integer,intent(in),optional,target :: kg_fft_k(:,:),kg_fft_kp(:,:)
2453 : : real(dp),intent(out),target :: gsc(:,:)
2454 : : real(dp),intent(inout) :: cwavef(:,:)
2455 : : real(dp),intent(out) :: ghc(:,:),gvnlxc(:,:)
2456 : : type(pawcprj_type),intent(inout),target :: cwaveprj(:,:)
2457 : :
2458 : : !Local variables-------------------------------
2459 : : !scalars
2460 : : logical :: filter_dilatmx_loc_
2461 : : integer :: firstelt, firstprj, lastelt, lastprj,usegvnlxc,usegsc
2462 : : integer :: nthreads,fftalga
2463 : : integer :: ithread
2464 : : integer :: chunk
2465 : : integer :: residuchunk
2466 : : integer :: firstband
2467 : : integer :: lastband
2468 : : integer :: spacedim, spacedim_prj
2469 : : logical :: fftw3_use_lib_threads_sav
2470 : : integer :: select_k_default
2471 : : ! *************************************************************************
2472 : :
2473 [ - + ]: 1151684 : select_k_default = 1; if (present(select_k)) select_k_default = select_k
2474 [ + + ]: 1151684 : filter_dilatmx_loc_ = .true.; if ( present(filter_dilatmx_loc) ) filter_dilatmx_loc_ = filter_dilatmx_loc
2475 : :
2476 : 1151684 : spacedim = size(cwavef ,dim=2)/ndat
2477 : 1151684 : spacedim_prj = size(cwaveprj,dim=2)/ndat
2478 : :
2479 : 1151684 : nthreads = xomp_get_num_threads(open_parallel=.True.)
2480 : 1151684 : fftalga = gs_ham%ngfft(7)/100
2481 [ - + - - ]: 1151684 : if (fftalga==FFT_SG.and.nthreads>1.and.ndat>1) then
2482 : 0 : ABI_ERROR("fftalg=1XX is not thread-safe, so it cannot be used in multi-threaded hamiltonian with nthreads>1 and ndat>1.")
2483 : : end if
2484 : :
2485 : : ! Disabling multithreading for GPU variants (getghc_ompgpu is not thread-safe for now)
2486 : : !$omp parallel default (none) &
2487 : : !$omp& private(ithread,nthreads,chunk,firstband,lastband,residuchunk,firstelt,lastelt), &
2488 : : !$omp& private(firstprj,lastprj,usegvnlxc,usegsc,fftw3_use_lib_threads_sav), &
2489 : : !$omp& shared(cwavef,ghc,gsc, gvnlxc,spacedim,spacedim_prj,ndat,kg_fft_k,kg_fft_kp,gs_ham,cwaveprj,mpi_enreg), &
2490 : : !$omp& shared(gemm_nonlop_use_gemm), &
2491 : : !$omp& firstprivate(cpopt,lambda,prtvol,sij_opt,tim_getghc,type_calc,select_k_default,filter_dilatmx_loc_) &
2492 : : !$omp& IF(gs_ham%gpu_option==ABI_GPU_DISABLED .and. .not. gemm_nonlop_use_gemm)
2493 : 1151684 : ithread = 0
2494 : 1151684 : nthreads = 1
2495 : 1151684 : fftw3_use_lib_threads_sav = .false.
2496 : : if(gs_ham%gpu_option==ABI_GPU_DISABLED .and. .not. gemm_nonlop_use_gemm) then
2497 : : #ifdef HAVE_OPENMP
2498 : : ithread = omp_get_thread_num()
2499 : : nthreads = omp_get_num_threads()
2500 : : !Ensure that libs are used without threads (mkl, openblas, fftw3, ...)
2501 : : #ifdef HAVE_LINALG_MKL_THREADS
2502 : : call mkl_set_num_threads(1)
2503 : : #endif
2504 : : !LB-23/07/24: OpenBLAS detects parallel sections automatically. To comment this line improves performances for some cases.
2505 : : !#ifdef HAVE_LINALG_OPENBLAS_THREADS
2506 : : ! call openblas_set_num_threads(1)
2507 : : !#endif
2508 : : #ifdef HAVE_LINALG_NVPL_THREADS
2509 : : call nvpl_blas_set_num_threads(1)
2510 : : #endif
2511 : : #ifdef HAVE_FFTW3_THREADS
2512 : : fftw3_use_lib_threads_sav=(.not.fftw3_spawn_threads_here(nthreads,nthreads))
2513 : : call fftw3_use_lib_threads(.false.)
2514 : : #endif
2515 : : #endif
2516 : : end if
2517 : 1151684 : chunk = ndat/nthreads ! Divide by 2 to construct chunk of even number of bands
2518 : 1151684 : residuchunk = ndat - nthreads*chunk
2519 : : if ( ithread < nthreads-residuchunk ) then
2520 : : firstband = ithread*chunk+1
2521 : : lastband = (ithread+1)*chunk
2522 : : else
2523 : : firstband = (nthreads-residuchunk)*chunk + ( ithread -(nthreads-residuchunk) )*(chunk+1) +1
2524 : : lastband = firstband+chunk
2525 : : end if
2526 : 1151684 : usegvnlxc=1
2527 [ + + + + ]: 3455052 : if (size(gvnlxc)<=1) usegvnlxc=0
2528 : 1151684 : usegsc=0
2529 [ + + + + ]: 1151684 : if (gs_ham%usepaw==1.and.sij_opt==1) usegsc=1
2530 : :
2531 [ + - ]: 1151684 : if ( lastband /= 0 ) then
2532 : 1151684 : firstelt = (firstband-1)*spacedim+1
2533 : 1151684 : firstprj = (firstband-1)*spacedim_prj+1
2534 : 1151684 : lastelt = lastband*spacedim
2535 : 1151684 : lastprj = lastband*spacedim_prj
2536 : : ! Don't know how to manage optional arguments .... :(
2537 [ + + ]: 1151684 : if ( present(kg_fft_k) ) then
2538 [ - + ]: 3274 : if (present(kg_fft_kp)) then
2539 : : call getghc(cpopt,cwavef(:,firstelt:lastelt),cwaveprj(:,firstprj:lastprj),&
2540 : : & ghc(:,firstelt:lastelt),gsc(:,firstelt:lastelt*usegsc),&
2541 : : & gs_ham,gvnlxc(:,firstelt:lastelt*usegvnlxc),lambda, mpi_enreg,lastband-firstband+1,&
2542 : : & prtvol,sij_opt,tim_getghc,type_calc,&
2543 : 0 : & select_k=select_k_default,kg_fft_k=kg_fft_k,kg_fft_kp=kg_fft_kp,filter_dilatmx_loc=filter_dilatmx_loc_)
2544 : : else
2545 : : call getghc(cpopt,cwavef(:,firstelt:lastelt),cwaveprj(:,firstprj:lastprj),&
2546 : : & ghc(:,firstelt:lastelt),gsc(:,firstelt:lastelt*usegsc),&
2547 : : & gs_ham,gvnlxc(:,firstelt:lastelt*usegvnlxc),lambda, mpi_enreg,lastband-firstband+1,&
2548 : : & prtvol,sij_opt,tim_getghc,type_calc,&
2549 : 3274 : & select_k=select_k_default,kg_fft_k=kg_fft_k,filter_dilatmx_loc=filter_dilatmx_loc_)
2550 : : end if
2551 : : else
2552 [ - + ]: 1148410 : if (present(kg_fft_kp)) then
2553 : : call getghc(cpopt,cwavef(:,firstelt:lastelt),cwaveprj(:,firstprj:lastprj),&
2554 : : & ghc(:,firstelt:lastelt),gsc(:,firstelt:lastelt*usegsc),&
2555 : : & gs_ham,gvnlxc(:,firstelt:lastelt*usegvnlxc),lambda, mpi_enreg,lastband-firstband+1,&
2556 : : & prtvol,sij_opt,tim_getghc,type_calc,&
2557 : 0 : & select_k=select_k_default,kg_fft_kp=kg_fft_kp,filter_dilatmx_loc=filter_dilatmx_loc_)
2558 : : else
2559 : : call getghc(cpopt,cwavef(:,firstelt:lastelt),cwaveprj(:,firstprj:lastprj),&
2560 : : & ghc(:,firstelt:lastelt),gsc(:,firstelt:lastelt*usegsc),&
2561 : : & gs_ham,gvnlxc(:,firstelt:lastelt*usegvnlxc),lambda, mpi_enreg,lastband-firstband+1,&
2562 : : & prtvol,sij_opt,tim_getghc,type_calc,&
2563 : 1148410 : & select_K=select_k_default,filter_dilatmx_loc=filter_dilatmx_loc_)
2564 : : end if
2565 : : end if
2566 : : end if
2567 : : if(gs_ham%gpu_option==ABI_GPU_DISABLED .and. .not. gemm_nonlop_use_gemm) then
2568 : : #ifdef HAVE_OPENMP
2569 : : !Restore libs behavior (mkl, openblas, fftw3, ...)
2570 : : #ifdef HAVE_LINALG_MKL_THREADS
2571 : : call mkl_set_num_threads(nthreads)
2572 : : #endif
2573 : : !LB-23/07/24: OpenBLAS detects parallel sections automatically. To comment this line improves performances for some cases.
2574 : : !#ifdef HAVE_LINALG_OPENBLAS_THREADS
2575 : : ! call openblas_set_num_threads(nthreads)
2576 : : !#endif
2577 : : #ifdef HAVE_LINALG_NVPL_THREADS
2578 : : call nvpl_blas_set_num_threads(nthreads)
2579 : : #endif
2580 : : #ifdef HAVE_FFTW3_THREADS
2581 : : call fftw3_use_lib_threads(fftw3_use_lib_threads_sav)
2582 : : #endif
2583 : : #endif
2584 : : end if
2585 : : !$omp end parallel
2586 : :
2587 : 1151684 : end subroutine multithreaded_getghc
2588 : : !!***
2589 : :
2590 : : end module m_getghc
2591 : : !!***
|
