Branch data 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 : : !!***
|
