Branch data Line data Source code
1 : : !!****m* ABINIT/m_xctk
2 : : !! NAME
3 : : !! m_xctk
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, DRH)
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 : : module m_xctk
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_errors
27 : :
28 : : use defs_abitypes, only : MPI_type
29 : : use m_time, only : timab
30 : : use m_mpinfo, only : ptabs_fourdp
31 : : use m_fft_mesh, only : phase
32 : : use m_fft, only : fourdp
33 : :
34 : : implicit none
35 : :
36 : : private
37 : : !!***
38 : :
39 : : public :: xcden
40 : : public :: xcpot
41 : : public :: xcpotdq
42 : : !!***
43 : :
44 : : contains
45 : : !!***
46 : :
47 : : !!****f* ABINIT/xcden
48 : : !! NAME
49 : : !! xcden
50 : : !!
51 : : !! FUNCTION
52 : : !! Prepare density data before calling local or semi-local xc evaluators.
53 : : !!
54 : : !! NOTES
55 : : !! Can take into account a shift of the grid, for purpose of better accuracy.
56 : : !! Can also compute the gradient of the density, for use with GGAs.
57 : : !! Also eliminate eventual negative densities.
58 : : !!
59 : : !! INPUTS
60 : : !! cplex=if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
61 : : !! gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1)
62 : : !! ishift : if ==0, do not shift the xc grid (usual case); if ==1, shift the xc grid
63 : : !! nfft=(effective) number of FFT grid points (for this processor)
64 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
65 : : !! ngrad : =1, only compute the density ; =2 also compute the
66 : : !! gradient of the density. Note : ngrad**2 is also used to dimension rhonow
67 : : !! nspden=number of spin-density components
68 : : !! qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
69 : : !! rhor(cplex*nfft,nspden)=real space electron density in electrons/bohr**3, on the
70 : : !! unshifted grid (total in first half and spin-up in second half if nspden=2)
71 : : !!
72 : : !! OUTPUT
73 : : !! rhonow(cplex*nfft,nspden,ngrad*ngrad)=electron (spin)-density in real space and
74 : : !! eventually its gradient, either on the unshifted grid (if ishift==0,
75 : : !! then equal to rhor),or on the shifted grid
76 : : !! rhonow(:,:,1)=electron density in electrons/bohr**3
77 : : !! if ngrad==2 : rhonow(:,:,2:4)=gradient of electron density in electrons/bohr**4
78 : : !! OPTIONAL OUTPUT
79 : : !! d2rhonow(cplex*nfft,nspden,6)=2nd derivatives of the electron (spin)-density in real space
80 : : !! in electrons/bohr**5 (in case of meta GGA) (Voigt notation)
81 : : !! lrhonow(cplex*nfft,nspden)=Laplacian of the electron (spin)-density in real space
82 : : !! in electrons/bohr**5 (in case of meta GGA)
83 : : !!
84 : : !! SOURCE
85 : :
86 : 136282 : subroutine xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden,qphon,rhor,rhonow, & !Mandatory arguments
87 : 272 : & d2rhonow,lrhonow) !Optional arguments
88 : :
89 : : !Arguments ------------------------------------
90 : : !scalars
91 : : integer,intent(in) :: cplex,ishift,nfft,ngrad,nspden
92 : : type(MPI_type),intent(in) :: mpi_enreg
93 : : !arrays
94 : : integer,intent(in) :: ngfft(18)
95 : : real(dp),intent(in) :: gprimd(3,3),qphon(3),rhor(cplex*nfft,nspden)
96 : : real(dp),intent(out) :: rhonow(cplex*nfft,nspden,ngrad*ngrad)
97 : : real(dp),intent(out),optional :: d2rhonow(cplex*nfft,nspden,6),lrhonow(cplex*nfft,nspden)
98 : :
99 : : !Local variables-------------------------------
100 : : !scalars
101 : : integer,parameter :: voigt1(6)=[1,2,3,3,3,2],voigt2(6)=[1,2,3,2,1,1]
102 : : integer :: i1,i2,i3,id1,id2,id3,idir,ifft,ig1,ig2,ig3
103 : : integer :: ispden,ivoigt,jdir,ndir,n1,n2,n3,qeq0
104 : : logical :: need_derivative2,need_laplacian
105 : : real(dp) :: gc23_idir,gc23_jdir,gcart_idir,gcart_jdir
106 : : real(dp) :: ph123i,ph123r,ph1i,ph1r,ph23i,ph23r,ph2i,ph2r,ph3i,ph3r
107 : : real(dp) :: work_im,work_re
108 : : character(len=500) :: message
109 : : !arrays
110 : 68141 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
111 : 68141 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
112 : : real(dp) :: tsec(2)
113 : 68141 : real(dp),allocatable :: gcart1(:),gcart2(:),gcart3(:)
114 : 68141 : real(dp),allocatable :: g2cart1(:),g2cart2(:),g2cart3(:)
115 : 68141 : real(dp),allocatable :: ph1(:),ph2(:),ph3(:)
116 : 68141 : real(dp),allocatable :: wkcmpx(:,:),work(:),workgr(:,:),workgr2(:,:)
117 : : ! *************************************************************************
118 : :
119 : : !DEBUG
120 : : !write(std_out,*)' xcden : enter '
121 : : !ENDDEBUG
122 : :
123 : 68141 : if (ishift/=0 .and. ishift/=1) then
124 : 0 : write(message, '(a,i0)' )'ishift must be 0 or 1 ; input was',ishift
125 : 0 : ABI_BUG(message)
126 : : end if
127 : :
128 [ - + ]: 68141 : if (ngrad/=1 .and. ngrad/=2) then
129 : 0 : write(message, '(a,i0)' )'ngrad must be 1 or 2 ; input was',ngrad
130 : 0 : ABI_BUG(message)
131 : : end if
132 : :
133 : 68141 : need_laplacian = present(lrhonow)
134 : 68141 : need_derivative2 = present(d2rhonow)
135 : :
136 : : !Keep local copy of fft dimensions
137 : 68141 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
138 : :
139 : : !Initialize computation of G in cartesian coordinates
140 : 68141 : id1=n1/2+2 ; id2=n2/2+2 ; id3=n3/2+2
141 : :
142 : : !Get the distrib associated with this fft_grid
143 : 68141 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
144 : :
145 : : !Check whether q=0
146 : 68141 : qeq0=0
147 [ + + ]: 68141 : if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1
148 : :
149 [ + + ]: 68141 : if(ishift==0)then
150 : :
151 : : ! Copy the input rhor in the new location. Will check on negative values later
152 : :
153 [ + + ]: 137989 : do ispden=1,nspden
154 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,nfft,rhonow,rhor)
155 [ + + ]: 1950359729 : do ifft=1,cplex*nfft
156 : 1950295791 : rhonow(ifft,ispden,1)=rhor(ifft,ispden)
157 : : end do
158 : : end do
159 : :
160 : : end if
161 : :
162 [ + + + + ]: 68141 : if(ishift==1 .or. ngrad==2)then
163 : :
164 [ - + + - : 73977 : ABI_MALLOC(work,(cplex*nfft))
+ - + - -
+ - + ]
165 [ - + + - : 73977 : ABI_MALLOC(wkcmpx,(2,nfft))
+ - + - -
+ - + ]
166 [ + + ]: 24659 : if(ngrad==2)then
167 [ - + + - : 41158 : ABI_MALLOC(workgr,(2,nfft))
- + - + ]
168 [ + + + + : 7155887 : if (need_laplacian) lrhonow(:,:)=zero
+ + ]
169 [ + + ]: 20579 : if (need_laplacian.or.need_derivative2) then
170 [ - + + - : 274 : ABI_MALLOC(workgr2,(2,nfft))
- + - + ]
171 : : end if
172 [ - + + - : 61737 : ABI_MALLOC(gcart1,(n1))
+ - + - -
+ - + ]
173 [ - + + - : 61737 : ABI_MALLOC(gcart2,(n2))
+ - + - -
+ - + ]
174 [ - + + - : 61737 : ABI_MALLOC(gcart3,(n3))
+ - + - -
+ - + ]
175 [ + + ]: 20579 : if (need_derivative2) then
176 [ - + + - : 270 : ABI_MALLOC(g2cart1,(n1))
- + - + ]
177 [ - + + - : 270 : ABI_MALLOC(g2cart2,(n2))
- + - + ]
178 [ - + + - : 270 : ABI_MALLOC(g2cart3,(n3))
- + - + ]
179 : : end if
180 : : end if
181 : :
182 [ + + ]: 24659 : if(ishift==1)then
183 : : ! Precompute phases (The phases correspond to a shift of density on real space
184 : : ! grid from center at 0 0 0 to (1/2)*(1/n1,1/n2,1/n3).)
185 [ + - - + : 12609 : ABI_MALLOC(ph1,(2*n1))
+ - + - +
- - + -
+ ]
186 [ + - - + : 12609 : ABI_MALLOC(ph2,(2*n2))
+ - + - +
- - + -
+ ]
187 [ + - - + : 12609 : ABI_MALLOC(ph3,(2*n3))
+ - + - +
- - + -
+ ]
188 : 4203 : call phase(n1,ph1)
189 : 4203 : call phase(n2,ph2)
190 : 4203 : call phase(n3,ph3)
191 : : end if
192 : :
193 [ + + ]: 54191 : do ispden=1,nspden
194 : :
195 : : ! Obtain rho(G) in wkcmpx from input rho(r)
196 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,nfft,rhor,work)
197 [ + + ]: 911391312 : do ifft=1,cplex*nfft
198 : 911391312 : work(ifft)=rhor(ifft,ispden)
199 : : end do
200 : :
201 : 29532 : call timab(82,1,tsec)
202 : 29532 : call fourdp(cplex,wkcmpx,work,-1,mpi_enreg,nfft,1,ngfft,0)
203 : 29532 : call timab(82,2,tsec)
204 : :
205 : : ! If shift is required, multiply now rho(G) by phase, then generate rho(r+delta)
206 [ + + ]: 29532 : if(ishift==1)then
207 : : !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,ph1i,ph1r,ph123i,ph123r,ph2i,ph2r,ph23i,ph23r,ph3i,ph3r,work_im,work_re) &
208 : : !$OMP&SHARED(n1,n2,n3,ph1,ph2,ph3,wkcmpx,mpi_enreg,fftn2_distrib)
209 [ + + ]: 96060 : do i3=1,n3
210 : 91054 : ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft)
211 : 91054 : ph3r=ph3(2*i3-1)
212 : 91054 : ph3i=ph3(2*i3 )
213 [ + + ]: 1895384 : do i2=1,n2
214 : 1799324 : ph2r=ph2(2*i2-1)
215 : 1799324 : ph2i=ph2(2*i2 )
216 : 1799324 : ph23r=ph2r*ph3r-ph2i*ph3i
217 : 1799324 : ph23i=ph2i*ph3r+ph2r*ph3i
218 [ + - ]: 1890378 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
219 [ + + ]: 54104780 : do i1=1,n1
220 : 52305456 : ifft=ifft+1
221 : 52305456 : ph1r=ph1(2*i1-1)
222 : 52305456 : ph1i=ph1(2*i1 )
223 : 52305456 : ph123r=ph1r*ph23r-ph1i*ph23i
224 : 52305456 : ph123i=ph1i*ph23r+ph1r*ph23i
225 : : ! Must use intermediate variables !
226 : 52305456 : work_re=ph123r*wkcmpx(1,ifft)-ph123i*wkcmpx(2,ifft)
227 : 52305456 : work_im=ph123i*wkcmpx(1,ifft)+ph123r*wkcmpx(2,ifft)
228 : 52305456 : wkcmpx(1,ifft)=work_re
229 : 54104780 : wkcmpx(2,ifft)=work_im
230 : : end do
231 : : end if
232 : : end do
233 : : end do
234 : 5006 : call timab(82,1,tsec)
235 : 5006 : call fourdp(cplex,wkcmpx,work,1,mpi_enreg,nfft,1,ngfft,0)
236 : 5006 : call timab(82,2,tsec)
237 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,rhonow,work)
238 [ + + ]: 52310462 : do ifft=1,cplex*nfft
239 : 52310462 : rhonow(ifft,ispden,1)=work(ifft)
240 : : end do
241 : : end if
242 : :
243 : : ! If gradient of the density is required, take care of the three components now
244 : : ! Note : this operation is applied on the eventually shifted rho(G)
245 [ + + ]: 54191 : if(ngrad==2)then
246 : :
247 : : ! Need 3 derivatives for the gradient and the Laplacian
248 : : ! Need 6 for the 2nd derivatives
249 [ + + ]: 24652 : ndir=3; if (need_derivative2) ndir=6
250 [ + + ]: 99040 : do ivoigt=1,ndir
251 : 74388 : idir=voigt1(ivoigt) ; jdir=voigt2(ivoigt)
252 : :
253 [ + + + + ]: 7683708288 : workgr=zero
254 [ + + + + : 124254165 : if (need_laplacian.or.need_derivative2) workgr2=zero
+ + ]
255 : :
256 [ + + ]: 1981542 : do i1=1,n1
257 : 1907154 : ig1=i1-(i1/id1)*n1-1
258 : 1981542 : gcart1(i1)=gprimd(idir,1)*two_pi*(dble(ig1)+qphon(1))
259 : : end do
260 : : ! Note that the G <-> -G symmetry must be maintained
261 [ + + + + ]: 74388 : if(mod(n1,2)==0 .and. qeq0==1)gcart1(n1/2+1)=zero
262 [ + + ]: 1942938 : do i2=1,n2
263 : 1868550 : ig2=i2-(i2/id2)*n2-1
264 : 1942938 : gcart2(i2)=gprimd(idir,2)*two_pi*(dble(ig2)+qphon(2))
265 : : end do
266 [ + + + + ]: 74388 : if(mod(n2,2)==0 .and. qeq0==1)gcart2(n2/2+1)=zero
267 [ + + ]: 1990698 : do i3=1,n3
268 : 1916310 : ig3=i3-(i3/id3)*n3-1
269 : 1990698 : gcart3(i3)=gprimd(idir,3)*two_pi*(dble(ig3)+qphon(3))
270 : : end do
271 [ + + + + ]: 74388 : if(mod(n3,2)==0 .and. qeq0==1)gcart3(n3/2+1)=zero
272 : :
273 : : !Need a second g-vector component for some 2nd derivatives
274 [ + + ]: 74388 : if (idir/=jdir) then
275 [ + + ]: 14538 : do i1=1,n1
276 : 14106 : ig1=i1-(i1/id1)*n1-1
277 : 14538 : g2cart1(i1)=gprimd(jdir,1)*two_pi*(dble(ig1)+qphon(1))
278 : : end do
279 : : ! Note that the G <-> -G symmetry must be maintained
280 [ + + + - ]: 432 : if(mod(n1,2)==0 .and. qeq0==1)g2cart1(n1/2+1)=zero
281 [ + + ]: 14538 : do i2=1,n2
282 : 14106 : ig2=i2-(i2/id2)*n2-1
283 : 14538 : g2cart2(i2)=gprimd(jdir,2)*two_pi*(dble(ig2)+qphon(2))
284 : : end do
285 [ + + + - ]: 432 : if(mod(n2,2)==0 .and. qeq0==1)g2cart2(n2/2+1)=zero
286 [ + + ]: 14934 : do i3=1,n3
287 : 14502 : ig3=i3-(i3/id3)*n3-1
288 : 14934 : g2cart3(i3)=gprimd(jdir,3)*two_pi*(dble(ig3)+qphon(3))
289 : : end do
290 [ + - + - ]: 432 : if(mod(n3,2)==0 .and. qeq0==1)g2cart3(n3/2+1)=zero
291 : : end if
292 : :
293 : : ! MG: Be careful here with OMP due to ifft. Disabled for the time being.
294 : : ! !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,gcart_idir,gcart_jdir,gc23_idir,gc23_jdir) &
295 : : ! !$OMP&SHARED(gcart1,gcart2,gcart3,g2cart1,g2cart2,g2cart3,n1,n2,n3,wkcmpx,workgr,workgr2)
296 : : ifft = 0
297 [ + + ]: 1990698 : do i3=1,n3
298 [ + + ]: 62980446 : do i2=1,n2
299 : 60989748 : gc23_idir=gcart2(i2)+gcart3(i3) ; gc23_jdir=gc23_idir
300 [ + + ]: 60989748 : if (idir/=jdir) gc23_jdir=g2cart2(i2)+g2cart3(i3)
301 [ + + ]: 62906058 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
302 [ + + ]: 2620407288 : do i1=1,n1
303 : 2561211300 : ifft=ifft+1
304 : 2561211300 : gcart_idir=gc23_idir+gcart1(i1) ; gcart_jdir=gcart_idir
305 [ + + ]: 2561211300 : if (idir/=jdir) gcart_jdir=gc23_jdir+g2cart1(i1)
306 : : ! Multiply by i 2pi G(idir)
307 : 2561211300 : workgr(2,ifft)= gcart_idir*wkcmpx(1,ifft)
308 : 2561211300 : workgr(1,ifft)=-gcart_idir*wkcmpx(2,ifft)
309 : : ! Do the same to the gradient in order to get the laplacian or the 2nd derivatives
310 [ + + ]: 2620407288 : if (need_laplacian.or.need_derivative2) then
311 : 41392968 : workgr2(2,ifft)= gcart_jdir*workgr(1,ifft)
312 : 41392968 : workgr2(1,ifft)=-gcart_jdir*workgr(2,ifft)
313 : : end if
314 : : end do
315 : : end if
316 : : end do
317 : : end do
318 : :
319 : : ! Store gradient of density
320 [ + + ]: 74388 : if (ivoigt<=3) then
321 : 73956 : call timab(82,1,tsec)
322 : 73956 : call fourdp(cplex,workgr,work,1,mpi_enreg,nfft,1,ngfft,0)
323 : 73956 : call timab(82,2,tsec)
324 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(idir,ispden,cplex,nfft,rhonow,work)
325 [ + + ]: 2584097280 : do ifft=1,cplex*nfft
326 : 2584097280 : rhonow(ifft,ispden,1+idir)=work(ifft)
327 : : end do
328 : : end if
329 : :
330 : : ! Store/accumulate 2nd derivative or Laplacian of density
331 [ + + ]: 99040 : if (need_laplacian.or.need_derivative2) then
332 : 873 : call timab(82,1,tsec)
333 : 873 : call fourdp(cplex,workgr2,work,1,mpi_enreg,nfft,1,ngfft,0)
334 : 873 : call timab(82,2,tsec)
335 [ + + ]: 873 : if (need_laplacian.and.ivoigt<=3) then
336 [ + + ]: 21405513 : do ifft=1,cplex*nfft
337 : 21405513 : lrhonow(ifft,ispden)=lrhonow(ifft,ispden)+work(ifft)
338 : : end do
339 : : end if
340 [ + + ]: 873 : if (need_derivative2) then
341 [ + + ]: 39976656 : do ifft=1,cplex*nfft
342 : 39976656 : d2rhonow(ifft,ispden,ivoigt)=work(ifft)
343 : : end do
344 : : end if
345 : : end if
346 : :
347 : : end do
348 : : end if
349 : :
350 : : end do ! End loop on spins
351 : :
352 : : ! Release memory
353 [ - + ]: 24659 : ABI_FREE(work)
354 [ - + ]: 24659 : ABI_FREE(wkcmpx)
355 [ + + ]: 24659 : if (allocated(workgr)) then
356 : 20579 : ABI_FREE(workgr)
357 : : end if
358 [ + + ]: 24659 : if (allocated(workgr2)) then
359 : 137 : ABI_FREE(workgr2)
360 : : end if
361 [ + + ]: 24659 : if(ishift==1) then
362 [ - + ]: 4203 : ABI_FREE(ph1)
363 [ - + ]: 4203 : ABI_FREE(ph2)
364 [ - + ]: 4203 : ABI_FREE(ph3)
365 : : end if
366 [ + + ]: 24659 : if(ngrad==2) then
367 [ - + ]: 20579 : ABI_FREE(gcart1)
368 [ - + ]: 20579 : ABI_FREE(gcart2)
369 [ - + ]: 20579 : ABI_FREE(gcart3)
370 [ + + ]: 20579 : if (need_derivative2) then
371 [ - + ]: 135 : ABI_FREE(g2cart1)
372 [ - + ]: 135 : ABI_FREE(g2cart2)
373 [ - + ]: 135 : ABI_FREE(g2cart3)
374 : : end if
375 : : end if
376 : :
377 : : end if ! End condition on ishift and ngrad
378 : :
379 [ + + + + : 136554 : end subroutine xcden
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
380 : : !!***
381 : :
382 : : !!****f* ABINIT/xcpot
383 : : !! NAME
384 : : !! xcpot
385 : : !!
386 : : !! FUNCTION
387 : : !! Process data after calling local or semi-local xc evaluators
388 : : !! The derivative of Exc with respect to the density, gradient of density
389 : : !! in case of GGAs, and Laplacian of density in case of meta-GGA
390 : : !! are contained in depsxc(:,:).
391 : : !! In case of GGAs (and meta-GGAs) the gradient of the density contained
392 : : !! in rhonow(:,:,2:4) is already multiplied by the local partial derivative
393 : : !! of the XC functional.
394 : : !! Can take into account a shift of the grid, for purpose of better accuracy
395 : : !!
396 : : !! INPUTS
397 : : !! cplex=if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
398 : : !! [depsxc(cplex*nfft,nspgrad)]=derivative of Exc with respect to the (spin-)density,
399 : : !! or to the norm of the gradient of the (spin-)density,
400 : : !! further divided by the norm of the gradient of the (spin-)density
401 : : !! The different components of depsxc will be
402 : : !! for nspden=1, depsxc(:,1)=d(rho.exc)/d(rho)
403 : : !! and if ngrad=2, depsxc(:,2)=1/2*1/|grad rho_up|*d(n.exc)/d(|grad rho_up|)
404 : : !! +1/|grad rho|*d(rho.exc)/d(|grad rho|)
405 : : !! and if use_laplacian=1, depsxc(:,3)=d(rho.exc)/d(lapl rho)
406 : : !! for nspden>=2, depsxc(:,1)=d(rho.exc)/d(rho_up)
407 : : !! depsxc(:,2)=d(rho.exc)/d(rho_down)
408 : : !! and if ngrad=2, depsxc(:,3)=1/|grad rho_up|*d(rho.exc)/d(|grad rho_up|)
409 : : !! depsxc(:,4)=1/|grad rho_down|*d(rho.exc)/d(|grad rho_down|)
410 : : !! depsxc(:,5)=1/|grad rho|*d(rho.exc)/d(|grad rho|)
411 : : !! and if use_laplacian=1, depsxc(:,6)=d(rho.exc)/d(lapl rho_up)
412 : : !! depsxc(:,7)=d(rho.exc)/d(lapl rho_down)
413 : : !! gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1)
414 : : !! ishift : if ==0, do not shift the xc grid (usual case);
415 : : !! if ==1, shift the xc grid
416 : : !! use_laplacian : 1 if we use a functional depending on the laplacian of the density
417 : : !! nfft=(effective) number of FFT grid points (for this processor)
418 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
419 : : !! ngrad : =1, only take into account derivative wrt the density ;
420 : : !! =2, also take into account derivative wrt the gradient of the density.
421 : : !! nspden=number of spin-density components
422 : : !! nspgrad=number of spin-density and spin-density-gradient components
423 : : !! qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
424 : : !! [rhonow(cplex*nfft,nspden,ngrad*ngrad)]=electron (spin)-density in real space and
425 : : !! eventually its gradient already multiplied by the local partial derivative
426 : : !! of the XC functional, either on the unshifted grid (if ishift==0,
427 : : !! then equal to rhor), or on the shifted grid
428 : : !! rhonow(:,:,1)=electron density in electrons/bohr**3
429 : : !! if ngrad==2 : rhonow(:,:,2:4)=gradient of electron density in el./bohr**4,
430 : : !! times local partial derivative of the functional, as required by the GGA
431 : : !! In this routine, rhonow is used only in the GGA case (ngrad=2).
432 : : !!
433 : : !! OUTPUT
434 : : !! (see side effects)
435 : : !!
436 : : !! SIDE EFFECTS
437 : : !! Input/Output (all optional:
438 : : !! [vxc(cplex*nfft,nspden)]=xc potential (spin up in first half and spin down in
439 : : !! second half if nspden>=2). Contribution from the present shifted
440 : : !! or unshifted grid is ADDED to the input vxc data.
441 : : !! [vxctau(cplex*nfft,nspden,4)]=derivative of XC energy density with respect to
442 : : !! kinetic energy density (depsxcdtau). The arrays vxctau(nfft,nspden,4) contains also
443 : : !! the gradient of vxctau (gvxctau) which will be computed here in vxctau(:,:,2:4).
444 : : !!
445 : : !! SOURCE
446 : :
447 : 66565 : subroutine xcpot (cplex,gprimd,ishift,use_laplacian,mpi_enreg,nfft,ngfft,ngrad,nspden,&
448 : : & nspgrad,qphon,&
449 : 200361 : & depsxc,rhonow,vxc,vxctau) ! optional argument
450 : :
451 : : !Arguments ------------------------------------
452 : : !scalars
453 : : integer,intent(in) :: cplex,ishift,nfft,ngrad,nspden,nspgrad,use_laplacian
454 : : type(MPI_type),intent(in) :: mpi_enreg
455 : : !arrays
456 : : integer,intent(in) :: ngfft(18)
457 : : real(dp),intent(in),optional :: rhonow(cplex*nfft,nspden,ngrad*ngrad)
458 : : real(dp),intent(in),optional :: depsxc(cplex*nfft,nspgrad),gprimd(3,3),qphon(3)
459 : : real(dp),intent(inout),optional :: vxc(cplex*nfft,nspden)
460 : : real(dp),intent(inout),optional :: vxctau(cplex*nfft,nspden,4)
461 : :
462 : : !Local variables-------------------------------
463 : : !scalars
464 : : integer :: i1,i2,i3,id1,id2,id3,idir,ifft,ig1,ig2,ig3,ispden,n1,n2,n3,qeq0
465 : : real(dp),parameter :: lowden=1.d-14,precis=1.d-15
466 : : real(dp) :: gc23_idir,gcart_idir,ph123i,ph123r,ph1i,ph1r,ph23i,ph23r,ph2i,ph2r
467 : : real(dp) :: ph3i,ph3r,work_im,work_re
468 : : character(len=500) :: message
469 : : !arrays
470 : 66565 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
471 : 66565 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
472 : : logical :: with_vxc,with_vxctau
473 : : real(dp) :: tsec(2)
474 : 66565 : real(dp),allocatable :: gcart1(:),gcart2(:),gcart3(:),ph1(:),ph2(:),ph3(:)
475 : 66565 : real(dp),allocatable :: wkcmpx(:,:),wkcmpxtau(:,:)
476 : 66565 : real(dp),allocatable :: work(:),workgr(:,:),worklp(:,:),worktau(:,:)
477 : : ! *************************************************************************
478 : :
479 : 66565 : if (ishift/=0 .and. ishift/=1) then
480 : 0 : write(message, '(a,i0)' )' ishift must be 0 or 1 ; input was',ishift
481 : 0 : ABI_BUG(message)
482 : : end if
483 : :
484 [ - + ]: 66565 : if (ngrad/=1 .and. ngrad/=2 ) then
485 : 0 : write(message, '(a,i0)' )' ngrad must be 1 or 2 ; input was',ngrad
486 : 0 : ABI_BUG(message)
487 : : end if
488 : :
489 : 66565 : with_vxc=present(vxc) ; with_vxctau=present(vxctau)
490 [ + - ]: 66565 : if (with_vxc) then
491 [ - + ]: 66565 : if ((.not.present(rhonow)).or.(.not.present(depsxc))) then
492 : 0 : message='need rhonow or depsxc!'
493 : 0 : ABI_BUG(message)
494 : : end if
495 : : end if
496 : :
497 : : !Keep local copy of fft dimensions
498 : 66565 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
499 : :
500 : : !Initialize computation of G in cartesian coordinates
501 : 66565 : id1=n1/2+2 ; id2=n2/2+2 ; id3=n3/2+2
502 : :
503 : : !Get the distrib associated with this fft_grid
504 : 66565 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
505 : :
506 : : !Check whether q=0
507 [ + + ]: 66565 : qeq0=0;if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1
508 : :
509 [ + - + + ]: 66565 : if(with_vxc.and.ishift==0)then ! Add the value of depsxc to vxc
510 [ + + ]: 133583 : do ispden=1,min(nspden,2)
511 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,depsxc,nfft,vxc,ispden)
512 [ + + ]: 1899843394 : do ifft=1,cplex*nfft
513 : 1899781032 : vxc(ifft,ispden)=vxc(ifft,ispden)+depsxc(ifft,ispden)
514 : : end do
515 : : end do
516 : : end if
517 : :
518 : : !If the grid is shifted, or if gradient corrections are present, there must be FFTs.
519 [ + + + + ]: 66565 : if(ishift==1 .or. ngrad==2)then
520 : :
521 [ + - ]: 23912 : if(with_vxc.or.with_vxctau) then
522 [ - + + - : 71736 : ABI_MALLOC(work,(cplex*nfft))
+ - + - -
+ - + ]
523 : : end if
524 [ + - ]: 23912 : if (with_vxc) then
525 [ - + + - : 71736 : ABI_MALLOC(wkcmpx,(2,nfft))
+ - + - -
+ - + ]
526 : : end if
527 : :
528 [ + + ]: 23912 : if(ishift==1)then
529 [ + - - + : 12609 : ABI_MALLOC(ph1,(2*n1))
+ - + - +
- - + -
+ ]
530 [ + - - + : 12609 : ABI_MALLOC(ph2,(2*n2))
+ - + - +
- - + -
+ ]
531 [ + - - + : 12609 : ABI_MALLOC(ph3,(2*n3))
+ - + - +
- - + -
+ ]
532 : : ! Precompute phases (The phases correspond to a shift of density on real space
533 : : ! grid from center at 0 0 0 to (1/2)*(1/n1,1/n2,1/n3).)
534 : 4203 : call phase(n1,ph1)
535 : 4203 : call phase(n2,ph2)
536 : 4203 : call phase(n3,ph3)
537 : : end if
538 : :
539 [ + + ]: 51468 : do ispden=1,min(nspden,2)
540 : :
541 : : ! Initialize wkcmpx either to 0 or to the shifted vxc value
542 [ + - ]: 27556 : if (with_vxc) then
543 [ + + ]: 27556 : if(ishift==0)then
544 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(nfft,wkcmpx)
545 [ + + ]: 826096069 : do ifft=1,nfft
546 [ + + ]: 2478243107 : wkcmpx(:,ifft)=zero
547 : : end do
548 : : else
549 : : ! Obtain depsxc(G)*phase in wkcmpx from input depsxc(r+delta)
550 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,depsxc,ispden,nfft,work)
551 [ + + ]: 52310462 : do ifft=1,cplex*nfft
552 : 52310462 : work(ifft)=depsxc(ifft,ispden)
553 : : end do
554 : 5006 : call timab(82,1,tsec)
555 : 5006 : call fourdp(cplex,wkcmpx,work,-1,mpi_enreg,nfft,1,ngfft,0)
556 : 5006 : call timab(82,2,tsec)
557 : : end if
558 : : end if
559 : :
560 : : ! If gradient correction is present, take care of the three components now
561 : : ! Note : this operation is done on the eventually shifted grid
562 [ + + ]: 27556 : if (ngrad==2) then
563 [ - + + - : 68028 : ABI_MALLOC(gcart1,(n1))
+ - + - -
+ - + ]
564 [ - + + - : 68028 : ABI_MALLOC(gcart2,(n2))
+ - + - -
+ - + ]
565 [ - + + - : 68028 : ABI_MALLOC(gcart3,(n3))
+ - + - -
+ - + ]
566 [ + - ]: 22676 : if (with_vxc) then
567 [ - + + - : 68028 : ABI_MALLOC(workgr,(2,nfft))
+ - + - -
+ - + -
+ ]
568 [ + + ]: 22676 : if (use_laplacian==1) then
569 [ - + + - : 288 : ABI_MALLOC(worklp,(2,nfft))
- + - + -
+ ]
570 : : end if
571 : : end if
572 [ + + ]: 22676 : if (with_vxctau) then
573 [ - + + - : 2073 : ABI_MALLOC(worktau,(2,nfft))
+ - + - -
+ - + -
+ ]
574 [ - + + - : 1382 : ABI_MALLOC(wkcmpxtau,(2,nfft))
- + - + -
+ ]
575 : : end if
576 : :
577 [ + + ]: 90704 : do idir=1,3
578 : :
579 [ + - ]: 68028 : if (with_vxc) then
580 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,idir,ispden,nfft,rhonow,work)
581 [ + + ]: 2527942857 : do ifft=1,cplex*nfft
582 : 2527942857 : work(ifft)=rhonow(ifft,ispden,1+idir)
583 : : end do
584 : 68028 : call timab(82,1,tsec)
585 : 68028 : call fourdp(cplex,workgr,work,-1,mpi_enreg,nfft,1,ngfft,0)
586 : 68028 : call timab(82,2,tsec)
587 : :
588 : : ! IF Meta-GGA then take care of the laplacian term involved.
589 : : ! Note : this operation is done on the eventually shifted grid
590 [ + + ]: 68028 : if(use_laplacian==1)then
591 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,idir,ispden,nspden,nfft,depsxc,work)
592 [ + + ]: 19988328 : do ifft=1,cplex*nfft
593 [ + + ]: 19988328 : if(nspden==1)then
594 : 17468472 : work(ifft)=depsxc(ifft,2+ispden)
595 [ + - ]: 2519424 : else if(nspden==2)then
596 : 2519424 : work(ifft)=depsxc(ifft,5+ispden)
597 : : end if
598 : : end do
599 : 432 : call timab(82,1,tsec)
600 : 432 : call fourdp(cplex,worklp,work,-1,mpi_enreg,nfft,1,ngfft,0)
601 : 432 : call timab(82,2,tsec)
602 : : end if
603 : : end if
604 : :
605 [ + + ]: 68028 : if(with_vxctau)then
606 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,vxctau,work)
607 [ + + ]: 113359377 : do ifft=1,cplex*nfft
608 : 113359377 : work(ifft)=vxctau(ifft,ispden,1)
609 : : end do
610 : 2073 : call timab(82,1,tsec)
611 : 2073 : call fourdp(cplex,worktau,work,-1,mpi_enreg,nfft,1,ngfft,0)
612 : 2073 : call timab(82,2,tsec)
613 : : end if ! present vxctau
614 : :
615 [ + + ]: 1843395 : do i1=1,n1
616 : 1775367 : ig1=i1-(i1/id1)*n1-1
617 : 1843395 : gcart1(i1)=gprimd(idir,1)*two_pi*(dble(ig1)+qphon(1))
618 : : end do
619 : : ! Note that the G <-> -G symmetry must be maintained
620 [ + + + + ]: 68028 : if(mod(n1,2)==0 .and. qeq0==1)gcart1(n1/2+1)=zero
621 [ + + ]: 1804803 : do i2=1,n2
622 : 1736775 : ig2=i2-(i2/id2)*n2-1
623 : 1804803 : gcart2(i2)=gprimd(idir,2)*two_pi*(dble(ig2)+qphon(2))
624 : : end do
625 [ + + + + ]: 68028 : if(mod(n2,2)==0 .and. qeq0==1)gcart2(n2/2+1)=zero
626 [ + + ]: 1851771 : do i3=1,n3
627 : 1783743 : ig3=i3-(i3/id3)*n3-1
628 : 1851771 : gcart3(i3)=gprimd(idir,3)*two_pi*(dble(ig3)+qphon(3))
629 : : end do
630 [ + + + + ]: 68028 : if(mod(n3,2)==0 .and. qeq0==1)gcart3(n3/2+1)=zero
631 : :
632 : : ! !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,gc23_idir,gcart_idir) &
633 : : ! !$OMP&SHARED(gcart1,gcart2,gcart3,n1,n2,n3,wkcmpx,workgr)
634 : : ifft = 0
635 [ + + ]: 1851771 : do i3=1,n3
636 [ + + ]: 59836398 : do i2=1,n2
637 : 57984627 : gc23_idir=gcart2(i2)+gcart3(i3)
638 [ + + ]: 59768370 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
639 [ + + ]: 2541265776 : do i1=1,n1
640 : 2485074909 : ifft=ifft+1
641 : 2485074909 : gcart_idir=gc23_idir+gcart1(i1)
642 [ + - ]: 2485074909 : if(with_vxc)then
643 : : ! Multiply by - i 2pi G(idir) and accumulate in wkcmpx
644 : 2485074909 : wkcmpx(1,ifft)=wkcmpx(1,ifft)+gcart_idir*workgr(2,ifft)
645 : 2485074909 : wkcmpx(2,ifft)=wkcmpx(2,ifft)-gcart_idir*workgr(1,ifft)
646 [ + + ]: 2485074909 : if(use_laplacian==1)then
647 : : ! Multiply by - i 2pi G(idir) and accumulate in wkcmpx
648 : 19987896 : wkcmpx(1,ifft)=wkcmpx(1,ifft)-gcart_idir**2*worklp(1,ifft)
649 : 19987896 : wkcmpx(2,ifft)=wkcmpx(2,ifft)-gcart_idir**2*worklp(2,ifft)
650 : : end if
651 : : end if
652 [ + + ]: 2541265776 : if(with_vxctau)then
653 : 113357304 : wkcmpxtau(1,ifft)= gcart_idir*worktau(2,ifft)
654 : 113357304 : wkcmpxtau(2,ifft)=-gcart_idir*worktau(1,ifft)
655 : : end if
656 : : end do
657 : : end if
658 : : end do
659 : : end do
660 : :
661 [ + + ]: 90704 : if (with_vxctau) then
662 : 2073 : call timab(82,1,tsec)
663 : 2073 : call fourdp(cplex,wkcmpxtau,work,1,mpi_enreg,nfft,1,ngfft,0)
664 : 2073 : call timab(82,2,tsec)
665 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,vxctau,work)
666 [ + + ]: 113359377 : do ifft=1,cplex*nfft
667 : 113359377 : vxctau(ifft,ispden,1+idir)=work(ifft)
668 : : end do
669 : : end if
670 : :
671 : : end do ! enddo idir
672 : :
673 [ - + ]: 22676 : ABI_FREE(gcart1)
674 [ - + ]: 22676 : ABI_FREE(gcart2)
675 [ - + ]: 22676 : ABI_FREE(gcart3)
676 [ + - ]: 22676 : if (with_vxc) then
677 [ - + ]: 22676 : ABI_FREE(workgr)
678 [ + + ]: 22676 : if (use_laplacian==1) then
679 [ - + ]: 144 : ABI_FREE(worklp)
680 : : end if
681 : : end if
682 [ + + ]: 22676 : if (with_vxctau) then
683 [ - + ]: 691 : ABI_FREE(worktau)
684 [ - + ]: 691 : ABI_FREE(wkcmpxtau)
685 : : end if
686 : :
687 : : end if
688 : :
689 : : ! wkcmpx(:,:) contains now the full exchange-correlation potential, but
690 : : ! eventually for the shifted grid
691 : :
692 [ + - ]: 51468 : if (with_vxc) then
693 [ + + ]: 27556 : if(ishift==1)then
694 : : ! Take away the phase to get depsxc(G)
695 : : ifft=0
696 [ + + ]: 96060 : do i3=1,n3
697 : 91054 : ph3r=ph3(2*i3-1)
698 : 91054 : ph3i=ph3(2*i3 )
699 [ + + ]: 1895384 : do i2=1,n2
700 : 1799324 : ph2r=ph2(2*i2-1)
701 : 1799324 : ph2i=ph2(2*i2 )
702 : 1799324 : ph23r=ph2r*ph3r-ph2i*ph3i
703 : 1799324 : ph23i=ph2i*ph3r+ph2r*ph3i
704 [ + - ]: 1890378 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
705 [ + + ]: 54104780 : do i1=1,n1
706 : 52305456 : ifft=ifft+1
707 : 52305456 : ph1r=ph1(2*i1-1)
708 : 52305456 : ph1i=ph1(2*i1 )
709 : 52305456 : ph123r=ph1r*ph23r-ph1i*ph23i
710 : 52305456 : ph123i=ph1i*ph23r+ph1r*ph23i
711 : : ! Multiply by phase. Must use intermediate variables !
712 : 52305456 : work_re= ph123r*wkcmpx(1,ifft)+ph123i*wkcmpx(2,ifft)
713 : 52305456 : work_im=-ph123i*wkcmpx(1,ifft)+ph123r*wkcmpx(2,ifft)
714 : 52305456 : wkcmpx(1,ifft)=work_re
715 : 54104780 : wkcmpx(2,ifft)=work_im
716 : : end do
717 : : end if
718 : : end do
719 : : end do
720 : : end if
721 : :
722 : 27556 : call timab(82,1,tsec)
723 : 27556 : call fourdp(cplex,wkcmpx,work,1,mpi_enreg,nfft,1,ngfft,0)
724 : 27556 : call timab(82,2,tsec)
725 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,vxc,work)
726 [ + + ]: 892673171 : do ifft=1,cplex*nfft
727 : 892673171 : vxc(ifft,ispden)=vxc(ifft,ispden)+work(ifft)
728 : : end do
729 : : end if
730 : :
731 : : end do ! End loop on spins
732 : :
733 [ + + ]: 23912 : if(ishift==1) then
734 [ - + ]: 4203 : ABI_FREE(ph1)
735 [ - + ]: 4203 : ABI_FREE(ph2)
736 [ - + ]: 4203 : ABI_FREE(ph3)
737 : : end if
738 [ + - ]: 23912 : if(with_vxc) then
739 [ - + ]: 23912 : ABI_FREE(wkcmpx)
740 : : end if
741 [ + - ]: 23912 : if(with_vxc.or.with_vxctau) then
742 [ - + ]: 23912 : ABI_FREE(work)
743 : : end if
744 : :
745 : : end if ! End condition on ishift/ngrad
746 : :
747 [ + - + + : 333491 : end subroutine xcpot
+ - + - -
+ - + - +
- + - + -
+ - + - +
- + - + ]
748 : : !!***
749 : :
750 : : !!****f* ABINIT/xcpotdq
751 : : !! NAME
752 : : !! xcpotdq
753 : : !!
754 : : !! FUNCTION
755 : : !! Equivalent to xcpot for the q-derivative of the GGA xc kernel.
756 : : !!
757 : : !! INPUTS
758 : : !! agradn(cplex*nfft,nspgrad,3)=kxc(:,4)*gradrho(:,:)*gradrho(:,qdir)*rho1(*)
759 : : !! cplex=if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
760 : : !! gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1)
761 : : !! ishift : if ==0, do not shift the xc grid (usual case);
762 : : !! if ==1, shift the xc grid (not implemented)
763 : : !! nfft=(effective) number of FFT grid points (for this processor)
764 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
765 : : !! ngrad : =1, only take into account derivative wrt the density ;
766 : : !! =2, also take into account derivative wrt the gradient of the density.
767 : : !! nspden=number of spin-density components
768 : : !! nspgrad=number of spin-density and spin-density-gradient components
769 : : !!
770 : : !! OUTPUT
771 : : !! vxc(cplex*nfft,nspden)]=q-derivative of the GGA xc potential.
772 : : !! At input already includes three terms.
773 : : !!
774 : : !! SOURCE
775 : :
776 : 912 : subroutine xcpotdq (agradn,cplex,gprimd,ishift,mpi_enreg, &
777 : 912 : & nfft,ngfft,ngrad,nspden,nspgrad,vxc)
778 : :
779 : : !Arguments ------------------------------------
780 : : !scalars
781 : : integer,intent(in) :: cplex,ishift,nfft,ngrad,nspden,nspgrad
782 : : type(MPI_type),intent(in) :: mpi_enreg
783 : : !arrays
784 : : integer,intent(in) :: ngfft(18)
785 : : real(dp),intent(in) :: agradn(cplex*nfft,nspgrad,3)
786 : : real(dp),intent(in) :: gprimd(3,3)
787 : : real(dp),intent(inout) :: vxc(2*nfft,nspden)
788 : :
789 : : !Local variables-------------------------------
790 : : !scalars
791 : : integer :: i1,i2,i3,id1,id2,id3,idir,ifft,ig1,ig2,ig3,ispden,n1,n2,n3
792 : : real(dp),parameter :: lowden=1.d-14,precis=1.d-15
793 : : real(dp) :: gc23_idir,gcart_idir
794 : : character(len=500) :: message
795 : : !arrays
796 : 912 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
797 : 912 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
798 : : real(dp) :: tsec(2)
799 : 912 : real(dp),allocatable :: gcart1(:),gcart2(:),gcart3(:)
800 : 912 : real(dp),allocatable :: wkcmpx(:,:)
801 : 912 : real(dp),allocatable :: work(:),workgr(:,:)
802 : : ! *************************************************************************
803 : :
804 : 912 : if (ishift/=0) then
805 : 0 : write(message, '(a,i0)' )' ishift must be 0 ; input was',ishift
806 : 0 : ABI_BUG(message)
807 : : end if
808 : :
809 [ - + ]: 912 : if (ngrad/=2) then
810 : 0 : write(message, '(a,i0)' )' ngrad must be 2 ; input was',ngrad
811 : 0 : ABI_BUG(message)
812 : : end if
813 : :
814 : : !Keep local copy of fft dimensions
815 : 912 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
816 : :
817 : : !Initialize computation of G in cartesian coordinates
818 : 912 : id1=n1/2+2 ; id2=n2/2+2 ; id3=n3/2+2
819 : :
820 : : !Get the distrib associated with this fft_grid
821 : 912 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
822 : :
823 : : !Compute the real-space gradient of de second term
824 [ - + + - : 2736 : ABI_MALLOC(work,(cplex*nfft))
+ - + - -
+ - + ]
825 [ - + + - : 2736 : ABI_MALLOC(wkcmpx,(2,nfft))
+ - + - -
+ - + ]
826 [ - + + - : 1824 : ABI_MALLOC(workgr,(2,nfft))
- + + - ]
827 : :
828 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(nfft,wkcmpx)
829 [ + + ]: 3182736 : do ifft=1,nfft
830 [ + + ]: 9546384 : wkcmpx(:,ifft)=zero
831 : : end do
832 : :
833 : : ! Obtain agradn(G)*phase in wkcmpx from input agradn(r)
834 : 912 : ispden=1
835 [ - + + - : 2736 : ABI_MALLOC(gcart1,(n1))
+ - + - -
+ - + ]
836 [ - + + - : 2736 : ABI_MALLOC(gcart2,(n2))
+ - + - -
+ - + ]
837 [ - + + - : 2736 : ABI_MALLOC(gcart3,(n3))
+ - + - -
+ + - ]
838 [ + + ]: 3648 : do idir=1, 3
839 : :
840 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,idir,agradn,ispden,nfft,work)
841 [ + + ]: 9548208 : do ifft=1,cplex*nfft
842 : 9548208 : work(ifft)=agradn(ifft,ispden,idir)
843 : : end do
844 : 2736 : call timab(82,1,tsec)
845 : 2736 : call fourdp(cplex,workgr,work,-1,mpi_enreg,nfft,1,ngfft,0)
846 : 2736 : call timab(82,2,tsec)
847 : :
848 [ + + ]: 44208 : do i1=1,n1
849 : 41472 : ig1=i1-(i1/id1)*n1-1
850 : 44208 : gcart1(i1)=gprimd(idir,1)*two_pi*dble(ig1)
851 : : end do
852 : : !Note that the G <-> -G symmetry must be maintained
853 [ + + ]: 2736 : if(mod(n1,2)==0) gcart1(n1/2+1)=zero
854 [ + + ]: 44208 : do i2=1,n2
855 : 41472 : ig2=i2-(i2/id2)*n2-1
856 : 44208 : gcart2(i2)=gprimd(idir,2)*two_pi*dble(ig2)
857 : : end do
858 [ + + ]: 2736 : if(mod(n2,2)==0) gcart2(n2/2+1)=zero
859 [ + + ]: 44208 : do i3=1,n3
860 : 41472 : ig3=i3-(i3/id3)*n3-1
861 : 44208 : gcart3(i3)=gprimd(idir,3)*two_pi*dble(ig3)
862 : : end do
863 [ + + ]: 2736 : if(mod(n3,2)==0) gcart3(n3/2+1)=zero
864 : :
865 : : ! !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,gc23_idir,gcart_idir) &
866 : : ! !$OMP&SHARED(gcart1,gcart2,gcart3,n1,n2,n3,wkcmpx,workgr)
867 : : ifft = 0
868 [ + + ]: 45120 : do i3=1,n3
869 [ + + ]: 673200 : do i2=1,n2
870 : 628992 : gc23_idir=gcart2(i2)+gcart3(i3)
871 [ + - ]: 670464 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
872 [ + + ]: 10174464 : do i1=1,n1
873 : 9545472 : ifft=ifft+1
874 : 9545472 : gcart_idir=gc23_idir+gcart1(i1)
875 : : ! Multiply by -i 2pi G(idir) and accumulate in wkcmpx
876 : 9545472 : wkcmpx(1,ifft)=wkcmpx(1,ifft)+gcart_idir*workgr(2,ifft)
877 : 10174464 : wkcmpx(2,ifft)=wkcmpx(2,ifft)-gcart_idir*workgr(1,ifft)
878 : : end do
879 : : end if
880 : : end do
881 : : end do
882 : :
883 : : end do
884 : :
885 [ - + ]: 912 : ABI_FREE(gcart1)
886 [ - + ]: 912 : ABI_FREE(gcart2)
887 [ - + ]: 912 : ABI_FREE(gcart3)
888 [ - + ]: 912 : ABI_FREE(workgr)
889 : :
890 : 912 : call timab(82,1,tsec)
891 : 912 : call fourdp(cplex,wkcmpx,work,1,mpi_enreg,nfft,1,ngfft,0)
892 : 912 : call timab(82,2,tsec)
893 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,vxc,work)
894 [ + + ]: 3182736 : do ifft=1,nfft
895 : 3181824 : vxc(2*ifft,ispden)=vxc(2*ifft,ispden)+work(ifft)
896 : : !Apply here the two pi factor
897 : 3182736 : vxc(2*ifft,ispden)=vxc(2*ifft,ispden)*two_pi
898 : : end do
899 [ - + ]: 912 : ABI_FREE(wkcmpx)
900 [ - + ]: 912 : ABI_FREE(work)
901 : :
902 [ - + ]: 912 : end subroutine xcpotdq
903 : : !!***
904 : :
905 : : end module m_xctk
906 : : !!***
|
