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