Branch data Line data Source code
1 : : !!****m* ABINIT/m_dfpt_mkrho
2 : : !! NAME
3 : : !! m_dfpt_mkrho
4 : : !!
5 : : !! FUNCTION
6 : : !! Compute RF charge density rho1(r) and rho1(G) in electrons/bohr**3
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, LSI, AR, MB, MT, SPr)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : ! nvtx related macro definition
23 : : #include "nvtx_macros.h"
24 : :
25 : : module m_dfpt_mkrho
26 : :
27 : : use defs_basis
28 : : use m_abicore
29 : : use m_errors
30 : : use m_cgtools
31 : : use m_xmpi
32 : :
33 : : use defs_abitypes, only : MPI_type
34 : : use m_time, only : timab
35 : : use m_io_tools, only : get_unit, iomode_from_fname
36 : : use m_fftcore, only : sphereboundary
37 : : use m_fft, only : fftpac, fourwf
38 : : use m_spacepar, only : symrhg
39 : : use m_hamiltonian, only : gs_hamiltonian_type
40 : : use m_pawrhoij, only : pawrhoij_type
41 : : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free
42 : : use m_paw_occupancies, only : pawaccrhoij
43 : : use m_paral_atom, only : get_my_atmtab
44 : : use m_mpinfo, only : proc_distrb_cycle
45 : : use m_cgprj, only : getcprj
46 : :
47 : : #if defined(HAVE_GPU_MARKERS)
48 : : use m_nvtx_data
49 : : #endif
50 : :
51 : : implicit none
52 : :
53 : : private
54 : : !!***
55 : :
56 : : public :: dfpt_mkrho
57 : : public :: dfpt_accrho
58 : : !!***
59 : :
60 : : contains
61 : : !!***
62 : :
63 : : !!****f* ABINIT/dfpt_mkrho
64 : : !! NAME
65 : : !! dfpt_mkrho
66 : : !!
67 : : !! FUNCTION
68 : : !! Compute RF charge density rho1(r) and rho1(G) in electrons/bohr**3
69 : : !! from input RF and GS wavefunctions, band occupations, and k point weights.
70 : : !!
71 : : !! INPUTS
72 : : !! cg(2,mpw*nspinor*mband*mkmem*nsppol)=wf in G space
73 : : !! cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=first-order wf in G space
74 : : !! cplex=1 if rhor1 is real, 2 if rhor1 is complex
75 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
76 : : !! irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
77 : : !! istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs
78 : : !! kg(3,mpw*mkmem)=reduced planewave coordinates, GS data.
79 : : !! kg1(3,mpw1*mkmem1)=reduced planewave coordinates, RF data.
80 : : !! mband=maximum number of bands
81 : : !! mgfft=maximum size of 1D FFTs
82 : : !! mkmem=Number of k points treated by this node (GS data)
83 : : !! mk1mem=Number of k points treated by this node (RF data)
84 : : !! mpi_enreg=information about MPI parallelization
85 : : !! mpw=maximum allowed value for npw (GS wfs)
86 : : !! mpw1=maximum allowed value for npw1 (RF data)
87 : : !! nband_rbz(nkpt_rbz*nsppol)=number of bands to be included in summation
88 : : !! at each k point for each spin channel.
89 : : !! nfft=(effective) number of FFT grid points (for this processor)
90 : : !! ngfft(18)=contain all needed information about 3D FFT,
91 : : !! see ~abinit/doc/variables/vargs.htm#ngfft
92 : : !! nkpt_rbz=number of k points in the reduced Brillouin zone
93 : : !! npwarr(nkpt_rbz)=number of planewaves and boundary planewaves at k points
94 : : !! npwar1(nkpt_rbz)=number of planewaves and boundary planewaves at k+q points
95 : : !! nspden=number of spin-density components
96 : : !! nspinor=number of spinorial components of the wavefunctions
97 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
98 : : !! nsym=number of symmetry elements in group (at least 1 for identity)
99 : : !! occ_rbz(mband*nkpt_rbz*nsppol)=occupation numbers for each band
100 : : !! (usually 2.0) at each k point of the reduced Brillouin zone
101 : : !! phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
102 : : !! rprimd(3,3)=dimensional real space primitive translations
103 : : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
104 : : !! symrel(3,3,nsym)=symmetry matrices in real space (integers)
105 : : !! ucvol=unit cell volume (Bohr**3)
106 : : !! wtk_rbz(nkpt_rbz)=k point weights (they sum to 1.0).
107 : : !!
108 : : !! OUTPUT
109 : : !! rhog1(2,nfft)=total electron density in G space
110 : : !! rhor1(cplex*nfft,nspden)=electron density in r space
111 : : !! (if spin polarized, array contains total density in first half and
112 : : !! spin-up density in second half)
113 : : !!
114 : : !! SOURCE
115 : :
116 : 8 : subroutine dfpt_mkrho(cg,cg1,cplex,gprimd,irrzon,istwfk_rbz,&
117 : 8 : & kg,kg1,mband,mband_mem,mgfft,mkmem,mk1mem,mpi_enreg,mpw,mpw1,nband_rbz,&
118 : 8 : & nfft,ngfft,nkpt_rbz,npwarr,npwar1,nspden,nspinor,nsppol,nsym,&
119 : 8 : & occ_rbz,phnons,rhog1,rhor1,rprimd,symafm,symrel,tnons,ucvol,wtk_rbz)
120 : :
121 : : !Arguments ------------------------------------
122 : : !scalars
123 : : integer,intent(in) :: cplex,mband,mband_mem,mgfft,mk1mem,mkmem,mpw,mpw1,nfft,nkpt_rbz
124 : : integer,intent(in) :: nspden,nspinor,nsppol,nsym
125 : : real(dp),intent(in) :: ucvol
126 : : type(MPI_type),intent(in) :: mpi_enreg
127 : : !arrays
128 : : integer,intent(in) :: irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))
129 : : integer,intent(in) :: istwfk_rbz(nkpt_rbz),kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem)
130 : : integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),ngfft(18),npwar1(nkpt_rbz)
131 : : integer,intent(in) :: npwarr(nkpt_rbz),symafm(nsym),symrel(3,3,nsym)
132 : : real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem*nsppol)
133 : : real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol),gprimd(3,3)
134 : : real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol)
135 : : real(dp),intent(in) :: phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))
136 : : real(dp),intent(in) :: rprimd(3,3),tnons(3,nsym)
137 : : real(dp),intent(in) :: wtk_rbz(nkpt_rbz)
138 : : real(dp),intent(out) :: rhog1(2,nfft),rhor1(cplex*nfft,nspden)
139 : :
140 : : !Local variables-------------------------------
141 : : !scalars
142 : : integer,parameter :: tim_fourwf7=7,tim_rwwf15=15
143 : : integer,save :: nskip=0
144 : : integer :: bdtot_index,i1,i2,i3,iband,icg,icg1,ierr,ifft,ikg,ptr
145 : : integer :: iband_me
146 : : integer :: ikg1,ikpt,ispden,ispinor,isppol,istwf_k,ptr1,ptr2
147 : : integer :: me,n1,n2,n3,n4,n5,n6,nband_k,npw1_k
148 : : integer :: npw_k,spaceworld
149 : : real(dp) :: im0,im1,re0,re1,weight
150 : : real(dp) :: im0_up,im1_up,re0_up,re1_up,im0_down,im1_down,re0_down,re1_down
151 : : character(len=500) :: message
152 : : !arrays
153 : 8 : integer,allocatable :: gbound(:,:),gbound1(:,:),kg1_k(:,:)
154 : 8 : integer,allocatable :: kg_k(:,:)
155 : : real(dp) :: tsec(2)
156 : 8 : real(dp),allocatable,target :: cwavef(:,:),cwavef1(:,:)
157 : 8 : real(dp),allocatable :: dummy(:,:),rhoaug(:,:,:,:)
158 : 8 : real(dp),allocatable :: rhoaug1(:,:,:,:),wfraug(:,:,:,:),wfraug1(:,:,:,:)
159 : 8 : real(dp),allocatable :: wfraug1_up(:,:,:,:),wfraug1_down(:,:,:,:)
160 : 8 : real(dp),allocatable :: wfraug_up(:,:,:,:),wfraug_down(:,:,:,:)
161 : 8 : real(dp),allocatable :: cwave0_up(:,:),cwave0_down(:,:),cwave1_up(:,:),cwave1_down(:,:)
162 : :
163 : : ! *************************************************************************
164 : :
165 : : !DBG_ENTER("COLL")
166 : : ABI_NVTX_START_RANGE(NVTX_DFPT_MKRHO)
167 : :
168 [ - + ]: 8 : if(nspden==4)then
169 : : ! NOTE: see mkrho for the modifications needed for non-collinear treatment
170 : : write(message, '(3a)' )&
171 : 0 : ' Linear-response calculations are under construction with nspden=4',ch10,&
172 : 0 : ' Action: modify value of nspden in input file unless you know what you are doing.'
173 : 0 : ABI_WARNING(message)
174 : : end if
175 : :
176 : : !Init spaceworld
177 : 8 : spaceworld=mpi_enreg%comm_cell
178 : 8 : me=mpi_enreg%me_kpt
179 : :
180 : : !zero the charge density array in real space
181 : : !$OMP PARALLEL DO
182 [ + + ]: 16 : do ispden=1,nspden
183 [ + + ]: 28976 : do ifft=1,cplex*nfft
184 : 28968 : rhor1(ifft,ispden)=zero
185 : : end do
186 : : end do
187 : :
188 : : !start loop over spin and k points
189 : 8 : bdtot_index=0; icg=0; icg1=0
190 : :
191 : 8 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
192 : 8 : n4=ngfft(4); n5=ngfft(5); n6=ngfft(6) !n4,n5,n6 are FFT dimensions, modified to avoid cache trashing
193 : :
194 : : !Note that the dimensioning of cwavef and cwavef1 does not include nspinor
195 [ - + + - : 24 : ABI_MALLOC(cwavef,(2,mpw))
+ - + - -
+ - + ]
196 [ - + + - : 24 : ABI_MALLOC(cwavef1,(2,mpw1))
+ - + - -
+ - + ]
197 : : !Actually, rhoaug is not needed, except for strong dimensioning requirement
198 [ - + ]: 8 : ABI_MALLOC(dummy,(2,1))
199 [ - + + - : 48 : ABI_MALLOC(rhoaug,(n4,n5,n6,nspinor**2))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
200 [ - + + - : 48 : ABI_MALLOC(rhoaug1,(cplex*n4,n5,n6,nspinor**2))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
201 [ - + + - : 40 : ABI_MALLOC(wfraug,(2,n4,n5,n6))
- + + - -
+ + - + -
+ - - + -
+ ]
202 [ - + + - : 32 : ABI_MALLOC(wfraug1,(2,n4,n5,n6))
- + + - -
+ + - - +
- + ]
203 : :
204 : : ! EB FR Separate collinear and non-collinear magnetism
205 [ + - ]: 8 : if (nspden /= 4) then ! EB FR nspden check
206 [ + + ]: 16 : do isppol=1,nsppol
207 : :
208 : 16 : ikg=0; ikg1=0
209 : :
210 [ + + + + : 34544 : rhoaug1(:,:,:,:)=zero
+ + + + ]
211 : :
212 [ + + ]: 370 : do ikpt=1,nkpt_rbz
213 : :
214 : 362 : nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
215 : 362 : istwf_k=istwfk_rbz(ikpt)
216 : 362 : npw_k=npwarr(ikpt)
217 : 362 : npw1_k=npwar1(ikpt)
218 : :
219 [ + + ]: 362 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
220 : 216 : bdtot_index=bdtot_index+nband_k
221 : 216 : cycle
222 : : end if
223 : :
224 [ + - - + : 584 : ABI_MALLOC(gbound,(2*mgfft+8,2))
+ - + - +
- + - - +
- + ]
225 [ - + + - : 438 : ABI_MALLOC(kg_k,(3,npw_k))
+ - + - -
+ - + ]
226 [ - + + - : 292 : ABI_MALLOC(gbound1,(2*mgfft+8,2))
- + - + ]
227 [ - + + - : 438 : ABI_MALLOC(kg1_k,(3,npw1_k))
+ - + - -
+ - + ]
228 : :
229 [ + + + + ]: 73582 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
230 : 146 : call sphereboundary(gbound,istwf_k,kg_k,mgfft,npw_k)
231 : :
232 [ + + + + ]: 73582 : kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
233 : 146 : call sphereboundary(gbound1,istwf_k,kg1_k,mgfft,npw1_k)
234 : :
235 : : ! Loop over bands to fft and square for rho(r)
236 : 146 : iband_me = 0
237 [ + + ]: 1306 : do iband=1,nband_k
238 [ + - ]: 1160 : if (mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me) cycle
239 : 1160 : iband_me = iband_me + 1
240 : : ! Only treat occupied states
241 [ + + ]: 1306 : if (abs(occ_rbz(iband+bdtot_index))>tol8) then
242 : : ! Treat separately the two spinor components
243 [ + + ]: 1168 : do ispinor=1,nspinor
244 : : ! Obtain Fourier transform in fft box and accumulate the density
245 : 584 : ptr = 1 + (ispinor-1)*npw_k + (iband_me-1)*npw_k*nspinor + icg
246 : 584 : call cg_zcopy(npw_k, cg(1,ptr), cwavef)
247 : :
248 : : ! In these two calls, rhoaug, rhoaug1 and weight are dummy variables, and are not modified
249 : : call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,&
250 : 584 : & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,0,tim_fourwf7,weight,weight)
251 : :
252 : : ! TODO: here ispinor should be ispinorp to get full matrix and nspden 4
253 : 584 : ptr = 1 + (ispinor-1)*npw1_k + (iband_me-1)*npw1_k*nspinor + icg1
254 : 584 : call cg_zcopy(npw1_k, cg1(1,ptr), cwavef1)
255 : :
256 : : call fourwf(cplex,rhoaug1,cwavef1,dummy,wfraug1,gbound1,gbound1,&
257 : : & istwf_k,kg1_k,kg1_k,mgfft,mpi_enreg,1,ngfft,npw1_k,1,n4,n5,n6,0,&
258 : 584 : & tim_fourwf7,weight,weight)
259 : :
260 : : ! Compute the weight, note that the factor 2 is
261 : : ! not the spin factor (see Eq.44 of PRB55,10337 (1997) [[cite:Gonze1997]])
262 : 584 : weight=two*occ_rbz(iband+bdtot_index)*wtk_rbz(ikpt)/ucvol
263 : :
264 : : ! Accumulate density
265 [ - + ]: 1168 : if(cplex==2)then
266 : : !$OMP PARALLEL DO PRIVATE(im0,im1,re0,re1)
267 [ # # ]: 0 : do i3=1,n3
268 [ # # ]: 0 : do i2=1,n2
269 [ # # ]: 0 : do i1=1,n1
270 : 0 : re0=wfraug(1,i1,i2,i3) ; im0=wfraug(2,i1,i2,i3)
271 : 0 : re1=wfraug1(1,i1,i2,i3); im1=wfraug1(2,i1,i2,i3)
272 : 0 : rhoaug1(2*i1-1,i2,i3,1)=rhoaug1(2*i1-1,i2,i3,1)+weight*(re0*re1+im0*im1)
273 : 0 : rhoaug1(2*i1 ,i2,i3,1)=rhoaug1(2*i1 ,i2,i3,1)+weight*(re0*im1-im0*re1)
274 : : end do
275 : : end do
276 : : end do
277 : : else
278 : : !$OMP PARALLEL DO
279 [ + + ]: 8200 : do i3=1,n3
280 [ + + ]: 121992 : do i2=1,n2
281 [ + + ]: 1981248 : do i1=1,n1
282 : : rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1)+&
283 : 1973632 : & weight*( wfraug(1,i1,i2,i3)*wfraug1(1,i1,i2,i3) + wfraug(2,i1,i2,i3)*wfraug1(2,i1,i2,i3) )
284 : : end do
285 : : end do
286 : : end do
287 : : end if ! cplex
288 : : end do ! ispinor
289 : : else !abs(occ_rbz(iband+bdtot_index))>tol8
290 : 576 : nskip=nskip+1 ! if the state is not occupied. Accumulate the number of one-way 3D ffts skipped
291 : : end if ! abs(occ_rbz(iband+bdtot_index))>tol8
292 : :
293 : : end do ! iband
294 : :
295 [ - + ]: 146 : ABI_FREE(gbound)
296 [ - + ]: 146 : ABI_FREE(kg_k)
297 [ - + ]: 146 : ABI_FREE(gbound1)
298 [ - + ]: 146 : ABI_FREE(kg1_k)
299 : :
300 : 146 : bdtot_index=bdtot_index+nband_k
301 : :
302 : : ! only increase indices for my bands on my proc
303 : 146 : icg=icg+npw_k*mband_mem*nspinor
304 : 146 : ikg=ikg+npw_k
305 : :
306 : 146 : icg1=icg1+npw1_k*mband_mem*nspinor
307 : 154 : ikg1=ikg1+npw1_k
308 : :
309 : : end do ! ikpt
310 : :
311 : : if (xmpi_paral==0) then ! Write the number of one-way 3D ffts skipped until now
312 : : write(message,'(a,i8)')' mkrho3 : number of one-way 3D ffts skipped in mkrho3 until now =',nskip
313 : : call wrtout(std_out,message,'PERS')
314 : : end if
315 : :
316 : : ! Transfer density on augmented fft grid to normal fft grid in real space
317 : : ! Take also into account the spin, to place it correctly in rhor1.
318 : : ! Note the use of cplex
319 : 16 : call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,rhor1,rhoaug1,1)
320 : :
321 : : end do ! loop over isppol spins
322 : :
323 : : else ! nspden = 4
324 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
325 : : ! Part added for the non collinear magnetism
326 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
327 : : ! The same lines of code are in 72_response/accrho3.F90
328 : : ! TODO: merge these lines in a single routine??!!
329 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
330 : :
331 : 0 : ikg=0; ikg1=0
332 : :
333 [ # # # # : 0 : rhoaug1(:,:,:,:)=zero
# # # # ]
334 : :
335 [ # # ]: 0 : do ikpt=1,nkpt_rbz
336 : :
337 : 0 : nband_k=nband_rbz(ikpt)
338 : 0 : istwf_k=istwfk_rbz(ikpt)
339 : 0 : npw_k=npwarr(ikpt)
340 : 0 : npw1_k=npwar1(ikpt)
341 : :
342 [ # # ]: 0 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,me)) then
343 : 0 : bdtot_index=bdtot_index+nband_k
344 : 0 : cycle
345 : : end if
346 : :
347 [ # # # # : 0 : ABI_MALLOC(gbound,(2*mgfft+8,2))
# # # # #
# # # # #
# # ]
348 [ # # # # : 0 : ABI_MALLOC(kg_k,(3,npw_k))
# # # # #
# # # ]
349 [ # # # # : 0 : ABI_MALLOC(gbound1,(2*mgfft+8,2))
# # # # ]
350 [ # # # # : 0 : ABI_MALLOC(kg1_k,(3,npw1_k))
# # # # #
# # # ]
351 : :
352 [ # # # # ]: 0 : kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
353 : 0 : call sphereboundary(gbound,istwf_k,kg_k,mgfft,npw_k)
354 : :
355 [ # # # # ]: 0 : kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
356 : 0 : call sphereboundary(gbound1,istwf_k,kg1_k,mgfft,npw1_k)
357 : :
358 : : ! Loop over bands to fft and square for rho(r)
359 : 0 : iband_me = 0
360 [ # # ]: 0 : do iband=1,nband_k
361 : :
362 [ # # ]: 0 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,1,me)) cycle
363 : 0 : iband_me = iband_me + 1
364 : :
365 : : ! Only treat occupied states
366 [ # # ]: 0 : if (abs(occ_rbz(iband+bdtot_index))>tol8) then
367 : :
368 : : ! Build the four components of rho. We use only norm quantities and, so fourwf.
369 : :
370 [ # # # # : 0 : ABI_MALLOC(wfraug_up,(2,n4,n5,n6))
# # # # #
# # # # #
# # ]
371 [ # # # # : 0 : ABI_MALLOC(wfraug_down,(2,n4,n5,n6))
# # # # #
# # # # #
# # ]
372 [ # # # # : 0 : ABI_MALLOC(wfraug1_up,(2,n4,n5,n6))
# # # # #
# # # # #
# # ]
373 [ # # # # : 0 : ABI_MALLOC(wfraug1_down,(2,n4,n5,n6))
# # # # #
# # # # #
# # ]
374 [ # # # # : 0 : ABI_MALLOC(cwave0_up,(2,npw_k))
# # # # #
# # # ]
375 [ # # # # : 0 : ABI_MALLOC(cwave0_down,(2,npw_k))
# # # # ]
376 [ # # # # : 0 : ABI_MALLOC(cwave1_up,(2,npw1_k))
# # # # #
# # # ]
377 [ # # # # : 0 : ABI_MALLOC(cwave1_down,(2,npw1_k))
# # # # ]
378 : :
379 : : ! EB FR build spinorial wavefunctions
380 : : ! Obtain Fourier transform in fft box and accumulate the density
381 : : ! zero order up and down spins
382 : 0 : ptr1 = 1 + (iband_me-1)*npw_k*nspinor + icg
383 : 0 : call cg_zcopy(npw_k, cg(1,ptr1), cwave0_up)
384 : 0 : ptr2 = 1 + npw_k + (iband_me-1)*npw_k*nspinor + icg
385 : 0 : call cg_zcopy(npw_k, cg(1,ptr2), cwave0_down)
386 : : ! first order up and down spins
387 : 0 : ptr1 = 1 + (iband_me-1)*npw1_k*nspinor + icg1
388 : 0 : call cg_zcopy(npw1_k, cg1(1,ptr1), cwave1_up)
389 : 0 : ptr2 = 1 + npw1_k + (iband_me-1)*npw1_k*nspinor + icg1
390 : 0 : call cg_zcopy(npw1_k, cg1(1,ptr2), cwave1_down)
391 : :
392 : : !density components
393 : : !GS wfk Fourrier Tranform
394 : : call fourwf(1,rhoaug(:,:,:,2),cwave0_up,dummy,wfraug_up,gbound,gbound,istwf_k,kg_k,kg_k,&
395 : : & mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,&
396 : 0 : & 0,tim_fourwf7,weight,weight)
397 : : call fourwf(1,rhoaug(:,:,:,2),cwave0_down,dummy,wfraug_down,gbound,gbound,istwf_k,kg_k,kg_k,&
398 : : & mgfft,mpi_enreg,1,ngfft,npw_k,1,n4,n5,n6,&
399 : 0 : & 0,tim_fourwf7,weight,weight)
400 : : !1st order wfk Fourrier Transform
401 : : call fourwf(cplex,rhoaug(:,:,:,1),cwave1_up,dummy,wfraug1_up,gbound1,gbound1,istwf_k,kg1_k,kg1_k,&
402 : 0 : & mgfft,mpi_enreg,1,ngfft,npw1_k,1,n4,n5,n6,0,tim_fourwf7,weight,weight)
403 : : call fourwf(cplex,rhoaug(:,:,:,1),cwave1_down,dummy,wfraug1_down,gbound1,gbound1,istwf_k,kg1_k,kg1_k,&
404 : 0 : & mgfft,mpi_enreg,1,ngfft,npw1_k,1,n4,n5,n6,0,tim_fourwf7,weight,weight)
405 : :
406 : : ! The factor 2 is not the spin factor (see Eq.44 of PRB55,10337 (1997) ?? [[cite:Gonze1997]])
407 : 0 : weight=two*occ_rbz(iband+bdtot_index)*wtk_rbz(ikpt)/ucvol
408 : :
409 : : ! Accumulate 1st-order density (x component)
410 [ # # ]: 0 : if (cplex==2) then
411 [ # # ]: 0 : do i3=1,n3
412 [ # # ]: 0 : do i2=1,n2
413 [ # # ]: 0 : do i1=1,n1
414 : 0 : re0_up=wfraug_up(1,i1,i2,i3) ; im0_up=wfraug_up(2,i1,i2,i3)
415 : 0 : re1_up=wfraug1_up(1,i1,i2,i3) ; im1_up=wfraug1_up(2,i1,i2,i3)
416 : 0 : re0_down=wfraug_down(1,i1,i2,i3) ; im0_down=wfraug_down(2,i1,i2,i3)
417 : 0 : re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3)
418 : 0 : rhoaug1(2*i1-1,i2,i3,1)=rhoaug1(2*i1-1,i2,i3,1)+weight*(re0_up*re1_up+im0_up*im1_up) !n_upup
419 : 0 : rhoaug1(2*i1 ,i2,i3,1)=rhoaug1(2*i1 ,i2,i3,1)+weight*(re0_up*im1_up-im0_up*re1_up)
420 : 0 : rhoaug1(2*i1-1,i2,i3,4)=rhoaug1(2*i1-1,i2,i3,4)+weight*(re0_down*re1_down+im0_down*im1_down) ! n_dndn
421 : 0 : rhoaug1(2*i1 ,i2,i3,4)=rhoaug1(2*i1 ,i2,i3,4)+weight*(re0_down*im1_down-im0_down*re1_down)
422 : :
423 : : rhoaug1(2*i1-1,i2,i3,2)=rhoaug1(2*i1-1,i2,i3,2)+weight*(re1_up*re0_down+im1_up*im0_down)& !Re[m1x]
424 : 0 : & +weight*(re1_down*re0_up+im1_down*im0_up)
425 : : rhoaug1(2*i1 ,i2,i3,2)=rhoaug1(2*i1 ,i2,i3,2)+weight*(-re1_up*im0_down+im1_up*re0_down)& !Im[m1x]
426 : 0 : & +weight*(-re1_down*im0_up+im1_down*re0_up)
427 : :
428 : : rhoaug1(2*i1-1,i2,i3,3)=rhoaug1(2*i1-1,i2,i3,3)+weight*(+re1_up*im0_down-im1_up*re0_down)& !Re[m1y]
429 : 0 : & +weight*(-re1_down*im0_up+im1_down*re0_up)
430 : : rhoaug1(2*i1 ,i2,i3,3)=rhoaug1(2*i1 ,i2,i3,3)+weight*(+re1_up*re0_down+im1_up*im0_down)& !Im[m1y]
431 : 0 : & +weight*(-re1_down*re0_up-im1_down*im0_up)
432 : : end do
433 : : end do
434 : : end do
435 : : else
436 : 0 : re0_up=zero;im0_up=zero;re1_up=zero;im1_up=zero;re0_down=zero;im0_down=zero
437 : 0 : re1_down=zero;im1_down=zero
438 [ # # ]: 0 : do i3=1,n3
439 [ # # ]: 0 : do i2=1,n2
440 [ # # ]: 0 : do i1=1,n1
441 : 0 : re0_up=wfraug_up(1,i1,i2,i3) ; im0_up=wfraug_up(2,i1,i2,i3)
442 : 0 : re1_up=wfraug1_up(1,i1,i2,i3) ; im1_up=wfraug1_up(2,i1,i2,i3)
443 : 0 : re0_down=wfraug_down(1,i1,i2,i3) ; im0_down=wfraug_down(2,i1,i2,i3)
444 : 0 : re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3)
445 : 0 : rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1)+weight*(re0_up*re1_up+im0_up*im1_up) ! n_upup
446 : 0 : rhoaug1(i1,i2,i3,4)=rhoaug1(i1,i2,i3,4)+weight*(re0_down*re1_down+im0_down*im1_down) ! n_dndn
447 : : rhoaug1(i1,i2,i3,2)=rhoaug1(i1,i2,i3,2)+weight*(re1_up*re0_down+re0_up*re1_down &
448 : 0 : & +im0_up*im1_down+im0_down*im1_up) !mx; the factor two is inside weight
449 : : rhoaug1(i1,i2,i3,3)=rhoaug1(i1,i2,i3,3)+weight*(re1_up*im0_down-im1_up*re0_down &
450 : 0 : & +re0_up*im1_down-im0_up*re1_down) !my; the factor two is inside weight
451 : : end do
452 : : end do
453 : : end do
454 : : end if
455 [ # # ]: 0 : ABI_FREE(wfraug_up)
456 [ # # ]: 0 : ABI_FREE(wfraug_down)
457 [ # # ]: 0 : ABI_FREE(wfraug1_up)
458 [ # # ]: 0 : ABI_FREE(wfraug1_down)
459 [ # # ]: 0 : ABI_FREE(cwave0_up)
460 [ # # ]: 0 : ABI_FREE(cwave0_down)
461 [ # # ]: 0 : ABI_FREE(cwave1_up)
462 [ # # ]: 0 : ABI_FREE(cwave1_down)
463 : :
464 : : end if ! occupied states
465 : : end do ! End loop on iband
466 : :
467 [ # # ]: 0 : ABI_FREE(gbound)
468 [ # # ]: 0 : ABI_FREE(kg_k)
469 [ # # ]: 0 : ABI_FREE(gbound1)
470 [ # # ]: 0 : ABI_FREE(kg1_k)
471 : :
472 : 0 : bdtot_index=bdtot_index+nband_k
473 : :
474 : : ! only increase indices for my bands on my proc
475 : 0 : icg=icg+npw_k*mband_mem*nspinor
476 : 0 : ikg=ikg+npw_k
477 : :
478 : 0 : icg1=icg1+npw1_k*mband_mem*nspinor
479 : 0 : ikg1=ikg1+npw1_k
480 : :
481 : : end do ! End loop on ikpt
482 : :
483 : :
484 : : if (xmpi_paral==0) then ! Write the number of one-way 3D ffts skipped until now
485 : : write(message,'(a,i8)')' dfpt_mkrho : number of one-way 3D ffts skipped in mkrho3 until now =',nskip
486 : : call wrtout(std_out,message,'PERS')
487 : : end if
488 : :
489 : : ! Transfer density on augmented fft grid to normal fft grid in real space
490 : : ! Take also into account the spin, to place it correctly in rhor1.
491 : 0 : call fftpac(1,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,rhor1,rhoaug1(:,:,:,1),1)
492 [ # # ]: 0 : do ispden=2,4
493 : 0 : call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,rhor1,rhoaug1(:,:,:,ispden),1)
494 : : end do
495 : :
496 : : end if ! nspden /= 4
497 : :
498 : : !if (xmpi_paral==1) then
499 : : !call timab(63,1,tsec)
500 : : !call wrtout(std_out,'dfpt_mkrho: loop on k-points and spins done in parallel','COLL')
501 : : !call xmpi_barrier(spaceworld)
502 : : !call timab(63,2,tsec)
503 : : !end if
504 : :
505 [ - + ]: 8 : ABI_FREE(cwavef)
506 [ - + ]: 8 : ABI_FREE(cwavef1)
507 [ - + ]: 8 : ABI_FREE(dummy)
508 [ - + ]: 8 : ABI_FREE(rhoaug)
509 [ - + ]: 8 : ABI_FREE(rhoaug1)
510 [ - + ]: 8 : ABI_FREE(wfraug)
511 [ - + ]: 8 : ABI_FREE(wfraug1)
512 : :
513 : : !Recreate full rhor1 on all proc.
514 : : !TODO : check this sums correctly on bands as well as k
515 : 8 : call timab(48,1,tsec)
516 : 8 : call timab(71,1,tsec)
517 : 8 : call xmpi_sum(rhor1,spaceworld,ierr)
518 : 8 : call timab(71,2,tsec)
519 : 8 : call timab(48,2,tsec)
520 : :
521 : : !This rotation is managed by symrhg
522 [ - + ]: 8 : if (nspden==4) then
523 [ # # ]: 0 : rhor1(:,2) = rhor1(:,2) + (rhor1(:,1) + rhor1(:,4))
524 [ # # ]: 0 : rhor1(:,3) = rhor1(:,3) + (rhor1(:,1) + rhor1(:,4))
525 : : end if
526 : :
527 : : call symrhg(cplex,gprimd,irrzon,mpi_enreg,nfft,nfft,ngfft,nspden,nsppol,nsym,phnons,&
528 : 8 : rhog1,rhor1,rprimd,symafm,symrel,tnons)
529 : :
530 : : !We now have both rho(r) and rho(G), symmetrized, and if nsppol=2
531 : : !we also have the spin-up density, symmetrized, in rhor1(:,2).
532 : :
533 : : ABI_NVTX_END_RANGE()
534 : : !DBG_EXIT("COLL")
535 : :
536 : 8 : end subroutine dfpt_mkrho
537 : : !!***
538 : :
539 : : !!****f* ABINIT/dfpt_accrho
540 : : !! NAME
541 : : !! dfpt_accrho
542 : : !!
543 : : !! FUNCTION
544 : : !! Response function calculation only:
545 : : !! Accumulate contribution to first-order density due to current (k,band)
546 : : !! Also accumulate zero-order potential part of the 2nd-order total energy (if needed)
547 : : !!
548 : : !! INPUTS
549 : : !! cplex=1 if 1st-order density is real, 2 if 1st-order density is complex
550 : : !! cwave0(2,npw*nspinor)=GS wavefunction at k, in reciprocal space
551 : : !! cwave1(2,npw1*nspinor)=1st-order wavefunction at k,q, in reciprocal space
552 : : !! cwavef(2,npw1*nspinor)=1st-order wavefunction at k,q, in reciprocal space, without correction due to occupation change
553 : : !! cwaveprj0(natom,nspinor*usecprj)= GS wave function at k projected with nl projectors
554 : : !! cwaveprj1(natom,nspinor*usecprj)= 1st-order wave function at k,q projected with nl projectors
555 : : !! gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
556 : : !! iband=index of current band
557 : : !! idir=direction of the current perturbation
558 : : !! ipert=type of the perturbation
559 : : !! isppol=1 index of current spin component
560 : : !! kptopt=option for the generation of k points
561 : : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
562 : : !! comm_atom=--optional-- MPI communicator over atoms
563 : : !! mpi_enreg=information about MPI parallelization
564 : : !! natom=number of atoms in cell
565 : : !! nband_k=number of bands at this k point for that spin polarization
566 : : !! ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
567 : : !! npw_k=number of planewaves in basis sphere at k
568 : : !! npw1_k=number of planewaves in basis sphere at k+q
569 : : !! nspinor=number of spinorial components of the wavefunctions
570 : : !! occ_k(nband_k)=occupation number for each band (usually 2) for each k.
571 : : !! option= 1: accumulate 1st-order density,
572 : : !! 2: accumulate 0-order potential part of the 2nd-order total energy
573 : : !! 3: accumulate both
574 : : !! tim_fourwf= timing code for fourwf (5 from dfpt_vtowfk, 18 from dfpt_nstwf)
575 : : !! wf_corrected=flag put to 1 if cwave1 is different from cwavef (if there is a contribution from occ. change)
576 : : !! wtk_k=weight assigned to the k point.
577 : : !!
578 : : !! OUTPUT
579 : : !! ====== if option=2 or option=3 =====
580 : : !! eloc0_k=zero-order local contribution to 2nd-order total energy for current band and k
581 : : !!
582 : : !! SIDE EFFECTS
583 : : !! ====== if option=1 or option=3 =====
584 : : !! rhoaug1(cplex*n4,n5,n6,nvloc)= density in electrons/bohr**3,
585 : : !! ==== if gs_hamkq%usepaw=1 =====
586 : : !! pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data
587 : : !! (cumulative, so input as well as output)
588 : : !!
589 : : !! NOTES
590 : : !! In this part of the treatment of one band, one has to
591 : : !! perform Fourier transforms, and to treat separately the
592 : : !! two spinorial components of the wavefunction.
593 : : !! Was part of dfpt_vtowfk before.
594 : : !!
595 : : !! SOURCE
596 : :
597 : 7275230 : subroutine dfpt_accrho(cplex,cwave0,cwave1,cwavef,cwaveprj0,cwaveprj1,&
598 : 3637615 : & eloc0_k,gs_hamkq,iband,idir,ipert,isppol,kptopt,&
599 : 3637615 : & mpi_enreg,ndat,natom,nband_k,ncpgr,npw_k,npw1_k,nspinor,occ_k,&
600 [ - + ]: 3637615 : & option,pawrhoij1,rhoaug1,tim_fourwf,wf_corrected,&
601 [ - + - - : 3637615 : & wtk_k,comm_atom,mpi_atmtab)
- - ]
602 : :
603 : : !Arguments ------------------------------------
604 : : !scalars
605 : : integer,intent(in) :: cplex,iband,idir,ipert,isppol,kptopt,natom,nband_k,ndat
606 : : integer,intent(in) :: ncpgr,npw_k,npw1_k,nspinor,option,tim_fourwf,wf_corrected
607 : : integer,optional,intent(in) :: comm_atom
608 : : integer,optional,target,intent(in) :: mpi_atmtab(:)
609 : : real(dp),intent(in) :: wtk_k
610 : : real(dp),intent(out) :: eloc0_k(ndat)
611 : : type(gs_hamiltonian_type),intent(inout),target :: gs_hamkq
612 : : type(MPI_type),intent(in) :: mpi_enreg
613 : : !arrays
614 : : real(dp),intent(in),target :: cwave0(2,npw_k*nspinor*ndat),cwave1(2,npw1_k*nspinor*ndat),cwavef(2,npw1_k*nspinor*ndat)
615 : : real(dp),intent(in) :: occ_k(nband_k)
616 : : real(dp),intent(inout) :: rhoaug1(cplex*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,gs_hamkq%nvloc)
617 : : type(pawcprj_type),intent(in) :: cwaveprj0(natom,nspinor*ndat*gs_hamkq%usecprj)
618 : : type(pawcprj_type),intent(in) :: cwaveprj1(natom,nspinor*ndat*gs_hamkq%usepaw)
619 : : type(pawrhoij_type),intent(inout) :: pawrhoij1(:)
620 : :
621 : : !Local variables-------------------------------
622 : : !scalars
623 : : integer,parameter :: level=14
624 : : integer :: choice,cplex_cprj,i1,i2,i3,idat,ispinor,my_comm_atom,my_natom,n1,n2,n3,option_rhoij,gpu_option
625 : : logical :: my_atmtab_allocated,paral_atom
626 : : logical :: use_timerev,use_zeromag
627 : : real(dp) :: valuer,diag,offdiag,weight
628 : : real(dp) :: im0_up,im1_up,re0_up,re1_up,im0_down,im1_down,re0_down,re1_down
629 : : !arrays
630 : 3637615 : integer,pointer :: my_atmtab(:)
631 : : real(dp) :: dummy(2,1)
632 : 3637615 : real(dp),allocatable :: rhoaug(:,:,:,:),wfraug(:,:,:,:),wfraug1(:,:,:,:)
633 : 3637615 : real(dp),allocatable :: wfraug1_up(:,:,:,:),wfraug1_down(:,:,:,:)
634 : 3637615 : real(dp),allocatable :: wfraug_up(:,:,:,:),wfraug_down(:,:,:,:)
635 : 3637615 : real(dp),pointer :: cwavef_sp(:,:),cwavef_up(:,:),cwavef_down(:,:)
636 : 3637615 : real(dp),pointer :: cwave0_up(:,:),cwave0_down(:,:),cwave1_up(:,:),cwave1_down(:,:)
637 : : real(dp), contiguous, pointer :: vlocal(:,:,:,:)=>null()
638 : 3637615 : type(pawcprj_type),allocatable :: cwaveprj_tmp(:,:)
639 : :
640 : : ! *********************************************************************
641 : : DBG_ENTER("COLL")
642 : : ABI_NVTX_START_RANGE(NVTX_DFPT_ACCRHO)
643 : :
644 [ - + ]: 102656 : if (gs_hamkq%nvloc==4 .and. ndat>1) then
645 : 0 : ABI_ERROR("nvloc==4 isn't supported with ndat~bandpp > 1 !")
646 : : end if
647 : :
648 [ + - ]: 3637615 : if (option/=1.and.option/=2.and.option/=3) return
649 : 3637615 : gpu_option=gs_hamkq%gpu_option
650 : :
651 : : !Initializations
652 [ - + + - : 21825690 : ABI_MALLOC(rhoaug,(cplex*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,gs_hamkq%nvloc))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
653 : 3637615 : n1=gs_hamkq%ngfft(1);n2=gs_hamkq%ngfft(2);n3=gs_hamkq%ngfft(3)
654 [ + + + + ]: 7015452 : if (option==2.or.option==3) eloc0_k=zero
655 [ + + ]: 3637615 : if (option==2.or.option==3) then
656 : : ! XG20241028 This coding confuses the gnu 8.5 compiler, wrt the CONTIGUOUS character of the pointer.
657 : : ! So, make a simple copy. This is not optimal, though...
658 : : ! Original ---
659 : : ! vlocal => gs_hamkq%vlocal
660 : : ! New ---
661 [ - + + - : 20267022 : ABI_MALLOC(vlocal, (size(gs_hamkq%vlocal,1),size(gs_hamkq%vlocal,2),size(gs_hamkq%vlocal,3),size(gs_hamkq%vlocal,4)) )
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
662 [ + + + + : 17425826288 : vlocal=gs_hamkq%vlocal
+ + + + ]
663 : : ! END XG20241028
664 : : endif
665 : :
666 : : !Loop on spinorial components
667 : : ! TODO : double loop on spinors for full rhoaug1 matrix if nspden =4
668 [ + + ]: 3637615 : if (gs_hamkq%nvloc/=4) then ! see later EB FR
669 [ - + + - : 17674795 : ABI_MALLOC(wfraug1,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6*ndat))
- + + - -
+ + - + -
+ - - + -
+ ]
670 : : #ifdef HAVE_OPENMP_OFFLOAD
671 : : !$OMP TARGET ENTER DATA MAP(alloc:wfraug1) IF(gpu_option==ABI_GPU_OPENMP)
672 : : !$OMP TARGET ENTER DATA MAP(to:rhoaug1) IF(gpu_option==ABI_GPU_OPENMP)
673 : : #endif
674 : :
675 [ + + ]: 7106456 : do ispinor=1,nspinor
676 : :
677 : : ! Part devoted to the accumulation of the 0-order potential part of the 2nd-order total energy
678 : : ! --------------------------------------------------------------------------------------------
679 : :
680 : : ! Fourier transform of cwavef. Here, rhoaug is a dummy variable.
681 [ + + + + : 3571497 : if (wf_corrected==0.or.option==2.or.option==3) then
+ - ]
682 [ + + ]: 3571497 : if (ispinor==1) then
683 : 3534959 : cwavef_sp => cwavef(:,1:npw1_k*ndat)
684 : : else
685 : 36538 : cwavef_sp => cwavef(:,1+npw1_k*ndat:2*npw1_k*ndat)
686 : : end if
687 : : !make an inverse FFT from cwavef_sp to wfraug1
688 : : call fourwf(cplex,rhoaug,cwavef_sp,dummy,wfraug1,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,&
689 : : & gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,ndat,gs_hamkq%ngfft,&
690 : : & gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,&
691 [ + - - - : 3571497 : & weight,weight,gpu_option=gpu_option)
- - - - +
- - - - -
- - + - -
- - - - -
+ - - - -
- - - - +
- + - + -
+ - + ]
692 : 3571497 : nullify(cwavef_sp)
693 : :
694 : : end if
695 : :
696 : : ! Compute contribution of bands in ndat to zero-order potential part of the 2nd-order total energy
697 : : ! NB: this is spinor diagonal
698 [ + + ]: 3571497 : if (option==2.or.option==3) then
699 : : call update_potential_contrib(eloc0_k,wfraug1,vlocal,&
700 : 3311719 : & n1,n2,n3,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,gs_hamkq%nvloc,gs_hamkq%nfft,ndat,gpu_option)
701 : : end if ! option
702 : :
703 : : ! Part devoted to the accumulation of the 1st-order density
704 : : ! ---------------------------------------------------------
705 [ + + ]: 7106456 : if (option==1.or.option==3) then
706 : :
707 : : ! Compute 1st-order WF in real space
708 : : ! One needs the Fourier transform of cwave1. However, only the one of
709 : : ! cwavef is available. If cwavef and cwave1 differs, this Fourier
710 : : ! transform must be computed. In both case the result is in wfraug1.
711 [ + + ]: 2786812 : if (wf_corrected==1) then
712 [ + - ]: 566975 : if (ispinor==1) then
713 : 566975 : cwavef_sp => cwave1(:,1:npw1_k*ndat)
714 : : else
715 : 0 : cwavef_sp => cwave1(:,1+npw1_k*ndat:2*npw1_k*ndat)
716 : : end if
717 : : call fourwf(cplex,rhoaug,cwavef_sp,dummy,wfraug1,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,&
718 : : & gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,ndat,gs_hamkq%ngfft,&
719 : : & gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,&
720 [ + - - - : 566975 : & weight,weight,gpu_option=gpu_option)
- - - - +
- - - - -
- - + - -
- - - - -
+ - - - -
- - - - +
- + - + -
+ - + ]
721 : 566975 : nullify(cwavef_sp)
722 : : end if
723 : :
724 : : ! Compute 0-order WF in real space
725 : : ! TODO: add loop over ispinor_prime here
726 [ - + + - : 13934060 : ABI_MALLOC(wfraug,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6*ndat))
- + + - -
+ + - + -
+ - - + -
+ ]
727 : : #ifdef HAVE_OPENMP_OFFLOAD
728 : : !$OMP TARGET ENTER DATA MAP(alloc:wfraug) IF(gpu_option==ABI_GPU_OPENMP)
729 : : #endif
730 [ + + ]: 2786812 : if (ispinor==1) then
731 : 2779986 : cwavef_sp => cwave0(:,1:npw_k*ndat)
732 : : else
733 : 6826 : cwavef_sp => cwave0(:,1+npw_k*ndat:2*npw_k*ndat)
734 : : end if
735 : : call fourwf(1,rhoaug,cwavef_sp,dummy,wfraug,gs_hamkq%gbound_k,gs_hamkq%gbound_k,&
736 : : & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_k,gs_hamkq%mgfft,mpi_enreg,ndat,gs_hamkq%ngfft,&
737 : : & gs_hamkq%npw_k,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,&
738 [ + - - - : 2786812 : & weight,weight,gpu_option=gpu_option)
- - - - +
- - - - -
- - - + -
+ - + ]
739 : 2786812 : nullify(cwavef_sp)
740 : :
741 : : ! The factor 2 is not the spin factor (see Eq.44 of PRB55,10337 (1997) [[cite:Gonze1997]])
742 : : ! Accumulate 1st-order density
743 : :
744 : : call accumulate_1st_order_density(rhoaug1,wfraug,wfraug1,&
745 : : & cplex,n1,n2,n3,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,gs_hamkq%nvloc,&
746 : 2786812 : & ndat,nband_k,iband,gs_hamkq%ucvol,wtk_k,occ_k,gpu_option)
747 : :
748 : : #ifdef HAVE_OPENMP_OFFLOAD
749 : : !$OMP TARGET EXIT DATA MAP(delete:wfraug) IF(gpu_option==ABI_GPU_OPENMP)
750 : : #endif
751 [ - + ]: 2786812 : ABI_FREE(wfraug)
752 : : end if ! option
753 : :
754 : :
755 : : end do ! Loop on spinorial components if nspden=1 or 2
756 : :
757 : : #ifdef HAVE_OPENMP_OFFLOAD
758 : : !$OMP TARGET EXIT DATA MAP(delete:wfraug1) IF(gpu_option==ABI_GPU_OPENMP)
759 : : !$OMP TARGET EXIT DATA MAP(from:rhoaug1) IF(gpu_option==ABI_GPU_OPENMP)
760 : : #endif
761 [ - + ]: 3534959 : ABI_FREE(wfraug1)
762 : : else ! nvloc = 4
763 : : ! The same lines of code are in 72_response/dfpt_mkrho.F90
764 : : ! TODO merge these lines in a single routine??!!
765 [ - + + - : 513280 : ABI_MALLOC(wfraug1_up,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
- + + - -
+ + - + -
+ - - + -
+ ]
766 [ - + + - : 410624 : ABI_MALLOC(wfraug1_down,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
- + + - -
+ + - - +
- + ]
767 : :
768 [ + + + + : 1517385680 : wfraug1_up(:,:,:,:)=zero
+ + + + ]
769 [ + + + + : 1517385680 : wfraug1_down(:,:,:,:)=zero
+ + + + ]
770 : :
771 : : ! Part devoted to the accumulation of the 0-order potential part of the 2nd-order total energy
772 : : ! --------------------------------------------------------------------------------------------
773 : :
774 : : ! Fourier transform of cwavef. Here, rhoaug is a dummy variable.
775 [ + + + - : 102656 : if (wf_corrected==0.or.option==2.or.option==3) then
+ - ]
776 : 102656 : cwavef_up => cwavef(:,1:npw1_k) ! wfs up spin-polarized
777 : : call fourwf(cplex,rhoaug(:,:,:,1),cwavef_up,dummy,wfraug1_up,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,&
778 : : & gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
779 : : & gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,&
780 [ + - - - : 102656 : & weight,weight,gpu_option=gpu_option)
- - - - +
- - - - -
- - + - -
- - - - -
+ - - - -
- - - - +
- + - + -
+ - + ]
781 : 102656 : nullify(cwavef_up)
782 : :
783 : 102656 : cwavef_down => cwavef(:,1+npw1_k:2*npw1_k) ! wfs down spin-polarized
784 : : call fourwf(cplex,rhoaug(:,:,:,1),cwavef_down,dummy,wfraug1_down,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,&
785 : : & gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
786 : : & gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,&
787 [ + - - - : 102656 : & weight,weight,gpu_option=gpu_option)
- - - - +
- - - - -
- - + - -
- - - - -
+ - - - -
- - - - +
- + - + -
+ - + ]
788 : 102656 : nullify(cwavef_down)
789 : : end if
790 [ + - ]: 102656 : if (option==2.or.option==3) then
791 : : valuer=zero
792 : 1855760 : diag=zero
793 : 1855760 : offdiag=zero
794 : : ! EB FR 2nd term in Eq. 91 PRB52,1096 [[cite:Gonze1995]] for non-collinear magnetism
795 [ + + ]: 1855760 : do i3=1,n3
796 [ + + ]: 30152912 : do i2=1,n2
797 [ + + ]: 497562384 : do i1=1,n1
798 : : diag=vlocal(i1,i2,i3,1)*(wfraug1_up(1,i1,i2,i3)**2+wfraug1_up(2,i1,i2,i3)**2)&
799 : 467512128 : & +vlocal(i1,i2,i3,2)*(wfraug1_down(1,i1,i2,i3)**2+wfraug1_down(2,i1,i2,i3)**2)
800 : : offdiag=(two*vlocal(i1,i2,i3,3)*((wfraug1_up(1,i1,i2,i3)*wfraug1_down(1,i1,i2,i3))+&
801 : : & (wfraug1_up(2,i1,i2,i3)*wfraug1_down(2,i1,i2,i3))))+&
802 : : & (two*vlocal(i1,i2,i3,4)*((-wfraug1_down(2,i1,i2,i3)*wfraug1_up(1,i1,i2,i3))+&
803 : 467512128 : & (wfraug1_down(1,i1,i2,i3)*wfraug1_up(2,i1,i2,i3))))
804 : 495809280 : valuer=valuer+diag+offdiag
805 : : end do
806 : : end do
807 : : end do
808 : : ! Local potential energy of this band
809 [ + + ]: 205312 : eloc0_k=eloc0_k+two*valuer/dble(gs_hamkq%nfft)
810 : : end if ! option
811 : :
812 : : ! Part devoted to the accumulation of the 1st-order density
813 : : ! ---------------------------------------------------------
814 : :
815 : : ! first order
816 : : !
817 [ + + ]: 102656 : if (option==1.or.option==3) then
818 : :
819 : : !SPr: condition on wf_corrected not to do FFTs of the same Bloch functions
820 [ + + ]: 86528 : if (wf_corrected==1) then
821 : 6392 : cwave1_up => cwave1(:,1:npw1_k)
822 : 6392 : cwave1_down => cwave1(:,1+npw1_k:2*npw1_k)
823 [ + + + + : 291635000 : wfraug1_up(:,:,:,:)=zero
+ + + + ]
824 [ + + + + : 291635000 : wfraug1_down(:,:,:,:)=zero
+ + + + ]
825 : :
826 : : call fourwf(cplex,rhoaug(:,:,:,1),cwave1_up,dummy,wfraug1_up,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,&
827 : : & gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
828 : : & gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,&
829 [ + - - - : 6392 : & weight,weight,gpu_option=gpu_option)
- - - - +
- - - - -
- - + - -
- - - - -
+ - - - -
- - - - +
- + - + -
+ - + ]
830 : 6392 : nullify(cwave1_up)
831 : :
832 : : call fourwf(cplex,rhoaug(:,:,:,1),cwave1_down,dummy,wfraug1_down,gs_hamkq%gbound_kp,gs_hamkq%gbound_kp,&
833 : : & gs_hamkq%istwf_k,gs_hamkq%kg_kp,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
834 : : & gs_hamkq%npw_kp,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,&
835 [ + - - - : 6392 : & weight,weight,gpu_option=gpu_option)
- - - - +
- - - - -
- - + - -
- - - - -
+ - - - -
- - - - +
- + - + -
+ - + ]
836 : 6392 : nullify(cwave1_down)
837 : : end if
838 : :
839 : :
840 : : ! EB FR build spinorial wavefunctions
841 : : ! zero order
842 : 86528 : cwave0_up => cwave0(:,1:npw_k*ndat)
843 : 86528 : cwave0_down => cwave0(:,1+npw_k*ndat:2*npw_k*ndat)
844 [ - + + - : 432640 : ABI_MALLOC(wfraug_up,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
- + + - -
+ + - + -
+ - - + -
+ ]
845 [ - + + - : 346112 : ABI_MALLOC(wfraug_down,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
- + + - -
+ + - - +
- + ]
846 [ + + + + : 1275782192 : wfraug_up(:,:,:,:)=zero
+ + + + ]
847 [ + + + + : 1275782192 : wfraug_down(:,:,:,:)=zero
+ + + + ]
848 : : !
849 : : !density components
850 : : !GS wfk Fourrier Tranform
851 : : ! EB FR in the fourwf calls rhoaug(:,:,:,2) is a dummy argument
852 : : call fourwf(1,rhoaug(:,:,:,2),cwave0_up,dummy,wfraug_up,gs_hamkq%gbound_k,gs_hamkq%gbound_k,&
853 : : & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_k,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
854 : : & gs_hamkq%npw_k,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,&
855 [ + - - - : 86528 : & weight,weight,gpu_option=gpu_option)
- - - - +
- - - - -
- - - + -
+ - + ]
856 : 86528 : nullify(cwave0_up)
857 : : call fourwf(1,rhoaug(:,:,:,2),cwave0_down,dummy,wfraug_down,gs_hamkq%gbound_k,gs_hamkq%gbound_k,&
858 : : & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_k,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
859 : : & gs_hamkq%npw_k,1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,0,tim_fourwf,&
860 [ + - - - : 86528 : & weight,weight,gpu_option=gpu_option)
- - - - +
- - - - -
- - - + -
+ - + ]
861 : 86528 : nullify(cwave0_down)
862 : : ! Accumulate 1st-order density (x component)
863 : 86528 : re0_up=zero;im0_up=zero;re1_up=zero;im1_up=zero;re0_down=zero;im0_down=zero
864 : 86528 : re1_down=zero;im1_down=zero
865 : : ! The factor 2 is not the spin factor (see Eq.44 of PRB55,10337 (1997) [[cite:Gonze1997]])
866 : : ! SPr: the following treatment with factor=2 is ok for perturbations not breaking the
867 : : ! time reversal symmetry of the Hamiltonian (due to Kramer's degeneracy) hence
868 : : ! not applicable for magnetic field perturbation (for phonons with SOC, H^(0) has
869 : : ! time reversal symmetry though). The formulas below are rectified in dfpt_scfcv
870 : : ! in case of broken time-reversal upon reconstructing rhor1_pq and rhor1_mq.
871 : 86528 : weight=two*occ_k(iband)*wtk_k/gs_hamkq%ucvol
872 [ + + ]: 86528 : if (cplex==2) then
873 [ + + ]: 690008 : do i3=1,n3
874 [ + + ]: 12113240 : do i2=1,n2
875 [ + + ]: 224286528 : do i1=1,n1
876 : 212226048 : re0_up=wfraug_up(1,i1,i2,i3) ; im0_up=wfraug_up(2,i1,i2,i3)
877 : 212226048 : re1_up=wfraug1_up(1,i1,i2,i3) ; im1_up=wfraug1_up(2,i1,i2,i3)
878 : 212226048 : re0_down=wfraug_down(1,i1,i2,i3) ; im0_down=wfraug_down(2,i1,i2,i3)
879 : 212226048 : re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3)
880 : : !SPr: in case of +q/-q calculation, the factor will be corrected later from dfpt_scfcv level
881 : : ! along with the reconstruction of correct rhor1_{+q} and rhor1_{-q}
882 : : ! here, rhoaug1_{sigma,sigma'} = \sum_{n,k} u1_{sigma} u0*_{sigma'} independent of the sign of q
883 : 212226048 : rhoaug1(2*i1-1,i2,i3,1)=rhoaug1(2*i1-1,i2,i3,1)+weight*(re0_up*re1_up+im0_up*im1_up) !n_upup
884 : 212226048 : rhoaug1(2*i1 ,i2,i3,1)=rhoaug1(2*i1 ,i2,i3,1)+weight*(re0_up*im1_up-im0_up*re1_up)
885 : 212226048 : rhoaug1(2*i1-1,i2,i3,4)=rhoaug1(2*i1-1,i2,i3,4)+weight*(re0_down*re1_down+im0_down*im1_down) ! n_dndn
886 : 212226048 : rhoaug1(2*i1 ,i2,i3,4)=rhoaug1(2*i1 ,i2,i3,4)+weight*(re0_down*im1_down-im0_down*re1_down)
887 : :
888 : : rhoaug1(2*i1-1,i2,i3,2)=rhoaug1(2*i1-1,i2,i3,2)+weight*(re1_up*re0_down+im1_up*im0_down)& !Re[m1x]
889 : 212226048 : & +weight*(re1_down*re0_up+im1_down*im0_up)
890 : : rhoaug1(2*i1 ,i2,i3,2)=rhoaug1(2*i1 ,i2,i3,2)+weight*(-re1_up*im0_down+im1_up*re0_down)& !Im[m1x]
891 : 212226048 : & +weight*(-re1_down*im0_up+im1_down*re0_up)
892 : :
893 : : rhoaug1(2*i1-1,i2,i3,3)=rhoaug1(2*i1-1,i2,i3,3)+weight*(+re1_up*im0_down-im1_up*re0_down)& !Re[m1y]
894 : 212226048 : & +weight*(-re1_down*im0_up+im1_down*re0_up)
895 : : rhoaug1(2*i1 ,i2,i3,3)=rhoaug1(2*i1 ,i2,i3,3)+weight*(+re1_up*re0_down+im1_up*im0_down)& !Im[m1y]
896 : 223649280 : & +weight*(-re1_down*re0_up-im1_down*im0_up)
897 : : end do
898 : : end do
899 : : end do
900 : : else !cplex
901 : 816984 : re0_up=zero;im0_up=zero;re1_up=zero;im1_up=zero;re0_down=zero;im0_down=zero
902 : 816984 : re1_down=zero;im1_down=zero
903 [ + + ]: 816984 : do i3=1,n3
904 [ + + ]: 12640824 : do i2=1,n2
905 [ + + ]: 191174256 : do i1=1,n1
906 : 178567200 : re0_up=wfraug_up(1,i1,i2,i3) ; im0_up=wfraug_up(2,i1,i2,i3)
907 : 178567200 : re1_up=wfraug1_up(1,i1,i2,i3) ; im1_up=wfraug1_up(2,i1,i2,i3)
908 : 178567200 : re0_down=wfraug_down(1,i1,i2,i3) ; im0_down=wfraug_down(2,i1,i2,i3)
909 : 178567200 : re1_down=wfraug1_down(1,i1,i2,i3) ; im1_down=wfraug1_down(2,i1,i2,i3)
910 : :
911 : 178567200 : rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1)+weight*(re0_up*re1_up+im0_up*im1_up) ! n_upup
912 : 178567200 : rhoaug1(i1,i2,i3,4)=rhoaug1(i1,i2,i3,4)+weight*(re0_down*re1_down+im0_down*im1_down) ! n_dndn
913 : : rhoaug1(i1,i2,i3,2)=rhoaug1(i1,i2,i3,2)+weight*(re1_up*re0_down+re0_up*re1_down &
914 : 178567200 : & +im0_up*im1_down+im0_down*im1_up) !mx; the factor two is inside weight
915 : : rhoaug1(i1,i2,i3,3)=rhoaug1(i1,i2,i3,3)+weight*(re1_up*im0_down-im1_up*re0_down &
916 : 190391040 : & +re0_up*im1_down-im0_up*re1_down) !my; the factor two is inside weight
917 : : end do
918 : : end do
919 : : end do
920 : : end if !cplex
921 : :
922 [ - + ]: 86528 : ABI_FREE(wfraug_up)
923 [ - + ]: 86528 : ABI_FREE(wfraug_down)
924 : :
925 : : end if ! option
926 : :
927 [ - + ]: 102656 : ABI_FREE(wfraug1_up)
928 [ - + ]: 102656 : ABI_FREE(wfraug1_down)
929 : :
930 : : end if ! nvloc /= 4
931 : :
932 [ - + ]: 3637615 : ABI_FREE(rhoaug)
933 : :
934 : : !XG20241028 See above
935 : : !New ---
936 [ + + ]: 3637615 : if (option==2.or.option==3) then
937 [ - + ]: 3377837 : ABI_FREE(vlocal)
938 : : endif
939 : : !END XG20241028
940 : :
941 : : !Part devoted to the accumulation of the 1st-order occupation matrix in PAW case
942 : : ! TODO: parse for more nspden 4 dependencies on spinors
943 : : ! EB FR CHECK: to be modified for non-collinear?????
944 : : !-------------------------------------------------------------------------------
945 : :
946 [ + + + + : 10912845 : if ((option==1.or.option==3).and.gs_hamkq%usepaw==1) then
+ + - + -
- - - -
- ]
947 : :
948 : : ! Set up parallelism over atoms
949 : 635829 : my_natom=natom; if(gs_hamkq%usepaw==1) my_natom=size(pawrhoij1)
950 [ - + - - ]: 635829 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
951 [ - + ]: 635829 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
952 [ - + ]: 635829 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
953 : 635829 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
954 : :
955 [ - + ]: 635829 : cplex_cprj=2;if (gs_hamkq%istwf_k>1) cplex_cprj=1
956 : 635829 : option_rhoij=2
957 : 635829 : use_timerev=(kptopt>0.and.kptopt<3)
958 [ + - - + : 1271658 : use_zeromag=.false.;if (my_natom>0) use_zeromag=(pawrhoij1(1)%nspden==4.and.gs_hamkq%nvloc==1)
- - ]
959 : :
960 [ + - ]: 635829 : if (gs_hamkq%usecprj==1) then
961 [ + + ]: 1271658 : do idat=1,ndat
962 : : call pawaccrhoij(gs_hamkq%atindx,cplex_cprj,cwaveprj0(:,1+(idat-1)*nspinor:idat*nspinor),&
963 : : & cwaveprj1(:,1+(idat-1)*nspinor:idat*nspinor),ipert,isppol,my_natom,&
964 : : & natom,nspinor,occ_k(iband+idat-1),option_rhoij,pawrhoij1,use_timerev,use_zeromag,wtk_k,&
965 [ - + - - : 1271658 : & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
- - - + -
- - - - -
- - - - -
- - - ]
966 : : end do
967 : : else
968 : 0 : ABI_BUG("toto")
969 [ # # # # : 0 : ABI_MALLOC(cwaveprj_tmp,(natom,nspinor*ndat))
# # # # #
# # # # #
# # # # #
# # # ]
970 : 0 : call pawcprj_alloc(cwaveprj_tmp,ncpgr,gs_hamkq%dimcprj)
971 : 0 : choice=2
972 : : call getcprj(choice,0,cwave0,cwaveprj_tmp,&
973 : : & gs_hamkq%ffnl_k,idir,gs_hamkq%indlmn,gs_hamkq%istwf_k,&
974 : : & gs_hamkq%kg_k,gs_hamkq%kpg_k,gs_hamkq%kpt_k,gs_hamkq%lmnmax,&
975 : : & gs_hamkq%mgfft,mpi_enreg,ndat,gs_hamkq%natom,gs_hamkq%nattyp,gs_hamkq%ngfft,&
976 : : & gs_hamkq%nloalg,gs_hamkq%npw_k,gs_hamkq%nspinor,gs_hamkq%ntypat,gs_hamkq%phkxred,&
977 [ # # # # : 0 : & gs_hamkq%ph1d,gs_hamkq%ph3d_k,gs_hamkq%ucvol,gs_hamkq%useylm)
# # # # #
# ]
978 : : call pawaccrhoij(gs_hamkq%atindx,cplex_cprj,cwaveprj_tmp,cwaveprj1,ipert,isppol,my_natom,&
979 : : & gs_hamkq%natom,nspinor,occ_k(iband),option_rhoij,pawrhoij1,use_timerev,use_zeromag,wtk_k, &
980 [ # # # # : 0 : & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
# # # # #
# # # # #
# # # # #
# # # ]
981 : 0 : call pawcprj_free(cwaveprj_tmp)
982 [ # # # # : 0 : ABI_FREE(cwaveprj_tmp)
# # # # #
# ]
983 : : end if
984 : :
985 : : end if
986 : :
987 : : ABI_NVTX_END_RANGE()
988 : : DBG_EXIT("COLL")
989 : :
990 : : contains
991 : :
992 : : ! Compute contribution of bands in ndat to zero-order potential part of the 2nd-order total energy
993 : : ! This routine was separated from main code to accomodate OpenMP offloading with NVHPC
994 : 3311719 : subroutine update_potential_contrib(eloc0_k,wfraug1,vlocal,n1,n2,n3,n4,n5,n6,nvloc,nfft,ndat,gpu_option)
995 : : integer,intent(in) :: n1,n2,n3,n4,n5,n6,nvloc,ndat,nfft,gpu_option
996 : : real(dp),intent(in) :: wfraug1(2,n4,n5,n6*ndat),vlocal(n4,n5,n6,nvloc)
997 : : real(dp),intent(inout) :: eloc0_k(ndat)
998 : : integer :: idat,i1,i2,i3
999 : : real(dp) :: valuer,nfft_r
1000 : :
1001 : 3311719 : nfft_r=dble(nfft)
1002 : : #ifdef HAVE_OPENMP_OFFLOAD
1003 : : !$OMP TARGET TEAMS DISTRIBUTE &
1004 : : !$OMP& MAP(to:wfraug1,vlocal) MAP(tofrom:eloc0_k) PRIVATE(idat,valuer) &
1005 : : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1006 : : #endif
1007 [ + + ]: 6623438 : do idat=1,ndat
1008 : 3311719 : valuer=zero
1009 : : !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(i1,i2,i3) REDUCTION(+:valuer)
1010 [ + + ]: 54708635 : do i3=1,n3
1011 [ + + ]: 838348545 : do i2=1,n2
1012 [ + + ]: 14241614306 : do i1=1,n1
1013 : 14190217390 : valuer=valuer+vlocal(i1,i2,i3,1)*(wfraug1(1,i1,i2,i3+n3*(idat-1))**2+wfraug1(2,i1,i2,i3+n3*(idat-1))**2)
1014 : : end do
1015 : : end do
1016 : : end do
1017 : : ! Local potential energy of this band
1018 : 6623438 : eloc0_k(idat)=eloc0_k(idat)+two*valuer/nfft_r
1019 : : end do
1020 : :
1021 : : #ifndef HAVE_OPENMP_OFFLOAD
1022 : : ! Make testfarm happy
1023 : : ABI_UNUSED((/gpu_option/))
1024 : : #endif
1025 : 3311719 : end subroutine update_potential_contrib
1026 : :
1027 : : ! Accumulate 1st-order density
1028 : : ! This routine was separated from main code to accomodate OpenMP offloading with NVHPC
1029 : 2786812 : subroutine accumulate_1st_order_density(rhoaug1,wfraug,wfraug1,&
1030 : 2786812 : & cplex,n1,n2,n3,n4,n5,n6,nvloc,ndat,nband_k,iband,ucvol,wtk_k,occ_k,gpu_option)
1031 : : integer,intent(in) :: cplex,n1,n2,n3,n4,n5,n6,nvloc,ndat,nband_k,iband,gpu_option
1032 : : real(dp),intent(in) :: ucvol,wtk_k
1033 : : real(dp),intent(inout) :: rhoaug1(cplex*n4,n5,n6,nvloc)
1034 : : real(dp),intent(in) :: wfraug(2,n4,n5,n6*ndat),wfraug1(2,n4,n5,n6*ndat),occ_k(nband_k)
1035 : :
1036 : : integer :: i1,i2,i3,idat
1037 : : real(dp) :: weight
1038 : :
1039 [ + + ]: 2786812 : if (cplex==2) then
1040 [ + + ]: 3126744 : do idat=1,ndat
1041 : 1563372 : weight=two*occ_k(iband+idat-1)*wtk_k/ucvol
1042 : : #ifdef HAVE_OPENMP_OFFLOAD
1043 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(3) &
1044 : : !$OMP& MAP(to:wfraug1,wfraug,rhoaug1) PRIVATE(i3,i2,i1) &
1045 : : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1046 : : #endif
1047 [ + + ]: 26921180 : do i3=1,n3
1048 [ + + ]: 373375940 : do i2=1,n2
1049 [ + + ]: 5922995620 : do i1=1,n1
1050 : : ! TODO: check which terms (ispinor ispinorp) enter a given element of rhoaug1
1051 : : rhoaug1(2*i1-1,i2,i3,1)=rhoaug1(2*i1-1,i2,i3,1) &
1052 : : & +weight*(wfraug(1,i1,i2,i3+n6*(idat-1))*wfraug1(1,i1,i2,i3+n6*(idat-1)) &
1053 : 5551183052 : & +wfraug(2,i1,i2,i3+n6*(idat-1))*wfraug1(2,i1,i2,i3+n6*(idat-1)))
1054 : : rhoaug1(2*i1 ,i2,i3,1)=rhoaug1(2*i1 ,i2,i3,1) &
1055 : : & +weight*(wfraug(1,i1,i2,i3+n6*(idat-1))*wfraug1(2,i1,i2,i3+n6*(idat-1)) &
1056 : 5899201184 : & -wfraug(2,i1,i2,i3+n6*(idat-1))*wfraug1(1,i1,i2,i3+n6*(idat-1)))
1057 : : end do
1058 : : end do
1059 : : end do
1060 : : end do
1061 : : else
1062 [ + + ]: 2446880 : do idat=1,ndat
1063 : 1223440 : weight=two*occ_k(iband+idat-1)*wtk_k/ucvol
1064 : : #ifdef HAVE_OPENMP_OFFLOAD
1065 : : !$OMP TARGET TEAMS DISTRIBUTE PARALLEL DO COLLAPSE(3) &
1066 : : !$OMP& MAP(to:wfraug1,wfraug,rhoaug1) PRIVATE(i3,i2,i1) &
1067 : : !$OMP& IF(gpu_option==ABI_GPU_OPENMP)
1068 : : #endif
1069 [ + + ]: 22272349 : do i3=1,n3
1070 [ + + ]: 332062328 : do i2=1,n2
1071 [ + + ]: 5799978031 : do i1=1,n1
1072 : : rhoaug1(i1,i2,i3,1)=rhoaug1(i1,i2,i3,1) &
1073 : : & +weight*(wfraug(1,i1,i2,i3+n6*(idat-1))*wfraug1(1,i1,i2,i3+n6*(idat-1)) &
1074 : 5780152562 : & +wfraug(2,i1,i2,i3+n6*(idat-1))*wfraug1(2,i1,i2,i3+n6*(idat-1)))
1075 : : end do
1076 : : end do
1077 : : end do
1078 : : end do
1079 : : end if
1080 : :
1081 : : #ifndef HAVE_OPENMP_OFFLOAD
1082 : : ! Make testfarm happy
1083 : : ABI_UNUSED((/gpu_option/))
1084 : : #endif
1085 : 2786812 : end subroutine accumulate_1st_order_density
1086 : :
1087 : : end subroutine dfpt_accrho
1088 : : !!***
1089 : :
1090 : : end module m_dfpt_mkrho
1091 : : !!***
|
