TABLE OF CONTENTS


ABINIT/m_xctk [ Modules ]

[ Top ] [ Modules ]

NAME

  m_xctk

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, DRH)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_xctk
28 
29  use defs_basis
30  use defs_abitypes
31  use m_abicore
32  use m_errors
33 
34  use m_time,     only : timab
35  use m_mpinfo,   only : ptabs_fourdp
36  use m_fft_mesh, only : phase
37  use m_fft,      only : fourdp
38 
39  implicit none
40 
41  private

ABINIT/xcden [ Functions ]

[ Top ] [ Functions ]

NAME

 xcden

FUNCTION

 Prepare density data before calling local or semi-local xc evaluators.

NOTES

 Can take into account a shift of the grid, for purpose of better accuracy.
 Can also compute the gradient of the density, for use with GGAs.
 Also eliminate eventual negative densities.

INPUTS

  cplex=if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
  gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1)
  ishift : if ==0, do not shift the xc grid (usual case); if ==1, shift the xc grid
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  ngrad : =1, only compute the density ; =2 also compute the
      gradient of the density. Note : ngrad**2 is also used to dimension rhonow
  nspden=number of spin-density components
  qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
  rhor(cplex*nfft,nspden)=real space electron density in electrons/bohr**3, on the
   unshifted grid (total in first half and spin-up in second half if nspden=2)

OUTPUT

  rhonow(cplex*nfft,nspden,ngrad*ngrad)=electron (spin)-density in real space and
     eventually its gradient, either on the unshifted grid (if ishift==0,
     then equal to rhor),or on the shifted grid
    rhonow(:,:,1)=electron density in electrons/bohr**3
    if ngrad==2 : rhonow(:,:,2:4)=gradient of electron density in electrons/bohr**4
  OPTIONAL OUTPUT
  lrhonow(cplex*nfft,nspden)=Laplacian of the electron (spin)-density in real space
    in electrons/bohr**5 (in case of meta GGA)

PARENTS

      afterscfloop,dfpt_mkvxcgga,dfpt_mkvxcstrgga,gammapositron_fft
      rhohxcpositron,rhotoxc

CHILDREN

      fourdp,phase,ptabs_fourdp,timab

SOURCE

 95 subroutine xcden (cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden,paral_kgb,qphon,rhor,rhonow, & !Mandatory arguments
 96 &  lrhonow)              !Optional arguments
 97 
 98 
 99 !This section has been created automatically by the script Abilint (TD).
100 !Do not modify the following lines by hand.
101 #undef ABI_FUNC
102 #define ABI_FUNC 'xcden'
103 !End of the abilint section
104 
105  implicit none
106 
107 !Arguments ------------------------------------
108 !scalars
109  integer,intent(in) :: cplex,ishift,nfft,ngrad,nspden,paral_kgb
110  type(MPI_type),intent(in) :: mpi_enreg
111 !arrays
112  integer,intent(in) :: ngfft(18)
113  real(dp),intent(in) :: gprimd(3,3),qphon(3),rhor(cplex*nfft,nspden)
114  real(dp),intent(out) :: rhonow(cplex*nfft,nspden,ngrad*ngrad)
115  real(dp),intent(out),optional :: lrhonow(cplex*nfft,nspden)
116 
117 !Local variables-------------------------------
118 !scalars
119  integer :: i1,i2,i3,id1,id2,id3,idir,ifft,ig1,ig2,ig3
120  integer :: ispden,n1,n2,n3,qeq0
121  real(dp) :: gc23_idir,gcart_idir,ph123i,ph123r,ph1i,ph1r,ph23i,ph23r,ph2i,ph2r
122  real(dp) :: ph3i,ph3r,work_im,work_re
123  character(len=500) :: message
124 !arrays
125  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
126  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
127  real(dp) :: tsec(2)
128  real(dp),allocatable :: gcart1(:),gcart2(:),gcart3(:),ph1(:),ph2(:),ph3(:)
129  real(dp),allocatable :: wkcmpx(:,:),work(:),workgr(:,:),worklp(:,:)
130 
131 ! *************************************************************************
132 
133 !DEBUG
134 !write(std_out,*)' xcden : enter '
135 !ENDDEBUG
136 
137  if (ishift/=0 .and. ishift/=1) then
138    write(message, '(a,i0)' )'ishift must be 0 or 1 ; input was',ishift
139    MSG_BUG(message)
140  end if
141 
142  if (ngrad/=1 .and. ngrad/=2) then
143    write(message, '(a,i0)' )'ngrad must be 1 or 2 ; input was',ngrad
144    MSG_BUG(message)
145  end if
146 
147 !Keep local copy of fft dimensions
148  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
149 
150 !Initialize computation of G in cartesian coordinates
151  id1=n1/2+2  ; id2=n2/2+2  ; id3=n3/2+2
152 
153 !Get the distrib associated with this fft_grid
154  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
155 
156 !Check whether q=0
157  qeq0=0
158  if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1
159 
160  if(ishift==0)then
161 
162 !  Copy the input rhor in the new location. Will check on negative values later
163 
164    do ispden=1,nspden
165 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,nfft,rhonow,rhor)
166      do ifft=1,cplex*nfft
167        rhonow(ifft,ispden,1)=rhor(ifft,ispden)
168      end do
169    end do
170 
171  end if
172 
173  if(ishift==1 .or. ngrad==2)then
174 
175    ABI_ALLOCATE(wkcmpx,(2,nfft))
176    ABI_ALLOCATE(work,(cplex*nfft))
177 
178    if(ishift==1)then
179 !    Precompute phases (The phases correspond to a shift of density on real space
180 !    grid from center at 0 0 0 to (1/2)*(1/n1,1/n2,1/n3).)
181      ABI_ALLOCATE(ph1,(2*n1))
182      ABI_ALLOCATE(ph2,(2*n2))
183      ABI_ALLOCATE(ph3,(2*n3))
184      call phase(n1,ph1)
185      call phase(n2,ph2)
186      call phase(n3,ph3)
187    end if
188 
189    do ispden=1,nspden
190 
191 !    Obtain rho(G) in wkcmpx from input rho(r)
192 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,nfft,rhor,work)
193      do ifft=1,cplex*nfft
194        work(ifft)=rhor(ifft,ispden)
195      end do
196 
197      call timab(82,1,tsec)
198      call fourdp(cplex,wkcmpx,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
199      call timab(82,2,tsec)
200 
201 !    If shift is required, multiply now rho(G) by phase, then generate rho(r+delta)
202      if(ishift==1)then
203 !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,ph1i,ph1r,ph123i,ph123r,ph2i,ph2r,ph23i,ph23r,ph3i,ph3r,work_im,work_re) &
204 !$OMP&SHARED(n1,n2,n3,ph1,ph2,ph3,wkcmpx,mpi_enreg,fftn2_distrib)
205        do i3=1,n3
206          ifft=(i3-1)*n1*(n2/mpi_enreg%nproc_fft)
207          ph3r=ph3(2*i3-1)
208          ph3i=ph3(2*i3  )
209          do i2=1,n2
210            ph2r=ph2(2*i2-1)
211            ph2i=ph2(2*i2  )
212            ph23r=ph2r*ph3r-ph2i*ph3i
213            ph23i=ph2i*ph3r+ph2r*ph3i
214            if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
215              do i1=1,n1
216                ifft=ifft+1
217                ph1r=ph1(2*i1-1)
218                ph1i=ph1(2*i1  )
219                ph123r=ph1r*ph23r-ph1i*ph23i
220                ph123i=ph1i*ph23r+ph1r*ph23i
221 !              Must use intermediate variables !
222                work_re=ph123r*wkcmpx(1,ifft)-ph123i*wkcmpx(2,ifft)
223                work_im=ph123i*wkcmpx(1,ifft)+ph123r*wkcmpx(2,ifft)
224                wkcmpx(1,ifft)=work_re
225                wkcmpx(2,ifft)=work_im
226              end do
227            end if
228          end do
229        end do
230        call timab(82,1,tsec)
231        call fourdp(cplex,wkcmpx,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
232        call timab(82,2,tsec)
233 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(ispden,nfft,rhonow,work)
234        do ifft=1,cplex*nfft
235          rhonow(ifft,ispden,1)=work(ifft)
236        end do
237      end if
238 
239 !    If gradient of the density is required, take care of the three components now
240 !    Note : this operation is applied on the eventually shifted rho(G)
241      if(ngrad==2)then
242        ABI_ALLOCATE(gcart1,(n1))
243        ABI_ALLOCATE(gcart2,(n2))
244        ABI_ALLOCATE(gcart3,(n3))
245        ABI_ALLOCATE(workgr,(2,nfft))
246        if (present(lrhonow)) then
247          ABI_ALLOCATE(worklp,(2,nfft))
248          lrhonow(:,ispden)=zero
249        end if
250        do idir=1,3
251 
252          do i1=1,n1
253            ig1=i1-(i1/id1)*n1-1
254            gcart1(i1)=gprimd(idir,1)*two_pi*(dble(ig1)+qphon(1))
255          end do
256 !        Note that the G <-> -G symmetry must be maintained
257          if(mod(n1,2)==0 .and. qeq0==1)gcart1(n1/2+1)=zero
258          do i2=1,n2
259            ig2=i2-(i2/id2)*n2-1
260            gcart2(i2)=gprimd(idir,2)*two_pi*(dble(ig2)+qphon(2))
261          end do
262          if(mod(n2,2)==0 .and. qeq0==1)gcart2(n2/2+1)=zero
263          do i3=1,n3
264            ig3=i3-(i3/id3)*n3-1
265            gcart3(i3)=gprimd(idir,3)*two_pi*(dble(ig3)+qphon(3))
266          end do
267          if(mod(n3,2)==0 .and. qeq0==1)gcart3(n3/2+1)=zero
268 
269 !        MG: Be careful here with OMP due to ifft. Disabled for the time being.
270 !        !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,gcart_idir,gc23_idir) &
271 !        !$OMP&SHARED(gcart1,gcart2,gcart3,n1,n2,n3,wkcmpx,workgr)
272          ifft = 0
273          do i3=1,n3
274            do i2=1,n2
275              gc23_idir=gcart2(i2)+gcart3(i3)
276              if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
277                do i1=1,n1
278                  ifft=ifft+1
279                  gcart_idir=gc23_idir+gcart1(i1)
280 !                Multiply by i 2pi G(idir)
281                  workgr(2,ifft)= gcart_idir*wkcmpx(1,ifft)
282                  workgr(1,ifft)=-gcart_idir*wkcmpx(2,ifft)
283 !                Do the same to the gradient in order to get the laplacian
284                  if (present(lrhonow)) then
285                    worklp(2,ifft)= gcart_idir*workgr(1,ifft)
286                    worklp(1,ifft)=-gcart_idir*workgr(2,ifft)
287                  end if
288                end do
289              end if
290            end do
291          end do
292          call timab(82,1,tsec)
293          call fourdp(cplex,workgr,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
294          call timab(82,2,tsec)
295 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(idir,ispden,cplex,nfft,rhonow,work)
296          do ifft=1,cplex*nfft
297            rhonow(ifft,ispden,1+idir)=work(ifft)
298          end do
299 
300          if (present(lrhonow)) then
301            call fourdp(cplex,worklp,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
302            do ifft=1,cplex*nfft
303              lrhonow(ifft,ispden)=lrhonow(ifft,ispden)+work(ifft)
304            end do
305          end if
306 
307        end do
308        ABI_DEALLOCATE(gcart1)
309        ABI_DEALLOCATE(gcart2)
310        ABI_DEALLOCATE(gcart3)
311        ABI_DEALLOCATE(workgr)
312        if (allocated(worklp))  then
313          ABI_DEALLOCATE(worklp)
314        end if
315      end if
316 
317    end do  ! End loop on spins
318 
319    ABI_DEALLOCATE(wkcmpx)
320    ABI_DEALLOCATE(work)
321    if(ishift==1) then
322      ABI_DEALLOCATE(ph1)
323      ABI_DEALLOCATE(ph2)
324      ABI_DEALLOCATE(ph3)
325    end if
326 
327  end if  ! End condition on ishift and ngrad
328 
329 end subroutine xcden

ABINIT/xcpot [ Functions ]

[ Top ] [ Functions ]

NAME

 xcpot

FUNCTION

 Process data after calling local or semi-local xc evaluators
 The derivative of Exc with respect to the density, gradient of density
 in case of GGAs, and Laplacian of density in case of meta-GGA
 are contained in depsxc(:,:).
 In case of GGAs (and meta-GGAs) the gradient of the density contained
 in rhonow(:,:,2:4) is already multiplied by the local partial derivative
 of the XC functional.
 Can take into account a shift of the grid, for purpose of better accuracy

INPUTS

  cplex=if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
  depsxc(cplex*nfft,nspgrad)=derivative of Exc with respect to the (spin-)density,
    or to the norm of the gradient of the (spin-)density,
    further divided by the norm of the gradient of the (spin-)density
   The different components of depsxc will be
   for nspden=1,         depsxc(:,1)=d(rho.exc)/d(rho)
     and if ngrad=2,     depsxc(:,2)=1/2*1/|grad rho_up|*d(n.exc)/d(|grad rho_up|)
                                     +1/|grad rho|*d(rho.exc)/d(|grad rho|)
     and if mgga=1,      depsxc(:,3)=d(rho.exc)/d(lapl rho)
   for nspden>=2,        depsxc(:,1)=d(rho.exc)/d(rho_up)
                         depsxc(:,2)=d(rho.exc)/d(rho_down)
     and if ngrad=2,     depsxc(:,3)=1/|grad rho_up|*d(rho.exc)/d(|grad rho_up|)
                         depsxc(:,4)=1/|grad rho_down|*d(rho.exc)/d(|grad rho_down|)
                         depsxc(:,5)=1/|grad rho|*d(rho.exc)/d(|grad rho|)
     and if mgga=1,      depsxc(:,6)=d(rho.exc)/d(lapl rho_up)
                         depsxc(:,7)=d(rho.exc)/d(lapl rho_down)
  gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1)
  ishift : if ==0, do not shift the xc grid (usual case);
           if ==1, shift the xc grid
  mgga : 1 if we use a meta-GGA functional.
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  ngrad : =1, only take into account derivative wrt the density ;
          =2, also take into account derivative wrt the gradient of the density.
  nspden=number of spin-density components
  nspgrad=number of spin-density and spin-density-gradient components
  qphon(3)=reduced coordinates for the phonon wavelength (needed if cplex==2).
  rhonow(cplex*nfft,nspden,ngrad*ngrad)=electron (spin)-density in real space and
     eventually its gradient already multiplied by the local partial derivative
     of the XC functional, either on the unshifted grid (if ishift==0,
     then equal to rhor), or on the shifted grid
    rhonow(:,:,1)=electron density in electrons/bohr**3
    if ngrad==2 : rhonow(:,:,2:4)=gradient of electron density in el./bohr**4,
     times local partial derivative of the functional, as required by the GGA
    In this routine, rhonow is used only in the GGA case (ngrad=2).

OUTPUT

  (see side effects)

SIDE EFFECTS

 Input/Output:
  vxc(cplex*nfft,nspden)=xc potential (spin up in first half and spin down in
   second half if nspden>=2). Contribution from the present shifted
   or unshifted grid is ADDED to the input vxc data.
  vxctau(cplex*nfft,nspden,4)=derivative of XC energy density with respect to
   kinetic energy density (depsxcdtau). The arrays vxctau(nfft,nspden,4) contains also
   the gradient of vxctau (gvxctau) which will be computed here in vxctau(:,:,2:4).

PARENTS

      dfpt_mkvxcgga,dfpt_mkvxcstrgga,rhotoxc

CHILDREN

      fourdp,phase,ptabs_fourdp,timab

SOURCE

403 subroutine xcpot (cplex,depsxc,gprimd,ishift,mgga,mpi_enreg,nfft,ngfft,ngrad,nspden,&
404 & nspgrad,paral_kgb,qphon,rhonow,vxc,&
405 & vxctau) ! optional argument
406 
407 
408 !This section has been created automatically by the script Abilint (TD).
409 !Do not modify the following lines by hand.
410 #undef ABI_FUNC
411 #define ABI_FUNC 'xcpot'
412 !End of the abilint section
413 
414  implicit none
415 
416 !Arguments ------------------------------------
417 !scalars
418  integer,intent(in) :: cplex,ishift,mgga,nfft,ngrad,nspden,nspgrad,paral_kgb
419  type(MPI_type),intent(in) :: mpi_enreg
420 !arrays
421  integer,intent(in) :: ngfft(18)
422  real(dp),intent(in) :: depsxc(cplex*nfft,nspgrad),gprimd(3,3),qphon(3)
423  real(dp),intent(in) :: rhonow(cplex*nfft,nspden,ngrad*ngrad)
424  real(dp),intent(inout) :: vxc(cplex*nfft,nspden) !vz_i
425  real(dp),intent(inout),optional :: vxctau(cplex*nfft,nspden,4)
426 
427 !Local variables-------------------------------
428 !scalars
429  integer :: i1,i2,i3,id1,id2,id3,idir,ifft,ig1,ig2,ig3,ispden,n1,n2,n3,qeq0
430  real(dp),parameter :: lowden=1.d-14,precis=1.d-15
431  real(dp) :: gc23_idir,gcart_idir,ph123i,ph123r,ph1i,ph1r,ph23i,ph23r,ph2i,ph2r
432  real(dp) :: ph3i,ph3r,work_im,work_re
433  character(len=500) :: message
434 !arrays
435  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
436  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
437  real(dp) :: tsec(2)
438  real(dp),allocatable :: gcart1(:),gcart2(:),gcart3(:),ph1(:),ph2(:),ph3(:)
439  real(dp),allocatable :: wkcmpx(:,:),wkcmpxtau(:,:)
440  real(dp),allocatable :: work(:),workgr(:,:),worklp(:,:),worktau(:,:)
441 
442 ! *************************************************************************
443 
444  if (ishift/=0 .and. ishift/=1) then
445    write(message, '(a,i0)' )' ishift must be 0 or 1 ; input was',ishift
446    MSG_BUG(message)
447  end if
448 
449  if (ngrad/=1 .and. ngrad/=2 ) then
450    write(message, '(a,i0)' )' ngrad must be 1 or 2 ; input was',ngrad
451    MSG_BUG(message)
452  end if
453 
454 !Keep local copy of fft dimensions
455  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
456 
457 !Initialize computation of G in cartesian coordinates
458  id1=n1/2+2  ; id2=n2/2+2  ; id3=n3/2+2
459 
460  !Get the distrib associated with this fft_grid
461  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
462 
463 !Check whether q=0
464  qeq0=0
465  if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) qeq0=1
466 
467  if(ishift==0)then ! Add the value of depsxc to vxc
468    do ispden=1,min(nspden,2)
469 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,depsxc,nfft,vxc,ispden)
470      do ifft=1,cplex*nfft
471        vxc(ifft,ispden)=vxc(ifft,ispden)+depsxc(ifft,ispden)
472      end do
473    end do
474 
475  end if
476 
477 !If the grid is shifted, or if gradient corrections are present, there must be FFTs.
478  if(ishift==1 .or. ngrad==2)then
479 
480    ABI_ALLOCATE( wkcmpx,(2,nfft))
481    ABI_ALLOCATE(work,(cplex*nfft))
482 
483    if(ishift==1)then
484      ABI_ALLOCATE( ph1,(2*n1))
485      ABI_ALLOCATE(ph2,(2*n2))
486      ABI_ALLOCATE(ph3,(2*n3))
487 !    Precompute phases (The phases correspond to a shift of density on real space
488 !    grid from center at 0 0 0 to (1/2)*(1/n1,1/n2,1/n3).)
489      call phase(n1,ph1)
490      call phase(n2,ph2)
491      call phase(n3,ph3)
492    end if
493 
494    do ispden=1,min(nspden,2)
495 
496 !    Initialize wkcmpx either to 0 or to the shifted vxc value
497      if(ishift==0)then
498 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(nfft,wkcmpx)
499        do ifft=1,nfft
500          wkcmpx(:,ifft)=zero
501        end do
502 
503      else
504 !      Obtain depsxc(G)*phase in wkcmpx from input depsxc(r+delta)
505 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,depsxc,ispden,nfft,work)
506        do ifft=1,cplex*nfft
507          work(ifft)=depsxc(ifft,ispden)
508        end do
509        call timab(82,1,tsec)
510        call fourdp(cplex,wkcmpx,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
511        call timab(82,2,tsec)
512      end if
513 
514 !    If gradient correction is present, take care of the three components now
515 !    Note : this operation is done on the eventually shifted grid
516      if(ngrad==2)then
517        ABI_ALLOCATE(gcart1,(n1))
518        ABI_ALLOCATE(gcart2,(n2))
519        ABI_ALLOCATE(gcart3,(n3))
520        ABI_ALLOCATE(workgr,(2,nfft))
521        if(mgga==1)  then
522          ABI_ALLOCATE(worklp,(2,nfft))
523        end if
524        if(present(vxctau))  then
525          ABI_ALLOCATE(worktau,(2,nfft))
526          ABI_ALLOCATE(wkcmpxtau,(2,nfft))
527        end if
528 
529        do idir=1,3
530 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,idir,ispden,nfft,rhonow,work)
531          do ifft=1,cplex*nfft
532            work(ifft)=rhonow(ifft,ispden,1+idir)
533          end do
534 
535          call timab(82,1,tsec)
536          call fourdp(cplex,workgr,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
537          call timab(82,2,tsec)
538 
539 !        IF Meta-GGA then take care of the laplacian term involved.
540 !        Note : this operation is done on the eventually shifted grid
541          if(mgga==1)then
542 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,idir,ispden,nspden,nfft,depsxc,work)
543            do ifft=1,cplex*nfft
544              if(nspden==1)then
545                work(ifft)=depsxc(ifft,2+ispden)
546              else if(nspden==2)then
547                work(ifft)=depsxc(ifft,5+ispden)
548              end if
549            end do
550 
551            call timab(82,1,tsec)
552            call fourdp(cplex,worklp,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
553            call timab(82,2,tsec)
554          end if
555 
556          if(present(vxctau))then
557 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,vxctau,work)
558            do ifft=1,cplex*nfft
559              work(ifft)=vxctau(ifft,ispden,1)
560            end do
561 
562            call timab(82,1,tsec)
563            call fourdp(cplex,worktau,work,-1,mpi_enreg,nfft,ngfft,paral_kgb,0)
564            call timab(82,2,tsec)
565          end if ! present vxctau
566 
567          do i1=1,n1
568            ig1=i1-(i1/id1)*n1-1
569            gcart1(i1)=gprimd(idir,1)*two_pi*(dble(ig1)+qphon(1))
570          end do
571 !        Note that the G <-> -G symmetry must be maintained
572          if(mod(n1,2)==0 .and. qeq0==1)gcart1(n1/2+1)=zero
573          do i2=1,n2
574            ig2=i2-(i2/id2)*n2-1
575            gcart2(i2)=gprimd(idir,2)*two_pi*(dble(ig2)+qphon(2))
576          end do
577          if(mod(n2,2)==0 .and. qeq0==1)gcart2(n2/2+1)=zero
578          do i3=1,n3
579            ig3=i3-(i3/id3)*n3-1
580            gcart3(i3)=gprimd(idir,3)*two_pi*(dble(ig3)+qphon(3))
581          end do
582          if(mod(n3,2)==0 .and. qeq0==1)gcart3(n3/2+1)=zero
583 
584 !        !$OMP PARALLEL DO PRIVATE(ifft,i1,i2,i3,gc23_idir,gcart_idir) &
585 !        !$OMP&SHARED(gcart1,gcart2,gcart3,n1,n2,n3,wkcmpx,workgr)
586          ifft = 0
587          do i3=1,n3
588            do i2=1,n2
589              gc23_idir=gcart2(i2)+gcart3(i3)
590              if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
591                do i1=1,n1
592                  ifft=ifft+1
593                  gcart_idir=gc23_idir+gcart1(i1)
594 !                Multiply by - i 2pi G(idir) and accumulate in wkcmpx
595                  wkcmpx(1,ifft)=wkcmpx(1,ifft)+gcart_idir*workgr(2,ifft)
596                  wkcmpx(2,ifft)=wkcmpx(2,ifft)-gcart_idir*workgr(1,ifft)
597                  if(mgga==1)then
598 !                  Multiply by - i 2pi G(idir) and accumulate in wkcmpx
599                    wkcmpx(1,ifft)=wkcmpx(1,ifft)-gcart_idir**2*worklp(1,ifft)
600                    wkcmpx(2,ifft)=wkcmpx(2,ifft)-gcart_idir**2*worklp(2,ifft)
601                  end if
602                  if(present(vxctau))then
603                    wkcmpxtau(1,ifft)= gcart_idir*worktau(2,ifft)
604                    wkcmpxtau(2,ifft)=-gcart_idir*worktau(1,ifft)
605                  end if
606                end do
607              end if
608            end do
609          end do
610 
611          if(present(vxctau))then
612            call timab(82,1,tsec)
613            call fourdp(cplex,wkcmpxtau,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
614            call timab(82,2,tsec)
615 
616 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,vxctau,work)
617            do ifft=1,cplex*nfft
618              vxctau(ifft,ispden,1+idir)=work(ifft)
619            end do
620          end if ! present vxctau
621 
622        end do ! enddo idir
623        ABI_DEALLOCATE(gcart1)
624        ABI_DEALLOCATE(gcart2)
625        ABI_DEALLOCATE(gcart3)
626        ABI_DEALLOCATE(workgr)
627        if(mgga==1)then
628          ABI_DEALLOCATE(worklp)
629        end if
630        if(present(vxctau)) then
631          ABI_DEALLOCATE(worktau)
632          ABI_DEALLOCATE(wkcmpxtau)
633        end if
634      end if
635 
636 !    wkcmpx(:,:) contains now the full exchange-correlation potential, but
637 !    eventually for the shifted grid
638 
639      if(ishift==1)then
640 !      Take away the phase to get depsxc(G)
641        ifft=0
642        do i3=1,n3
643          ph3r=ph3(2*i3-1)
644          ph3i=ph3(2*i3  )
645          do i2=1,n2
646            ph2r=ph2(2*i2-1)
647            ph2i=ph2(2*i2  )
648            ph23r=ph2r*ph3r-ph2i*ph3i
649            ph23i=ph2i*ph3r+ph2r*ph3i
650            if (fftn2_distrib(i2)==mpi_enreg%me_fft) then
651              do i1=1,n1
652                ifft=ifft+1
653                ph1r=ph1(2*i1-1)
654                ph1i=ph1(2*i1  )
655                ph123r=ph1r*ph23r-ph1i*ph23i
656                ph123i=ph1i*ph23r+ph1r*ph23i
657 !              Multiply by phase.  Must use intermediate variables !
658                work_re= ph123r*wkcmpx(1,ifft)+ph123i*wkcmpx(2,ifft)
659                work_im=-ph123i*wkcmpx(1,ifft)+ph123r*wkcmpx(2,ifft)
660                wkcmpx(1,ifft)=work_re
661                wkcmpx(2,ifft)=work_im
662              end do
663            end if
664          end do
665        end do
666      end if
667 
668      call timab(82,1,tsec)
669      call fourdp(cplex,wkcmpx,work,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
670      call timab(82,2,tsec)
671 
672 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,vxc,work)
673      do ifft=1,cplex*nfft
674        vxc(ifft,ispden)=vxc(ifft,ispden)+work(ifft)
675      end do
676 
677    end do ! End loop on spins
678 
679    if(ishift==1)  then
680      ABI_DEALLOCATE(ph1)
681      ABI_DEALLOCATE(ph2)
682      ABI_DEALLOCATE(ph3)
683    end if
684    ABI_DEALLOCATE(wkcmpx)
685    ABI_DEALLOCATE(work)
686 
687  end if ! End condition on ishift
688 
689 end subroutine xcpot