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 134864 : 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 67432 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
111 67432 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
112 : real(dp) :: tsec(2)
113 67432 : real(dp),allocatable :: gcart1(:),gcart2(:),gcart3(:)
114 67432 : real(dp),allocatable :: g2cart1(:),g2cart2(:),g2cart3(:)
115 67432 : real(dp),allocatable :: ph1(:),ph2(:),ph3(:)
116 67432 : real(dp),allocatable :: wkcmpx(:,:),work(:),workgr(:,:),workgr2(:,:)
117 : ! *************************************************************************
118 :
119 : !DEBUG
120 : !write(std_out,*)' xcden : enter '
121 : !ENDDEBUG
122 :
123 67432 : 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 67432 : 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 67432 : need_laplacian = present(lrhonow)
134 67432 : need_derivative2 = present(d2rhonow)
135 :
136 : !Keep local copy of fft dimensions
137 67432 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
138 :
139 : !Initialize computation of G in cartesian coordinates
140 67432 : id1=n1/2+2 ; id2=n2/2+2 ; id3=n3/2+2
141 :
142 : !Get the distrib associated with this fft_grid
143 67432 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
144 :
145 : !Check whether q=0
146 67432 : qeq0=0
147 67432 : if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1
148 :
149 67432 : if(ishift==0)then
150 :
151 : ! Copy the input rhor in the new location. Will check on negative values later
152 :
153 136397 : do ispden=1,nspden
154 : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,nfft,rhonow,rhor)
155 1722545665 : do ifft=1,cplex*nfft
156 1722482436 : rhonow(ifft,ispden,1)=rhor(ifft,ispden)
157 : end do
158 : end do
159 :
160 : end if
161 :
162 67432 : if(ishift==1 .or. ngrad==2)then
163 :
164 72273 : ABI_MALLOC(work,(cplex*nfft))
165 72273 : ABI_MALLOC(wkcmpx,(2,nfft))
166 24091 : if(ngrad==2)then
167 40022 : ABI_MALLOC(workgr,(2,nfft))
168 7155319 : if (need_laplacian) lrhonow(:,:)=zero
169 20011 : if (need_laplacian.or.need_derivative2) then
170 274 : ABI_MALLOC(workgr2,(2,nfft))
171 : end if
172 60033 : ABI_MALLOC(gcart1,(n1))
173 60033 : ABI_MALLOC(gcart2,(n2))
174 60033 : ABI_MALLOC(gcart3,(n3))
175 20011 : 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 24091 : 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 52956 : 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 695456843 : do ifft=1,cplex*nfft
198 695456843 : work(ifft)=rhor(ifft,ispden)
199 : end do
200 :
201 28865 : call timab(82,1,tsec)
202 28865 : call fourdp(cplex,wkcmpx,work,-1,mpi_enreg,nfft,1,ngfft,0)
203 28865 : call timab(82,2,tsec)
204 :
205 : ! If shift is required, multiply now rho(G) by phase, then generate rho(r+delta)
206 28865 : 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 52956 : if(ngrad==2)then
246 :
247 : ! Need 3 derivatives for the gradient and the Laplacian
248 : ! Need 6 for the 2nd derivatives
249 23985 : ndir=3; if (need_derivative2) ndir=6
250 96372 : do ivoigt=1,ndir
251 72387 : idir=voigt1(ivoigt) ; jdir=voigt2(ivoigt)
252 :
253 5740302069 : workgr=zero
254 124252164 : if (need_laplacian.or.need_derivative2) workgr2=zero
255 :
256 1871967 : do i1=1,n1
257 1799580 : ig1=i1-(i1/id1)*n1-1
258 1871967 : 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 72387 : if(mod(n1,2)==0 .and. qeq0==1)gcart1(n1/2+1)=zero
262 1833363 : do i2=1,n2
263 1760976 : ig2=i2-(i2/id2)*n2-1
264 1833363 : gcart2(i2)=gprimd(idir,2)*two_pi*(dble(ig2)+qphon(2))
265 : end do
266 72387 : if(mod(n2,2)==0 .and. qeq0==1)gcart2(n2/2+1)=zero
267 1881123 : do i3=1,n3
268 1808736 : ig3=i3-(i3/id3)*n3-1
269 1881123 : gcart3(i3)=gprimd(idir,3)*two_pi*(dble(ig3)+qphon(3))
270 : end do
271 72387 : 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 72387 : 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 1881123 : do i3=1,n3
298 54347253 : do i2=1,n2
299 52466130 : gc23_idir=gcart2(i2)+gcart3(i3) ; gc23_jdir=gc23_idir
300 52466130 : if (idir/=jdir) gc23_jdir=g2cart2(i2)+g2cart3(i3)
301 54274866 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
302 1965582264 : do i1=1,n1
303 1913409894 : ifft=ifft+1
304 1913409894 : gcart_idir=gc23_idir+gcart1(i1) ; gcart_jdir=gcart_idir
305 1913409894 : if (idir/=jdir) gcart_jdir=gc23_jdir+g2cart1(i1)
306 : ! Multiply by i 2pi G(idir)
307 1913409894 : workgr(2,ifft)= gcart_idir*wkcmpx(1,ifft)
308 1913409894 : 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 1965582264 : 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 72387 : if (ivoigt<=3) then
321 71955 : call timab(82,1,tsec)
322 71955 : call fourdp(cplex,workgr,work,1,mpi_enreg,nfft,1,ngfft,0)
323 71955 : call timab(82,2,tsec)
324 : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(idir,ispden,cplex,nfft,rhonow,work)
325 1936293873 : do ifft=1,cplex*nfft
326 1936293873 : rhonow(ifft,ispden,1+idir)=work(ifft)
327 : end do
328 : end if
329 :
330 : ! Store/accumulate 2nd derivative or Laplacian of density
331 96372 : 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 24091 : ABI_FREE(work)
354 24091 : ABI_FREE(wkcmpx)
355 24091 : if (allocated(workgr)) then
356 20011 : ABI_FREE(workgr)
357 : end if
358 24091 : if (allocated(workgr2)) then
359 137 : ABI_FREE(workgr2)
360 : end if
361 24091 : if(ishift==1) then
362 4203 : ABI_FREE(ph1)
363 4203 : ABI_FREE(ph2)
364 4203 : ABI_FREE(ph3)
365 : end if
366 24091 : if(ngrad==2) then
367 20011 : ABI_FREE(gcart1)
368 20011 : ABI_FREE(gcart2)
369 20011 : ABI_FREE(gcart3)
370 20011 : 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 135136 : 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 65856 : subroutine xcpot (cplex,gprimd,ishift,use_laplacian,mpi_enreg,nfft,ngfft,ngrad,nspden,&
448 : & nspgrad,qphon,&
449 198234 : & 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 65856 : integer, contiguous, pointer :: fftn2_distrib(:),ffti2_local(:)
471 65856 : integer, contiguous, pointer :: fftn3_distrib(:),ffti3_local(:)
472 : logical :: with_vxc,with_vxctau
473 : real(dp) :: tsec(2)
474 65856 : real(dp),allocatable :: gcart1(:),gcart2(:),gcart3(:),ph1(:),ph2(:),ph3(:)
475 65856 : real(dp),allocatable :: wkcmpx(:,:),wkcmpxtau(:,:)
476 65856 : real(dp),allocatable :: work(:),workgr(:,:),worklp(:,:),worktau(:,:)
477 : ! *************************************************************************
478 :
479 65856 : 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 65856 : 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 65856 : with_vxc=present(vxc) ; with_vxctau=present(vxctau)
490 65856 : if (with_vxc) then
491 65856 : 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 65856 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
499 :
500 : !Initialize computation of G in cartesian coordinates
501 65856 : id1=n1/2+2 ; id2=n2/2+2 ; id3=n3/2+2
502 :
503 : !Get the distrib associated with this fft_grid
504 65856 : call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
505 :
506 : !Check whether q=0
507 65856 : qeq0=0;if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1
508 :
509 65856 : if(with_vxc.and.ishift==0)then ! Add the value of depsxc to vxc
510 132057 : do ispden=1,min(nspden,2)
511 : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,depsxc,nfft,vxc,ispden)
512 1673252232 : do ifft=1,cplex*nfft
513 1673190579 : 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 65856 : if(ishift==1 .or. ngrad==2)then
520 :
521 23344 : if(with_vxc.or.with_vxctau) then
522 70032 : ABI_MALLOC(work,(cplex*nfft))
523 : end if
524 23344 : if (with_vxc) then
525 70032 : ABI_MALLOC(wkcmpx,(2,nfft))
526 : end if
527 :
528 23344 : 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 50299 : do ispden=1,min(nspden,2)
540 :
541 : ! Initialize wkcmpx either to 0 or to the shifted vxc value
542 26955 : if (with_vxc) then
543 26955 : if(ishift==0)then
544 : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(nfft,wkcmpx)
545 611384502 : do ifft=1,nfft
546 1834109608 : 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 26955 : if (ngrad==2) then
563 66225 : ABI_MALLOC(gcart1,(n1))
564 66225 : ABI_MALLOC(gcart2,(n2))
565 66225 : ABI_MALLOC(gcart3,(n3))
566 22075 : if (with_vxc) then
567 66225 : ABI_MALLOC(workgr,(2,nfft))
568 22075 : if (use_laplacian==1) then
569 288 : ABI_MALLOC(worklp,(2,nfft))
570 : end if
571 : end if
572 22075 : if (with_vxctau) then
573 2073 : ABI_MALLOC(worktau,(2,nfft))
574 1382 : ABI_MALLOC(wkcmpxtau,(2,nfft))
575 : end if
576 :
577 88300 : do idir=1,3
578 :
579 66225 : if (with_vxc) then
580 : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,idir,ispden,nfft,rhonow,work)
581 1883808156 : do ifft=1,cplex*nfft
582 1883808156 : work(ifft)=rhonow(ifft,ispden,1+idir)
583 : end do
584 66225 : call timab(82,1,tsec)
585 66225 : call fourdp(cplex,workgr,work,-1,mpi_enreg,nfft,1,ngfft,0)
586 66225 : 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 66225 : 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 66225 : 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 1739310 : do i1=1,n1
616 1673085 : ig1=i1-(i1/id1)*n1-1
617 1739310 : 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 66225 : if(mod(n1,2)==0 .and. qeq0==1)gcart1(n1/2+1)=zero
621 1700718 : do i2=1,n2
622 1634493 : ig2=i2-(i2/id2)*n2-1
623 1700718 : gcart2(i2)=gprimd(idir,2)*two_pi*(dble(ig2)+qphon(2))
624 : end do
625 66225 : if(mod(n2,2)==0 .and. qeq0==1)gcart2(n2/2+1)=zero
626 1747686 : do i3=1,n3
627 1681461 : ig3=i3-(i3/id3)*n3-1
628 1747686 : gcart3(i3)=gprimd(idir,3)*two_pi*(dble(ig3)+qphon(3))
629 : end do
630 66225 : 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 1747686 : do i3=1,n3
636 51348699 : do i2=1,n2
637 49601013 : gc23_idir=gcart2(i2)+gcart3(i3)
638 51282474 : if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
639 1890249264 : do i1=1,n1
640 1840942011 : ifft=ifft+1
641 1840942011 : gcart_idir=gc23_idir+gcart1(i1)
642 1840942011 : if(with_vxc)then
643 : ! Multiply by - i 2pi G(idir) and accumulate in wkcmpx
644 1840942011 : wkcmpx(1,ifft)=wkcmpx(1,ifft)+gcart_idir*workgr(2,ifft)
645 1840942011 : wkcmpx(2,ifft)=wkcmpx(2,ifft)-gcart_idir*workgr(1,ifft)
646 1840942011 : 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 1890249264 : 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 88300 : 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 22075 : ABI_FREE(gcart1)
674 22075 : ABI_FREE(gcart2)
675 22075 : ABI_FREE(gcart3)
676 22075 : if (with_vxc) then
677 22075 : ABI_FREE(workgr)
678 22075 : if (use_laplacian==1) then
679 144 : ABI_FREE(worklp)
680 : end if
681 : end if
682 22075 : 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 50299 : if (with_vxc) then
693 26955 : 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 26955 : call timab(82,1,tsec)
723 26955 : call fourdp(cplex,wkcmpx,work,1,mpi_enreg,nfft,1,ngfft,0)
724 26955 : call timab(82,2,tsec)
725 : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,vxc,work)
726 677961604 : do ifft=1,cplex*nfft
727 677961604 : vxc(ifft,ispden)=vxc(ifft,ispden)+work(ifft)
728 : end do
729 : end if
730 :
731 : end do ! End loop on spins
732 :
733 23344 : if(ishift==1) then
734 4203 : ABI_FREE(ph1)
735 4203 : ABI_FREE(ph2)
736 4203 : ABI_FREE(ph3)
737 : end if
738 23344 : if(with_vxc) then
739 23344 : ABI_FREE(wkcmpx)
740 : end if
741 23344 : if(with_vxc.or.with_vxctau) then
742 23344 : ABI_FREE(work)
743 : end if
744 :
745 : end if ! End condition on ishift/ngrad
746 :
747 329946 : 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 : !!***
|
