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 7238222 : subroutine dfpt_accrho(cplex,cwave0,cwave1,cwavef,cwaveprj0,cwaveprj1,&
598 3619111 : & eloc0_k,gs_hamkq,iband,idir,ipert,isppol,kptopt,&
599 3619111 : & mpi_enreg,ndat,natom,nband_k,ncpgr,npw_k,npw1_k,nspinor,occ_k,&
600 3619111 : & option,pawrhoij1,rhoaug1,tim_fourwf,wf_corrected,&
601 3619111 : & 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 3619111 : integer,pointer :: my_atmtab(:)
631 : real(dp) :: dummy(2,1)
632 3619111 : real(dp),allocatable :: rhoaug(:,:,:,:),wfraug(:,:,:,:),wfraug1(:,:,:,:)
633 3619111 : real(dp),allocatable :: wfraug1_up(:,:,:,:),wfraug1_down(:,:,:,:)
634 3619111 : real(dp),allocatable :: wfraug_up(:,:,:,:),wfraug_down(:,:,:,:)
635 3619111 : real(dp),pointer :: cwavef_sp(:,:),cwavef_up(:,:),cwavef_down(:,:)
636 3619111 : real(dp),pointer :: cwave0_up(:,:),cwave0_down(:,:),cwave1_up(:,:),cwave1_down(:,:)
637 : real(dp), ABI_CONTIGUOUS pointer :: vlocal(:,:,:,:)=>null()
638 3619111 : 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 3619111 : if (option/=1.and.option/=2.and.option/=3) return
649 3619111 : gpu_option=gs_hamkq%gpu_option
650 :
651 : !Initializations
652 21714666 : ABI_MALLOC(rhoaug,(cplex*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,gs_hamkq%nvloc))
653 3619111 : n1=gs_hamkq%ngfft(1);n2=gs_hamkq%ngfft(2);n3=gs_hamkq%ngfft(3)
654 6985356 : if (option==2.or.option==3) eloc0_k=zero
655 3619111 : 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 20197470 : ABI_MALLOC(vlocal, (size(gs_hamkq%vlocal,1),size(gs_hamkq%vlocal,2),size(gs_hamkq%vlocal,3),size(gs_hamkq%vlocal,4)) )
662 16983811232 : 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 3619111 : if (gs_hamkq%nvloc/=4) then ! see later EB FR
669 17582275 : 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 7069448 : 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 3552993 : if (wf_corrected==0.or.option==2.or.option==3) then
682 3552993 : if (ispinor==1) then
683 3516455 : 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 3552993 : & weight,weight,gpu_option=gpu_option)
692 3552993 : 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 3552993 : if (option==2.or.option==3) then
699 : call update_potential_contrib(eloc0_k,wfraug1,vlocal,&
700 3300127 : & 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 7069448 : 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 2768380 : if (wf_corrected==1) then
712 555455 : if (ispinor==1) then
713 555455 : 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 555455 : & weight,weight,gpu_option=gpu_option)
721 555455 : nullify(cwavef_sp)
722 : end if
723 :
724 : ! Compute 0-order WF in real space
725 : ! TODO: add loop over ispinor_prime here
726 13841900 : 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 2768380 : if (ispinor==1) then
731 2761554 : 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 2768380 : & weight,weight,gpu_option=gpu_option)
739 2768380 : 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 2768380 : & 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 2768380 : 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 3516455 : 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 3619111 : ABI_FREE(rhoaug)
933 :
934 : !XG20241028 See above
935 : !New ---
936 3619111 : if (option==2.or.option==3) then
937 3366245 : 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 10857333 : if ((option==1.or.option==3).and.gs_hamkq%usepaw==1) then
947 :
948 : ! Set up parallelism over atoms
949 617397 : my_natom=natom; if(gs_hamkq%usepaw==1) my_natom=size(pawrhoij1)
950 617397 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
951 617397 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
952 617397 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
953 617397 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
954 :
955 617397 : cplex_cprj=2;if (gs_hamkq%istwf_k>1) cplex_cprj=1
956 617397 : option_rhoij=2
957 617397 : use_timerev=(kptopt>0.and.kptopt<3)
958 1234794 : use_zeromag=.false.;if (my_natom>0) use_zeromag=(pawrhoij1(1)%nspden==4.and.gs_hamkq%nvloc==1)
959 :
960 617397 : if (gs_hamkq%usecprj==1) then
961 1234794 : 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 1234794 : & 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 3300127 : 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 3300127 : 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 6600254 : do idat=1,ndat
1008 3300127 : valuer=zero
1009 : !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(i1,i2,i3) REDUCTION(+:valuer)
1010 54323219 : do i3=1,n3
1011 825793401 : do i2=1,n2
1012 13824709538 : do i1=1,n1
1013 13773686446 : 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 6600254 : 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 3300127 : 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 2768380 : subroutine accumulate_1st_order_density(rhoaug1,wfraug,wfraug1,&
1030 2768380 : & 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 2768380 : 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 2410016 : do idat=1,ndat
1063 1205008 : 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 21645661 : do i3=1,n3
1070 312579704 : do i2=1,n2
1071 5176534063 : 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 5157298418 : & +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 2768380 : end subroutine accumulate_1st_order_density
1086 :
1087 : end subroutine dfpt_accrho
1088 : !!***
1089 :
1090 : end module m_dfpt_mkrho
1091 : !!***
|
