TABLE OF CONTENTS


ABINIT/atm2fft [ Functions ]

[ Top ] [ Functions ]

NAME

 atm2fft

FUNCTION

 This routine sums atomic functions (density or potential) defined
 (in rec. space) on a radial grid to get global quantities on the
 fine FFT grid. It can also compute contribution to energy derivatives
 of these atomic functions.

 Possible options:
   optn=1: compute a sum of local atomic densities or contrib. to energy derivatives
   optv=1: compute a sum of local atomic potentials or contrib. to energy derivatives

   optatm =1: computes sum of atomic potentials/densities
   optgr  =1: computes contribution of atomic pot./dens. to forces
   optstr =1: computes contribution of atomic pot./dens. to stress tensor
   optdyfr=1: computes contribution of atomic pot./dens. to frozen part of dyn. matrix

COPYRIGHT

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

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx
  distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism
  gauss(2,ntypat)= params for gaussian atm density (optn2=3) for each atom type
  gmet(3,3)=reciprocal space metric
  gprimd(3,3)=reciprocal space dimensional primitive translations
  gsqcut=cutoff on |G|^2: see setup1 for definition (doubled sphere)
  kxc(2,nfft)=exchange and correlation kernel
  me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present)
  mgfft=maximum size of 1D FFTs
  comm_fft=--optional-- MPI communicator over FFT components
  mqgrid=number of grid pts in q array for f(q) spline.
  natom=number of atoms in unit cell.
  nattyp(ntypat)=number of atoms of each type in cell.
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT
  ntypat=number of types of atoms.
  optatm,optdyfr,optgr,optn,optn2,optstr,optv= (see NOTES below)
  paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present)
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phase information
  qgrid(mqgrid)=q grid for spline from 0 to qmax.
  qprtrb(3)= integer wavevector of possible perturbing potential
             in basis of reciprocal lattice translations
  rhog(2,nfft)=electron density rho(G) in reciprocal space
               (used only if optv=1 and (optgr=1 or optstr=1 or optdyfr=1))
  ucvol=unit cell volume
  usepaw= 0 for non paw calculation; =1 for paw calculation
  vspl(mqgrid,2,ntypat)=q^2 v(q) spline of an atomic potential
                        (used only if optv=1)
  vprtrb(2)=complex amplitude of possible perturbing potential; if nonzero,
            perturbing potential is added of the form V(G)=(vprtrb(1)+I*vprtrb(2))/2
            at the values G=qprtrb and (vprtrb(1)-I*vprtrb(2))/2 at G=-qprtrb
  vg(2,nfft)= potential V(G) in reciprocal space
              (used only if optn=1 and (optgr=1 or optstr=1 or optdyfr=1))
  vg1(2,nfft)= 1st-order potential V(G) in reciprocal space
               (used only if opteltfr==1)
  vg1_core(2,nfft)= 1-order potential V(G) in reciprocal space with only core contribution
                    (used only if opteltfr==1)

OUTPUT

  ======= if optv==1 =======
  ============================
   --- if optatm==1
    atmvloc(nfft)=sum of local atomic potentials in real space
   --- if optgr==1
    grv(3,natom)=contribution of atomic potentials to forces
   --- if optstr==1
    strv(6)=contribution of atomic potentials to stress tensor
            cart. coordinates, symmetric tensor, 6 comp. in order 11,22,33,32,31,21
   --- if optdyfr==1
    dyfrv(3,3,natom)=contribution of atomic potentials to frozen part of dyn. matrix

  ======= if optn==1 =======
  ============================
   --- if optatm==1
    atmrho(nfft)=sum of atomic densities in real space
   --- if optgr==1
    grn(3,natom)=contribution of atomic densities to forces
   --- if optstr==1
    strn(6)=contribution of atomic densities to stress tensor
            cart. coordinates, symmetric tensor, 6 comp. in order 11,22,33,32,31,21
   --- if optdyfr==1
    dyfrn(3,3,natom)=contribution of atomic densities to frozen part of dyn. matrix
   --- if opteltfr==1
    eltfrn(6+3*natom,6)=contribution of atomic density to frozen part of stress tensor

NOTES

 Details on possible options:
 ============================
 optv: controls the computation of a local potential as sum of atomic potentials
          Vloc(r)=Sum_R[V^AT(r-R)]
       or its contributions to energy derivatives, i.e. derivatives of Int[Vloc(r).rho(r).dr]
          rho(r) is stored in reciprocal space in array rhog()
          V^AT is stored in reciprocal space in array vspl (in practice vspl(q)=q^2.V^AT(q))

 optn: controls the computation of a density as sum of atomic densities
          n(r)=Sum_R[n^AT(r-R)]
       or its contributions to energy derivatives, i.e. derivatives of Int[n(r).V(r).dr]
          V(r) is stored in reciprocal space in array vg()
          n^AT is stored in reciprocal space:
          if optn2=1: n^AT is the atomic PAW PS core density stored in array pawtab%tcorespl()
                   2: n^AT is the atomic PAW PS valence density stored in array pawtab%tvalespl()
                   3: n^AT is a gaussian density: n(g)=gauss(1,ityp)*exp[-(gauss(2,ityp)*G)^2]
 Note: optv and optn can be activated together

 Options controlling which contrib. to Etot derivatives are computed:
   optatm  =1: computes Vloc(r) or n(r) as sum of atomic potentials/densities
   optgr   =1: computes contribution of atomic Vloc(r) or n(r) to forces
   optstr  =1: computes contribution of atomic Vloc(r) or n(r) to stress tensor
   optdyfr =1: computes contribution of atomic Vloc(r) or n(r) to fr part of dyn. matrix
   opteltfr=1: computes contribution of atomic Vloc(r) or n(r) to elastic tensor
 Note: optatm, optgr, optstr, optelfr and optdyfr can be activated together

 Typical uses:
 =============
 Computation of:
  - local potential: optv=1, optatm=1
  - contrib. of local potential to Etot derivatives: optv=1, rhog=total valence density
                                                     optgr=1 or optstr=1 or optdyfr=1
  - PS core density: optn=1, optn2=1, optatm=1
  - contrib. of NLCC to Etot derivatives: optn=1, optn2=1, vg=XC potential
                                          optgr=1 or optstr=1 or optdyfr=1
  - sum of atomic valence densities: optn=1, optn2=2 or 3, optatm=1
  - correction of forces due to potential residual: optn=1, optn2=2 or 3, optgr=1
                                                    vg=potential residual
    etc...

PARENTS

      dfpt_dyfro,dfpt_eltfrxc,extraprho,forces,fresidrsp,prcref,prcref_PMA
      respfn,setvtr,stress

CHILDREN

      destroy_distribfft,fourdp,init_distribfft_seq,initmpi_seq
      set_mpi_enreg_fft,timab,unset_mpi_enreg_fft,wrtout,xmpi_sum,zerosym

SOURCE

145 #if defined HAVE_CONFIG_H
146 #include "config.h"
147 #endif
148 
149 #include "abi_common.h"
150 
151 subroutine atm2fft(atindx1,atmrho,atmvloc,dyfrn,dyfrv,eltfrn,gauss,gmet,gprimd,&
152 &                  grn,grv,gsqcut,mgfft,mqgrid,natom,nattyp,nfft,ngfft,ntypat,&
153 &                  optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,&
154 &                  psps,pawtab,ph1d,qgrid,qprtrb,rhog,strn,strv,ucvol,usepaw,vg,vg1,vg1_core,vprtrb,vspl,&
155 &                  is2_in,comm_fft,me_g0,paral_kgb,distribfft) ! optional arguments
156 
157  use m_profiling_abi
158  use m_errors
159  use defs_basis
160  use defs_abitypes
161  use m_errors
162  use m_xmpi
163 
164  use defs_datatypes,only : pseudopotential_type
165  use m_pawtab,      only : pawtab_type
166  use m_distribfft,  only : distribfft_type
167  use m_fft,         only : zerosym
168  use m_mpinfo,      only : set_mpi_enreg_fft, unset_mpi_enreg_fft
169 
170 !This section has been created automatically by the script Abilint (TD).
171 !Do not modify the following lines by hand.
172 #undef ABI_FUNC
173 #define ABI_FUNC 'atm2fft'
174  use interfaces_14_hidewrite
175  use interfaces_18_timing
176  use interfaces_51_manage_mpi
177  use interfaces_53_ffts
178 !End of the abilint section
179 
180  implicit none
181 
182 !Arguments ------------------------------------
183 !scalars
184  integer,intent(in) :: mgfft,mqgrid,natom,nfft,ntypat,optatm,optdyfr,opteltfr
185  integer,intent(in) :: optgr,optn,optn2,optstr,optv,usepaw
186  integer,optional,intent(in) :: is2_in,me_g0,comm_fft,paral_kgb
187  real(dp),intent(in) :: gsqcut,ucvol
188  type(pseudopotential_type),target,intent(in) :: psps
189  type(distribfft_type),optional,intent(in),target :: distribfft
190 !arrays
191  integer,intent(in) :: atindx1(natom),nattyp(ntypat),ngfft(18),qprtrb(3)
192  real(dp),intent(in) :: gauss(2,ntypat*(optn2/3)),gmet(3,3),gprimd(3,3)
193  real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid)
194  real(dp),intent(in) :: rhog(2,nfft*optv*max(optgr,optstr,optdyfr,opteltfr))
195  real(dp),intent(in) :: vg(2,nfft*optn*max(optgr,optstr,optdyfr,opteltfr))
196  real(dp),intent(in) :: vg1(2,nfft*optn*opteltfr),vg1_core(2,nfft*optn*opteltfr)
197  real(dp),intent(in) :: vprtrb(2),vspl(mqgrid,2,ntypat*optv)
198  real(dp),intent(out) :: atmrho(nfft*optn)
199  real(dp),intent(inout) :: atmvloc(nfft*optv)
200  real(dp),intent(out) :: dyfrn(3,3,natom*optn*optdyfr),dyfrv(3,3,natom*optv*optdyfr)
201  real(dp),intent(out) :: eltfrn(6+3*natom,6)
202  real(dp),intent(inout) :: grn(3,natom*optn*optgr)
203  real(dp),intent(out) :: grv(3,natom*optv*optgr),strn(6*optn*optstr)
204  real(dp),intent(out) :: strv(6*optv*optstr)
205  type(pawtab_type),target,intent(in) :: pawtab(ntypat*usepaw)
206 
207 !Local variables ------------------------------
208 !scalars
209  integer,parameter :: im=2,re=1
210  integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ierr,ig1,ig1_,ig2,ig2_,ig3,ig3_,ii,is1,is2
211  integer :: itypat,jj,js,ka,kb,kd,kg,me_fft,my_comm_fft,ndir,n1,n2,n3,nproc_fft,paral_kgb_fft
212  integer :: shift1,shift2,shift3
213  logical :: have_g0
214  real(dp),parameter :: tolfix=1.0000001_dp
215  real(dp) :: aa,alf2pi2,bb,cc,cutoff,dbl_ig1,dbl_ig2,dbl_ig3,dd,dg1,dg2,d2g,diff
216  real(dp) :: dn_at,d2n_at,d2n_at2,dq,dq2div6,dqdiv6,dqm1,dv_at,ee,ff,gauss1,gauss2,gg,gmag,gsquar,n_at
217  real(dp) :: ph12i,ph12r,ph1i,ph1r,ph2i,ph2r,ph3i,ph3r,sfi,sfr,term,term1,term2,tmpni,tmpnr
218  real(dp) :: tmpvi,tmpvr,v_at,xnorm
219  character(len=500) :: message
220  type(distribfft_type),pointer :: my_distribfft
221  type(mpi_type) :: mpi_enreg_fft
222 !arrays
223  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
224  real(dp), ABI_CONTIGUOUS pointer :: tvalespl(:,:),tcorespl(:,:)
225  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
226  integer  :: delta(6)=(/1,1,1,0,0,0/)
227  real(dp) :: dgm(3,3,6),d2gm(3,3,6,6),gcart(3),tsec(2)
228  real(dp),allocatable :: dyfrn_indx(:,:,:),dyfrv_indx(:,:,:),grn_indx(:,:)
229  real(dp),allocatable :: grv_indx(:,:),phim_igia(:),phre_igia(:),workn(:,:)
230  real(dp),allocatable :: workv(:,:)
231 
232 ! *************************************************************************
233 
234  DBG_ENTER("COLL")
235 
236 !Check optional arguments
237  if (present(comm_fft)) then
238    if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then
239      MSG_BUG(' Need paral_kgb and me_g0 with comm_fft !')
240    end if
241  end if
242 
243  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
244  
245  me_fft=ngfft(11)
246  nproc_fft=ngfft(10)
247 
248 
249 !Get the distrib associated with this fft_grid
250  if (present(distribfft)) then
251    my_distribfft => distribfft
252  else
253    ABI_DATATYPE_ALLOCATE(my_distribfft,)
254    call init_distribfft_seq(my_distribfft,'f',n2,n3,'fourdp')
255  end if
256  if (n2==my_distribfft%n2_coarse) then
257    fftn2_distrib => my_distribfft%tab_fftdp2_distrib
258    ffti2_local => my_distribfft%tab_fftdp2_local
259  else if (n2 == my_distribfft%n2_fine) then
260    fftn2_distrib => my_distribfft%tab_fftdp2dg_distrib
261    ffti2_local => my_distribfft%tab_fftdp2dg_local
262  else
263    MSG_BUG("Unable to find an allocated distrib for this fft grid")
264  end if
265 
266  if (present(is2_in)) then
267    if(is2_in<1.or.is2_in>6) then
268      MSG_BUG("is2_in must be between 1 and 6")
269    else
270      ndir = 1
271    end if
272    ndir = -1
273  end if
274 
275 !Zero out arrays to permit accumulation over atom types
276  if (optv==1.and.optatm==1) then
277    ABI_ALLOCATE(workv,(2,nfft))
278    workv(:,:)=zero
279  end if
280  if (optn==1.and.optatm==1) then
281    ABI_ALLOCATE(workn,(2,nfft))
282    workn(:,:)=zero
283  end if
284  if (optv==1.and.optgr==1) then
285    ABI_ALLOCATE(grv_indx,(3,natom))
286    grv_indx(:,:)=zero
287  end if
288  if (optn==1.and.optgr==1) then
289    ABI_ALLOCATE(grn_indx,(3,natom))
290    grn_indx(:,:)=zero
291  end if
292  if (optv==1.and.optdyfr==1) then
293    ABI_ALLOCATE(dyfrv_indx,(3,3,natom))
294    dyfrv_indx(:,:,:)=zero
295  end if
296  if (optn==1.and.optdyfr==1) then
297    ABI_ALLOCATE(dyfrn_indx,(3,3,natom))
298    dyfrn_indx(:,:,:)=zero
299  end if
300  if (optv==1.and.optstr==1) strv(:)=zero
301  if (optn==1.and.optstr==1) strn(:)=zero
302  if (opteltfr==1) eltfrn(:,:) = zero
303 
304 !Compute 1st and 2nd derivatives of metric tensor wrt all strain components
305 !and store for use in inner loop below for elastic tensor.
306  if (opteltfr==1) then
307    dgm(:,:,:)=zero
308    d2gm(:,:,:,:)=zero
309 !  Loop over 2nd strain index
310    do is2=1,6
311      kg=idx(2*is2-1);kd=idx(2*is2)
312      do jj = 1,3
313        dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj))
314      end do
315 
316 !    Loop over 1st strain index, upper triangle only
317      do is1=1,is2
318        ka=idx(2*is1-1);kb=idx(2*is1)
319        do jj = 1,3
320          if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
321 &         +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj)
322          if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
323 &         +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj)
324          if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
325 &         +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj)
326          if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
327 &         +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj)
328        end do
329      end do !is1
330    end do !is2
331  end if
332 
333 
334  dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
335  dqm1=1.0_dp/dq
336  dqdiv6=dq/6.0_dp
337  dq2div6=dq**2/6.0_dp
338  cutoff=gsqcut*tolfix
339  id1=n1/2+2
340  id2=n2/2+2
341  id3=n3/2+2
342 
343  ABI_ALLOCATE(phre_igia,(natom))
344  ABI_ALLOCATE(phim_igia,(natom))
345 
346  ia1=1
347  do itypat=1,ntypat
348 !  ia1,ia2 sets range of loop over atoms:
349    ia2=ia1+nattyp(itypat)-1
350    ii=0
351 
352    if (optn2==3)then
353      gauss1=gauss(1,itypat)
354      gauss2=gauss(2,itypat)
355      alf2pi2=(two_pi*gauss2)**2
356    end if
357 
358    if (usepaw == 1) then
359      tcorespl => pawtab(itypat)%tcorespl
360      tvalespl => pawtab(itypat)%tvalespl
361    else
362      tcorespl => psps%nctab(itypat)%tcorespl
363      tvalespl => psps%nctab(itypat)%tvalespl
364    end if
365    
366    do i3=1,n3
367      ig3=i3-(i3/id3)*n3-1
368      ig3_=ig3;if (ig3_==(n3/2+1)) ig3_=0
369      do i2=1,n2
370        ig2=i2-(i2/id2)*n2-1
371        ig2_=ig2;if (ig2_==(n2/2+1)) ig2_=0
372        if(fftn2_distrib(i2)==me_fft) then
373          do i1=1,n1
374            ig1=i1-(i1/id1)*n1-1
375            ig1_=ig1;if (ig1_==(n1/2+1)) ig1_=0
376            ii=ii+1
377            gsquar=gsq_atm(ig1,ig2,ig3)
378 
379 !          Skip G**2 outside cutoff:
380            if (gsquar<=cutoff) then
381 
382              gmag=sqrt(gsquar)
383              have_g0=(ig1==0.and.ig2==0.and.ig3==0)
384 
385              jj=1+int(gmag*dqm1)
386              diff=gmag-qgrid(jj)
387 
388 !            Compute structure factor for all atoms of given type:
389              do ia=ia1,ia2
390                shift1=1+n1+(ia-1)*(2*n1+1)
391                shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)
392                shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)
393                ph1r=ph1d(1,ig1+shift1);ph1i=ph1d(2,ig1+shift1)
394                ph2r=ph1d(1,ig2+shift2);ph2i=ph1d(2,ig2+shift2)
395                ph3r=ph1d(1,ig3+shift3);ph3i=ph1d(2,ig3+shift3)
396                ph12r=ph1r*ph2r-ph1i*ph2i
397                ph12i=ph1r*ph2i+ph1i*ph2r
398                phre_igia(ia)=ph12r*ph3r-ph12i*ph3i
399                phim_igia(ia)=ph12r*ph3i+ph12i*ph3r
400              end do
401 
402 !            Assemble structure factors for this type of atom= sum[exp(-i.piG.R)]
403              if (optatm==1.or.optstr==1.or.opteltfr==1) then
404                sfr=zero;sfi=zero
405                do ia=ia1,ia2
406                  sfr=sfr+phre_igia(ia)
407                  sfi=sfi-phim_igia(ia)
408                end do
409              end if
410 
411 !            Compute V^AT(G) and/or n^AT(G) for given type of atom
412 !            Evaluate spline fit: p. 86 Numerical Recipes, Press et al;
413 !            NOTE: error in book for sign of "aa" term in derivative;
414 !            !           also see splfit routine.
415              if (optv==1.or.optn2/=3) then
416                bb = diff*dqm1
417                aa = 1.0_dp-bb
418                cc = aa*(aa**2-1.0_dp)*dq2div6
419                dd = bb*(bb**2-1.0_dp)*dq2div6
420              end if
421              if (optv==1) then
422                if (have_g0) then
423                  v_at=zero
424                else
425                  v_at=(aa*vspl(jj,1,itypat)+bb*vspl(jj+1,1,itypat)+&
426 &                 cc*vspl(jj,2,itypat)+dd*vspl(jj+1,2,itypat)) &
427 &                 /gsquar
428                end if
429              end if
430              if (optn==1) then
431                if (optn2==1) then
432                  n_at=(aa*tcorespl(jj,1)+bb*tcorespl(jj+1,1)+cc*tcorespl(jj,2)+dd*tcorespl(jj+1,2))
433                else if (optn2==2) then
434                  n_at=(aa*tvalespl(jj,1)+bb*tvalespl(jj+1,1)+cc*tvalespl(jj,2)+dd*tvalespl(jj+1,2))
435                else if (optn2==3) then
436                  n_at=gauss1*exp(-gsquar*alf2pi2)
437                else
438                  n_at=zero
439                end if
440              end if
441 
442 !            Compute sum of local atomic potentials or densities
443 !            ---------------------------------------------------
444              if(optatm==1) then
445 !              Accumulate V^AT(G)*SF(G) or n^AT(G)*SF(G)
446                if (optv==1) then
447                  workv(re,ii)=workv(re,ii)+sfr*v_at
448                  workv(im,ii)=workv(im,ii)+sfi*v_at
449                end if
450                if (optn==1) then
451                  workn(re,ii)=workn(re,ii)+sfr*n_at
452                  workn(im,ii)=workn(im,ii)+sfi*n_at
453                end if
454 
455 !            Compute contrib. to forces and/or frozen part of dyn. matrix
456 !            -------------------------------------------------------------
457              else if (optgr==1.or.optdyfr==1) then
458                dbl_ig1=dble(ig1_);dbl_ig2=dble(ig2_);dbl_ig3=dble(ig3_)
459 !              Compute (2Pi)*V^AT(G)*rho(G) or (2Pi)*n^AT(G)*V(G)
460                if (optv==1) then
461                  tmpvr=(two_pi*v_at)*rhog(re,ii)
462                  tmpvi=(two_pi*v_at)*rhog(im,ii)
463                end if
464                if (optn==1) then
465                  tmpnr=(two_pi*n_at)*vg(re,ii)
466                  tmpni=(two_pi*n_at)*vg(im,ii)
467                end if
468 !              === contrib. to forces
469                if (optgr==1) then
470 !                Accumulate -(2Pi.G)*V^AT(G)*rho(G)*SF(G)
471 !                or -(2Pi)*n^AT(G)*V(G)*SF(G) into forces
472                  if (optv==1) then
473                    do ia=ia1,ia2
474                      term=tmpvi*phre_igia(ia)+tmpvr*phim_igia(ia)
475                      grv_indx(1,ia)=grv_indx(1,ia)-dbl_ig1*term
476                      grv_indx(2,ia)=grv_indx(2,ia)-dbl_ig2*term
477                      grv_indx(3,ia)=grv_indx(3,ia)-dbl_ig3*term
478                    end do
479                  end if
480                  if (optn==1) then
481                    do ia=ia1,ia2
482                      term=tmpni*phre_igia(ia)+tmpnr*phim_igia(ia)
483                      grn_indx(1,ia)=grn_indx(1,ia)-dbl_ig1*term
484                      grn_indx(2,ia)=grn_indx(2,ia)-dbl_ig2*term
485                      grn_indx(3,ia)=grn_indx(3,ia)-dbl_ig3*term
486                    end do
487                  end if
488                end if
489 !              === contrib. to frozen part of dyn. matrix
490                if (optdyfr==1) then
491 !                Accumulate -(2Pi^2.Gi.Gj)*V^AT(G)*rho(G)*SF(G)
492 !                or -(2Pi^2.Gi.Gj)*n^AT(G)*V(G)*SF(G) into dyn. matrix
493                  if (optv==1) then
494                    do ia=ia1,ia2
495                      term=two_pi*(tmpvr*phre_igia(ia)-tmpvi*phim_igia(ia))
496                      dyfrv_indx(1,1,ia)=dyfrv_indx(1,1,ia)-dbl_ig1*dbl_ig1*term
497                      dyfrv_indx(1,2,ia)=dyfrv_indx(1,2,ia)-dbl_ig1*dbl_ig2*term
498                      dyfrv_indx(1,3,ia)=dyfrv_indx(1,3,ia)-dbl_ig1*dbl_ig3*term
499                      dyfrv_indx(2,2,ia)=dyfrv_indx(2,2,ia)-dbl_ig2*dbl_ig2*term
500                      dyfrv_indx(2,3,ia)=dyfrv_indx(2,3,ia)-dbl_ig2*dbl_ig3*term
501                      dyfrv_indx(3,3,ia)=dyfrv_indx(3,3,ia)-dbl_ig3*dbl_ig3*term
502                    end do
503                  end if
504                  if (optn==1) then
505                    do ia=ia1,ia2
506                      term=two_pi*(tmpnr*phre_igia(ia)-tmpni*phim_igia(ia))
507                      dyfrn_indx(1,1,ia)=dyfrn_indx(1,1,ia)-dbl_ig1*dbl_ig1*term
508                      dyfrn_indx(1,2,ia)=dyfrn_indx(1,2,ia)-dbl_ig1*dbl_ig2*term
509                      dyfrn_indx(1,3,ia)=dyfrn_indx(1,3,ia)-dbl_ig1*dbl_ig3*term
510                      dyfrn_indx(2,2,ia)=dyfrn_indx(2,2,ia)-dbl_ig2*dbl_ig2*term
511                      dyfrn_indx(2,3,ia)=dyfrn_indx(2,3,ia)-dbl_ig2*dbl_ig3*term
512                      dyfrn_indx(3,3,ia)=dyfrn_indx(3,3,ia)-dbl_ig3*dbl_ig3*term
513                    end do
514                  end if
515                end if
516              end if
517 
518 !            Compute (dV^AT(q)/dq)/q and/or (dn^AT(q)/dq)/q
519 !            For stress tensor and/or elastic tensor
520 !            ---------------------------------
521              if (optstr==1.or.opteltfr==1) then
522 !              Note: correction of Numerical Recipes sign error before (3._dp*aa**2-1._dp)
523 !              ee*dqm1 + ff*dqdiv6 is the best estimate of dV(q)/dq from splines
524                if (optv==1) then
525                  if (have_g0) then
526                    dv_at=zero
527                  else
528                    ee=vspl(jj+1,1,itypat)-vspl(jj,1,itypat)
529                    ff=(3._dp*bb**2-1._dp)*vspl(jj+1,2,itypat)&
530 &                   -(3._dp*aa**2-1._dp)*vspl(jj  ,2,itypat)
531                    dv_at=((ee*dqm1+ff*dqdiv6)/gmag-2.0_dp*v_at)/gsquar
532                  end if
533                end if
534                if (optn==1) then
535                  if (have_g0) then
536                    if (optn2==1) then
537                      if (usepaw ==1) then
538                        dn_at=pawtab(itypat)%dncdq0
539                      else
540                        dn_at=psps%nctab(itypat)%dncdq0
541                      end if
542                    else if (optn2==2) then
543                      if (usepaw == 1) then
544                        dn_at=pawtab(itypat)%dnvdq0
545                      else
546                        dn_at=psps%nctab(itypat)%dnvdq0
547                      end if
548                    else if (optn2==3) then
549                      dn_at=-two*gauss1*alf2pi2
550                    end if
551                    if (opteltfr==1) then
552                      d2n_at = 0
553                    end if
554                  else
555                    if (optn2==1) then
556                      ee=tcorespl(jj+1,1)-tcorespl(jj,1)
557                      ff=(3._dp*bb**2-1._dp)*tcorespl(jj+1,2) &
558 &                     -(3._dp*aa**2-1._dp)*tcorespl(jj,2)
559 !                    Also get nc''(q)
560                      if (opteltfr==1) then
561                        gg=aa*tcorespl(jj,2)+bb*tcorespl(jj+1,2)
562                      end if
563                    else if (optn2==2) then
564                      ee=tvalespl(jj+1,1)-tvalespl(jj,1)
565                      ff=(3._dp*bb**2-1._dp)*tvalespl(jj+1,2) &
566 &                     -(3._dp*aa**2-1._dp)*tvalespl(jj,2)
567 !                    Also get nc''(q)
568                      if (opteltfr==1) then
569                        gg=aa*tvalespl(jj,2)+bb*tvalespl(jj+1,2)
570                      end if
571                    else if (optn2==3) then
572                      dn_at=-two*gauss1*alf2pi2*exp(-gsquar*alf2pi2)
573                    else
574                    end if
575                    dn_at  = (ee*dqm1+ff*dqdiv6)/gmag
576                    if (opteltfr==1) then
577                      d2n_at = (gg-dn_at)/gsquar
578                      d2n_at2 = gg/gmag**3
579                    end if
580                  end if
581                end if
582 
583 !              Compute G in cartesian coordinates
584                gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+&
585 &               gprimd(1,3)*dble(ig3)
586                gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+&
587 &               gprimd(2,3)*dble(ig3)
588                gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+&
589 &               gprimd(3,3)*dble(ig3)
590              end if
591 
592 !            Compute contrib. to stress tensor
593 !            ---------------------------------
594              if (optstr==1)then
595 !              Accumulate -dV^AT/dG*rho(G)*SF(G)*Gi.Gj/G
596 !              or -dn^AT/dG*V(G)*SF(G)*Gi.Gj/G
597 !              into stress tensor
598                if (optv==1) then
599                  term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi)
600                  strv(1)=strv(1)-term*(dv_at*gcart(1)*gcart(1)+v_at)
601                  strv(2)=strv(2)-term*(dv_at*gcart(2)*gcart(2)+v_at)
602                  strv(3)=strv(3)-term*(dv_at*gcart(3)*gcart(3)+v_at)
603                  strv(4)=strv(4)-term*dv_at*gcart(3)*gcart(2)
604                  strv(5)=strv(5)-term*dv_at*gcart(3)*gcart(1)
605                  strv(6)=strv(6)-term*dv_at*gcart(2)*gcart(1)
606                end if
607                if (optn==1) then
608                  term=(vg(re,ii)*sfr+vg(im,ii)*sfi)*dn_at
609                  strn(1)=strn(1)-term*gcart(1)*gcart(1)
610                  strn(2)=strn(2)-term*gcart(2)*gcart(2)
611                  strn(3)=strn(3)-term*gcart(3)*gcart(3)
612                  strn(4)=strn(4)-term*gcart(3)*gcart(2)
613                  strn(5)=strn(5)-term*gcart(3)*gcart(1)
614                  strn(6)=strn(6)-term*gcart(2)*gcart(1)
615                end if
616              end if
617 
618 !            Compute contrib. to elastic tensor
619 !            ---------------------------------
620              if (opteltfr==1) then
621                dbl_ig1=dble(ig1_);dbl_ig2=dble(ig2_);dbl_ig3=dble(ig3_)
622 !             if (ig1==0 .and. ig2==0 .and. ig3==0) cycle
623 !              Compute G*dG/Deps_{\gamme\delta}
624                dg2=0.5_dp*dgsqds_atm(ig1,ig2,ig3,is2_in)
625 
626                term  = vg (re,ii)*sfr + vg (im,ii)*sfi
627                term1 = vg1(re,ii)*sfr + vg1(im,ii)*sfi
628                term2 = vg1_core(re,ii)*sfr + vg1_core(im,ii)*sfi
629                
630                do is1=1,6
631 
632 !                Compute G*dG/Deps_{\alpha\beta}
633                  dg1=0.5_dp*dgsqds_atm(ig1,ig2,ig3,is1)
634 !                Compute G^2*d2G/Deps_{alphabetagammadelta}
635                  d2g=(0.25_dp*d2gsqds_atm(ig1,ig2,ig3,is1,is2_in))
636 
637                  eltfrn(is1,is2_in) = eltfrn(is1,is2_in) + (term*(d2n_at*dg1*dg2 + d2g*dn_at))
638                  eltfrn(is1,is2_in) = eltfrn(is1,is2_in) + 0.5*(term1*dn_at*dg1)
639                  eltfrn(is2_in,is1) = eltfrn(is2_in,is1) + 0.5*(term1*dn_at*dg1)
640 
641                  if(is2_in<=3)then
642                    eltfrn(is1,is2_in) = eltfrn(is1,is2_in) - term*dn_at*dg1
643                  end if
644                  if(is1<=3)then
645                    eltfrn(is1,is2_in) = eltfrn(is1,is2_in) - term*dn_at*dg2
646                  end if
647                  if(is2_in<=3.and.is1<=3)then
648                    eltfrn(is1,is2_in) = eltfrn(is1,is2_in) - (term1-term)*n_at
649                  end if
650                end do
651 
652 !              internal strain
653                do ia=ia1,ia2
654                  js=7+3*(ia-1)
655 !                Compute -2pi*G*i*vxcis_core(G)*nat(G)*(exp(-iGr))
656                  term=(vg1_core(im,ii)*phre_igia(ia)+vg1_core(re,ii)*phim_igia(ia))*n_at*two_pi
657                  eltfrn(js  ,is2_in) = eltfrn(js  ,is2_in) - dbl_ig1*term
658                  eltfrn(js+1,is2_in) = eltfrn(js+1,is2_in) - dbl_ig2*term
659                  eltfrn(js+2,is2_in) = eltfrn(js+2,is2_in) - dbl_ig3*term
660 
661 !                Compute -2pi*G*i*vxc(G)*(dnat*dG/deps-delta*nat(G))*(exp(-iGr))
662                  term=(vg(im,ii)*phre_igia(ia)+vg(re,ii)*phim_igia(ia))*&
663 &                 (dn_at*dg2-delta(is2_in)*n_at)*two_pi
664                  eltfrn(js  ,is2_in) = eltfrn(js  ,is2_in) - dbl_ig1*term 
665                  eltfrn(js+1,is2_in) = eltfrn(js+1,is2_in) - dbl_ig2*term 
666                  eltfrn(js+2,is2_in) = eltfrn(js+2,is2_in) - dbl_ig3*term 
667                end do
668              end if
669 !            End skip G**2 outside cutoff:
670            end if
671 !          End loop on n1, n2, n3
672          end do
673        end if ! this plane is for me_fft
674      end do
675    end do
676 
677 !  Symmetrize the dynamical matrix with respect to indices
678    if (optdyfr==1) then
679      if (optv==1) then
680        do ia=ia1,ia2
681          dyfrv_indx(2,1,ia)=dyfrv_indx(1,2,ia)
682          dyfrv_indx(3,1,ia)=dyfrv_indx(1,3,ia)
683          dyfrv_indx(3,2,ia)=dyfrv_indx(2,3,ia)
684        end do
685      end if
686      if (optn==1) then
687        do ia=ia1,ia2
688          dyfrn_indx(2,1,ia)=dyfrn_indx(1,2,ia)
689          dyfrn_indx(3,1,ia)=dyfrn_indx(1,3,ia)
690          dyfrn_indx(3,2,ia)=dyfrn_indx(2,3,ia)
691        end do
692      end if
693    end if
694 
695    ia1=ia2+1
696 
697 !  End loop on type of atoms
698  end do
699 
700  ABI_DEALLOCATE(phre_igia)
701  ABI_DEALLOCATE(phim_igia)
702 
703 !Get local potential or density back to real space
704  if(optatm==1)then
705 !  Allow for the addition of a perturbing potential
706    if ((optv==1).and.(vprtrb(1)**2+vprtrb(2)**2) > 1.d-30) then
707 !    Find the linear indices which correspond with the input wavevector qprtrb
708 !    The double modulus handles both i>=n and i<0, mapping into [0,n-1];
709 !    then add 1 to get range [1,n] for each
710      i3=1+mod(n3+mod(qprtrb(3),n3),n3)
711      i2=1+mod(n2+mod(qprtrb(2),n2),n2)
712      i1=1+mod(n1+mod(qprtrb(1),n1),n1)
713 !    Compute the linear index in the 3 dimensional array
714      ii=i1+n1*((ffti2_local(i2)-1)+(n2/nproc_fft)*(i3-1))
715 !    Add in the perturbation at G=qprtrb
716      workv(re,ii)=workv(re,ii)+0.5_dp*vprtrb(1)
717      workv(im,ii)=workv(im,ii)+0.5_dp*vprtrb(2)
718 !    Same thing for G=-qprtrb
719      i3=1+mod(n3+mod(-qprtrb(3),n3),n3)
720      i2=1+mod(n2+mod(-qprtrb(2),n2),n2)
721      i1=1+mod(n1+mod(-qprtrb(1),n1),n1)
722 !    ii=i1+n1*((i2-1)+n2*(i3-1))
723      workv(re,ii)=workv(re,ii)+0.5_dp*vprtrb(1)
724      workv(im,ii)=workv(im,ii)-0.5_dp*vprtrb(2)
725      write(message, '(a,1p,2e12.4,a,0p,3i4,a)' )&
726 &     ' atm2fft: perturbation of vprtrb=', vprtrb,&
727 &     ' and q=',qprtrb,' has been added'
728      call wrtout(std_out,message,'COLL')
729    end if
730 
731    if (optv==1.or.optn==1) then
732 !    Create fake mpi_enreg to wrap fourdp
733      call initmpi_seq(mpi_enreg_fft)
734      ABI_DATATYPE_DEALLOCATE(mpi_enreg_fft%distribfft)
735      if (present(comm_fft)) then
736        call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb)
737        my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb
738      else
739        my_comm_fft=xmpi_comm_self;paral_kgb_fft=0;
740        mpi_enreg_fft%distribfft => my_distribfft
741      end if
742 !    Non-symetrized non-zero elements have to be nullified
743 !    Transform back to real space; divide by unit cell volume
744      xnorm=one/ucvol
745      if (optv==1) then
746        call zerosym(workv,2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft)
747        call fourdp(1,workv,atmvloc,1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0)
748        atmvloc(:)=atmvloc(:)*xnorm
749        ABI_DEALLOCATE(workv)
750      end if
751      if (optn==1) then
752        call zerosym(workn,2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft)
753        call fourdp(1,workn,atmrho,1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0)
754        atmrho(:)=atmrho(:)*xnorm
755        ABI_DEALLOCATE(workn)
756      end if
757 !    Destroy fake mpi_enreg
758      call unset_mpi_enreg_fft(mpi_enreg_fft)
759    end if
760 
761  end if
762 
763 !Additional treatment in case of parallelization
764  if (present(comm_fft)) then
765    if ((xmpi_comm_size(comm_fft)>1).and.&
766 &   (optgr==1.or.optstr==1.or.optdyfr==1.or.opteltfr==1)) then
767      call timab(48,1,tsec)
768      if (optv==1) then
769        if (optgr==1)then
770          call xmpi_sum(grv_indx,comm_fft,ierr)
771        end if
772        if (optstr==1)then
773          call xmpi_sum(strv,comm_fft,ierr)
774        end if
775        if (optdyfr==1)then
776          call xmpi_sum(dyfrv_indx,comm_fft,ierr)
777        end if
778      end if
779      if (optn==1) then
780        if (optgr==1)then
781          call xmpi_sum(grn_indx,comm_fft,ierr)
782        end if
783        if (optstr==1)then
784          call xmpi_sum(strn,comm_fft,ierr)
785        end if
786        if (optdyfr==1)then
787          call xmpi_sum(dyfrn_indx,comm_fft,ierr)
788        end if
789      end if
790      call timab(48,2,tsec)
791    end if
792  end if
793 
794 !Forces: re-order atoms
795  if (optgr==1) then
796    if (optv==1) then
797      do ia=1,natom
798        grv(1:3,atindx1(ia))=grv_indx(1:3,ia)
799      end do
800      ABI_DEALLOCATE(grv_indx)
801    end if
802    if (optn==1) then
803      do ia=1,natom
804        grn(1:3,atindx1(ia))=grn_indx(1:3,ia)
805      end do
806      ABI_DEALLOCATE(grn_indx)
807    end if
808  end if
809 
810 !Elastic tensor: Fill in lower triangle
811 ! if (opteltfr==1) then
812 !   do is2=2,6
813 !     do is1=1,is2-1
814 !       eltfrn(is2,is1)=eltfrn(is1,is2)
815 !     end do
816 !   end do
817 ! end if
818 
819 !Normalize stress tensor: 
820  if (optstr==1) then
821    if (optv==1) then
822      strv(:)=strv(:)/ucvol
823    end if
824    if (optn==1) strn(:)=strn(:)/ucvol
825  end if
826 
827 !Dynamical matrix: re-order atoms
828  if (optdyfr==1) then
829    if (optv==1) then
830      do ia=1,natom
831        dyfrv(1:3,1:3,atindx1(ia))=dyfrv_indx(1:3,1:3,ia)
832      end do
833      ABI_DEALLOCATE(dyfrv_indx)
834    end if
835    if (optn==1) then
836      do ia=1,natom
837        dyfrn(1:3,1:3,atindx1(ia))=dyfrn_indx(1:3,1:3,ia)
838      end do
839      ABI_DEALLOCATE(dyfrn_indx)
840    end if
841  end if
842 
843  if (.not.present(distribfft)) then
844    call destroy_distribfft(my_distribfft)
845    ABI_DATATYPE_DEALLOCATE(my_distribfft)
846  end if
847 
848  DBG_EXIT("COLL")
849 
850  contains
851 
852    function gsq_atm(i1,i2,i3)
853 
854 
855 !This section has been created automatically by the script Abilint (TD).
856 !Do not modify the following lines by hand.
857 #undef ABI_FUNC
858 #define ABI_FUNC 'gsq_atm'
859 !End of the abilint section
860 
861    real(dp) :: gsq_atm
862    integer,intent(in) :: i1,i2,i3
863    gsq_atm=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+dble(i3*i3)*gmet(3,3) &
864 &   +two*(dble(i1*i2)*gmet(1,2)+dble(i2*i3)*gmet(2,3)+dble(i3*i1)*gmet(3,1))
865  end function gsq_atm
866 
867    function dgsqds_atm(i1,i2,i3,is)
868 !Define dG^2/ds based on G space metric derivative
869 
870 !This section has been created automatically by the script Abilint (TD).
871 !Do not modify the following lines by hand.
872 #undef ABI_FUNC
873 #define ABI_FUNC 'dgsqds_atm'
874 !End of the abilint section
875 
876    real(dp) :: dgsqds_atm
877    integer,intent(in) :: i1,i2,i3,is
878    dgsqds_atm=dble(i1*i1)*dgm(1,1,is)+dble(i2*i2)*dgm(2,2,is)+&
879 &   dble(i3*i3)*dgm(3,3,is)+&
880 &   dble(i1*i2)*(dgm(1,2,is)+dgm(2,1,is))+&
881 &   dble(i1*i3)*(dgm(1,3,is)+dgm(3,1,is))+&
882 &   dble(i2*i3)*(dgm(2,3,is)+dgm(3,2,is))
883  end function dgsqds_atm
884 
885    function d2gsqds_atm(i1,i2,i3,is1,is2)
886 !  Define 2dG^2/ds1ds2  based on G space metric derivative
887 
888 !This section has been created automatically by the script Abilint (TD).
889 !Do not modify the following lines by hand.
890 #undef ABI_FUNC
891 #define ABI_FUNC 'd2gsqds_atm'
892 !End of the abilint section
893 
894    real(dp) :: d2gsqds_atm
895    integer,intent(in) :: i1,i2,i3,is1,is2
896    d2gsqds_atm=dble(i1*i1)*d2gm(1,1,is1,is2)+&
897 &   dble(i2*i2)*d2gm(2,2,is1,is2)+dble(i3*i3)*d2gm(3,3,is1,is2)+&
898 &   dble(i1*i2)*(d2gm(1,2,is1,is2)+d2gm(2,1,is1,is2))+&
899 &   dble(i1*i3)*(d2gm(1,3,is1,is2)+d2gm(3,1,is1,is2))+&
900 &   dble(i2*i3)*(d2gm(2,3,is1,is2)+d2gm(3,2,is1,is2))
901  end function d2gsqds_atm
902 
903 end subroutine atm2fft