TABLE OF CONTENTS


ABINIT/dfpt_atm2fft [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_atm2fft

FUNCTION

 This routine sums 1st-order atomic functions (density or potential)
 defined (in rec. space) on a radial grid to get global 1st-order
 quantities on the fine FFT grid.

 Possible options:
   optn=1: compute a sum of local 1st-order atomic densities
   optv=1: compute a sum of local 1st-order atomic potentials

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (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

  atindx(natom)=index table for atoms ordered by type
  cplex: if 1, real space 1-order functions on FFT grid
  distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism
  eei=local pseudopotential part of total energy
  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)
  idir=direction of atomic displacement (in case of phonons perturb.)
       used only if ndir=1 (see below)
  ipert=nindex of perturbation
  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.
  ndir=number of directions of atomic displacement (in case of phonon):
       can be 1 (idir direction in then used) or 3 (all directions)
       6 cartesian strain component 11,22,33,32,31,21 (in case strain perturbation)
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT
  nattyp(ntypat)=array describing how many atoms of each type in cell 
  ntypat=number of types of atoms.
  optn,optn2,optv= (see NOTES below)
  paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present)
  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.
  qphon(3)=wavevector of the phonon
  typat(natom)=type of each atom
  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)
  xred(3,natom)=reduced atomic coordinates
  psps <type(pseudopotential_type)>=variables related to pseudopotentials

OUTPUT

  ======= if optv==1 =======
    atmvlocr1(cplex*nfft,ndir)=sum of local 1st-order atomic potentials in real space
    atmvlocg1(2,nfft,ndir)=sum of local 1st-order atomic potentials in G-space
  ======= if optn==1 =======
   --- if optatm==1
    atmrhor1(cplex*nfft,ndir)=sum of 1st-order atomic densities in real space
    atmrhog1(2,nfft,ndir)=sum of 1st-order atomic densities in G-space

NOTES

 Details on possible options:
 ============================
 optv: controls the computation of a local 1st-order potential as sum of atomic potentials
          Vloc(r)=Sum_R[V1^AT(r-R)]
 optn: controls the computation of a 1st-order density as sum of atomic densities
          n(r)=Sum_R[n1^AT(r-R)]
          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

 Typical uses:
 =============
 Computation of:
  - 1st-order local potential: optv=1
  - 1st-order PS core density: optn=1, optn2=1

PARENTS

      dfpt_dyxc1,dfpt_eltfrxc,dfpt_looppert,dfpt_nstpaw,pawgrnl

CHILDREN

      destroy_distribfft,fourdp,init_distribfft_seq,initmpi_seq
      set_mpi_enreg_fft,unset_mpi_enreg_fft,zerosym

SOURCE

 96 #if defined HAVE_CONFIG_H
 97 #include "config.h"
 98 #endif
 99 
100 #include "abi_common.h"
101 
102 subroutine dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,&
103 &                   mgfft,mqgrid,natom,ndir,nfft,ngfft,ntypat,&
104 &                   ph1d,qgrid,qphon,typat,ucvol,usepaw,xred,psps,pawtab,&
105 &                   atmrhor1,atmrhog1,atmvlocr1,atmvlocg1,distribfft,gauss,comm_fft,me_g0,optn_in,&
106 &                   optn2_in,optv_in,paral_kgb,vspl) ! optional arguments
107 
108  use defs_basis
109  use defs_abitypes
110  use m_profiling_abi
111  use m_errors
112 
113  use m_xmpi,         only : xmpi_comm_self, xmpi_comm_size
114  use defs_datatypes, only : pseudopotential_type
115  use m_pawtab,       only : pawtab_type
116  use m_distribfft,   only : distribfft_type
117  use m_fft,          only : zerosym
118  use m_mpinfo,       only : set_mpi_enreg_fft, unset_mpi_enreg_fft
119 
120 !This section has been created automatically by the script Abilint (TD).
121 !Do not modify the following lines by hand.
122 #undef ABI_FUNC
123 #define ABI_FUNC 'dfpt_atm2fft'
124  use interfaces_51_manage_mpi
125  use interfaces_53_ffts
126 !End of the abilint section
127 
128  implicit none
129 
130 !Arguments ------------------------------------
131 !scalars
132  integer,intent(in) :: cplex,idir,ipert,mgfft,mqgrid,natom,ndir,nfft,ntypat,usepaw
133  integer,optional,intent(in) :: optn_in,optn2_in,optv_in
134  integer,optional,intent(in) :: me_g0,comm_fft,paral_kgb
135  real(dp),intent(in) :: gsqcut,ucvol
136  type(pseudopotential_type),target,intent(in) :: psps
137  type(distribfft_type),optional,intent(in),target :: distribfft
138 !arrays
139  integer,intent(in) :: atindx(natom),ngfft(18),typat(natom)
140  real(dp),intent(in) :: gmet(3,3),gprimd(3,3)
141  real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid),qphon(3)
142  real(dp),intent(in) :: xred(3,natom)
143  real(dp),optional,intent(in)  :: gauss(2,ntypat),vspl(mqgrid,2,ntypat)
144  real(dp),optional,intent(out) :: atmrhor1(cplex*nfft,ndir)
145  real(dp),optional,intent(out) :: atmrhog1(2,nfft,ndir)
146  real(dp),optional,intent(out) :: atmvlocr1(cplex*nfft,ndir)
147  real(dp),optional,intent(out) :: atmvlocg1(2,nfft,ndir)
148  type(pawtab_type),target,intent(in) :: pawtab(ntypat*usepaw)
149 
150 !Local variables ------------------------------
151 !scalars
152  integer,parameter :: im=2,re=1
153  integer :: i1,i2,i3,ia,ia1,ia2,iatm,iatom,id,id1,id2,id3
154  integer :: ig1,ig1max,ig1min,ig2,ig2max,ig2min,ig3,ig3max,ig3min
155  integer :: ii,itypat,jj,me_fft,my_comm_fft,n1,n2,n3,nattyp,nproc_fft,ntype,paral_kgb_fft
156  integer :: optn,optv,optn2,shift1,shift2,shift3,type1,type2
157  logical :: have_g0,qeq0,qeq05
158  real(dp),parameter :: tolfix=1.0000001_dp
159  real(dp) :: aa,alf2pi2,bb,cc,cutoff,dd,diff,dq,dq2div6,dqdiv6,dqm1,ee,ff
160  real(dp) :: gauss1,gauss2,gmag,gq1,gq2,gq3,gsquar,n_at,dn_at,ph12i,ph12r,ph1i
161  real(dp) :: ph1r,ph2i,ph2r,ph3i,ph3r,phqim,phqre,qxred2pi
162  real(dp) :: sfi,sfqi,sfqr,sfr,term_n,term_v,v_at,dv_at,xnorm
163  type(distribfft_type),pointer :: my_distribfft
164  type(mpi_type) :: mpi_enreg_fft
165 !arrays
166  integer :: eps1(6)=(/1,2,3,2,3,1/),eps2(6)=(/1,2,3,3,1,2/),jdir(ndir)
167  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:)
168  real(dp), ABI_CONTIGUOUS pointer :: tvalespl(:,:),tcorespl(:,:)
169  real(dp) ::  gq(6),gcart(3)
170  real(dp),allocatable :: phim_igia(:),phre_igia(:),workn(:,:,:),workv(:,:,:)
171 
172 !no_abirules
173 !Define G^2 based on G space metric gmet.
174 ! gsq(g1,g2,g3)=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+g3*g3*gmet(3,3) &
175 ! &       +two*(g1*g2*gmet(1,2)+g2*g3*gmet(2,3)+g3*g1*gmet(3,1))
176 ! *************************************************************************
177 
178  DBG_ENTER("COLL")
179 
180 !  Check optional arguments
181  if (present(comm_fft)) then
182    if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then
183      MSG_BUG('Need paral_kgb and me_g0 with comm_fft !')
184    end if
185  end if
186 
187  if (present(gauss))then
188    if (.not.present(optn2_in)) then   
189      optn2 = 3
190    else
191      if(optn2_in/=3)then
192        MSG_BUG('optn2 must be set to 3!')
193      else
194        optn2 = optn2_in
195      end if
196    end if     
197  end if
198 
199  if (present(atmrhor1))then
200    if(.not.present(optn_in))then
201      optn = 1
202    else
203      optn = optn_in   
204    end if
205    if (.not.present(optn2_in)) then   
206      MSG_BUG('rho1 calculation need optn2 !')
207    else
208      optn2 = optn2_in
209    end if
210  else
211    optn  = 0
212    optn2 = 0
213  end if
214  
215  if (present(atmvlocr1))then
216    if(.not.present(optv_in))then
217      optv = 1
218    else
219      if(optv_in/=1)then
220        MSG_BUG('optv_in must be set to 1!')
221      else
222        optv = optv_in
223      end if
224    end if
225    if(.not.present(vspl))then
226      MSG_BUG('vloc1 calculation need vspl!')
227    end if
228  else
229    optv  = 0 
230  end if
231  
232  if(ipert==natom+1.or.ipert==natom+2.or.ipert==natom+10.or.ipert==natom+11) then
233 
234 !  (In case of d/dk or an electric/magnetic field)
235    if (optn==1) then
236      atmrhor1(1:cplex*nfft,1:ndir)=zero
237      if (present(atmrhog1)) atmrhog1 = zero
238    end if
239    if (optv==1) then 
240      atmvlocr1(1:cplex*nfft,1:ndir)=zero
241      if (present(atmvlocg1)) atmvlocg1 = zero
242    end if
243 
244  else
245 
246 !   Useful quantities
247    if (ipert/=natom+3.and.ipert/=natom+4) then
248      iatom=ipert;iatm=atindx(iatom)
249      itypat=typat(iatom)
250    else
251       !sum of all (strain pertubation)
252      iatom  = 1
253      iatm   = 1
254      itypat = 1
255    end if
256 
257    n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
258    me_fft=ngfft(11)
259    nproc_fft=ngfft(10)
260    if (ndir==1) then
261      jdir(1)=idir
262    else
263      do id=1,ndir
264        jdir(id)=id
265      end do
266    end if
267 
268 !   Get the distrib associated with this fft_grid
269    if (present(distribfft)) then
270      my_distribfft => distribfft
271    else
272      ABI_DATATYPE_ALLOCATE(my_distribfft,)
273      call init_distribfft_seq(my_distribfft,'f',n2,n3,'fourdp')
274    end if
275    if (n2==my_distribfft%n2_coarse) then
276      fftn2_distrib => my_distribfft%tab_fftdp2_distrib
277    else if (n2 == my_distribfft%n2_fine) then
278      fftn2_distrib => my_distribfft%tab_fftdp2dg_distrib
279    else
280      MSG_BUG("Unable to find an allocated distrib for this fft grid")
281    end if
282 
283    qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)
284    qeq05=(abs(abs(qphon(1))-half)<tol12.or. &
285 &   abs(abs(qphon(2))-half)<tol12.or. &
286 &   abs(abs(qphon(3))-half)<tol12)
287 
288    if (nproc_fft>1.and.qeq05) then
289      MSG_ERROR('not compatible with FFT parallelism')
290    end if
291 
292    if (optn2==3)then
293      gauss1=gauss(1,itypat)
294      gauss2=gauss(2,itypat)
295      alf2pi2=(two_pi*gauss2)**2
296    end if
297 
298    dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
299    dqm1=one/dq
300    dqdiv6=dq/six
301    dq2div6=dq**2/six
302    cutoff=gsqcut*tolfix
303    id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
304    ig1max=-1;ig2max=-1;ig3max=-1
305    ig1min=n1;ig2min=n2;ig3min=n3
306 
307 !   Determination of phase qxred*
308    qxred2pi=two_pi*(qphon(1)*xred(1,iatom)+ &
309 &   qphon(2)*xred(2,iatom)+ &
310 &   qphon(3)*xred(3,iatom) )
311    phqre=cos(qxred2pi)
312    phqim=sin(qxred2pi)
313 
314 !   Zero out temporary arrays
315    if (optv==1) then
316      ABI_ALLOCATE(workv,(2,nfft,ndir))
317      workv(:,:,:)=zero
318    end if
319    if (optn==1) then
320      ABI_ALLOCATE(workn,(2,nfft,ndir))
321      workn(:,:,:)=zero
322    end if
323 
324    if (ipert==natom+3.or.ipert==natom+4) then
325      ntype = ntypat
326      ia1=1
327      type1 = 1
328      type2 = ntype
329      ABI_ALLOCATE(phre_igia,(natom))
330      ABI_ALLOCATE(phim_igia,(natom))
331    else
332      type1 = itypat
333      type2 = itypat
334      ntype = 1
335      ABI_ALLOCATE(phre_igia,(iatm:iatm))
336      ABI_ALLOCATE(phim_igia,(iatm:iatm))
337    end if
338 
339 
340    do itypat=type1,type2
341 !    ia1,ia2 sets range of loop over atoms:
342      if (ipert==natom+3.or.ipert==natom+4) then
343        nattyp = count(typat(:)==itypat)
344        ia2=ia1+nattyp-1
345      else
346        ia1 = iatm
347        ia2 = iatm
348      end if
349 
350      if (usepaw == 1) then
351        tcorespl => pawtab(itypat)%tcorespl
352        tvalespl => pawtab(itypat)%tvalespl
353      else
354        tcorespl => psps%nctab(itypat)%tcorespl
355        tvalespl => psps%nctab(itypat)%tvalespl
356      end if
357 
358      ii=0
359      do i3=1,n3
360        ig3=i3-(i3/id3)*n3-1
361        gq3=dble(ig3)+qphon(3)
362        gq(3)=gq3
363 
364        do i2=1,n2
365          if (fftn2_distrib(i2)==me_fft) then
366            ig2=i2-(i2/id2)*n2-1
367            gq2=dble(ig2)+qphon(2)
368            gq(2)=gq2
369            
370            do i1=1,n1
371              ig1=i1-(i1/id1)*n1-1
372              gq1=dble(ig1)+qphon(1)
373              gq(1)=gq1
374 
375              ii=ii+1
376 !            gsquar=gsq(gq1,gq2,gq3)
377              gsquar=gq1*gq1*gmet(1,1)+gq2*gq2*gmet(2,2)+gq3*gq3*gmet(3,3) &
378 &             +two*(gq1*gq2*gmet(1,2)+gq2*gq3*gmet(2,3)+gq3*gq1*gmet(3,1))
379 
380              
381 !             Skip G**2 outside cutoff:
382              if (gsquar<=cutoff) then
383                
384 !               Identify min/max indexes (to cancel unbalanced contributions later)
385                if (qeq05) then
386                  ig1max=max(ig1max,ig1);ig1min=min(ig1min,ig1)
387                  ig2max=max(ig2max,ig2);ig2min=min(ig2min,ig2)
388                  ig3max=max(ig3max,ig3);ig3min=min(ig3min,ig3)
389                end if
390                
391                gmag=sqrt(gsquar)
392                have_g0=(ig1==0.and.ig2==0.and.ig3==0.and.qeq0)
393 
394                jj=1+int(gmag*dqm1)
395                diff=gmag-qgrid(jj)
396                
397 !               Compute structure factor
398                phre_igia(:) = zero
399                phim_igia(:) = zero
400                
401                do ia=ia1,ia2
402                  shift1=1+n1+(ia-1)*(2*n1+1)
403                  shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)
404                  shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)
405                  ph1r=ph1d(re,ig1+shift1);ph1i=ph1d(im,ig1+shift1)
406                  ph2r=ph1d(re,ig2+shift2);ph2i=ph1d(im,ig2+shift2)
407                  ph3r=ph1d(re,ig3+shift3);ph3i=ph1d(im,ig3+shift3)
408                  ph12r=ph1r*ph2r-ph1i*ph2i
409                  ph12i=ph1r*ph2i+ph1i*ph2r
410                  phre_igia(ia)=ph12r*ph3r-ph12i*ph3i
411                  phim_igia(ia)=ph12r*ph3i+ph12i*ph3r
412                end do
413 !              Compute V^AT(g+q) and/or n^AT(g+q) for given type of atom
414 !              Evaluate spline fit: p. 86 Numerical Recipes, Press et al;
415 !              Note the error in book for sign of "aa" term in derivative.
416                if (optv==1.or.optn2/=3) then
417                  bb = diff*dqm1
418                  aa = one-bb
419                  cc = aa*(aa**2-one)*dq2div6
420                  dd = bb*(bb**2-one)*dq2div6
421                end if
422 
423                if (optv==1) then
424                  if (have_g0) then
425                    v_at=zero
426                    dv_at=zero
427                  else
428                    v_at=(aa*vspl(jj,1,itypat)+bb*vspl(jj+1,1,itypat)+&
429 &                   cc*vspl(jj,2,itypat)+dd*vspl(jj+1,2,itypat))/gsquar
430 
431 !                   Also get (dV(q)/dq)/q:
432 !                   (note correction of Numerical Recipes sign error
433 !                   before (3._dp*aa**2-1._dp)
434                    ee= vspl(jj+1,1,itypat)-vspl(jj,1,itypat)
435                    ff=  (3._dp*bb**2-1._dp)*vspl(jj+1,2,itypat) &
436 &                   - (3._dp*aa**2-1._dp)*vspl(jj,2,itypat)
437                    dv_at = ( ( ee*dqm1 + ff*dqdiv6 )/gmag&
438 &                   - 2.0_dp*v_at) / gsquar                   
439                  end if
440                end if ! end optv
441 
442                if (optn==1) then
443                  if (optn2==1) then
444                    n_at=(aa*tcorespl(jj,1)+bb*tcorespl(jj+1,1)+cc*tcorespl(jj,2)+dd*tcorespl(jj+1,2))
445                  else if (optn2==2) then
446                    n_at=(aa*tvalespl(jj,1)+bb*tvalespl(jj+1,1)+cc*tvalespl(jj,2)+dd*tvalespl(jj+1,2))
447                  else if (optn2==3) then
448                    n_at=gauss1*exp(-gsquar*alf2pi2)
449                  else
450                    n_at=zero
451                  end if
452 
453 !                Also get (dn^AT(q)/dq)/q:
454                  if (have_g0) then
455                    if (optn2==1) then
456                      if (usepaw == 1) then
457                        dn_at=pawtab(itypat)%dncdq0
458                      else
459                        dn_at=psps%nctab(itypat)%dncdq0
460                      end if
461                    else if (optn2==2) then
462                      if (usepaw == 1) then
463                        dn_at=pawtab(itypat)%dnvdq0
464                      else
465                        dn_at=psps%nctab(itypat)%dnvdq0
466                      end if
467                    else if (optn2==3) then
468                      dn_at=-two*gauss1*alf2pi2
469                    end if
470                  else
471                    if (optn2==1) then
472                      ee=tcorespl(jj+1,1)-tcorespl(jj,1)
473                      ff=(3._dp*bb**2-1._dp)*tcorespl(jj+1,2) &
474 &                     -(3._dp*aa**2-1._dp)*tcorespl(jj,2)
475                    else if (optn2==2) then
476                      ee=tvalespl(jj+1,1)-tvalespl(jj,1)
477                      ff=(3._dp*bb**2-1._dp)*tvalespl(jj+1,2) &
478 &                     -(3._dp*aa**2-1._dp)*tvalespl(jj,2)
479                    else if (optn2==3) then
480                      dn_at=-two*gauss1*alf2pi2*exp(-gsquar*alf2pi2)
481                    else
482                    end if
483                    dn_at=(ee*dqm1+ff*dqdiv6)/gmag
484                  end if
485                end if ! end optn
486                
487                do id=1,ndir
488                  sfr=zero;sfi=zero                
489                  do ia=ia1,ia2
490                    if (ipert==natom+3.or.ipert==natom+4) then
491 !                    sum[Exp(-i.2pi.g.xred)]  
492                      sfr=sfr+phre_igia(ia)
493                      sfi=sfi-phim_igia(ia)
494                    else
495 !                    Exp(-i.2pi.g.xred)  * -i.2pi.(g+q)
496                      sfr=-(two_pi*gq(jdir(id))*phim_igia(ia))
497                      sfi=-(two_pi*gq(jdir(id))*phre_igia(ia))
498                    end if
499                  end do
500 
501                  if (ipert/=natom+3.and.ipert/=natom+4) then
502 !                  Phonons case
503 
504 !                  Exp(-i.2pi.q.xred)       => -i.2pi.(g+q).Exp(-i.2pi.(g+q).xred)
505                    sfqr= sfr*phqre+sfi*phqim
506                    sfqi=-sfr*phqim+sfi*phqre
507 
508                    if (optv == 1) then
509                      workv(re,ii,id) = workv(re,ii,id) + sfqr*v_at
510                      workv(im,ii,id) = workv(im,ii,id) + sfqi*v_at
511                    end if
512                    if (optn == 1) then
513                      workn(re,ii,id) = workn(re,ii,id) + sfqr*n_at
514                      workn(im,ii,id) = workn(im,ii,id) + sfqi*n_at
515                    end if
516 
517                  else
518 !                  Strain case
519                    
520 !                  Compute G in cartesian coordinates
521                    gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+&
522 &                   gprimd(1,3)*dble(ig3)
523                    gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+&
524 &                   gprimd(2,3)*dble(ig3)
525                    gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+&
526 &                   gprimd(3,3)*dble(ig3)
527 
528 !                  Accumulate -dV^AT/dG*rho(G)*SF(G)*Gi.Gj/G
529 !                  or -dn^AT/dG*V(G)*SF(G)*Gi.Gj/G              
530                    if (optv==1) then
531                      if(jdir(id)<=3) then
532                        term_v = dv_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) + v_at
533                      else
534                        term_v = dv_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id)))
535                      end if
536                      workv(re,ii,id) = workv(re,ii,id) - (sfr*term_v)
537                      workv(im,ii,id) = workv(im,ii,id) - (sfi*term_v)
538                    end if
539 
540                    if (optn==1) then
541                      if(jdir(id)<=3) then
542                        term_n = dn_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) + n_at
543                      else
544                        term_n = dn_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id)))
545                      end if
546 
547                      workn(re,ii,id) = workn(re,ii,id) - (sfr*term_n)
548                      workn(im,ii,id) = workn(im,ii,id) - (sfi*term_n)
549                    end if
550 
551                  end if
552 !                End loop on ndir
553 
554                end do
555 !              End skip G**2 outside cutoff
556              end if
557 !            End loop on n1, n2, n3
558            end do
559          end if ! this plane is selected
560        end do
561      end do
562      ia1=ia2+1
563    end do ! end loop itype
564 
565    ABI_DEALLOCATE(phre_igia)
566    ABI_DEALLOCATE(phim_igia)
567 
568    if(ipert==natom+3.or.ipert==natom+4) then
569 !    Set Vloc(G=0)=0:
570      if (optv==1) then
571        workv(re,1,:)=zero
572        workv(im,1,:)=zero
573      end if
574    end if
575    
576 !  Identify unbalanced g-vectors
577    if (qeq05) then  !This doesn't work in parallel
578      ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2
579      ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2
580      ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2
581      if (abs(abs(qphon(1))-half)<tol12) then
582        if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max)
583        if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min)
584      end if
585      if (abs(abs(qphon(2))-half)<tol12) then
586        if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max)
587        if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min)
588      end if
589      if (abs(abs(qphon(3))-half)<tol12) then
590        if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max)
591        if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min)
592      end if
593    end if
594 
595 
596 !  Get 1st-order potential/density back to real space
597 !  Non-symetrized non-zero elements have to be nullified
598 !  Divide by unit cell volume
599 
600    if (optv==1.or.optn==1) then
601 
602      xnorm=one/ucvol
603 !    Create fake mpi_enreg to wrap fourdp
604      call initmpi_seq(mpi_enreg_fft)
605      ABI_DATATYPE_DEALLOCATE(mpi_enreg_fft%distribfft)
606      if (present(comm_fft)) then
607        call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb)
608        my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb
609      else
610        my_comm_fft=xmpi_comm_self;paral_kgb_fft=0;
611        mpi_enreg_fft%distribfft => my_distribfft
612      end if
613      
614      if (optv==1) then
615        do id=1,ndir
616 !        Eliminate unbalanced g-vectors
617          if (qeq0) then       !q=0
618            call zerosym(workv(:,:,id),2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft)
619          else if (qeq05) then !q=1/2; this doesn't work in parallel
620            call zerosym(workv(:,:,id),2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3)
621          end if
622          call fourdp(cplex,workv(:,:,id),atmvlocr1(:,id),1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0)  
623          atmvlocr1(:,id)=atmvlocr1(:,id)*xnorm
624        end do
625 
626        !if (present(atmvlocg1)) atmvlocg1 = workv
627        ABI_DEALLOCATE(workv)
628      end if
629      
630      if (optn==1) then
631        do id=1,ndir
632 !        Eliminate unbalanced g-vectors
633          if (qeq0) then       !q=0
634            call zerosym(workn(:,:,id),2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft)
635          else if (qeq05) then !q=1/2; this doesn't work in parallel
636            call zerosym(workn(:,:,id),2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3)
637          end if
638          call fourdp(cplex,workn(:,:,id),atmrhor1(:,id),1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0)
639          atmrhor1(:,id)=atmrhor1(:,id)*xnorm
640        end do
641        !if (present(atmrhog1)) atmrhog1 = workn
642        ABI_DEALLOCATE(workn)
643      end if
644 
645 !    Destroy fake mpi_enreg
646      call unset_mpi_enreg_fft(mpi_enreg_fft)
647    end if
648    
649    if (.not.present(distribfft)) then
650      call destroy_distribfft(my_distribfft)
651      ABI_DATATYPE_DEALLOCATE(my_distribfft)
652    end if
653 
654 !  End the condition of non-electric-field
655  end if
656 
657  DBG_EXIT("COLL")
658  
659 end subroutine dfpt_atm2fft