Line data Source code
1 : !!****m* ABINIT/m_fock_getghc
2 : !! NAME
3 : !!
4 : !! FUNCTION
5 : !!
6 : !! COPYRIGHT
7 : !! Copyright (C) 2013-2026 ABINIT group (CMartins, FJ, MT, XG)
8 : !! This file is distributed under the terms of the
9 : !! GNU General Public License, see ~abinit/COPYING
10 : !! or http://www.gnu.org/copyleft/gpl.txt .
11 : !!
12 : !! SOURCE
13 :
14 : #if defined HAVE_CONFIG_H
15 : #include "config.h"
16 : #endif
17 :
18 : #include "abi_common.h"
19 :
20 : module m_fock_getghc
21 :
22 : use, intrinsic :: iso_c_binding, only: c_size_t, c_loc
23 : use defs_basis
24 : use m_abicore
25 : use m_errors
26 : use m_xmpi
27 : use m_fock
28 : use m_pawcprj
29 : !use m_cgtools
30 :
31 : use defs_abitypes, only : mpi_type
32 : use defs_datatypes, only : pseudopotential_type
33 : use m_time, only : timab, time_accu
34 : use m_matrix, only : matr3inv
35 : use m_cgtools, only : dotprod_g,dotprod_g_batch_half,dotprod_g_batch_full
36 : use m_kg, only : mkkpg
37 : use m_fftcore, only : sphereboundary
38 : use m_fft, only : fftpac, fourwf, fourdp
39 : use m_fstrings, only : sjoin, itoa
40 : use m_hamiltonian, only : gs_hamiltonian_type, K_H_KPRIME
41 : use m_paw_nhat, only : pawmknhat_psipsi_ndat, pawdijhat_ndat
42 : use m_spacepar, only : hartre
43 : use m_nonlop, only : nonlop
44 : use m_bandfft_kpt, only : bandfft_kpt, bandfft_kpt_type, bandfft_kpt_savetabs,bandfft_kpt_restoretabs, &
45 : prep_bandfft_tabs
46 : use m_pawtab, only : pawtab_type
47 : use m_paw_ij, only : paw_ij_type
48 : use m_mkffnl, only : mkffnl
49 : use m_mpinfo, only : proc_distrb_cycle
50 : use m_gputk
51 : use m_abi_linalg
52 :
53 : #if defined(HAVE_GPU)
54 : use m_gpu_toolbox
55 : #endif
56 :
57 : implicit none
58 :
59 : private
60 : !!***
61 :
62 : public :: fock_getghc
63 : public :: fock2ACE
64 : public :: fock_ACE_getghc
65 : !!***
66 :
67 : contains
68 : !!***
69 :
70 0 : subroutine select_ndat_occ_for_gpu(ndat_occ,nband_k,ndat,npw,cplex_fock,nfftf,ngfft,&
71 0 : n4,n5,n6,natom,nspinor,lmn2_size,usepaw,cprj,ieigen,need_ghc,optfor,optstr,nfgd_max)
72 :
73 : !Arguments ------------------------------------
74 : ! Scalars
75 : integer,intent(in) :: nband_k,ndat,npw,cplex_fock,nfftf,n4,n5,n6,ngfft(18)
76 : integer,intent(in) :: natom,nspinor,lmn2_size,usepaw,ieigen,nfgd_max
77 : logical,intent(in) :: optfor,optstr,need_ghc
78 : integer,intent(out) :: ndat_occ
79 : type(pawcprj_type),intent(in) :: cprj(natom,nspinor*ndat)
80 :
81 : !Local variables-------------------------------
82 : integer :: i,ider,nprojs
83 : integer,target :: t_fft(3)
84 : integer(kind=c_size_t) :: sum_mem,free_mem,fourdp_smem
85 :
86 : ! *************************************************************************
87 :
88 0 : ider=0
89 0 : if (usepaw==1) then
90 : nprojs=0
91 0 : do i = 1,natom
92 0 : nprojs = nprojs + cprj(i, 1)%nlmn
93 : end do
94 : end if
95 :
96 0 : t_fft(1) = ngfft(3);
97 0 : t_fft(2) = ngfft(2);
98 0 : t_fft(3) = ngfft(1);
99 :
100 : #ifdef HAVE_GPU
101 : call gpu_get_max_mem(free_mem)
102 : free_mem = 0.85 * free_mem ! Cutting 15% out to be safe
103 : #endif
104 :
105 0 : do i=1,nband_k
106 0 : if(modulo(nband_k,i)/=0) cycle
107 0 : ndat_occ=nband_k/i
108 0 : sum_mem = 0
109 :
110 : ! cwavef
111 0 : sum_mem = sum_mem + INT(2,c_size_t)*n4*n5*n6*ndat
112 0 : if(need_ghc) then
113 : ! ghc
114 0 : sum_mem = sum_mem + INT(2,c_size_t)*npw*ndat
115 : ! ghc1
116 0 : sum_mem = sum_mem + INT(2,c_size_t)*npw*ndat
117 : ! ghc2
118 0 : sum_mem = sum_mem + INT(2,c_size_t)*npw*ndat
119 : end if
120 :
121 : ! cwavef_r
122 0 : sum_mem = sum_mem + INT(2,c_size_t)*n4*n5*n6*ndat
123 : ! vlocpsi_r
124 0 : sum_mem = sum_mem + INT(cplex_fock,c_size_t)*nfftf*ndat
125 : ! work (ompgpu_fourwf internal array)
126 0 : sum_mem = sum_mem + INT(2,c_size_t)*n4*n5*n6*ndat*ndat_occ
127 : #ifdef HAVE_GPU
128 : call gpu_fft_get_estimate_work_size(3, c_loc(t_fft), FFT_Z2Z, ndat*ndat_occ, fourdp_smem);
129 : #endif
130 0 : sum_mem = sum_mem + fourdp_smem
131 :
132 : ! rhor_munu
133 0 : sum_mem = sum_mem + INT(cplex_fock,c_size_t)*nfftf*ndat_occ*ndat
134 : ! rhog_munu
135 0 : sum_mem = sum_mem + INT(2,c_size_t)*nfftf*ndat_occ*ndat
136 : ! vfock
137 0 : sum_mem = sum_mem + INT(cplex_fock,c_size_t)*nfftf*ndat_occ*ndat
138 : ! occ
139 0 : sum_mem = sum_mem + INT(ndat_occ,c_size_t)
140 :
141 : !*Additional arrays in case of paw
142 0 : if (usepaw==1) then
143 0 : if ((optfor).and.(ieigen/=0)) then
144 0 : ider=3
145 : ! forout
146 0 : sum_mem = sum_mem + INT(2,c_size_t)*npw*nspinor*ndat_occ*ndat
147 :
148 : ! dprojs (gemm_nonlop_ompgpu internal work array)
149 0 : sum_mem = sum_mem + INT(2,c_size_t)*nprojs*npw*3
150 : end if
151 :
152 0 : if (optstr.and.(ieigen/=0)) then
153 0 : ider=3
154 : ! strout
155 0 : sum_mem = sum_mem + INT(2,c_size_t)*npw*nspinor*ndat_occ*ndat
156 :
157 : ! dprojs (gemm_nonlop_ompgpu internal work array)
158 0 : sum_mem = sum_mem + INT(2,c_size_t)*nprojs*npw*6
159 : end if
160 : ! grnhat_12
161 0 : ider=ider*2 ! Overestimate this buffer to ensure it fits as we don't manage OpenMP pool of GPU memory
162 0 : sum_mem = sum_mem + INT(2,c_size_t)*nfgd_max*nspinor**2*3*natom*(ider/3)*ndat_occ*ndat
163 : ! gvnlxc
164 0 : sum_mem = sum_mem + INT(2,c_size_t)*npw*nspinor*ndat_occ*ndat
165 : ! rho12 (=nhat12 output of pawmknhat_psipsi_ndat)
166 0 : sum_mem = sum_mem + INT(2,c_size_t)*nfgd_max*nspinor**2*ndat_occ*ndat*natom
167 : ! cprj1 (paw_psipsi internal work array)
168 0 : sum_mem = sum_mem + INT(2,c_size_t)*nprojs*nspinor*ndat
169 : ! cprj2 (paw_psipsi internal work array)
170 0 : sum_mem = sum_mem + INT(2,c_size_t)*nprojs*nspinor*ndat
171 : ! cpf (paw_psipsi internal work array)
172 0 : sum_mem = sum_mem + INT(2,c_size_t)*lmn2_size*nspinor*ndat*ndat_occ*natom
173 :
174 : ! dijhat
175 0 : sum_mem = sum_mem + INT(2,c_size_t)*lmn2_size*nspinor*ndat*ndat_occ*natom
176 :
177 : ! projs (gemm_nonlop_ompgpu internal work array)
178 0 : sum_mem = sum_mem + INT(2,c_size_t)*nprojs*npw
179 : end if
180 :
181 : ! cwaveocc_r
182 0 : sum_mem = sum_mem + INT(2,c_size_t)*n4*n5*n6*ndat_occ
183 :
184 0 : sum_mem = sum_mem*dp
185 :
186 0 : if(sum_mem < free_mem) exit
187 : end do
188 :
189 0 : if(sum_mem > free_mem) then
190 0 : ABI_WARNING("Test case doesn't fit in GPU memory. Try to lower bandpp.")
191 : end if
192 : !#ifdef DEBUG_VERBOSE
193 0 : write(std_out,*) "-----------DEBUG fock_getghc%select_ndat_occ_for_gpu : "
194 0 : write(std_out,'(A,F10.3,1x,A)') "Considered free GPU memory : ", real(free_mem,dp)/(1024*1024), "MiB"
195 0 : write(std_out,'(A,I4)') "selected ndat_occ : ", ndat_occ
196 0 : write(std_out,'(A,F10.3,1x,A)') "Forecasted consumed GPU memory : ", real(sum_mem,dp)/(1024*1024), "MiB"
197 0 : write(std_out,*) "-----------END DEBUG fock_getghc%select_ndat_occ_for_gpu : "
198 : !#endif
199 :
200 0 : end subroutine select_ndat_occ_for_gpu
201 :
202 : !!****f* ABINIT/fock_getghc
203 : !! NAME
204 : !! fock_getghc
205 : !!
206 : !! FUNCTION
207 : !! Compute the matrix elements <G|Vx|psi> of the Fock operator.
208 : !!
209 : !! INPUTS
210 : !! cwavef(2,npw*nspinor*ndat)= planewave coefficients of wavefunctions on which Fock operator is applied.
211 : !! cwaveprj <type(pawcprj_type> = <cwavevf|proj>
212 : !! gs_ham <type(gs_hamiltonian_type)>=all data for the Hamiltonian to be applied
213 : !! mpi_enreg= information about MPI parallelization
214 : !!
215 : !! SIDE EFFECTS
216 : !! ghc(2,npw*ndat)= matrix elements <G|H|C> or <G|H-lambda.S|C> (if sij_opt>=0 or =-1 in getghc)
217 : !! contains the fock exchange term for cwavef at the end.
218 : !!
219 : !! NOTES
220 : !! The current version assumes that:
221 : !! * nspinor = 1
222 : !! * no "my_nspinor"
223 : !! * no restriction to the value of istwfk_bz (but must be tested in all case)
224 : !! * all the data for the occupied states (cgocc_bz) are the same as those for the current states (cg)
225 : !!
226 : !! SOURCE
227 :
228 29003 : subroutine fock_getghc(cwavef,cwaveprj,ghc,gs_ham,mpi_enreg,ndat)
229 :
230 : !Arguments ------------------------------------
231 : ! Scalars
232 : integer,intent(in) :: ndat
233 : type(MPI_type),intent(in) :: mpi_enreg
234 : type(gs_hamiltonian_type),target,intent(inout) :: gs_ham
235 : ! Arrays
236 : type(pawcprj_type),intent(inout) :: cwaveprj(:,:)
237 : real(dp),intent(inout) :: cwavef(:,:)!,ghc(2,gs_ham%npw_k)
238 : real(dp),intent(inout) :: ghc(:,:)
239 :
240 : !Local variables-------------------------------
241 : ! Scalars
242 : integer,parameter :: tim_fourwf_fock_getghc=10,tim_fourdp_fock_getghc=10
243 : integer :: bdtot_jindex,choice,cplex_fock,cplex_dij,cpopt,i1,i2,i3,ia,iatom,idat,idat_occ,idat_tot,iatm
244 : integer :: iband_cprj,ider,idir,idir1,ier,ii,ind,ipw,ieigen,ifft,itypat,izero,jband,jbg,jcg,jkg
245 : integer :: jkpt,my_jsppol,jstwfk,lmn2_size,mgfftf,mpw,n1,n2,n3,n4,n5,n6,ndat_occ,ndat_tot
246 : integer :: n1f,n2f,n3f,n4f,n5f,n6f,natom,nband_k,ndij,nfft,nfftf,nfftotf,nhat12_grdim,nnlout
247 : integer :: nfgd_max
248 : integer :: npw,npwj,nspden_fock,nspinor,nkpg,paw_opt,signs,tim_nonlop,gpu_option
249 : integer, save :: ncount=0
250 : logical :: need_ghc,qeq0
251 : real(dp),parameter :: weight1=one
252 : real(dp) :: doti,eigen,imcwf,imcwocc,imvloc,invucvol,recwf,recwocc,revloc,wtk,esum,esumi
253 : complex(dp) :: cinvucvol,cucvol
254 : type(fock_common_type),pointer :: fockcommon
255 : type(fock_BZ_type),pointer :: fockbz
256 : ! Arrays
257 : integer :: ngfft(18),ngfftf(18)
258 29003 : integer,pointer :: gboundf(:,:),kg_occ(:,:),gbound_kp(:,:)
259 : real(dp) :: fockstr(6),qphon(3),qvec_j(3),tsec(2),gsc_dum(2,0),rhodum(2,1)
260 : real(dp) :: rhodum0(0,1,1)
261 29003 : real(dp), allocatable :: dummytab(:,:),dijhat(:,:,:,:,:),dijhat_tmp(:,:,:),ffnl_kp_dum(:,:,:,:),kpg_kp(:,:),occ(:)
262 29003 : real(dp), allocatable, target :: gvnlxc(:,:),ghc1(:,:),ghc2(:,:),grnhat12(:,:,:,:,:,:),grnhat_12(:,:,:,:,:,:,:),forikpt(:,:,:)
263 29003 : real(dp), allocatable :: rho12(:,:,:,:,:,:),rhog_munu(:,:,:,:),rhor_munu(:,:,:,:),vlocpsi_r(:,:),strdat(:,:,:,:)
264 29003 : real(dp), allocatable :: vfock(:,:,:),psilocal(:,:,:),enlout_dum(:),vectin_dum(:,:),vqg(:),forout(:,:),strout(:,:),for1(:,:,:,:)
265 29003 : real(dp), allocatable,target ::cwavef_r(:,:,:,:),cwavef_rep(:,:),vdotr(:,:,:,:),vdoti(:),vfockstr(:,:,:)
266 29003 : real(dp), ABI_CONTIGUOUS pointer :: cwaveocc_r(:,:,:,:,:)
267 29003 : type(pawcprj_type),pointer :: cwaveocc_prj(:,:)
268 29003 : type(pawcprj_type),pointer :: cwaveocc_prj_rep(:,:)
269 :
270 : real(dp) :: rprimd(3,3),for12(3)
271 29003 : integer, ABI_CONTIGUOUS pointer :: atom_ifftsph(:,:),atom_nfgd(:)
272 29003 : real(dp), ABI_CONTIGUOUS pointer :: stress_ikpt(:,:),atom_rfgd(:,:,:)
273 :
274 :
275 : ! *************************************************************************
276 : !return
277 :
278 29003 : ncount=ncount+1
279 :
280 29003 : call timab(1504,1,tsec) ; call timab(1505,-1,tsec) ; call timab(1515,-1,tsec) ; call timab(1541,-1,tsec)
281 :
282 29003 : ABI_CHECK(associated(gs_ham%fockcommon),"fock_common must be associated!")
283 29003 : fockcommon => gs_ham%fockcommon
284 29003 : ABI_CHECK(associated(gs_ham%fockbz),"fock_bz must be associated!")
285 29003 : fockbz => gs_ham%fockbz
286 :
287 29003 : ABI_CHECK(gs_ham%nspinor==1,"only allowed for nspinor=1!")
288 29003 : ABI_CHECK(gs_ham%npw_k==gs_ham%npw_kp,"only allowed for npw_k=npw_kp (ground state)!")
289 29003 : if (fockcommon%usepaw==1) then
290 2061 : ABI_CHECK((size(cwaveprj,1)==gs_ham%natom.and.size(cwaveprj,2)==gs_ham%nspinor*ndat),"error on cwaveprj dims")
291 : end if
292 29003 : need_ghc=(size(ghc,2)>0)
293 :
294 : !Some constants
295 29003 : invucvol=1.d0/sqrt(gs_ham%ucvol)
296 29003 : cinvucvol=dcmplx(invucvol,0.0_dp)
297 29003 : cucvol=dcmplx(sqrt(gs_ham%ucvol),0.0_dp)
298 29003 : call matr3inv(gs_ham%gprimd,rprimd)
299 29003 : cplex_fock=2;nspden_fock=1
300 29003 : natom=fockcommon%natom
301 29003 : nspinor=gs_ham%nspinor
302 : mpw=maxval(fockbz%npwarr)
303 29003 : npw=gs_ham%npw_k
304 29003 : gpu_option=gs_ham%gpu_option
305 29003 : ider=0;izero=0
306 29003 : if (fockcommon%usepaw==1) then
307 39159 : nfft =fockcommon%pawfgr%nfftc ; ngfft =fockcommon%pawfgr%ngfftc
308 39159 : nfftf=fockcommon%pawfgr%nfft ; ngfftf=fockcommon%pawfgr%ngfft
309 2061 : mgfftf=fockcommon%pawfgr%mgfft
310 : else
311 26942 : nfft =gs_ham%nfft ; nfftf =nfft
312 511898 : ngfft=gs_ham%ngfft ; ngfftf=ngfft
313 26942 : mgfftf=gs_ham%mgfft
314 : end if
315 29003 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
316 29003 : n4=ngfft(4);n5=ngfft(5);n6=ngfft(6)
317 29003 : n1f=ngfftf(1);n2f=ngfftf(2);n3f=ngfftf(3)
318 29003 : n4f=ngfftf(4);n5f=ngfftf(5);n6f=ngfftf(6)
319 :
320 : !*Max number of fine-grid points in a PAW augmentation sphere (over all atoms).
321 : !*Used to size the small per-atom rho12/grnhat_12 sphere buffers instead of the full FFT grid.
322 29003 : nfgd_max=1
323 29003 : if (fockcommon%usepaw==1) then
324 6183 : nfgd_max=max(1,maxval(fockcommon%pawfgrtab(1:natom)%nfgd))
325 : end if
326 :
327 : ! ===========================
328 : ! === Initialize arrays ===
329 : ! ===========================
330 : ! transient optfor and optstress
331 : ! fockcommon%optfor=.false.
332 : ! fockcommon%optstr=.false.
333 : !*Initialization of local pointers
334 : !*Initialization of the array cwavef_r
335 : !*cwavef_r = current wavefunction in r-space
336 145015 : ABI_MALLOC(cwavef_r,(2,n4f,n5f,n6f*ndat))
337 : #ifdef HAVE_OPENMP_OFFLOAD
338 : !$OMP TARGET ENTER DATA MAP(alloc:cwavef_r) IF(gpu_option==ABI_GPU_OPENMP)
339 : #endif
340 : !*dummytab = variables for fourwf
341 87009 : ABI_MALLOC(dummytab,(2,nfft*ndat))
342 : !*vqg = 4pi/(G+q)**2
343 87009 : ABI_MALLOC(vqg,(nfftf))
344 :
345 29003 : if(need_ghc) then
346 : !*Initialization of the array ghc1
347 : !*ghc1 will contain the exact exchange contribution to the Hamiltonian
348 80853 : ABI_MALLOC(ghc1,(2,npw*ndat))
349 53902 : ABI_MALLOC(ghc2,(2,npw*ndat))
350 : #ifdef HAVE_OPENMP_OFFLOAD
351 : !$OMP TARGET ENTER DATA MAP(alloc:ghc1,ghc2,vqg) IF(gpu_option==ABI_GPU_OPENMP)
352 : #endif
353 26951 : if(gpu_option==ABI_GPU_DISABLED) then
354 20048969 : ghc1=zero
355 20048969 : ghc2=zero
356 0 : else if(gpu_option==ABI_GPU_OPENMP) then
357 0 : call gpu_set_to_zero(ghc1, int(2,c_size_t)*npw*ndat)
358 0 : call gpu_set_to_zero(ghc2, int(2,c_size_t)*npw*ndat)
359 : end if
360 : end if
361 : !*Initialization of the array vlocpsi_r
362 : !*vlocpsi_r = partial local Fock operator applied to cwavef in r-space and summed over all occupied (jkpt,mu)
363 116012 : ABI_MALLOC(vlocpsi_r,(cplex_fock*nfftf,ndat))
364 : #ifdef HAVE_OPENMP_OFFLOAD
365 : !$OMP TARGET ENTER DATA MAP(alloc:vlocpsi_r) IF(gpu_option==ABI_GPU_OPENMP)
366 : #endif
367 29003 : if(gpu_option==ABI_GPU_DISABLED) then
368 456957546 : vlocpsi_r=zero
369 0 : else if(gpu_option==ABI_GPU_OPENMP) then
370 0 : call gpu_set_to_zero(vlocpsi_r, int(cplex_fock,c_size_t)*nfftf*ndat)
371 : end if
372 :
373 : !*Additional arrays in case of paw
374 29003 : if (fockcommon%usepaw==1) then
375 2061 : nhat12_grdim=0
376 : end if
377 :
378 29003 : if (fockcommon%usepaw==1) then
379 2061 : if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then
380 4208 : ABI_MALLOC(forikpt,(3,natom,ndat))
381 10520 : forikpt=zero
382 : end if
383 : end if
384 29003 : if (fockcommon%usepaw==1.or.fockcommon%optstr) then
385 15596 : ABI_MALLOC(gboundf,(2*mgfftf+8,2))
386 3899 : call sphereboundary(gboundf,gs_ham%istwf_k,gs_ham%kg_k,mgfftf,npw)
387 : else
388 25104 : gboundf=>gs_ham%gbound_k
389 : end if
390 : ! ==========================================
391 : ! === Get cwavef in real space using FFT ===
392 : ! ==========================================
393 29003 : if(gpu_option==ABI_GPU_DISABLED) then
394 740383784 : cwavef_r=zero
395 0 : else if(gpu_option==ABI_GPU_OPENMP) then
396 0 : call gpu_set_to_zero(cwavef_r, int(2,c_size_t)*n4f*n5f*n6f*ndat)
397 : end if
398 29003 : call timab(1515,2,tsec) ; call timab(1541,-2,tsec) ; call timab(1512,-1,tsec)
399 : call fourwf(0,rhodum0,cwavef,rhodum,cwavef_r,gboundf,gboundf,gs_ham%istwf_k,gs_ham%kg_k,gs_ham%kg_k,&
400 : & mgfftf,mpi_enreg,ndat,ngfftf,npw,1,n4f,n5f,n6f,0,tim_fourwf_fock_getghc,weight1,weight1,&
401 29003 : & gpu_option=gs_ham%gpu_option)
402 29003 : call timab(1512,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1541,-1,tsec)
403 29003 : if(gpu_option==ABI_GPU_DISABLED) then
404 740383784 : cwavef_r=cwavef_r*invucvol
405 0 : else if(gpu_option==ABI_GPU_OPENMP) then
406 0 : call abi_xscal(n4f*n5f*n6f*ndat,cinvucvol,cwavef_r,1,x_cplx=2,gpu_option=gpu_option)
407 : end if
408 :
409 : ! =====================================================
410 : ! === Select the states in cgocc_bz with the same spin ===
411 : ! =====================================================
412 : !* Initialization of the indices/shifts, according to the value of isppol
413 : !* bdtot_jindex = shift to be applied on the location of data in the array occ_bz ?
414 29003 : bdtot_jindex=0
415 : !* jbg = shift to be applied on the location of data in the array cprj/occ
416 29003 : jbg=0;jcg=0
417 29003 : my_jsppol=fockcommon%isppol
418 29003 : if((fockcommon%isppol==2).and.(mpi_enreg%nproc_spkpt/=1)) my_jsppol=1
419 :
420 : !===================================
421 : !=== Loop on the k-points in IBZ ===
422 : !===================================
423 29003 : jkg=0
424 :
425 29003 : if (associated(gs_ham%ph3d_kp)) then
426 29003 : nullify (gs_ham%ph3d_kp)
427 : end if
428 :
429 29003 : call timab(1505,2,tsec) ; call timab(1506,-1,tsec) ; call timab(1541,-2,tsec)
430 :
431 231529 : do jkpt=1,fockbz%mkpt
432 :
433 202526 : if(fockbz%nbandocc_bz(jkpt,my_jsppol)==0) cycle
434 :
435 202526 : call timab(1521,1,tsec)
436 :
437 : !* nband_k = number of bands at point k_j
438 202526 : nband_k=fockbz%nbandocc_bz(jkpt,my_jsppol)
439 :
440 : ! Select ndat_occ value :
441 : ! GPU : check GPU memory available and compute maximum value for ndat_occ
442 : ! CPU : 4, because it seems to be optimal from my observations
443 202526 : if(gpu_option/=ABI_GPU_DISABLED) then
444 0 : lmn2_size=0
445 0 : if(fockcommon%usepaw==1) lmn2_size=fockcommon%pawtab(1)%lmn2_size
446 : call select_ndat_occ_for_gpu(ndat_occ,nband_k,ndat,npw,cplex_fock,&
447 : & nfftf,ngfftf,n4f,n5f,n6f,natom,nspinor,lmn2_size,&
448 0 : & fockcommon%usepaw,cwaveprj,fockcommon%ieigen,need_ghc,fockcommon%optfor,fockcommon%optstr,nfgd_max)
449 : else
450 202526 : ndat_occ=min(nband_k,16)
451 202526 : do ii=1,nband_k
452 202526 : if(modulo(nband_k,ndat_occ)==0) exit
453 202526 : ndat_occ=ndat_occ-1
454 : end do
455 : end if
456 : !* ndat_tot = combined batch size (ndat*ndat_occ), used to batch nonlop calls
457 202526 : ndat_tot=ndat*ndat_occ
458 :
459 : !* wtk = weight in BZ of this k point
460 202526 : wtk=fockbz%wtk_bz(jkpt) !*sqrt(gs_ham%ucvol)
461 : !* jstwfk= how is stored the wavefunction
462 202526 : jstwfk=fockbz%istwfk_bz(jkpt)
463 : !* npwj= number of plane wave in basis for the wavefunction
464 202526 : npwj=fockbz%npwarr(jkpt)
465 : !* Basis sphere of G vectors
466 202526 : if (allocated(fockbz%cgocc)) then
467 202526 : gbound_kp => fockbz%gbound_bz(:,:,jkpt)
468 202526 : kg_occ => fockbz%kg_bz(:,1+jkg:npwj+jkg)
469 : end if
470 :
471 : !*rhormunu = overlap matrix between cwavef and (jkpt,mu) in R-space
472 1215156 : ABI_MALLOC(rhor_munu,(cplex_fock,nfftf,ndat_occ,ndat))
473 : #ifdef HAVE_OPENMP_OFFLOAD
474 : !$OMP TARGET ENTER DATA MAP(alloc:rhor_munu) IF(gpu_option==ABI_GPU_OPENMP)
475 : #endif
476 : !*rhogmunu = overlap matrix between cwavef and (jkpt,mu) in G-space
477 1012630 : ABI_MALLOC(rhog_munu,(2,nfftf,ndat_occ,ndat))
478 : #ifdef HAVE_OPENMP_OFFLOAD
479 : !$OMP TARGET ENTER DATA MAP(alloc:rhog_munu) IF(gpu_option==ABI_GPU_OPENMP)
480 : #endif
481 : !*vfock = Fock potential
482 1012630 : ABI_MALLOC(vfock,(cplex_fock*nfftf,ndat_occ,ndat))
483 : #ifdef HAVE_OPENMP_OFFLOAD
484 : !$OMP TARGET ENTER DATA MAP(alloc:vfock) IF(gpu_option==ABI_GPU_OPENMP)
485 : #endif
486 607578 : ABI_MALLOC(occ,(ndat_occ))
487 : #ifdef HAVE_OPENMP_OFFLOAD
488 : !$OMP TARGET ENTER DATA MAP(alloc:occ) IF(gpu_option==ABI_GPU_OPENMP)
489 : #endif
490 :
491 : !*Additional arrays in case of paw
492 202526 : if (fockcommon%usepaw==1) then
493 9678 : if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then
494 2804 : ider=3
495 8412 : ABI_MALLOC(forout,(2,npw*nspinor*ndat_tot))
496 : #ifdef HAVE_OPENMP_OFFLOAD
497 : !$OMP TARGET ENTER DATA MAP(alloc:forout) IF(gpu_option==ABI_GPU_OPENMP)
498 : #endif
499 : end if
500 :
501 9678 : if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then
502 1870 : ider=3
503 5610 : ABI_MALLOC(strout,(2,npw*nspinor*ndat_tot))
504 : #ifdef HAVE_OPENMP_OFFLOAD
505 : !$OMP TARGET ENTER DATA MAP(alloc:strout) IF(gpu_option==ABI_GPU_OPENMP)
506 : #endif
507 : end if
508 77424 : ABI_MALLOC(grnhat_12,(2,nfgd_max,nspinor**2,3,natom*(ider/3),ndat_occ,ndat))
509 29034 : ABI_MALLOC(gvnlxc,(2,npw*nspinor*ndat_tot))
510 9678 : ABI_MALLOC(grnhat12,(1,1,1,1,1,1))
511 67746 : ABI_MALLOC(rho12,(2,nfgd_max,nspinor**2,ndat_occ,ndat,natom))
512 :
513 29034 : ABI_MALLOC(atom_nfgd, (natom))
514 29034 : do iatom=1,natom
515 29034 : atom_nfgd(iatom) = fockcommon%pawfgrtab(iatom)%nfgd
516 : end do
517 38712 : ABI_MALLOC(atom_ifftsph, (nfgd_max, natom))
518 38712 : ABI_MALLOC(atom_rfgd, (3, nfgd_max, natom))
519 29034 : do iatom=1,natom
520 15460232 : atom_ifftsph(1:atom_nfgd(iatom),iatom) = fockcommon%pawfgrtab(iatom)%ifftsph(1:atom_nfgd(iatom))
521 61792538 : atom_rfgd(:,1:atom_nfgd(iatom),iatom) = fockcommon%pawfgrtab(iatom)%rfgd(:,1:atom_nfgd(iatom))
522 : end do
523 :
524 : #ifdef HAVE_OPENMP_OFFLOAD
525 : !$OMP TARGET ENTER DATA MAP(alloc:grnhat_12) IF(gpu_option==ABI_GPU_OPENMP .and. ider==3)
526 : !$OMP TARGET ENTER DATA MAP(alloc:gvnlxc) IF(gpu_option==ABI_GPU_OPENMP)
527 : !$OMP TARGET ENTER DATA MAP(alloc:rho12) IF(gpu_option==ABI_GPU_OPENMP)
528 : !$OMP TARGET ENTER DATA MAP(alloc:atom_ifftsph,atom_nfgd,atom_rfgd) IF(gpu_option==ABI_GPU_OPENMP)
529 : !$OMP TARGET UPDATE TO(atom_ifftsph,atom_nfgd,atom_rfgd) IF(gpu_option==ABI_GPU_OPENMP)
530 : #endif
531 : end if
532 :
533 : !* Load k^prime hamiltonian in the gs_ham datastructure
534 : ! Note: ffnl_kp / ph3d_kp / gbound_kp are not used
535 :
536 202526 : if (.not. associated(gs_ham%ph3d_kp)) then
537 810104 : ABI_MALLOC(gs_ham%ph3d_kp,(2,npwj,gs_ham%matblk))
538 : end if
539 :
540 : call gs_ham%load_kprime(kpt_kp=fockbz%kptns_bz(:,jkpt),&
541 202526 : & istwf_kp=jstwfk,npw_kp=npwj,kg_kp=fockbz%kg_bz(:,1+jkg:npwj+jkg))
542 : !* Some temporary allocations needed for PAW
543 202526 : if (fockcommon%usepaw==1) then
544 29034 : ABI_MALLOC(enlout_dum,(ndat_tot))
545 29034 : ABI_MALLOC(vectin_dum,(2,npwj*nspinor*ndat_tot))
546 28480275 : vectin_dum=zero
547 48390 : ABI_MALLOC(ffnl_kp_dum,(npwj,1,gs_ham%lmnmax,gs_ham%ntypat))
548 9678 : nkpg=size(gs_ham%kpg_k,2)
549 38712 : ABI_MALLOC(kpg_kp,(npwj,nkpg))
550 9678 : if (nkpg>0) then
551 9678 : call mkkpg(gs_ham%kg_kp,kpg_kp,gs_ham%kpt_kp,nkpg,npwj)
552 : end if
553 9678 : call gs_ham%load_kprime(ffnl_kp=ffnl_kp_dum,kpg_kp=kpg_kp)
554 : end if
555 :
556 : ! ======================================
557 : ! === Calculate the vector q=k_i-k_j ===
558 : ! ======================================
559 : !* Evaluation of kpoint_j, the considered k-point in reduced coordinates
560 : ! kpoint_j(:)=fockbz%kptns_bz(:,jkpt)
561 : !* the vector qvec is expressed in reduced coordinates.
562 : ! qvec(:)=kpoint_i(:)-kpoint_j(:)
563 810104 : qvec_j(:)=gs_ham%kpt_k(:)-fockbz%kptns_bz(:,jkpt)
564 202526 : qeq0=(qvec_j(1)**2+qvec_j(2)**2+qvec_j(3)**2<1.d-15)
565 :
566 : ! Get the Coulomb interaction in reciprocal space
567 202526 : call bare_vqg(qvec_j,fockcommon,gs_ham%gmet,nfftf,fockbz%nkpt_bz,ngfftf,gs_ham%ucvol,vqg)
568 : #ifdef HAVE_OPENMP_OFFLOAD
569 : !$OMP TARGET UPDATE TO(vqg) IF(gpu_option==ABI_GPU_OPENMP)
570 : #endif
571 :
572 202526 : call timab(1521,2,tsec)
573 :
574 : ! =================================================
575 : ! === Loop on the band indices jband of cgocc_k ===
576 : ! =================================================
577 405052 : do jband=1,nband_k,ndat_occ
578 :
579 : !* occ = occupancy of jband at this k point
580 1115372 : occ(1:ndat_occ)=fockbz%occ_bz(jband+bdtot_jindex:jband+ndat_occ-1+bdtot_jindex,my_jsppol)
581 1317898 : if(maxval(occ(1:ndat_occ))<tol8) cycle
582 : #ifdef HAVE_OPENMP_OFFLOAD
583 : !$OMP TARGET UPDATE TO(occ) IF(gpu_option==ABI_GPU_OPENMP)
584 : #endif
585 :
586 : ! This timing is placed after the cycle ...
587 202526 : call timab(1522,1,tsec) ; call timab(1542,-1,tsec)
588 :
589 : ! ==============================================
590 : ! === Get cwaveocc_r in real space using FFT ===
591 : ! ==============================================
592 202526 : if (allocated(fockbz%cwaveocc_bz)) then
593 0 : cwaveocc_r => fockbz%cwaveocc_bz(:,:,:,:,jband+jbg:jband+jbg+ndat_occ-1,my_jsppol)
594 : else
595 1215156 : ABI_MALLOC(cwaveocc_r,(2,n4f,n5f,n6f,ndat_occ))
596 : #ifdef HAVE_OPENMP_OFFLOAD
597 : !$OMP TARGET ENTER DATA MAP(alloc:cwaveocc_r) IF(gpu_option==ABI_GPU_OPENMP)
598 : #endif
599 202526 : if(gpu_option==ABI_GPU_DISABLED) then
600 13556461613 : cwaveocc_r=zero
601 0 : else if(gpu_option==ABI_GPU_OPENMP) then
602 0 : call gpu_set_to_zero(cwaveocc_r, int(2,c_size_t)*n4f*n5f*n6f*ndat_occ)
603 : end if
604 202526 : call timab(1515,2,tsec) ; call timab(1512,-1,tsec) ; call timab(1542,-2,tsec)
605 : call fourwf(1,rhodum0,fockbz%cgocc(:,1+jcg+npwj*(jband-1):jcg+(jband+ndat_occ-1)*npwj,my_jsppol),rhodum,cwaveocc_r, &
606 : & gbound_kp,gbound_kp,jstwfk,kg_occ,kg_occ,mgfftf,mpi_enreg,ndat_occ,ngfftf,&
607 202526 : & npwj,1,n4f,n5f,n6f,0,tim_fourwf_fock_getghc,weight1,weight1,gpu_option=gs_ham%gpu_option)
608 202526 : call timab(1512,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1542,-1,tsec)
609 202526 : if(gpu_option==ABI_GPU_DISABLED) then
610 13556461613 : cwaveocc_r=cwaveocc_r*invucvol
611 0 : else if(gpu_option==ABI_GPU_OPENMP) then
612 0 : call abi_xscal(n4f*n5f*n6f*ndat_occ,cinvucvol,cwaveocc_r,1,x_cplx=2,gpu_option=gpu_option)
613 : end if
614 : end if
615 :
616 : ! ================================================
617 : ! === Get the overlap density matrix rhor_munu ===
618 : ! ================================================
619 : !* Calculate the overlap density matrix in real space = conj(cwaveocc_r)*cwavef_r
620 : !* rhor_munu will contain the overlap density matrix.
621 : ! vfock=-int{conj(cwaveocc_r)*cwavef_r*dr'/|r-r'|}
622 :
623 202526 : call timab(1522,2,tsec) ; call timab(1542,-2,tsec) ; call timab(1523,-1,tsec)
624 :
625 202526 : if(gpu_option==ABI_GPU_DISABLED) then
626 : !$OMP PARALLEL DO COLLAPSE(2) &
627 : !$OMP& PRIVATE(ind,imcwf,recwf,recwocc,imcwocc)
628 405052 : do idat=1,ndat
629 1317898 : do idat_occ=1,ndat_occ
630 15580767 : do i3=1,n3f
631 251108924 : do i2=1,n2f
632 4349064708 : do i1=1,n1f
633 4098868630 : ind=i1+(i2-1)*n1f+(i3-1)*n2f*n1f
634 4098868630 : recwf =cwavef_r(1,i1,i2,(idat-1)*n3f+i3)
635 4098868630 : imcwf =cwavef_r(2,i1,i2,(idat-1)*n3f+i3)
636 4098868630 : recwocc=cwaveocc_r(1,i1,i2,i3,idat_occ)
637 4098868630 : imcwocc=cwaveocc_r(2,i1,i2,i3,idat_occ)
638 4098868630 : rhor_munu(1,ind,idat_occ,idat)= recwocc*recwf+imcwocc*imcwf
639 4334599313 : rhor_munu(2,ind,idat_occ,idat)= recwocc*imcwf-imcwocc*recwf
640 : end do ! i1
641 : end do ! i2
642 : end do ! i3
643 : end do ! idat_occ
644 : end do ! idat
645 : else if(gpu_option==ABI_GPU_OPENMP) then
646 : #ifdef HAVE_OPENMP_OFFLOAD
647 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
648 : !$OMP& MAP(to:rhor_munu,cwavef_r,cwaveocc_r) PRIVATE(idat,idat_occ)
649 : do idat=1,ndat
650 : do idat_occ=1,ndat_occ
651 : !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(ind,i3,i2,i1)
652 : do i3=1,n3f
653 : do i2=1,n2f
654 : do i1=1,n1f
655 : ind=i1+(i2-1)*n1f+(i3-1)*n2f*n1f
656 : rhor_munu(1,ind,idat_occ,idat)= cwaveocc_r(1,i1,i2,i3,idat_occ)*cwavef_r(1,i1,i2,(idat-1)*n3f+i3)&
657 : +cwaveocc_r(2,i1,i2,i3,idat_occ)*cwavef_r(2,i1,i2,(idat-1)*n3f+i3)
658 : rhor_munu(2,ind,idat_occ,idat)= cwaveocc_r(1,i1,i2,i3,idat_occ)*cwavef_r(2,i1,i2,(idat-1)*n3f+i3)&
659 : -cwaveocc_r(2,i1,i2,i3,idat_occ)*cwavef_r(1,i1,i2,(idat-1)*n3f+i3)
660 : end do ! i1
661 : end do ! i2
662 : end do ! i3
663 : end do ! idat_occ
664 : end do ! idat
665 : #endif
666 : end if ! gpu_option
667 :
668 202526 : call timab(1523,2,tsec)
669 :
670 : ! =======================================================
671 : ! === Add compensation charge density in the PAW case ===
672 : ! =======================================================
673 :
674 202526 : call timab(1524,-1,tsec) ; call timab(1544,-1,tsec)
675 :
676 202526 : if (fockcommon%usepaw==1) then
677 :
678 9678 : iband_cprj=(my_jsppol-1)*fockbz%mkptband+jbg+jband
679 9678 : cwaveocc_prj=>fockbz%cwaveocc_prj(:,iband_cprj:iband_cprj+ndat_occ*nspinor-1)
680 :
681 : call pawmknhat_psipsi_ndat(cwaveprj(:,:),cwaveocc_prj(:,:),&
682 : & ider,izero,natom,natom,nfftf,ngfftf,&
683 : & nhat12_grdim,nspinor,fockcommon%ntypat,ndat,ndat_occ,fockbz%pawang,fockcommon%pawfgrtab,grnhat12,&
684 : & rho12,gs_ham%nattyp,&
685 : & fockcommon%pawtab,gprimd=gs_ham%gprimd,grnhat_12=grnhat_12,qphon=qvec_j,&
686 9678 : & xred=gs_ham%xred,atindx=gs_ham%atindx,gpu_option=gpu_option)
687 :
688 9678 : if(gpu_option==ABI_GPU_DISABLED) then
689 : !$OMP PARALLEL DO COLLAPSE(2) &
690 : !$OMP& PRIVATE(idat,idat_occ,iatom,ifft,ind)
691 19356 : do idat=1,ndat
692 57729 : do idat_occ=1,ndat_occ
693 124797 : do iatom=1,natom
694 60932723 : do ifft=1,atom_nfgd(iatom)
695 60817604 : ind = atom_ifftsph(ifft,iatom)
696 60817604 : rhor_munu(1,ind,idat_occ,idat) = rhor_munu(1,ind,idat_occ,idat) + rho12(1,ifft,nspinor,idat_occ,idat,iatom)
697 60894350 : rhor_munu(2,ind,idat_occ,idat) = rhor_munu(2,ind,idat_occ,idat) - rho12(2,ifft,nspinor,idat_occ,idat,iatom)
698 : end do ! ifft
699 : end do ! iatom
700 : end do ! idat_occ
701 : end do ! idat
702 : else if(gpu_option==ABI_GPU_OPENMP) then
703 : #ifdef HAVE_OPENMP_OFFLOAD
704 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
705 : !$OMP& MAP(to:rhor_munu,rho12,atom_nfgd,atom_ifftsph) PRIVATE(idat,idat_occ)
706 : do idat=1,ndat
707 : do idat_occ=1,ndat_occ
708 : !$OMP PARALLEL DO PRIVATE(iatom,ifft,ind)
709 : do iatom=1,natom
710 : do ifft=1,atom_nfgd(iatom)
711 : ind = atom_ifftsph(ifft,iatom)
712 : rhor_munu(1,ind,idat_occ,idat) = rhor_munu(1,ind,idat_occ,idat) + rho12(1,ifft,nspinor,idat_occ,idat,iatom)
713 : rhor_munu(2,ind,idat_occ,idat) = rhor_munu(2,ind,idat_occ,idat) - rho12(2,ifft,nspinor,idat_occ,idat,iatom)
714 : end do ! ifft
715 : end do ! iatom
716 : end do ! idat_occ
717 : end do ! idat
718 : #endif
719 : end if
720 : end if
721 :
722 :
723 202526 : call timab(1515,2,tsec) ; call timab(1513,-1,tsec) ; call timab(1544,-2,tsec)
724 : ! Perform an FFT using fourwf to get rhog_munu = FFT^-1(rhor_munu)
725 : call fourdp(cplex_fock,rhog_munu,rhor_munu,-1,mpi_enreg,nfftf,ndat*ndat_occ,&
726 202526 : & ngfftf,tim_fourdp_fock_getghc,gpu_option=gpu_option)
727 202526 : call timab(1513,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1544,-1,tsec)
728 :
729 202526 : if(fockcommon%optstr.and.(fockcommon%ieigen/=0)) then
730 51840 : ABI_MALLOC(vfockstr, (6,ndat_occ,ndat))
731 : call strfock(fockcommon,gs_ham%gprimd,vfockstr,&
732 : & mpi_enreg,nfftf,ngfftf,fockbz%nkpt_bz,ndat*ndat_occ,rhog_munu,gs_ham%ucvol,&
733 12960 : & qvec_j,gpu_option=gpu_option)
734 25920 : do idat=1,ndat
735 87080 : do idat_occ=1,ndat_occ
736 808040 : fockcommon%stress_ikpt(:,fockcommon%ieigen+idat-1)=fockcommon%stress_ikpt(:,fockcommon%ieigen+idat-1)+vfockstr(:,idat_occ,idat)*occ(idat_occ)*wtk
737 : end do ! idat_occ
738 : end do ! idat
739 12960 : ABI_FREE(vfockstr)
740 12960 : if (fockcommon%usepaw==0.and.(.not.need_ghc)) then
741 11090 : if (allocated(fockbz%cgocc)) then
742 : #ifdef HAVE_OPENMP_OFFLOAD
743 : !$OMP TARGET EXIT DATA MAP(delete:cwaveocc_r) IF(gpu_option==ABI_GPU_OPENMP)
744 : #endif
745 11090 : ABI_FREE(cwaveocc_r)
746 : end if
747 11090 : call timab(1524,2,tsec) ; call timab(1544,-2,tsec)
748 11090 : cycle
749 : end if
750 : end if
751 191436 : call timab(1524,2,tsec) ; call timab(1544,-2,tsec)
752 :
753 : ! ===================================================
754 : ! === Calculate the local potential vfockloc_munu ===
755 : ! ===================================================
756 : !* Apply the Poisson solver to "rhog_munu" while taking into account the effect of the vector "qvec"
757 : !* This is precisely what is done in the subroutine hartre, with option cplex=2.
758 : !* vfock will contain the local Fock potential, the result of hartre routine.
759 : !* vfock = FFT( rhog_munu/|g+qvec|^2 )
760 191436 : call timab(1525,-1,tsec) ; call timab(1545,-1,tsec)
761 : #if 0
762 :
763 : do idat=1,ndat
764 : do idat_occ=1,ndat_occ
765 : call timab(1515,-2,tsec) ; call timab(1513,-1,tsec)
766 : call hartre(cplex_fock,fockcommon%gsqcut,fockcommon%usepaw,mpi_enreg,nfftf,ngfftf,&
767 : & mpi_enreg%paral_kgb,rhog_munu(:,:,idat_occ,idat),rprimd,vfock(:,idat_occ,idat),divgq0=fock%divgq0,qpt=qvec_j)
768 : call timab(1513,2,tsec) ; call timab(1515,-1,tsec)
769 : end do ! idat_occ
770 : end do ! idat
771 :
772 : #else
773 191436 : if(gpu_option==ABI_GPU_DISABLED) then
774 : !$OMP PARALLEL DO COLLAPSE(2)
775 382872 : do idat=1,ndat
776 1242043 : do idat_occ=1,ndat_occ
777 3863349637 : do ifft=1,nfftf
778 3862299030 : rhog_munu(1,ifft,idat_occ,idat) = rhog_munu(1,ifft,idat_occ,idat) * vqg(ifft)
779 3863158201 : rhog_munu(2,ifft,idat_occ,idat) = rhog_munu(2,ifft,idat_occ,idat) * vqg(ifft)
780 : end do
781 : end do ! idat_occ
782 : end do ! idat
783 : else if(gpu_option==ABI_GPU_OPENMP) then
784 : #ifdef HAVE_OPENMP_OFFLOAD
785 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
786 : !$OMP& MAP(to:rhog_munu,vqg) PRIVATE(idat,idat_occ)
787 : do idat=1,ndat
788 : do idat_occ=1,ndat_occ
789 : !$OMP PARALLEL DO PRIVATE(ifft)
790 : do ifft=1,nfftf
791 : rhog_munu(1,ifft,idat_occ,idat) = rhog_munu(1,ifft,idat_occ,idat) * vqg(ifft)
792 : rhog_munu(2,ifft,idat_occ,idat) = rhog_munu(2,ifft,idat_occ,idat) * vqg(ifft)
793 : end do
794 : end do ! idat_occ
795 : end do ! idat
796 : #endif
797 : end if ! gpu_option
798 :
799 191436 : call timab(1515,2,tsec) ; call timab(1513,-1,tsec) ; call timab(1545,-2,tsec)
800 191436 : call fourdp(cplex_fock,rhog_munu,vfock,+1,mpi_enreg,nfftf,ndat*ndat_occ,ngfftf,tim_fourdp_fock_getghc,gpu_option=gpu_option)
801 191436 : call timab(1513,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1545,-1,tsec)
802 : #endif
803 191436 : call timab(1525,-2,tsec) ; call timab(1545,-2,tsec)
804 :
805 : !===============================================================
806 : !======== Calculate Dij_Fock_hat contribution in case of PAW ===
807 : !===============================================================
808 :
809 191436 : call timab(1526,-1,tsec) ; call timab(1546,-1,tsec)
810 :
811 191436 : if (fockcommon%usepaw==1) then
812 38712 : qphon=qvec_j;nfftotf=product(ngfftf(1:3))
813 9678 : ndij=nspden_fock
814 : ! dimekb1 is dimensioned as cplex_dij*lmnmax*(lmnmax+1)/2
815 9678 : cplex_dij=2*gs_ham%dimekb1/(gs_ham%lmnmax*(gs_ham%lmnmax+1))
816 67746 : ABI_MALLOC(dijhat,(gs_ham%dimekb1,natom,ndij,ndat_tot,cplex_fock))
817 5861730 : dijhat=zero
818 :
819 : #ifdef HAVE_OPENMP_OFFLOAD
820 : !!$OMP TARGET UPDATE FROM(vfock) IF(gpu_option==ABI_GPU_OPENMP)
821 : #endif
822 9678 : iatm=0
823 20070 : do itypat=1,gs_ham%ntypat
824 10392 : lmn2_size=fockcommon%pawtab(itypat)%lmn2_size
825 51960 : ABI_MALLOC(dijhat_tmp,(cplex_fock*cplex_dij*lmn2_size,ndij*ndat_occ*ndat,gs_ham%nattyp(itypat)))
826 : call pawdijhat_ndat(dijhat_tmp,cplex_dij,cplex_fock,gs_ham%gprimd,iatm,&
827 : & natom,ndij,nfftf,nfftotf,nspden_fock,nspden_fock,ndat_occ*ndat,&
828 : & gs_ham%nattyp(itypat),fockbz%pawang,fockcommon%pawfgrtab,&
829 : & fockcommon%pawtab(itypat),vfock,qphon,gs_ham%ucvol,gs_ham%xred,&
830 10392 : & gpu_option=gpu_option)
831 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(ia,idat,idat_occ,idat_tot,ii,iatom,ind)
832 29748 : do ia=1,gs_ham%nattyp(itypat)
833 49104 : do idat=1,ndat
834 115458 : do idat_occ=1,ndat_occ
835 76746 : idat_tot=(idat-1)*ndat_occ+idat_occ
836 76746 : iatom=iatm+ia
837 249594 : do ii=1,cplex_fock
838 153492 : ind=(ii-1)*lmn2_size*cplex_dij
839 : dijhat(1:cplex_dij*lmn2_size,iatom,:,idat_tot,ii)=&
840 : & dijhat_tmp(ind+1:ind+cplex_dij*lmn2_size,&
841 5761800 : 1+(idat_occ-1)*ndij+(idat-1)*ndat_occ*ndij:idat_occ*ndij+(idat-1)*ndat_occ*ndij,ia)
842 : end do
843 : end do ! idat_occ
844 : end do ! idat
845 : end do ! ia
846 10392 : ABI_FREE(dijhat_tmp)
847 20070 : iatm=iatm+gs_ham%nattyp(itypat)
848 : end do
849 :
850 : !*NOTE Build a replica of cwaveocc_prj covering the full ndat*ndat_occ batch: cwaveocc_prj
851 : !*itself does not depend on idat, but gemm_nonlop requires cprjin to be sized on the
852 : !*full batch passed as its "ndat" argument.
853 : !*This is a bit tedious but remains cheaper than calling nonlop ndat times
854 153831 : ABI_MALLOC(cwaveocc_prj_rep,(natom,nspinor*ndat_tot))
855 29034 : call pawcprj_alloc(cwaveocc_prj_rep,cwaveocc_prj(1,1)%ncpgr,cwaveocc_prj(:,1)%nlmn)
856 19356 : do idat=1,ndat
857 : call pawcprj_copy(cwaveocc_prj,&
858 19356 : & cwaveocc_prj_rep(:,(idat-1)*nspinor*ndat_occ+1:idat*nspinor*ndat_occ))
859 : end do
860 :
861 9678 : signs=2; cpopt=2;idir=0; paw_opt=1;nnlout=1;tim_nonlop=17
862 :
863 9678 : if(need_ghc) then
864 7712 : choice=1
865 7712 : call timab(1515,2,tsec) ; call timab(1514,-1,tsec) ; call timab(1546,-2,tsec)
866 : call nonlop(choice,cpopt,cwaveocc_prj_rep,enlout_dum,gs_ham,idir,(/zero/),&
867 : & mpi_enreg,ndat_tot,nnlout,paw_opt,signs,gsc_dum,tim_nonlop,vectin_dum,&
868 : & gvnlxc,enl_ndat=dijhat,&
869 7712 : & select_k=K_H_KPRIME)
870 :
871 7712 : if(gpu_option==ABI_GPU_DISABLED) then
872 15424 : do idat=1,ndat
873 45928 : do idat_occ=1,ndat_occ
874 30504 : idat_tot=(idat-1)*ndat_occ+idat_occ
875 : ghc2(:,1+(idat-1)*npw*nspinor:idat*npw*nspinor) = ghc2(:,1+(idat-1)*npw*nspinor:idat*npw*nspinor)&
876 24180404 : & - gvnlxc(:,1+(idat_tot-1)*npw*nspinor:idat_tot*npw*nspinor) * occ(idat_occ) * wtk
877 : end do ! idat_occ
878 : end do ! idat
879 : else if(gpu_option==ABI_GPU_OPENMP) then
880 : #ifdef HAVE_OPENMP_OFFLOAD
881 : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(2) &
882 : !$OMP& MAP(to:gvnlxc,ghc2,occ) PRIVATE(idat,ipw,idat_occ,idat_tot)
883 : do idat=1,ndat
884 : do ipw=1,npw
885 : do idat_occ=1,ndat_occ
886 : idat_tot=(idat-1)*ndat_occ+idat_occ
887 : ghc2(1,ipw+(idat-1)*npw*nspinor) = ghc2(1,ipw+(idat-1)*npw*nspinor)&
888 : & - gvnlxc(1,ipw+(idat_tot-1)*npw*nspinor) * occ(idat_occ) * wtk
889 : ghc2(2,ipw+(idat-1)*npw*nspinor) = ghc2(2,ipw+(idat-1)*npw*nspinor)&
890 : & - gvnlxc(2,ipw+(idat_tot-1)*npw*nspinor) * occ(idat_occ) * wtk
891 : end do
892 : end do
893 : end do ! idat
894 : #endif
895 : end if
896 7712 : call timab(1514,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1546,-1,tsec)
897 : end if
898 :
899 : ! Forces calculation
900 :
901 9678 : if (fockcommon%optfor.and.(fockcommon%ieigen/=0)) then
902 16824 : ABI_MALLOC(vdotr,(ndat_occ,ndat,3,natom))
903 8412 : ABI_MALLOC(vdoti,(ndat_tot))
904 14020 : ABI_MALLOC(for1,(ndat_occ,3,natom,ndat))
905 : #ifdef HAVE_OPENMP_OFFLOAD
906 : !$OMP TARGET ENTER DATA MAP(alloc:vdotr,vdoti,for1) IF(gpu_option==ABI_GPU_OPENMP)
907 : #endif
908 : !*NOTE:cwavef_rep replicates cwavef ndat_occ times (per idat block) so that single
909 : !*nonlop and dotprod_g_batch_full call can process the whole ndat*ndat_occ batch at once.
910 8412 : ABI_MALLOC(cwavef_rep,(2,npw*nspinor*ndat_tot))
911 : #ifdef HAVE_OPENMP_OFFLOAD
912 : !$OMP TARGET ENTER DATA MAP(alloc:cwavef_rep) IF(gpu_option==ABI_GPU_OPENMP)
913 : #endif
914 2804 : if(gpu_option==ABI_GPU_DISABLED) then
915 5608 : do idat=1,ndat
916 16480 : do idat_occ=1,ndat_occ
917 10872 : idat_tot=(idat-1)*ndat_occ+idat_occ
918 : cwavef_rep(:,1+(idat_tot-1)*npw*nspinor:idat_tot*npw*nspinor)=&
919 17315012 : & cwavef(:,1+(idat-1)*npw*nspinor:idat*npw*nspinor)
920 : end do ! idat_occ
921 : end do ! idat
922 : else if(gpu_option==ABI_GPU_OPENMP) then
923 : #ifdef HAVE_OPENMP_OFFLOAD
924 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:cwavef_rep,cwavef) PRIVATE(idat,idat_occ,idat_tot)
925 : do idat=1,ndat
926 : do idat_occ=1,ndat_occ
927 : idat_tot=(idat-1)*ndat_occ+idat_occ
928 : !$OMP PARALLEL DO PRIVATE(ipw)
929 : do ipw=1,npw*nspinor
930 : cwavef_rep(1,ipw+(idat_tot-1)*npw*nspinor)=cwavef(1,ipw+(idat-1)*npw*nspinor)
931 : cwavef_rep(2,ipw+(idat_tot-1)*npw*nspinor)=cwavef(2,ipw+(idat-1)*npw*nspinor)
932 : end do
933 : end do ! idat_occ
934 : end do ! idat
935 : #endif
936 : end if
937 107292 : choice=2; vdotr=zero;doti=zero;cpopt=4;tim_nonlop=17
938 11216 : do idir=1,3
939 28040 : do iatom=1,natom
940 16824 : call timab(1515,2,tsec) ; call timab(1514,-1,tsec) ; call timab(1546,-2,tsec)
941 : call nonlop(choice,cpopt,cwaveocc_prj_rep,enlout_dum,gs_ham,idir,(/zero/),mpi_enreg,&
942 : & ndat_tot,nnlout,paw_opt,signs,gsc_dum,tim_nonlop,vectin_dum,&
943 : & forout,enl_ndat=dijhat,iatom_only=iatom,&
944 16824 : & select_k=K_H_KPRIME)
945 16824 : call timab(1514,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1546,-1,tsec)
946 : call dotprod_g_batch_full(vdotr(:,:,idir,iatom),vdoti,gs_ham%istwf_k,npw,ndat_tot,2,&
947 25236 : & cwavef_rep,forout,mpi_enreg%me_g0,mpi_enreg%comm_fft,gpu_option=gpu_option)
948 : end do ! iatom
949 : end do ! idir
950 : #ifdef HAVE_OPENMP_OFFLOAD
951 : !$OMP TARGET EXIT DATA MAP(delete:cwavef_rep) IF(gpu_option==ABI_GPU_OPENMP)
952 : #endif
953 2804 : ABI_FREE(cwavef_rep)
954 :
955 2804 : if(gpu_option==ABI_GPU_DISABLED) then
956 5608 : do idat=1,ndat
957 11216 : do iatom=1,natom
958 25236 : do idir=1,3
959 87664 : do idat_occ=1,ndat_occ
960 65232 : esum=0
961 55569312 : do ifft=1,atom_nfgd(iatom)
962 55504080 : ind=atom_ifftsph(ifft,iatom)
963 : esum = esum &
964 : & + vfock(2*ind-1,idat_occ,idat)*grnhat_12(1,ifft,1,idir,iatom,idat_occ,idat) &
965 55569312 : & - vfock(2*ind,idat_occ,idat)*grnhat_12(2,ifft,1,idir,iatom,idat_occ,idat)
966 : end do
967 82056 : for1(idat_occ,idir,iatom,idat)=esum
968 : end do ! idat_occ
969 : end do ! idir
970 : end do ! iatom
971 : end do ! idat
972 : else if(gpu_option==ABI_GPU_OPENMP) then
973 : #ifdef HAVE_OPENMP_OFFLOAD
974 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
975 : !$OMP& MAP(to:vfock,grnhat_12,for1,atom_nfgd,atom_ifftsph) &
976 : !$OMP& PRIVATE(ifft,ind,iatom) PRIVATE(idat_occ,idir,esum)
977 : do idat=1,ndat
978 : do iatom=1,natom
979 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(esum,ind,idir,idat_occ,ifft)
980 : do idir=1,3
981 : do idat_occ=1,ndat_occ
982 : esum=0
983 : do ifft=1,atom_nfgd(iatom)
984 : ind=atom_ifftsph(ifft,iatom)
985 : esum = esum &
986 : & + vfock(2*ind-1,idat_occ,idat)*grnhat_12(1,ifft,1,idir,iatom,idat_occ,idat) &
987 : & - vfock(2*ind,idat_occ,idat)*grnhat_12(2,ifft,1,idir,iatom,idat_occ,idat)
988 : end do
989 : for1(idat_occ,idir,iatom,idat)=esum
990 : end do ! idat_occ
991 : end do ! idir
992 : end do ! iatom
993 : end do ! idat
994 : !$OMP TARGET UPDATE FROM(for1)
995 : #endif
996 : end if
997 :
998 : if(.true.) then
999 : !if(gpu_option==ABI_GPU_DISABLED) then
1000 : #ifdef HAVE_OPENMP_OFFLOAD
1001 : !$OMP TARGET UPDATE FROM(vdotr) IF(gpu_option==ABI_GPU_OPENMP)
1002 : #endif
1003 5608 : do idat=1,ndat
1004 11216 : do iatom=1,natom
1005 25236 : do idir=1,3
1006 87664 : do idat_occ=1,ndat_occ
1007 : for12(idir)=rprimd(1,idir)*for1(idat_occ,1,iatom,idat)&
1008 : & +rprimd(2,idir)*for1(idat_occ,2,iatom,idat)&
1009 65232 : & +rprimd(3,idir)*for1(idat_occ,3,iatom,idat)
1010 : forikpt(idir,iatom,idat)=forikpt(idir,iatom,idat)&
1011 : & - (for12(idir)*gs_ham%ucvol/nfftf+vdotr(idat_occ,idat,idir,iatom))&
1012 82056 : & * occ(idat_occ) * wtk
1013 : end do ! idat_occ
1014 : end do ! idir
1015 : end do ! iatom
1016 : end do ! idat
1017 : end if
1018 : #ifdef HAVE_OPENMP_OFFLOAD
1019 : !$OMP TARGET EXIT DATA MAP(delete:vdotr,vdoti,for1) IF(gpu_option==ABI_GPU_OPENMP)
1020 : #endif
1021 2804 : ABI_FREE(vdotr)
1022 2804 : ABI_FREE(vdoti)
1023 2804 : ABI_FREE(for1)
1024 : end if
1025 :
1026 : ! Stresses calculation
1027 9678 : if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then
1028 1870 : signs=2;choice=3;cpopt=4;tim_nonlop=17
1029 :
1030 : ! first contribution
1031 9350 : ABI_MALLOC(vdotr,(ndat_occ,ndat,6,1))
1032 5610 : ABI_MALLOC(vdoti,(ndat_tot))
1033 : #ifdef HAVE_OPENMP_OFFLOAD
1034 : !$OMP TARGET ENTER DATA MAP(alloc:vdotr,vdoti) IF(gpu_option==ABI_GPU_OPENMP)
1035 : #endif
1036 : !*cwavef_rep replicates cwavef ndat_occ times (per idat block) so that a single
1037 : !*dotprod_g_batch_full call can process the whole ndat*ndat_occ batch at once,
1038 : !*matching the merged nonlop call below.
1039 5610 : ABI_MALLOC(cwavef_rep,(2,npw*nspinor*ndat_tot))
1040 : #ifdef HAVE_OPENMP_OFFLOAD
1041 : !$OMP TARGET ENTER DATA MAP(alloc:cwavef_rep) IF(gpu_option==ABI_GPU_OPENMP)
1042 : #endif
1043 1870 : if(gpu_option==ABI_GPU_DISABLED) then
1044 3740 : do idat=1,ndat
1045 11225 : do idat_occ=1,ndat_occ
1046 7485 : idat_tot=(idat-1)*ndat_occ+idat_occ
1047 : cwavef_rep(:,1+(idat_tot-1)*npw*nspinor:idat_tot*npw*nspinor)=&
1048 4097380 : & cwavef(:,1+(idat-1)*npw*nspinor:idat*npw*nspinor)
1049 : end do ! idat_occ
1050 : end do ! idat
1051 : else if(gpu_option==ABI_GPU_OPENMP) then
1052 : #ifdef HAVE_OPENMP_OFFLOAD
1053 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) MAP(to:cwavef_rep,cwavef) PRIVATE(idat,idat_occ,idat_tot)
1054 : do idat=1,ndat
1055 : do idat_occ=1,ndat_occ
1056 : idat_tot=(idat-1)*ndat_occ+idat_occ
1057 : !$OMP PARALLEL DO PRIVATE(ipw)
1058 : do ipw=1,npw*nspinor
1059 : cwavef_rep(1,ipw+(idat_tot-1)*npw*nspinor)=cwavef(1,ipw+(idat-1)*npw*nspinor)
1060 : cwavef_rep(2,ipw+(idat_tot-1)*npw*nspinor)=cwavef(2,ipw+(idat-1)*npw*nspinor)
1061 : end do
1062 : end do ! idat_occ
1063 : end do ! idat
1064 : #endif
1065 : end if
1066 13090 : do idir=1,6
1067 11220 : call timab(1515,2,tsec) ; call timab(1514,-1,tsec) ; call timab(1546,-2,tsec)
1068 : call nonlop(choice,cpopt,cwaveocc_prj_rep,enlout_dum,gs_ham,idir,(/zero/),mpi_enreg,&
1069 : & ndat_tot,nnlout,paw_opt,signs,gsc_dum,tim_nonlop,vectin_dum,&
1070 11220 : & strout,enl_ndat=dijhat,select_k=K_H_KPRIME)
1071 11220 : call timab(1514,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1546,-1,tsec)
1072 : call dotprod_g_batch_full(vdotr(:,:,idir,1),vdoti,gs_ham%istwf_k,npw,ndat_tot,2,&
1073 : & cwavef_rep,strout,&
1074 11220 : & mpi_enreg%me_g0,mpi_enreg%comm_fft,gpu_option=gpu_option)
1075 : #ifdef HAVE_OPENMP_OFFLOAD
1076 : !$OMP TARGET UPDATE FROM(vdotr) IF(gpu_option==ABI_GPU_OPENMP)
1077 : #endif
1078 24310 : do idat=1,ndat
1079 67350 : do idat_occ=1,ndat_occ
1080 : fockcommon%stress_ikpt(idir,fockcommon%ieigen+idat-1)=fockcommon%stress_ikpt(idir,fockcommon%ieigen+idat-1)-&
1081 56130 : & vdotr(idat_occ,idat,idir,1)*occ(idat_occ)*wtk/gs_ham%ucvol
1082 : end do ! idat_occ
1083 : end do ! idat
1084 : end do ! idir
1085 : #ifdef HAVE_OPENMP_OFFLOAD
1086 : !$OMP TARGET EXIT DATA MAP(delete:vdotr,vdoti,cwavef_rep) IF(gpu_option==ABI_GPU_OPENMP)
1087 : #endif
1088 1870 : ABI_FREE(vdotr)
1089 1870 : ABI_FREE(vdoti)
1090 1870 : ABI_FREE(cwavef_rep)
1091 :
1092 1870 : stress_ikpt => fockcommon%stress_ikpt
1093 1870 : ieigen = fockcommon%ieigen
1094 :
1095 : ! second contribution
1096 : !if(.true.) then
1097 1870 : if(gpu_option==ABI_GPU_DISABLED) then
1098 : #ifdef HAVE_OPENMP_OFFLOAD
1099 : !$OMP TARGET UPDATE FROM(vfock,grnhat_12) IF(gpu_option==ABI_GPU_OPENMP)
1100 : #endif
1101 7480 : ABI_MALLOC(strdat, (3,3,ndat_occ,ndat))
1102 101045 : strdat=zero
1103 3740 : do idat=1,ndat
1104 11225 : do idat_occ=1,ndat_occ
1105 31810 : do idir=1,3
1106 97305 : do idir1=1,3
1107 67365 : esum=0
1108 202095 : do iatom=1,natom
1109 104707395 : do ifft=1,atom_nfgd(iatom)
1110 104505300 : ind = atom_ifftsph(ifft,iatom)
1111 : esum = esum + (vfock(2*ind-1,idat_occ,idat)*grnhat_12(1,ifft,1,idir,iatom,idat_occ,idat)&
1112 : & - vfock(2*ind,idat_occ,idat)*grnhat_12(2,ifft,1,idir,iatom,idat_occ,idat))&
1113 104640030 : & * atom_rfgd(idir1,ifft,iatom)
1114 : end do
1115 : end do
1116 89820 : strdat(idir,idir1,idat_occ,idat)=esum
1117 : end do
1118 : end do
1119 : end do ! idat_occ
1120 : end do ! idat
1121 3740 : do idat=1,ndat
1122 11225 : do idat_occ=1,ndat_occ
1123 29940 : do idir=1,3
1124 29940 : fockstr(idir)=strdat(idir,idir,idat_occ,idat)
1125 : end do
1126 7485 : fockstr(4)=(strdat(3,2,idat_occ,idat)+strdat(2,3,idat_occ,idat))*half
1127 7485 : fockstr(5)=(strdat(3,1,idat_occ,idat)+strdat(1,3,idat_occ,idat))*half
1128 7485 : fockstr(6)=(strdat(1,2,idat_occ,idat)+strdat(2,1,idat_occ,idat))*half
1129 54265 : do idir=1,6
1130 : fockcommon%stress_ikpt(idir,fockcommon%ieigen+idat-1)=fockcommon%stress_ikpt(idir,fockcommon%ieigen+idat-1)+&
1131 52395 : & fockstr(idir)/nfftf*occ(idat_occ)*wtk
1132 : end do
1133 : end do ! idat_occ
1134 : end do ! idat
1135 1870 : ABI_FREE(strdat)
1136 : else if(gpu_option==ABI_GPU_OPENMP) then
1137 : #ifdef HAVE_OPENMP_OFFLOAD
1138 : ABI_MALLOC(strdat, (3,3,ndat_occ,ndat))
1139 : strdat=zero
1140 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(4) MAP(tofrom:strdat) &
1141 : !$OMP& MAP(to:vfock,grnhat_12,atom_nfgd,atom_rfgd,atom_ifftsph) &
1142 : !$OMP& PRIVATE(idat,idat_occ,idir,idir1,esum,iatom)
1143 : do idat=1,ndat
1144 : do idat_occ=1,ndat_occ
1145 : do idir=1,3
1146 : do idir1=1,3
1147 : do iatom=1,natom
1148 : esum=0
1149 : !$OMP PARALLEL DO PRIVATE(ifft,ind) REDUCTION(+:esum)
1150 : do ifft=1,atom_nfgd(iatom)
1151 : ind=atom_ifftsph(ifft,iatom)
1152 : esum = esum + (vfock(2*ind-1,idat_occ,idat)*grnhat_12(1,ifft,1,idir,iatom,idat_occ,idat)&
1153 : & - vfock(2*ind,idat_occ,idat)*grnhat_12(2,ifft,1,idir,iatom,idat_occ,idat))&
1154 : & * atom_rfgd(idir1,ifft,iatom)
1155 : end do
1156 : strdat(idir,idir1,idat_occ,idat)=strdat(idir,idir1,idat_occ,idat)+esum
1157 : end do
1158 : end do
1159 : end do
1160 : end do ! idat_occ
1161 : end do ! idat
1162 : do idat=1,ndat
1163 : do idat_occ=1,ndat_occ
1164 : do idir=1,3
1165 : fockstr(idir)=strdat(idir,idir,idat_occ,idat)
1166 : end do
1167 : fockstr(4)=(strdat(3,2,idat_occ,idat)+strdat(2,3,idat_occ,idat))*half
1168 : fockstr(5)=(strdat(3,1,idat_occ,idat)+strdat(1,3,idat_occ,idat))*half
1169 : fockstr(6)=(strdat(1,2,idat_occ,idat)+strdat(2,1,idat_occ,idat))*half
1170 : do idir=1,6
1171 : stress_ikpt(idir,ieigen+idat-1)=stress_ikpt(idir,ieigen+idat-1)+&
1172 : & fockstr(idir)/nfftf*occ(idat_occ)*wtk
1173 : end do
1174 : end do ! idat_occ
1175 : end do ! idat
1176 : ABI_FREE(strdat)
1177 : #endif
1178 : end if
1179 :
1180 : ! third contribution
1181 1870 : if(gpu_option==ABI_GPU_DISABLED) then
1182 3740 : do idat=1,ndat
1183 11225 : do idat_occ=1,ndat_occ
1184 7485 : doti=zero
1185 22455 : do iatom=1,natom
1186 11634155 : do ifft=1,atom_nfgd(iatom)
1187 11611700 : ind = atom_ifftsph(ifft,iatom)
1188 : doti = doti &
1189 : & + vfock(2*ind-1,idat_occ,idat) * rho12(1,ifft,nspinor,idat_occ,idat,iatom)&
1190 11626670 : & - vfock(2*ind,idat_occ,idat) * rho12(2,ifft,nspinor,idat_occ,idat,iatom)
1191 : end do
1192 : end do
1193 31810 : fockcommon%stress_ikpt(1:3,fockcommon%ieigen+idat-1)=fockcommon%stress_ikpt(1:3,fockcommon%ieigen+idat-1)-doti/nfftf*occ(idat_occ)*wtk
1194 : end do ! idat_occ
1195 : end do ! idat
1196 : else if(gpu_option==ABI_GPU_OPENMP) then
1197 : #ifdef HAVE_OPENMP_OFFLOAD
1198 : ABI_MALLOC(strdat, (1,1,ndat_occ,ndat))
1199 : strdat=zero
1200 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(2) &
1201 : !$OMP& MAP(tofrom:strdat) MAP(to:vfock,rho12,atom_nfgd,atom_ifftsph) &
1202 : !$OMP& PRIVATE(idat,idat_occ,idir,idir1,doti)
1203 : do idat=1,ndat
1204 : do idat_occ=1,ndat_occ
1205 : doti=zero
1206 : !$OMP PARALLEL DO PRIVATE(iatom,ifft,ind) REDUCTION(+:doti)
1207 : do iatom=1,natom
1208 : do ifft=1,atom_nfgd(iatom)
1209 : ind = atom_ifftsph(ifft,iatom)
1210 : doti = doti &
1211 : & + vfock(2*ind-1,idat_occ,idat) * rho12(1,ifft,nspinor,idat_occ,idat,iatom)&
1212 : & - vfock(2*ind,idat_occ,idat) * rho12(2,ifft,nspinor,idat_occ,idat,iatom)
1213 : end do
1214 : end do
1215 : strdat(1,1,idat_occ,idat)=doti
1216 : end do ! idat_occ
1217 : end do ! idat
1218 : do idat=1,ndat
1219 : do idat_occ=1,ndat_occ
1220 : fockcommon%stress_ikpt(1:3,fockcommon%ieigen+idat-1)=fockcommon%stress_ikpt(1:3,fockcommon%ieigen+idat-1) &
1221 : -strdat(1,1,idat_occ,idat)/nfftf*occ(idat_occ)*wtk
1222 : end do ! idat_occ
1223 : end do ! idat
1224 : ABI_FREE(strdat)
1225 : #endif
1226 : end if
1227 : end if ! end stresses
1228 :
1229 9678 : call pawcprj_free(cwaveocc_prj_rep)
1230 9678 : ABI_FREE(cwaveocc_prj_rep)
1231 9678 : ABI_FREE(dijhat)
1232 : end if !end PAW
1233 191436 : call timab(1526,2,tsec) ; call timab(1546,-2,tsec)
1234 :
1235 : ! =============================================================
1236 : ! === Apply the local potential vfockloc_munu to cwaveocc_r ===
1237 : ! =============================================================
1238 191436 : call timab(1527,-1,tsec)
1239 191436 : if(gpu_option==ABI_GPU_DISABLED) then
1240 : ABI_UNUSED(esumi)
1241 : !$OMP PARALLEL DO &
1242 : !$OMP& PRIVATE(ind,recwocc,imcwocc,revloc,imvloc)
1243 382872 : do idat=1,ndat
1244 1242043 : do idat_occ=1,ndat_occ
1245 14648932 : do i3=1,ngfftf(3)
1246 235971433 : do i2=1,ngfftf(2)
1247 4097411292 : do i1=1,ngfftf(1)
1248 3862299030 : ind=i1+(i2-1)*ngfftf(1)+(i3-1)*ngfftf(2)*ngfftf(1)
1249 3862299030 : revloc=vfock(2*ind-1,idat_occ,idat) ; imvloc=vfock(2*ind,idat_occ,idat)
1250 3862299030 : recwocc=cwaveocc_r(1,i1,i2,i3,idat_occ)
1251 3862299030 : imcwocc=cwaveocc_r(2,i1,i2,i3,idat_occ)
1252 3862299030 : vlocpsi_r(2*ind-1,idat)=vlocpsi_r(2*ind-1,idat)-(revloc*recwocc-imvloc*imcwocc)*occ(idat_occ)*wtk
1253 4083812967 : vlocpsi_r(2*ind ,idat)=vlocpsi_r(2*ind ,idat)-(revloc*imcwocc+imvloc*recwocc)*occ(idat_occ)*wtk
1254 : end do
1255 : end do
1256 : end do
1257 : end do ! idat_occ
1258 : end do ! idat
1259 : else if(gpu_option==ABI_GPU_OPENMP) then
1260 : #ifdef HAVE_OPENMP_OFFLOAD
1261 : !$OMP TARGET TEAMS DISTRIBUTE COLLAPSE(3) &
1262 : !$OMP& MAP(to:vlocpsi_r,cwaveocc_r,occ,vfock) PRIVATE(idat,i3,i2)
1263 : do idat=1,ndat
1264 : do i3=1,ngfftf(3)
1265 : do i2=1,ngfftf(2)
1266 : !$OMP PARALLEL DO PRIVATE(ind,i1,idat_occ,recwocc,imcwocc,revloc,imvloc,esum,esumi)
1267 : do i1=1,ngfftf(1)
1268 : ind=i1+(i2-1)*ngfftf(1)+(i3-1)*ngfftf(2)*ngfftf(1)
1269 : esum=zero ; esumi=zero
1270 : do idat_occ=1,ndat_occ
1271 : revloc=vfock(2*ind-1,idat_occ,idat) ; imvloc=vfock(2*ind,idat_occ,idat)
1272 : recwocc=cwaveocc_r(1,i1,i2,i3,idat_occ)
1273 : imcwocc=cwaveocc_r(2,i1,i2,i3,idat_occ)
1274 : esum=esum+(revloc*recwocc-imvloc*imcwocc)*occ(idat_occ)
1275 : esumi=esumi+(revloc*imcwocc+imvloc*recwocc)*occ(idat_occ)
1276 : end do
1277 : vlocpsi_r(2*ind-1,idat)=vlocpsi_r(2*ind-1,idat)-esum*wtk
1278 : vlocpsi_r(2*ind ,idat)=vlocpsi_r(2*ind ,idat)-esumi*wtk
1279 : end do
1280 : end do
1281 : end do
1282 : end do ! idat
1283 : #endif
1284 : end if
1285 191436 : if (allocated(fockbz%cgocc)) then
1286 : #ifdef HAVE_OPENMP_OFFLOAD
1287 : !$OMP TARGET EXIT DATA MAP(delete:cwaveocc_r) IF(gpu_option==ABI_GPU_OPENMP)
1288 : #endif
1289 191436 : ABI_FREE(cwaveocc_r)
1290 : end if
1291 393962 : call timab(1527,2,tsec)
1292 : end do ! jband
1293 :
1294 : ! ========================================================
1295 : ! === End of loop : update of shifts and deallocations ===
1296 : ! ==============================:=========================
1297 : !* Update of the shifts to be applied (reminder : mkmem is not 0, nspinor=1)
1298 202526 : call timab(1528,1,tsec)
1299 202526 : jcg=jcg+npwj*nband_k
1300 202526 : jbg=jbg+nband_k
1301 202526 : bdtot_jindex=bdtot_jindex+nband_k
1302 202526 : jkg=jkg+npwj
1303 202526 : if (fockcommon%usepaw==1) then
1304 9678 : ABI_FREE(enlout_dum)
1305 9678 : ABI_FREE(vectin_dum)
1306 9678 : ABI_FREE(ffnl_kp_dum)
1307 9678 : ABI_FREE(kpg_kp)
1308 : end if
1309 202526 : if (associated(gs_ham%ph3d_kp)) then
1310 202526 : ABI_FREE(gs_ham%ph3d_kp)
1311 : end if
1312 : #ifdef HAVE_OPENMP_OFFLOAD
1313 : !$OMP TARGET EXIT DATA MAP(delete:rhor_munu,rhog_munu,vfock,occ,rho12,gvnlxc) IF(gpu_option==ABI_GPU_OPENMP)
1314 : #endif
1315 202526 : ABI_FREE(rhor_munu)
1316 202526 : ABI_FREE(rhog_munu)
1317 202526 : ABI_FREE(vfock)
1318 202526 : ABI_FREE(occ)
1319 : !*Additional arrays in case of paw
1320 202526 : if (fockcommon%usepaw==1) then
1321 9678 : if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then
1322 : #ifdef HAVE_OPENMP_OFFLOAD
1323 : !$OMP TARGET EXIT DATA MAP(delete:forout) IF(gpu_option==ABI_GPU_OPENMP)
1324 : #endif
1325 2804 : ABI_FREE(forout)
1326 : end if
1327 9678 : if (fockcommon%optstr.and.(fockcommon%ieigen/=0)) then
1328 : #ifdef HAVE_OPENMP_OFFLOAD
1329 : !$OMP TARGET EXIT DATA MAP(delete:strout) IF(gpu_option==ABI_GPU_OPENMP)
1330 : #endif
1331 1870 : ABI_FREE(strout)
1332 : end if
1333 : #ifdef HAVE_OPENMP_OFFLOAD
1334 : !$OMP TARGET EXIT DATA MAP(delete:grnhat_12) IF(gpu_option==ABI_GPU_OPENMP .and. ider==3)
1335 : !$OMP TARGET EXIT DATA MAP(delete:rho12,gvnlxc) IF(gpu_option==ABI_GPU_OPENMP)
1336 : !$OMP TARGET EXIT DATA MAP(delete:atom_ifftsph,atom_nfgd,atom_rfgd) IF(gpu_option==ABI_GPU_OPENMP)
1337 : #endif
1338 9678 : ABI_FREE(grnhat_12)
1339 9678 : ABI_FREE(gvnlxc)
1340 9678 : ABI_FREE(grnhat12)
1341 9678 : ABI_FREE(rho12)
1342 9678 : ABI_FREE(atom_ifftsph)
1343 9678 : ABI_FREE(atom_nfgd)
1344 9678 : ABI_FREE(atom_rfgd)
1345 : end if
1346 :
1347 231529 : call timab(1528,2,tsec)
1348 :
1349 : end do ! jkpt
1350 :
1351 : ! ========================================================
1352 : ! === After loop ===
1353 : ! ========================================================
1354 :
1355 29003 : call timab(1506,2,tsec) ; call timab(1507,1,tsec) ; call timab(1547,-1,tsec)
1356 :
1357 29003 : if (fockcommon%usepaw==1) then
1358 2061 : if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then
1359 1052 : call timab(1547,2,tsec) ; call timab(1548,-1,tsec)
1360 1052 : call xmpi_sum(forikpt,mpi_enreg%comm_hf,ier)
1361 1052 : call timab(1548,2,tsec) ; call timab(1547,-1,tsec)
1362 2104 : do idat=1,ndat
1363 4208 : do iatom=1,natom !Loop over atom
1364 2104 : ia=gs_ham%atindx(iatom)
1365 15780 : fockcommon%forces_ikpt(:,ia,fockcommon%ieigen+idat-1)=forikpt(:,iatom,idat)
1366 : end do
1367 : end do
1368 : end if
1369 : end if
1370 29003 : if(fockcommon%optstr.and.(fockcommon%ieigen/=0)) then
1371 2040 : call timab(1547,2,tsec) ; call timab(1548,-1,tsec)
1372 2040 : call xmpi_sum(fockcommon%stress_ikpt,mpi_enreg%comm_hf,ier)
1373 2040 : call timab(1548,2,tsec) ; call timab(1547,-1,tsec)
1374 : end if
1375 :
1376 29003 : if (.not.need_ghc) then
1377 :
1378 : ! ===============================
1379 : ! === Deallocate local arrays ===
1380 : ! ===============================
1381 : #ifdef HAVE_OPENMP_OFFLOAD
1382 : !$OMP TARGET EXIT DATA MAP(delete:cwavef_r,vqg,vlocpsi_r) IF(gpu_option==ABI_GPU_OPENMP)
1383 : #endif
1384 2052 : ABI_FREE(cwavef_r)
1385 2052 : ABI_FREE(vlocpsi_r)
1386 2052 : ABI_FREE(dummytab)
1387 2052 : ABI_FREE(vqg)
1388 2052 : if(fockcommon%usepaw==1.or.fockcommon%optstr) then
1389 2052 : ABI_FREE(gboundf)
1390 : end if
1391 2052 : if (fockcommon%usepaw==1) then
1392 214 : if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then
1393 136 : ABI_FREE(forikpt)
1394 : end if
1395 : end if
1396 : !*Restore gs_ham datastructure
1397 :
1398 2052 : if (associated(gs_ham%ph3d_kp)) then
1399 0 : ABI_MALLOC(gs_ham%ph3d_kp,(2,gs_ham%npw_k,gs_ham%matblk))
1400 : end if
1401 : call gs_ham%load_kprime(kpt_kp=gs_ham%kpt_k,istwf_kp=gs_ham%istwf_k,&
1402 2052 : & npw_kp=gs_ham%npw_k,kg_kp=gs_ham%kg_k,ffnl_kp=gs_ham%ffnl_k,ph3d_kp=gs_ham%ph3d_k)
1403 :
1404 : ! if (fockcommon%ieigen/=0) fockcommon%ieigen=0
1405 :
1406 : else
1407 :
1408 : ! *Restore gs_ham datastructure
1409 :
1410 26951 : if (associated(gs_ham%ph3d_kp)) then
1411 0 : ABI_MALLOC(gs_ham%ph3d_kp,(2,gs_ham%npw_k,gs_ham%matblk))
1412 : end if
1413 : call gs_ham%load_kprime(kpt_kp=gs_ham%kpt_k,istwf_kp=gs_ham%istwf_k,&
1414 26951 : & npw_kp=gs_ham%npw_k,kg_kp=gs_ham%kg_k,ffnl_kp=gs_ham%ffnl_k,ph3d_kp=gs_ham%ph3d_k)
1415 :
1416 : ! * Perform an FFT using fourwf to get ghc1 = FFT^-1(vlocpsi_r)
1417 : #ifdef HAVE_OPENMP_OFFLOAD
1418 : !$OMP TARGET UPDATE FROM(vlocpsi_r) IF(gpu_option==ABI_GPU_OPENMP)
1419 : #endif
1420 134755 : ABI_MALLOC(psilocal,(cplex_fock*n4f,n5f,n6f*ndat))
1421 53902 : do idat=1,ndat
1422 : call fftpac(1,mpi_enreg,nspden_fock,cplex_fock*n1f,n2f,n3f,&
1423 53902 : cplex_fock*n4f,n5f,n6f,ngfft,vlocpsi_r(:,idat),psilocal(:,:,1+(idat-1)*n6f:idat*n6f),2)
1424 : end do ! idat
1425 :
1426 26951 : call timab(1515,2,tsec) ; call timab(1512,-1,tsec) ; call timab(1547,-2,tsec)
1427 : call fourwf(0,rhodum0,rhodum,ghc1,psilocal,gboundf,gboundf,gs_ham%istwf_k,gs_ham%kg_k,gs_ham%kg_k,&
1428 : & mgfftf,mpi_enreg,ndat,ngfftf,1,npw,n4f,n5f,n6f,3,tim_fourwf_fock_getghc,weight1,weight1,&
1429 26951 : & gpu_option=gs_ham%gpu_option)
1430 26951 : call timab(1512,2,tsec) ; call timab(1515,-1,tsec) ; call timab(1547,-1,tsec)
1431 26951 : ABI_FREE(psilocal)
1432 :
1433 26951 : if(gpu_option==ABI_GPU_DISABLED) then
1434 20075920 : ghc1=ghc1*sqrt(gs_ham%ucvol)+ghc2
1435 0 : else if(gpu_option==ABI_GPU_OPENMP) then
1436 0 : call abi_xaxpy(npw*ndat,cucvol,ghc1,1,ghc2,1,x_cplx=2,gpu_option=gpu_option)
1437 : #ifdef HAVE_OPENMP_OFFLOAD
1438 : call gpu_copy(ghc1,ghc2,int(2,c_size_t)*npw*ndat)
1439 : !$OMP TARGET UPDATE FROM(ghc1)
1440 : #endif
1441 : end if
1442 :
1443 : ! * If the calculation is parallelized, perform an MPI_allreduce to sum all the contributions in the array ghc
1444 20048969 : ghc(:,:)=ghc(:,:)/mpi_enreg%nproc_hf + ghc1(:,:)
1445 :
1446 26951 : call timab(1547,2,tsec) ; call timab(1548,-1,tsec)
1447 26951 : call xmpi_sum(ghc,mpi_enreg%comm_hf,ier)
1448 26951 : call timab(1548,2,tsec) ; call timab(1547,-1,tsec)
1449 :
1450 : ! ===============================
1451 : ! === Deallocate local PAW arrays ===
1452 : ! ===============================
1453 26951 : if (fockcommon%usepaw==1) then
1454 1847 : if ((fockcommon%optfor).and.(fockcommon%ieigen/=0)) then
1455 916 : ABI_FREE(forikpt)
1456 : end if
1457 : end if
1458 26951 : if(fockcommon%usepaw==1.or.fockcommon%optstr) then
1459 1847 : ABI_FREE(gboundf)
1460 : end if
1461 : ! ============================================
1462 : ! === Calculate the contribution to energy ===
1463 : ! ============================================
1464 : ! * Only the contribution when cwavef=cgocc_bz are calculated, in order to cancel exactly the self-interaction
1465 : ! * at each convergence step. (consistent definition with the definition of hartree energy)
1466 26951 : if (fockcommon%ieigen/=0) then
1467 53902 : do idat=1,ndat
1468 : eigen=zero
1469 : ! * Dot product of cwavef and ghc
1470 : ! * inspired from the routine 54_spacepar/meanvalue_g but without the reference to parallelism and filtering
1471 26951 : if(gs_ham%istwf_k==2) then
1472 6 : eigen=half*cwavef(1,1+(idat-1)*npw)*ghc1(1,1+(idat-1)*npw)
1473 : else
1474 26945 : eigen=cwavef(1,1+(idat-1)*npw)*ghc1(1,1+(idat-1)*npw)+cwavef(2,1+(idat-1)*npw)*ghc1(2,1+(idat-1)*npw)
1475 : end if
1476 6674006 : do ipw=2,npw
1477 6674006 : eigen=eigen+cwavef(1,ipw+(idat-1)*npw)*ghc1(1,ipw+(idat-1)*npw)+cwavef(2,ipw+(idat-1)*npw)*ghc1(2,ipw+(idat-1)*npw)
1478 : end do
1479 26951 : if(gs_ham%istwf_k>=2) eigen=two*eigen
1480 26951 : call timab(1547,2,tsec) ; call timab(1548,-1,tsec)
1481 26951 : call xmpi_sum(eigen,mpi_enreg%comm_hf,ier)
1482 26951 : call timab(1548,2,tsec) ; call timab(1547,-1,tsec)
1483 26951 : fockcommon%eigen_ikpt(fockcommon%ieigen+idat-1)= eigen
1484 80853 : if(fockcommon%use_ACE==0) fockcommon%ieigen = 0
1485 : end do ! idat
1486 : end if
1487 :
1488 : ! ===============================
1489 : ! === Deallocate local arrays ===
1490 : ! ===============================
1491 : #ifdef HAVE_OPENMP_OFFLOAD
1492 : !$OMP TARGET EXIT DATA MAP(delete:cwavef_r,ghc1,ghc2,vqg,vlocpsi_r) IF(gpu_option==ABI_GPU_OPENMP)
1493 : #endif
1494 26951 : ABI_FREE(cwavef_r)
1495 26951 : ABI_FREE(ghc1)
1496 26951 : ABI_FREE(ghc2)
1497 26951 : ABI_FREE(vlocpsi_r)
1498 26951 : ABI_FREE(dummytab)
1499 26951 : ABI_FREE(vqg)
1500 :
1501 : endif
1502 :
1503 29003 : call timab(1504,2,tsec) ; call timab(1507,-2,tsec) ; call timab(1515,-2,tsec) ; call timab(1547,-2,tsec)
1504 :
1505 58006 : end subroutine fock_getghc
1506 : !!***
1507 :
1508 : !!****f* ABINIT/fock2ACE
1509 : !! NAME
1510 : !! fock2ACE
1511 : !!
1512 : !! FUNCTION
1513 : !! Compute nonlocal contribution to the Fock part of the hamiltonian in the ACE formalism.
1514 : !! optionally contribution to Fock forces
1515 : !!
1516 : !! INPUTS
1517 : !! cg(2,mcg)=wavefunctions (may be read from disk file)
1518 : !! cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each NL proj |p_lmn>
1519 : !! fock <type(fock_type)>= quantities to calculate Fock exact exchange
1520 : !! istwfk(nkpt)=input option parameter that describes the storage of wfs
1521 : !! kg(3,mpw*mkmem)=reduced coordinates (integers) of G vecs in basis
1522 : !! kpt(3,nkpt)=k points in reduced coordinates
1523 : !! mband=maximum number of bands
1524 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
1525 : !! mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
1526 : !! mgfft=maximum size of 1D FFTs
1527 : !! mkmem=number of k points treated by this node.
1528 : !! mpi_enreg=information about MPI parallelization
1529 : !! mpsang=
1530 : !! mpw= maximum number of plane waves
1531 : !! my_natom=number of atoms treated by current processor
1532 : !! natom=number of atoms in cell.
1533 : !! nband(nkpt)=number of bands at each k point
1534 : !! nfft=number of FFT grid points
1535 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1536 : !! nkpt=number of k points in Brillouin zone
1537 : !! nloalg(3)=governs the choice of the algorithm for non-local operator.
1538 : !! npwarr(nkpt)=number of planewaves in basis and boundary at each k
1539 : !! nspden=Number of spin Density components
1540 : !! nspinor=number of spinorial components of the wavefunctions
1541 : !! nsppol=1 for unpolarized, 2 for spin-polarized
1542 : !! ntypat=number of types of atoms
1543 : !! occ(mband*nkpt*nsppol)=occupation numbers for each band over all k points
1544 : !! optfor=1 if computation of forces is required
1545 : !! paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
1546 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
1547 : !! ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
1548 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
1549 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
1550 : !! typat(natom)=type of each atom
1551 : !! usecprj=1 if cprj datastructure has been allocated
1552 : !! gpu_option= GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
1553 : !! wtk(nkpt)=weight associated with each k point
1554 : !! xred(3,natom)=reduced dimensionless atomic coordinates
1555 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
1556 : !!
1557 : !! OUTPUT
1558 : !!
1559 : !! fock%fockACE(ikpt,isppol)%xi
1560 : !! if optfor=1, fock%fock_common%forces
1561 : !!
1562 : !! SOURCE
1563 :
1564 1238 : subroutine fock2ACE(cg,cprj,fock,istwfk,kg,kpt,mband,mcg,mcprj,mgfft,mkmem,mpi_enreg,mpsang,&
1565 1238 : & mpw,my_natom,natom,nband,nfft,ngfft,nkpt,nloalg,npwarr,nspden,nspinor,nsppol,&
1566 1238 : & ntypat,occ,optfor,paw_ij,pawtab,ph1d,psps,rprimd,typat,usecprj,gpu_option,wtk,xred,ylm)
1567 :
1568 : !Arguments ------------------------------------
1569 : !scalars
1570 : integer,intent(in) :: mband,mcg,mcprj,mgfft,mkmem,mpsang,mpw,my_natom,natom,nfft,nkpt
1571 : integer,intent(in) :: nspden,nsppol,nspinor,ntypat,optfor
1572 : integer,intent(in) :: usecprj,gpu_option
1573 : type(MPI_type),intent(inout) :: mpi_enreg
1574 : type(pseudopotential_type),intent(in) :: psps
1575 : !arrays
1576 : integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),nband(nkpt*nsppol)
1577 : integer,intent(in) :: ngfft(18),nloalg(3),npwarr(nkpt)
1578 : integer,intent(in) :: typat(natom)
1579 : real(dp),intent(in) :: cg(2,mcg)
1580 : real(dp),intent(in) :: kpt(3,nkpt)
1581 : real(dp),intent(in) :: occ(mband*nkpt*nsppol),ph1d(2,3*(2*mgfft+1)*natom)
1582 : real(dp),intent(in) :: rprimd(3,3),wtk(nkpt),xred(3,natom)
1583 : real(dp),intent(in) :: ylm(mpw*mkmem,mpsang*mpsang*psps%useylm)
1584 : type(pawcprj_type),intent(inout) :: cprj(natom,mcprj*usecprj)
1585 : type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw)
1586 : type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
1587 : type(fock_type),pointer, intent(inout) :: fock
1588 : !Local variables-------------------------------
1589 : !scalars
1590 : integer :: bandpp,bdtot_index,dimffnl,iband,iband_cprj,iband_last,ibg,icg,ider
1591 : integer :: ierr,info,idir,ikg,ikpt,ilm,ipw,isppol,istwf_k,kk,ll
1592 : integer :: mband_cprj,me_distrb,my_ikpt,my_nspinor,nband_k,nband_cprj_k,nkpg
1593 : integer :: npw_k,spaceComm
1594 : integer :: use_ACE_old
1595 : integer :: blocksize,iblock,jblock,iblocksize,jblocksize,nblockbd
1596 : !integer, save :: counter=0
1597 1238 : type(gs_hamiltonian_type) :: gs_hamk
1598 : logical :: compute_gbound
1599 : character(len=500) :: msg
1600 : type(fock_common_type),pointer :: fockcommon
1601 : !arrays
1602 1238 : integer,allocatable :: kg_k(:,:)
1603 : real(dp) :: kpoint(3),rmet(3,3),tsec(2)
1604 1238 : real(dp),allocatable :: bb(:,:,:),cwavef(:,:),cwavefk(:,:),ffnl_sav(:,:,:,:)
1605 1238 : real(dp),allocatable :: kpg_k(:,:),kpg_k_sav(:,:)
1606 1238 : real(dp),allocatable :: mkl(:,:,:),occblock(:),ph3d(:,:,:),ph3d_sav(:,:,:)
1607 1238 : real(dp),allocatable :: wi(:,:,:),weight(:),ylm_k(:,:),ylmgr_k(:,:,:)
1608 1238 : real(dp),allocatable,target :: ffnl(:,:,:,:)
1609 : type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null()
1610 1238 : type(pawcprj_type),target,allocatable :: cwaveprj(:,:)
1611 :
1612 : !*************************************************************************
1613 :
1614 1238 : call timab(1560,1,tsec)
1615 1238 : call timab(1561,1,tsec)
1616 :
1617 : !DEBUG
1618 : !if(counter>0)return
1619 : !counter=counter+1
1620 : !ENDDEBUG
1621 :
1622 : !Init mpicomm and me
1623 1238 : if(mpi_enreg%paral_kgb==1)then
1624 0 : spaceComm=mpi_enreg%comm_kpt
1625 0 : me_distrb=mpi_enreg%me_kpt
1626 : else
1627 : !* In case of HF calculation
1628 1238 : if (mpi_enreg%paral_hf==1) then
1629 0 : spaceComm=mpi_enreg%comm_kpt
1630 0 : me_distrb=mpi_enreg%me_kpt
1631 : else
1632 1238 : spaceComm=mpi_enreg%comm_cell
1633 1238 : me_distrb=mpi_enreg%me_cell
1634 : end if
1635 : end if
1636 :
1637 : !Some initializations
1638 1238 : my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
1639 1238 : compute_gbound=.true.
1640 1238 : fockcommon => fock%fock_common
1641 1238 : use_ACE_old=fockcommon%use_ACE
1642 1238 : fockcommon%use_ACE=0
1643 1238 : fockcommon%e_fock0=zero
1644 :
1645 : !Initialize Hamiltonian (k- and spin-independent terms)
1646 :
1647 : call gs_hamk%init(psps,pawtab,nspinor,nsppol,nspden,natom,&
1648 : & typat,xred,nfft,mgfft,ngfft,rprimd,nloalg,usecprj=usecprj,&
1649 : & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
1650 : & paw_ij=paw_ij,ph1d=ph1d,fock=fock,&
1651 1238 : & gpu_option=gpu_option)
1652 49520 : rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
1653 1238 : fockcommon%use_ACE=use_ACE_old
1654 :
1655 : !need to reorder cprj=<p_lmn|Cnk> (from unsorted to atom-sorted)
1656 1238 : if (psps%usepaw==1) then
1657 204 : call pawcprj_reorder(cprj,gs_hamk%atindx)
1658 : end if
1659 :
1660 : !LOOP OVER SPINS
1661 1238 : bdtot_index=0;ibg=0;icg=0
1662 1238 : call timab(1561,2,tsec) ; call timab(1562,-1,tsec)
1663 :
1664 2611 : do isppol=1,nsppol
1665 1373 : fockcommon%isppol=isppol
1666 : ! Continue to initialize the Hamiltonian (PAW DIJ coefficients)
1667 1373 : call gs_hamk%load_spin(isppol,with_nonlocal=.true.)
1668 :
1669 : ! Loop over k points
1670 1373 : ikg=0
1671 9909 : do ikpt=1,nkpt
1672 7298 : fockcommon%ikpt=ikpt
1673 7298 : nband_k=nband(ikpt+(isppol-1)*nkpt)
1674 7298 : npw_k=npwarr(ikpt)
1675 29192 : kpoint(:)=kpt(:,ikpt)
1676 7298 : istwf_k=istwfk(ikpt)
1677 7298 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
1678 1720 : bdtot_index=bdtot_index+nband_k
1679 1720 : cycle
1680 : end if
1681 :
1682 : ! Parallelism over FFT and/or bands: define sizes and tabs
1683 5578 : if (mpi_enreg%paral_kgb==1) then
1684 0 : my_ikpt=mpi_enreg%my_kpttab(ikpt)
1685 0 : nblockbd=nband_k/(mpi_enreg%nproc_band*mpi_enreg%bandpp)
1686 0 : bandpp=mpi_enreg%bandpp
1687 0 : my_bandfft_kpt => bandfft_kpt(my_ikpt)
1688 : else
1689 5578 : my_ikpt=ikpt
1690 5578 : bandpp=mpi_enreg%bandpp
1691 5578 : nblockbd=nband_k/bandpp
1692 : end if
1693 5578 : blocksize=nband_k/nblockbd
1694 5578 : mband_cprj=mband/mpi_enreg%nproc_band
1695 5578 : nband_cprj_k=nband_k/mpi_enreg%nproc_band
1696 :
1697 16734 : ABI_MALLOC(cwavef,(2,npw_k*my_nspinor*blocksize))
1698 5578 : if (psps%usepaw==1) then
1699 2912 : ABI_MALLOC(cwaveprj,(natom,my_nspinor*bandpp))
1700 416 : call pawcprj_alloc(cwaveprj,0,gs_hamk%dimcprj)
1701 : else
1702 5162 : ABI_MALLOC(cwaveprj,(0,0))
1703 : end if
1704 :
1705 16734 : ABI_MALLOC(kg_k,(3,mpw))
1706 : !$OMP PARALLEL DO
1707 1240915 : do ipw=1,npw_k
1708 4946926 : kg_k(:,ipw)=kg(:,ipw+ikg)
1709 : end do
1710 :
1711 22312 : ABI_MALLOC(ylm_k,(npw_k,mpsang*mpsang*psps%useylm))
1712 5578 : ABI_MALLOC(ylmgr_k,(0,0,0))
1713 5578 : if (psps%useylm==1) then
1714 : !$OMP PARALLEL DO COLLAPSE(2)
1715 2080 : do ilm=1,mpsang*mpsang
1716 1280044 : do ipw=1,npw_k
1717 1279628 : ylm_k(ipw,ilm)=ylm(ipw+ikg,ilm)
1718 : end do
1719 : end do
1720 : end if
1721 :
1722 : ! Compute (k+G) vectors
1723 5578 : nkpg=3*nloalg(3)
1724 22312 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
1725 5578 : if (nkpg>0) then
1726 416 : call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
1727 : end if
1728 :
1729 :
1730 : ! Compute nonlocal form factors ffnl at all (k+G)
1731 5578 : ider=0;idir=0;dimffnl=1
1732 27890 : ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat))
1733 : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,gs_hamk%gmet,gs_hamk%gprimd,&
1734 : & ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,&
1735 5578 : & nkpg,npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_k)
1736 :
1737 : ! Load k-dependent part in the Hamiltonian datastructure
1738 : ! - Compute 3D phase factors
1739 : ! - Prepare various tabs in case of band-FFT parallelism
1740 : ! - Load k-dependent quantities in the Hamiltonian
1741 :
1742 22312 : ABI_MALLOC(ph3d,(2,npw_k,gs_hamk%matblk))
1743 : call gs_hamk%load_k(kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,&
1744 5578 : & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,compute_gbound=compute_gbound,compute_ph3d=.true.)
1745 :
1746 : ! Load band-FFT tabs (transposed k-dependent arrays)
1747 5578 : if (mpi_enreg%paral_kgb==1) then
1748 0 : call bandfft_kpt_savetabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav)
1749 0 : call prep_bandfft_tabs(gs_hamk,ikpt,mkmem,mpi_enreg)
1750 : call gs_hamk%load_k(npw_fft_k=my_bandfft_kpt%ndatarecv, &
1751 : & kg_k =my_bandfft_kpt%kg_k_gather, &
1752 : & kpg_k =my_bandfft_kpt%kpg_k_gather, &
1753 : ffnl_k =my_bandfft_kpt%ffnl_gather, &
1754 0 : ph3d_k =my_bandfft_kpt%ph3d_gather,compute_gbound=compute_gbound)
1755 : end if
1756 :
1757 : ! The following is now wrong. In sequential, nblockbd=nband_k/bandpp
1758 : ! blocksize= bandpp (JB 2016/04/16)
1759 : ! Note that in sequential mode iblock=iband, nblockbd=nband_k and blocksize=1
1760 : !
1761 16734 : ABI_MALLOC(occblock,(blocksize))
1762 11156 : ABI_MALLOC(weight,(blocksize))
1763 16734 : occblock=zero;weight=zero
1764 :
1765 5578 : if (fockcommon%optfor) then
1766 13345 : fockcommon%forces_ikpt=zero
1767 : end if
1768 :
1769 22312 : ABI_MALLOC(wi,(2,npw_k*my_nspinor*blocksize,nblockbd))
1770 20054547 : wi=zero
1771 22312 : ABI_MALLOC(mkl,(2,nband_k,nband_k))
1772 541338 : mkl=zero
1773 : ! Calculate all the Wi for the current k-point
1774 :
1775 32529 : do iblock=1,nblockbd
1776 :
1777 26951 : iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband_k)
1778 26951 : iband_cprj=(iblock-1)*bandpp+1
1779 26951 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,me_distrb)) cycle
1780 :
1781 : ! Select occupied bandsddk
1782 53902 : occblock(:)=occ(1+(iblock-1)*blocksize+bdtot_index:iblock*blocksize+bdtot_index)
1783 53902 : weight(:)=wtk(ikpt)*occblock(:)
1784 :
1785 : ! Load contribution from n,k
1786 : cwavef(:,1:npw_k*my_nspinor*blocksize)=&
1787 20048969 : & cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg)
1788 26951 : if (psps%usepaw==1) then
1789 : call pawcprj_get(gs_hamk%atindx1,cwaveprj,cprj,natom,iband_cprj,ibg,ikpt,0,isppol,&
1790 : & mband_cprj,mkmem,natom,bandpp,nband_cprj_k,my_nspinor,nsppol,0,&
1791 1847 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
1792 : end if
1793 :
1794 26951 : if (mpi_enreg%paral_kgb==1) then
1795 0 : msg='fock2ACE: Paral_kgb is not yet implemented for fock calculations'
1796 0 : ABI_BUG(msg)
1797 : end if
1798 :
1799 26951 : fockcommon%ieigen=(iblock-1)*blocksize+1
1800 26951 : fockcommon%iband=(iblock-1)*blocksize+1
1801 26951 : call timab(1562,2,tsec) ; call timab(1563,-1,tsec)
1802 : call fock_getghc(cwavef,&
1803 : cwaveprj,&
1804 : wi(:,:,iblock),&
1805 26951 : gs_hamk,mpi_enreg,blocksize)
1806 :
1807 :
1808 26951 : call timab(1563,2,tsec) ; call timab(1562,-1,tsec)
1809 59480 : do iblocksize=1,blocksize
1810 26951 : mkl(1,fockcommon%ieigen+iblocksize-1,fockcommon%ieigen+iblocksize-1)=fockcommon%eigen_ikpt(fockcommon%ieigen+iblocksize-1)
1811 26951 : fockcommon%e_fock0=fockcommon%e_fock0+half*weight(iblocksize)*fockcommon%eigen_ikpt(fockcommon%ieigen+iblocksize-1)
1812 53902 : if (fockcommon%optfor) then
1813 12996 : fockcommon%forces(:,:)=fockcommon%forces(:,:)+weight(iblocksize)*fockcommon%forces_ikpt(:,:,fockcommon%ieigen+iblocksize-1)
1814 : end if
1815 : end do
1816 :
1817 :
1818 : end do ! End of loop on block of bands
1819 :
1820 : ! Calculate Mkl for the current k-point
1821 16734 : ABI_MALLOC(cwavefk,(2,npw_k*my_nspinor))
1822 32529 : do iblock=1,nblockbd
1823 : cwavef(:,1:npw_k*my_nspinor*blocksize)=&
1824 20048969 : & cg(:,1+(iblock-1)*npw_k*my_nspinor*blocksize+icg:iblock*npw_k*my_nspinor*blocksize+icg)
1825 59480 : do iblocksize=1,blocksize
1826 26951 : kk=(iblock-1)*blocksize+iblocksize
1827 20048969 : cwavefk(:,:)=cwavef(:,1+(iblocksize-1)*npw_k*my_nspinor:iblocksize*npw_k*my_nspinor)
1828 152179 : do jblock=1,iblock
1829 223505 : do jblocksize=1,blocksize
1830 98277 : ll=(jblock-1)*blocksize+jblocksize
1831 196554 : if (ll<kk) then
1832 : call dotprod_g(mkl(1,kk,ll),mkl(2,kk,ll),gs_hamk%istwf_k,npw_k,2,wi(:,1+(jblocksize-1)*npw_k*my_nspinor:&
1833 71326 : & jblocksize*npw_k*my_nspinor,jblock),cwavefk,mpi_enreg%me_g0,mpi_enreg%comm_fft)
1834 : end if
1835 : end do
1836 : end do
1837 : end do
1838 : end do ! End of loop on block of bands
1839 :
1840 5578 : ABI_FREE(cwavefk)
1841 541338 : mkl=-mkl
1842 :
1843 : ! Cholesky factorisation of -mkl=Lx(trans(L)*. On output mkl=L
1844 5578 : call abi_xpotrf("L",nband_k,mkl,nband_k,info,x_cplx=2)
1845 :
1846 : ! calculate trans(L-1)
1847 16734 : ABI_MALLOC(bb,(2,nband_k,nband_k))
1848 541338 : bb=zero
1849 32529 : do kk=1,nband_k
1850 32529 : bb(1,kk,kk)=one
1851 : end do
1852 5578 : call ztrtrs("L","T","N",nband_k,nband_k,mkl,nband_k,bb,nband_k,info)
1853 20054547 : fock%fockACE(ikpt,isppol)%xi=zero
1854 :
1855 : ! Calculate ksi
1856 32529 : do kk=1,nband_k
1857 202132 : do jblock=1,nblockbd
1858 366157 : do jblocksize=1,blocksize
1859 169603 : ll=(jblock-1)*blocksize+jblocksize
1860 : fock%fockACE(ikpt,isppol)%xi(1,:,kk)=fock%fockACE(ikpt,isppol)%xi(1,:,kk)+bb(1,ll,kk)*wi(1,1+(jblocksize-1)*&
1861 : & npw_k*my_nspinor:jblocksize*npw_k*my_nspinor,jblock)-&
1862 42655861 : & bb(2,ll,kk)*wi(2,1+(jblocksize-1)*npw_k*my_nspinor:jblocksize*npw_k*my_nspinor,jblock)
1863 : fock%fockACE(ikpt,isppol)%xi(2,:,kk)=fock%fockACE(ikpt,isppol)%xi(2,:,kk)+bb(1,ll,kk)*wi(2,1+(jblocksize-1)*&
1864 : npw_k*my_nspinor:jblocksize*npw_k*my_nspinor,jblock)+&
1865 42825464 : & bb(2,ll,kk)*wi(1,1+(jblocksize-1)*npw_k*my_nspinor:jblocksize*npw_k*my_nspinor,jblock)
1866 : end do
1867 : end do
1868 : end do
1869 :
1870 : ! DEBUG
1871 : ! fock%fockACE(ikpt,isppol)%xi=zero
1872 : ! ENDDEBUG
1873 :
1874 5578 : ABI_FREE(wi)
1875 5578 : ABI_FREE(mkl)
1876 :
1877 : ! Restore the bandfft tabs
1878 5578 : if (mpi_enreg%paral_kgb==1) then
1879 0 : call bandfft_kpt_restoretabs(my_bandfft_kpt,ffnl=ffnl_sav,ph3d=ph3d_sav,kpg=kpg_k_sav)
1880 : end if
1881 :
1882 : ! Increment indices
1883 5578 : bdtot_index=bdtot_index+nband_k
1884 5578 : if (mkmem/=0) then
1885 5578 : ibg=ibg+my_nspinor*nband_cprj_k
1886 5578 : icg=icg+npw_k*my_nspinor*nband_k
1887 5578 : ikg=ikg+npw_k
1888 : end if
1889 :
1890 5578 : if (psps%usepaw==1) then
1891 416 : call pawcprj_free(cwaveprj)
1892 : end if
1893 6410 : ABI_FREE(cwaveprj)
1894 5578 : ABI_FREE(cwavef)
1895 5578 : ABI_FREE(bb)
1896 5578 : ABI_FREE(occblock)
1897 5578 : ABI_FREE(weight)
1898 5578 : ABI_FREE(ffnl)
1899 5578 : ABI_FREE(kg_k)
1900 5578 : ABI_FREE(kpg_k)
1901 5578 : ABI_FREE(ylm_k)
1902 5578 : ABI_FREE(ylmgr_k)
1903 6951 : ABI_FREE(ph3d)
1904 : end do ! End k point loop
1905 : end do ! End loop over spins
1906 :
1907 1238 : call timab(1562,2,tsec)
1908 1238 : call timab(1565,1,tsec)
1909 :
1910 : !Parallel case: accumulate (n,k) contributions
1911 : if (xmpi_paral==1) then
1912 1238 : call xmpi_sum(fockcommon%e_fock0,spaceComm,ierr)
1913 : ! Forces
1914 1238 : if (optfor==1) then
1915 220 : if (psps%usepaw==1) then
1916 144 : call xmpi_sum(fockcommon%forces,spaceComm,ierr)
1917 : end if
1918 : end if
1919 : end if
1920 :
1921 : !need to reorder cprj=<p_lmn|Cnk> (from atom-sorted to unsorted)
1922 1238 : if (psps%usepaw==1) then
1923 204 : call pawcprj_reorder(cprj,gs_hamk%atindx1)
1924 : end if
1925 : !Deallocate temporary space
1926 1238 : call gs_hamk%free()
1927 :
1928 1238 : call timab(1565,2,tsec)
1929 1238 : call timab(1560,2,tsec)
1930 :
1931 2476 : end subroutine fock2ACE
1932 : !!***
1933 :
1934 : !!****f* ABINIT/fock_ACE_getghc
1935 : !! NAME
1936 : !! fock_ACE_getghc
1937 : !!
1938 : !! FUNCTION
1939 : !! Compute the matrix elements <G|Vx|psi> of the Fock operator in the ACE context.
1940 : !!
1941 : !! INPUTS
1942 : !! cwavef(2,npw*nspinor*ndat)= planewave coefficients of wavefunctions on which Fock operator is applied.
1943 : !! gs_ham <type(gs_hamiltonian_type)>=all data for the Hamiltonian to be applied
1944 : !! mpi_enreg= information about MPI parallelization
1945 : !!
1946 : !! SIDE EFFECTS
1947 : !! ghc(2,npw*ndat)= matrix elements <G|H|C> or <G|H-lambda.S|C> (if sij_opt>=0 or =-1 in getghc)
1948 : !! contains the fock exchange term for cwavef at the end.
1949 : !!
1950 : !! NOTES
1951 : !! The current version assumes that :
1952 : !! * nspinor = 1
1953 : !! * no "my_nspinor"
1954 : !! * no restriction to the value of istwfk_bz (but must be tested in all case)
1955 : !! * all the data for the occupied states (cgocc_bz) are the same as those for the current states (cg)
1956 : !!
1957 : !! SOURCE
1958 :
1959 142256 : subroutine fock_ACE_getghc(cwavef,ghc,gs_ham,mpi_enreg,ndat,gpu_option)
1960 :
1961 : !Arguments ------------------------------------
1962 : ! Scalars
1963 : integer :: ndat
1964 : integer,optional :: gpu_option
1965 : type(MPI_type),intent(in) :: mpi_enreg
1966 : type(gs_hamiltonian_type),target,intent(inout) :: gs_ham
1967 : ! Arrays
1968 : real(dp),target,intent(inout) :: cwavef(:,:)!,ghc(2,gs_ham%npw_k*ndat)
1969 : real(dp),target,intent(inout) :: ghc(:,:)
1970 :
1971 : !Local variables-------------------------------
1972 : ! Scalars
1973 : complex(dp), parameter :: cminusone = (-1._dp,0._dp)
1974 : integer :: iband,ikpt,ipw,my_nspinor,nband_k,npw,idat,gpu_option_
1975 : real(dp) :: eigen
1976 : type(fock_common_type),pointer :: fockcommon
1977 : ! Arrays
1978 : real(dp) :: tsec(2)
1979 142256 : real(dp), target, allocatable :: mat(:,:,:),ghc1(:,:)
1980 142256 : real(dp), allocatable :: vdotr(:),vdoti(:)
1981 142256 : real(dp), ABI_CONTIGUOUS pointer :: xi(:,:,:)
1982 :
1983 : ! *************************************************************************
1984 :
1985 142256 : call timab(1580,1,tsec)
1986 :
1987 142256 : ABI_CHECK(associated(gs_ham%fockcommon),"fock must be associated!")
1988 142256 : fockcommon => gs_ham%fockcommon
1989 :
1990 142256 : ABI_CHECK(gs_ham%nspinor==1,"only allowed for nspinor=1!")
1991 142256 : ABI_CHECK(gs_ham%npw_k==gs_ham%npw_kp,"only allowed for npw_k=npw_kp (ground state)!")
1992 :
1993 142256 : ikpt=fockcommon%ikpt
1994 142256 : npw=gs_ham%npw_k
1995 142256 : nband_k=fockcommon%nband(ikpt)
1996 142256 : my_nspinor=max(1,gs_ham%nspinor/mpi_enreg%nproc_spinor)
1997 142256 : gpu_option_=ABI_GPU_DISABLED; if(present(gpu_option)) gpu_option_ = gpu_option
1998 : !*Initialization of the array ghc1
1999 : !*ghc1 will contain the exact exchange contribution to the Hamiltonian
2000 426768 : ABI_MALLOC(ghc1,(2,npw*my_nspinor*ndat))
2001 107443100 : ghc1=zero
2002 :
2003 142256 : xi => gs_ham%fockACE_k%xi(:,:,:)
2004 :
2005 142256 : if(gpu_option_==ABI_GPU_DISABLED) then
2006 142256 : if(gs_ham%istwf_k==1) then
2007 568480 : ABI_MALLOC(mat,(2,nband_k,ndat))
2008 : call abi_zgemm_2r('C', 'N', nband_k, ndat, npw, cone, &
2009 : xi, npw, &
2010 : cwavef, npw, &
2011 : czero, &
2012 142120 : mat, nband_k)
2013 : call abi_zgemm_2r('N', 'N', npw, ndat, nband_k, cminusone, &
2014 : xi, npw, &
2015 : mat, nband_k, &
2016 : czero, &
2017 142120 : ghc1, npw)
2018 142120 : ABI_FREE(mat)
2019 : else
2020 408 : ABI_MALLOC(vdotr,(nband_k))
2021 272 : ABI_MALLOC(vdoti,(nband_k))
2022 272 : do idat=1,ndat
2023 : call dotprod_g_batch_half(vdotr,vdoti,gs_ham%istwf_k,npw*my_nspinor,nband_k,2,&
2024 136 : & cwavef(:,1+(idat-1)*npw:idat*npw),xi(:,:,:),mpi_enreg%me_g0,mpi_enreg%comm_fft)
2025 :
2026 408 : do iband=1, nband_k
2027 820624 : ghc1(1,1+(idat-1)*npw:idat*npw)=ghc1(1,1+(idat-1)*npw:idat*npw)-vdotr(iband)*xi(1,:,iband)
2028 820760 : ghc1(2,1+(idat-1)*npw:idat*npw)=ghc1(2,1+(idat-1)*npw:idat*npw)-vdotr(iband)*xi(2,:,iband)
2029 : end do
2030 : end do
2031 136 : ABI_FREE(vdotr)
2032 136 : ABI_FREE(vdoti)
2033 : end if
2034 :
2035 : !* If the calculation is parallelized, perform an MPI_allreduce to sum all the contributions in the array ghc
2036 : ! ghc(:,:)=ghc(:,:)/mpi_enreg%nproc_spkpt + ghc1(:,:)
2037 214886200 : ghc(:,:)=ghc(:,:) + ghc1(:,:)
2038 :
2039 : ! call xmpi_sum(ghc,mpi_enreg%comm_kpt,ier)
2040 :
2041 : else if(gpu_option_==ABI_GPU_OPENMP) then
2042 : #ifdef HAVE_OPENMP_OFFLOAD
2043 : !$OMP TARGET ENTER DATA MAP(alloc:xi,ghc1)
2044 : if(gs_ham%istwf_k==1) then
2045 : ABI_MALLOC(mat,(2,nband_k,ndat))
2046 : !$OMP TARGET ENTER DATA MAP(alloc:mat)
2047 : !$OMP TARGET UPDATE TO(xi)
2048 : call abi_zgemm_2r('C', 'N', nband_k, ndat, npw, cone, &
2049 : & xi, npw, &
2050 : & cwavef, npw, &
2051 : & czero, &
2052 : & mat, nband_k, &
2053 : & gpu_option=gpu_option_)
2054 : call abi_zgemm_2r('N', 'N', npw, ndat, nband_k, cminusone, &
2055 : & xi, npw, &
2056 : & mat, nband_k, &
2057 : & czero, &
2058 : & ghc1, npw, &
2059 : & gpu_option=gpu_option_)
2060 :
2061 : !$OMP TARGET EXIT DATA MAP(delete:mat)
2062 : ABI_FREE(mat)
2063 : else
2064 : ABI_MALLOC(vdotr,(nband_k))
2065 : ABI_MALLOC(vdoti,(nband_k))
2066 : !$OMP TARGET ENTER DATA MAP(alloc:vdotr,vdoti)
2067 : do idat=1,ndat
2068 : call dotprod_g_batch_half(vdotr,vdoti,gs_ham%istwf_k,npw*my_nspinor,nband_k,2,&
2069 : & cwavef(:,1+(idat-1)*npw:idat*npw),xi(:,:,:),mpi_enreg%me_g0,mpi_enreg%comm_fft,&
2070 : & gpu_option=gpu_option_)
2071 :
2072 : !$OMP TARGET TEAMS DISTRIBUTE PRIVATE(iband) MAP(ghc1,xi,vdotr)
2073 : do iband=1, nband_k
2074 : !$OMP PARALLEL DO PRIVATE(ipw)
2075 : do ipw=1,npw
2076 : ghc1(1,ipw+(idat-1)*npw)=ghc1(1,ipw+(idat-1)*npw)-vdotr(iband)*xi(1,ipw,iband)
2077 : ghc1(2,ipw+(idat-1)*npw)=ghc1(2,ipw+(idat-1)*npw)-vdotr(iband)*xi(2,ipw,iband)
2078 : end do
2079 : end do
2080 : end do
2081 : !$OMP TARGET EXIT DATA MAP(delete:vdotr,vdoti)
2082 : ABI_FREE(vdotr)
2083 : ABI_FREE(vdoti)
2084 : end if
2085 :
2086 : !* If the calculation is parallelized, perform an MPI_allreduce to sum all the contributions in the array ghc
2087 : ! ghc(:,:)=ghc(:,:)/mpi_enreg%nproc_spkpt + ghc1(:,:)
2088 :
2089 : call abi_xaxpy(npw*ndat,cone,ghc1,1,ghc,1,x_cplx=2,gpu_option=gpu_option)
2090 :
2091 : ! call xmpi_sum(ghc,mpi_enreg%comm_kpt,ier)
2092 : !$OMP TARGET UPDATE FROM(ghc1)
2093 :
2094 : !$OMP TARGET EXIT DATA MAP(delete:xi,ghc1)
2095 : #endif
2096 : end if
2097 :
2098 :
2099 : ! ============================================
2100 : ! === Calculate the contribution to energy ===
2101 : ! ============================================
2102 : !* Only the contribution when cwavef=cgocc_bz are calculated, in order to cancel exactly the self-interaction
2103 : !* at each convergence step. (consistent definition with the definition of hartree energy)
2104 142256 : if (fockcommon%ieigen/=0) then
2105 68518 : do idat=1,ndat
2106 34259 : eigen=zero
2107 : ! * Dot product of cwavef and ghc
2108 : ! * inspired from the routine 54_spacepar/meanvalue_g but without the reference to parallelism and filtering
2109 34259 : if(gs_ham%istwf_k==2) then
2110 29 : eigen=half*cwavef(1,1+(idat-1)*npw)*ghc1(1,1+(idat-1)*npw)
2111 : else
2112 34230 : eigen=cwavef(1,1+(idat-1)*npw)*ghc1(1,1+(idat-1)*npw)+cwavef(2,1+(idat-1)*npw)*ghc1(2,1+(idat-1)*npw)
2113 : end if
2114 8415726 : do ipw=2,npw
2115 8415726 : eigen=eigen+cwavef(1,ipw+(idat-1)*npw)*ghc1(1,ipw+(idat-1)*npw)+cwavef(2,ipw+(idat-1)*npw)*ghc1(2,ipw+(idat-1)*npw)
2116 : end do
2117 34259 : if(gs_ham%istwf_k>=2) eigen=two*eigen
2118 : ! call xmpi_sum(eigen,mpi_enreg%comm_kpt,ier)
2119 34259 : fockcommon%eigen_ikpt(fockcommon%ieigen+idat-1)= eigen
2120 68518 : fockcommon%ieigen = 0
2121 : end do ! idat
2122 : end if
2123 :
2124 : ! ===============================
2125 : ! === Deallocate local arrays ===
2126 : ! ===============================
2127 :
2128 142256 : ABI_FREE(ghc1)
2129 142256 : call timab(1580,2,tsec)
2130 :
2131 142256 : end subroutine fock_ACE_getghc
2132 : !!***
2133 :
2134 : end module m_fock_getghc
2135 : !!***
|
