TABLE OF CONTENTS


ABINIT/atm2fft [ Functions ]

[ Top ] [ Functions ]

NAME

 atm2fft

FUNCTION

 This routine sums atomic functions (density, kinetic 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 [kinetic] 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

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]
                   4: n^AT is the atomic PAW PS core kinetic density stored in array pawtab%ttaucorespl()
 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...

SOURCE

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

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

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

SOURCE

 995 subroutine dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,idir,ipert,&
 996 &                   mgfft,mqgrid,natom,ndir,nfft,ngfft,ntypat,&
 997 &                   ph1d,qgrid,qphon,typat,ucvol,usepaw,xred,psps,pawtab,&
 998 &                   atmrhor1,atmrhog1,atmvlocr1,atmvlocg1,distribfft,gauss,comm_fft,me_g0,optn_in,&
 999 &                   optn2_in,optv_in,paral_kgb,vspl) ! optional arguments
1000 
1001 !Arguments ------------------------------------
1002 !scalars
1003  integer,intent(in) :: cplex,idir,ipert,mgfft,mqgrid,natom,ndir,nfft,ntypat,usepaw
1004  integer,optional,intent(in) :: optn_in,optn2_in,optv_in
1005  integer,optional,intent(in) :: me_g0,comm_fft,paral_kgb
1006  real(dp),intent(in) :: gsqcut,ucvol
1007  type(pseudopotential_type),target,intent(in) :: psps
1008  type(distribfft_type),optional,intent(in),target :: distribfft
1009 !arrays
1010  integer,intent(in) :: atindx(natom),ngfft(18),typat(natom)
1011  real(dp),intent(in) :: gmet(3,3),gprimd(3,3)
1012  real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*natom),qgrid(mqgrid),qphon(3)
1013  real(dp),intent(in) :: xred(3,natom)
1014  real(dp),optional,intent(in)  :: gauss(2,ntypat),vspl(mqgrid,2,ntypat)
1015  real(dp),optional,intent(out) :: atmrhor1(cplex*nfft,ndir)
1016  real(dp),optional,intent(out) :: atmrhog1(2,nfft,ndir)
1017  real(dp),optional,intent(out) :: atmvlocr1(cplex*nfft,ndir)
1018  real(dp),optional,intent(out) :: atmvlocg1(2,nfft,ndir)
1019  type(pawtab_type),target,intent(in) :: pawtab(ntypat*usepaw)
1020 
1021 !Local variables ------------------------------
1022 !scalars
1023  integer,parameter :: im=2,re=1
1024  integer :: i1,i2,i3,ia,ia1,ia2,iatm,iatom,id,id1,id2,id3
1025  integer :: ig1,ig1max,ig1min,ig2,ig2max,ig2min,ig3,ig3max,ig3min
1026  integer :: ii,itypat,jj,me_fft,my_comm_fft,n1,n2,n3,nattyp,nproc_fft,ntype,paral_kgb_fft
1027  integer :: optn,optv,optn2,shift1,shift2,shift3,type1,type2
1028  logical :: have_g0,qeq0,qeq05
1029 #ifdef FC_NVHPC
1030 !Silly trick to prevent NVHPC optimization issue
1031  logical :: nothing=.false.
1032 #endif
1033  real(dp),parameter :: tolfix=1.0000001_dp
1034  real(dp) :: aa,alf2pi2,bb,cc,cutoff,dd,diff,dq,dq2div6,dqdiv6,dqm1,ee,ff
1035  real(dp) :: gauss1,gauss2,gmag,gq1,gq2,gq3,gsquar,n_at,dn_at,ph12i,ph12r,ph1i
1036  real(dp) :: ph1r,ph2i,ph2r,ph3i,ph3r,phqim,phqre,qxred2pi
1037  real(dp) :: sfi,sfqi,sfqr,sfr,term_n,term_v,v_at,dv_at,xnorm
1038  type(distribfft_type),pointer :: my_distribfft
1039  type(mpi_type) :: mpi_enreg_fft
1040 !arrays
1041  integer :: eps1(6)=(/1,2,3,2,3,1/),eps2(6)=(/1,2,3,3,1,2/),jdir(ndir)
1042  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:)
1043  real(dp), ABI_CONTIGUOUS pointer :: tvalespl(:,:),tcorespl(:,:)
1044  real(dp) ::  gq(6),gcart(3)
1045  real(dp),allocatable :: phim_igia(:),phre_igia(:),workn(:,:,:),workv(:,:,:)
1046 
1047 !no_abirules
1048 !Define G^2 based on G space metric gmet.
1049 ! gsq(g1,g2,g3)=g1*g1*gmet(1,1)+g2*g2*gmet(2,2)+g3*g3*gmet(3,3) &
1050 ! &       +two*(g1*g2*gmet(1,2)+g2*g3*gmet(2,3)+g3*g1*gmet(3,1))
1051 ! *************************************************************************
1052 
1053  DBG_ENTER("COLL")
1054 
1055 !  Check optional arguments
1056  if (present(comm_fft)) then
1057    if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then
1058      ABI_BUG('Need paral_kgb and me_g0 with comm_fft !')
1059    end if
1060  end if
1061 
1062  if (present(gauss))then
1063    if (.not.present(optn2_in)) then
1064      optn2 = 3
1065    else
1066      if(optn2_in/=3)then
1067        ABI_BUG('optn2 must be set to 3!')
1068      else
1069        optn2 = optn2_in
1070      end if
1071    end if
1072  end if
1073 
1074  if (present(atmrhor1))then
1075    if(.not.present(optn_in))then
1076      optn = 1
1077    else
1078      optn = optn_in
1079    end if
1080    if (.not.present(optn2_in)) then
1081      ABI_BUG('rho1 calculation need optn2 !')
1082    else
1083      optn2 = optn2_in
1084    end if
1085  else
1086    optn  = 0
1087    optn2 = 0
1088  end if
1089 
1090  if (present(atmvlocr1))then
1091    if(.not.present(optv_in))then
1092      optv = 1
1093    else
1094      if(optv_in/=1)then
1095        ABI_BUG('optv_in must be set to 1!')
1096      else
1097        optv = optv_in
1098      end if
1099    end if
1100    if(.not.present(vspl))then
1101      ABI_BUG('vloc1 calculation need vspl!')
1102    end if
1103  else
1104    optv  = 0
1105  end if
1106 
1107  if(ipert==natom+1.or.ipert==natom+2.or.ipert==natom+10.or.ipert==natom+11) then
1108 
1109 !  (In case of d/dk or an electric/magnetic field)
1110    if (optn==1) then
1111      atmrhor1(1:cplex*nfft,1:ndir)=zero
1112      if (present(atmrhog1)) atmrhog1 = zero
1113    end if
1114    if (optv==1) then
1115      atmvlocr1(1:cplex*nfft,1:ndir)=zero
1116      if (present(atmvlocg1)) atmvlocg1 = zero
1117    end if
1118 
1119  else
1120 
1121 !   Useful quantities
1122    if (ipert/=natom+3.and.ipert/=natom+4) then
1123      iatom=ipert;iatm=atindx(iatom)
1124      itypat=typat(iatom)
1125    else
1126       !sum of all (strain pertubation)
1127      iatom  = 1
1128      iatm   = 1
1129      itypat = 1
1130    end if
1131 
1132    n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1133    me_fft=ngfft(11)
1134    nproc_fft=ngfft(10)
1135    if (ndir==1) then
1136      jdir(1)=idir
1137    else
1138      do id=1,ndir
1139        jdir(id)=id
1140      end do
1141    end if
1142 
1143 !   Get the distrib associated with this fft_grid
1144    if (present(distribfft)) then
1145      my_distribfft => distribfft
1146    else
1147      ABI_MALLOC(my_distribfft,)
1148      call init_distribfft_seq(my_distribfft,'f',n2,n3,'fourdp')
1149    end if
1150    if (n2==my_distribfft%n2_coarse) then
1151      fftn2_distrib => my_distribfft%tab_fftdp2_distrib
1152    else if (n2 == my_distribfft%n2_fine) then
1153      fftn2_distrib => my_distribfft%tab_fftdp2dg_distrib
1154    else
1155      ABI_BUG("Unable to find an allocated distrib for this fft grid")
1156    end if
1157 
1158    qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)
1159    qeq05=(abs(abs(qphon(1))-half)<tol12.or. &
1160 &   abs(abs(qphon(2))-half)<tol12.or. &
1161 &   abs(abs(qphon(3))-half)<tol12)
1162 
1163    if (nproc_fft>1.and.qeq05) then
1164      ABI_ERROR('not compatible with FFT parallelism')
1165    end if
1166 
1167    if (optn2==3)then
1168      gauss1=gauss(1,itypat)
1169      gauss2=gauss(2,itypat)
1170      alf2pi2=(two_pi*gauss2)**2
1171    end if
1172 
1173    dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
1174    dqm1=one/dq
1175    dqdiv6=dq/six
1176    dq2div6=dq**2/six
1177    cutoff=gsqcut*tolfix
1178    id1=n1/2+2;id2=n2/2+2;id3=n3/2+2
1179    ig1max=-1;ig2max=-1;ig3max=-1
1180    ig1min=n1;ig2min=n2;ig3min=n3
1181 
1182 !   Determination of phase qxred*
1183    qxred2pi=two_pi*(qphon(1)*xred(1,iatom)+ &
1184 &   qphon(2)*xred(2,iatom)+ &
1185 &   qphon(3)*xred(3,iatom) )
1186    phqre=cos(qxred2pi)
1187    phqim=sin(qxred2pi)
1188 
1189 !   Zero out temporary arrays
1190    if (optv==1) then
1191      ABI_MALLOC(workv,(2,nfft,ndir))
1192      workv(:,:,:)=zero
1193    end if
1194    if (optn==1) then
1195      ABI_MALLOC(workn,(2,nfft,ndir))
1196      workn(:,:,:)=zero
1197    end if
1198 
1199    if (ipert==natom+3.or.ipert==natom+4) then
1200      ntype = ntypat
1201      ia1=1
1202      type1 = 1
1203      type2 = ntype
1204      ABI_MALLOC(phre_igia,(natom))
1205      ABI_MALLOC(phim_igia,(natom))
1206    else
1207      type1 = itypat
1208      type2 = itypat
1209      ntype = 1
1210      ABI_MALLOC(phre_igia,(iatm:iatm))
1211      ABI_MALLOC(phim_igia,(iatm:iatm))
1212    end if
1213 
1214 
1215    do itypat=type1,type2
1216 !    ia1,ia2 sets range of loop over atoms:
1217      if (ipert==natom+3.or.ipert==natom+4) then
1218        nattyp = count(typat(:)==itypat)
1219        ia2=ia1+nattyp-1
1220      else
1221        ia1 = iatm
1222        ia2 = iatm
1223      end if
1224 
1225      if (usepaw == 1) then
1226        tcorespl => pawtab(itypat)%tcorespl
1227        tvalespl => pawtab(itypat)%tvalespl
1228      else
1229        tcorespl => psps%nctab(itypat)%tcorespl
1230        tvalespl => psps%nctab(itypat)%tvalespl
1231      end if
1232 
1233      ii=0
1234      do i3=1,n3
1235        ig3=i3-(i3/id3)*n3-1
1236        gq3=dble(ig3)+qphon(3)
1237        gq(3)=gq3
1238 
1239        do i2=1,n2
1240          if (fftn2_distrib(i2)==me_fft) then
1241            ig2=i2-(i2/id2)*n2-1
1242            gq2=dble(ig2)+qphon(2)
1243            gq(2)=gq2
1244 
1245            do i1=1,n1
1246              ig1=i1-(i1/id1)*n1-1
1247              gq1=dble(ig1)+qphon(1)
1248              gq(1)=gq1
1249 
1250              ii=ii+1
1251 !            gsquar=gsq(gq1,gq2,gq3)
1252              gsquar=gq1*gq1*gmet(1,1)+gq2*gq2*gmet(2,2)+gq3*gq3*gmet(3,3) &
1253 &             +two*(gq1*gq2*gmet(1,2)+gq2*gq3*gmet(2,3)+gq3*gq1*gmet(3,1))
1254 
1255 
1256 !             Skip G**2 outside cutoff:
1257              if (gsquar<=cutoff) then
1258 
1259 !               Identify min/max indexes (to cancel unbalanced contributions later)
1260                if (qeq05) then
1261                  ig1max=max(ig1max,ig1);ig1min=min(ig1min,ig1)
1262                  ig2max=max(ig2max,ig2);ig2min=min(ig2min,ig2)
1263                  ig3max=max(ig3max,ig3);ig3min=min(ig3min,ig3)
1264                end if
1265 
1266                gmag=sqrt(gsquar)
1267                have_g0=(ig1==0.and.ig2==0.and.ig3==0.and.qeq0)
1268 
1269                jj=1+int(gmag*dqm1)
1270                diff=gmag-qgrid(jj)
1271 
1272 !               Compute structure factor
1273                phre_igia(:) = zero
1274                phim_igia(:) = zero
1275 
1276                do ia=ia1,ia2
1277 #ifdef FC_NVHPC
1278                  !Silly trick to prevent NVHPC optimization issue
1279                  if(nothing) write(100,*) shift1,shift2,shift3
1280 #endif                 
1281                  shift1=1+n1+(ia-1)*(2*n1+1)
1282                  shift2=1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1)
1283                  shift3=1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1)
1284                  ph1r=ph1d(re,ig1+shift1);ph1i=ph1d(im,ig1+shift1)
1285                  ph2r=ph1d(re,ig2+shift2);ph2i=ph1d(im,ig2+shift2)
1286                  ph3r=ph1d(re,ig3+shift3);ph3i=ph1d(im,ig3+shift3)
1287                  ph12r=ph1r*ph2r-ph1i*ph2i
1288                  ph12i=ph1r*ph2i+ph1i*ph2r
1289                  phre_igia(ia)=ph12r*ph3r-ph12i*ph3i
1290                  phim_igia(ia)=ph12r*ph3i+ph12i*ph3r
1291                end do
1292 !              Compute V^AT(g+q) and/or n^AT(g+q) for given type of atom
1293 !              Evaluate spline fit: p. 86 Numerical Recipes, Press et al;
1294 !              Note the error in book for sign of "aa" term in derivative.
1295                if (optv==1.or.optn2/=3) then
1296                  bb = diff*dqm1
1297                  aa = one-bb
1298                  cc = aa*(aa**2-one)*dq2div6
1299                  dd = bb*(bb**2-one)*dq2div6
1300                end if
1301 
1302                if (optv==1) then
1303                  if (have_g0) then
1304                    v_at=zero
1305                    dv_at=zero
1306                  else
1307                    v_at=(aa*vspl(jj,1,itypat)+bb*vspl(jj+1,1,itypat)+&
1308 &                   cc*vspl(jj,2,itypat)+dd*vspl(jj+1,2,itypat))/gsquar
1309 
1310 !                   Also get (dV(q)/dq)/q:
1311 !                   (note correction of Numerical Recipes sign error
1312 !                   before (3._dp*aa**2-1._dp)
1313                    ee= vspl(jj+1,1,itypat)-vspl(jj,1,itypat)
1314                    ff=  (3._dp*bb**2-1._dp)*vspl(jj+1,2,itypat) &
1315 &                   - (3._dp*aa**2-1._dp)*vspl(jj,2,itypat)
1316                    dv_at = ( ( ee*dqm1 + ff*dqdiv6 )/gmag&
1317 &                   - 2.0_dp*v_at) / gsquar
1318                  end if
1319                end if ! end optv
1320 
1321                if (optn==1) then
1322                  if (optn2==1) then
1323                    n_at=(aa*tcorespl(jj,1)+bb*tcorespl(jj+1,1)+cc*tcorespl(jj,2)+dd*tcorespl(jj+1,2))
1324                  else if (optn2==2) then
1325                    n_at=(aa*tvalespl(jj,1)+bb*tvalespl(jj+1,1)+cc*tvalespl(jj,2)+dd*tvalespl(jj+1,2))
1326                  else if (optn2==3) then
1327                    n_at=gauss1*exp(-gsquar*alf2pi2)
1328                  else
1329                    n_at=zero
1330                  end if
1331 
1332 !                Also get (dn^AT(q)/dq)/q:
1333                  if (have_g0) then
1334                    if (optn2==1) then
1335                      if (usepaw == 1) then
1336                        dn_at=pawtab(itypat)%dncdq0
1337                      else
1338                        dn_at=psps%nctab(itypat)%dncdq0
1339                      end if
1340                    else if (optn2==2) then
1341                      if (usepaw == 1) then
1342                        dn_at=pawtab(itypat)%dnvdq0
1343                      else
1344                        dn_at=psps%nctab(itypat)%dnvdq0
1345                      end if
1346                    else if (optn2==3) then
1347                      dn_at=-two*gauss1*alf2pi2
1348                    end if
1349                  else
1350                    if (optn2==1) then
1351                      ee=tcorespl(jj+1,1)-tcorespl(jj,1)
1352                      ff=(3._dp*bb**2-1._dp)*tcorespl(jj+1,2) &
1353 &                     -(3._dp*aa**2-1._dp)*tcorespl(jj,2)
1354                    else if (optn2==2) then
1355                      ee=tvalespl(jj+1,1)-tvalespl(jj,1)
1356                      ff=(3._dp*bb**2-1._dp)*tvalespl(jj+1,2) &
1357 &                     -(3._dp*aa**2-1._dp)*tvalespl(jj,2)
1358                    else if (optn2==3) then
1359                      dn_at=-two*gauss1*alf2pi2*exp(-gsquar*alf2pi2)
1360                    else
1361                    end if
1362                    dn_at=(ee*dqm1+ff*dqdiv6)/gmag
1363                  end if
1364                end if ! end optn
1365 
1366                do id=1,ndir
1367                  sfr=zero;sfi=zero
1368                  do ia=ia1,ia2
1369                    if (ipert==natom+3.or.ipert==natom+4) then
1370 !                    sum[Exp(-i.2pi.g.xred)]
1371                      sfr=sfr+phre_igia(ia)
1372                      sfi=sfi-phim_igia(ia)
1373                    else
1374 !                    Exp(-i.2pi.g.xred)  * -i.2pi.(g+q)
1375                      sfr=-(two_pi*gq(jdir(id))*phim_igia(ia))
1376                      sfi=-(two_pi*gq(jdir(id))*phre_igia(ia))
1377                    end if
1378                  end do
1379 
1380                  if (ipert/=natom+3.and.ipert/=natom+4) then
1381 !                  Phonons case
1382 
1383 !                  Exp(-i.2pi.q.xred)       => -i.2pi.(g+q).Exp(-i.2pi.(g+q).xred)
1384                    sfqr= sfr*phqre+sfi*phqim
1385                    sfqi=-sfr*phqim+sfi*phqre
1386 
1387                    if (optv == 1) then
1388                      workv(re,ii,id) = workv(re,ii,id) + sfqr*v_at
1389                      workv(im,ii,id) = workv(im,ii,id) + sfqi*v_at
1390                    end if
1391                    if (optn == 1) then
1392                      workn(re,ii,id) = workn(re,ii,id) + sfqr*n_at
1393                      workn(im,ii,id) = workn(im,ii,id) + sfqi*n_at
1394                    end if
1395 
1396                  else
1397 !                  Strain case
1398 
1399 !                  Compute G in cartesian coordinates
1400                    gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+&
1401 &                   gprimd(1,3)*dble(ig3)
1402                    gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+&
1403 &                   gprimd(2,3)*dble(ig3)
1404                    gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+&
1405 &                   gprimd(3,3)*dble(ig3)
1406 
1407 !                  Accumulate -dV^AT/dG*rho(G)*SF(G)*Gi.Gj/G
1408 !                  or -dn^AT/dG*V(G)*SF(G)*Gi.Gj/G
1409                    if (optv==1) then
1410                      if(jdir(id)<=3) then
1411                        term_v = dv_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) + v_at
1412                      else
1413                        term_v = dv_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id)))
1414                      end if
1415                      workv(re,ii,id) = workv(re,ii,id) - (sfr*term_v)
1416                      workv(im,ii,id) = workv(im,ii,id) - (sfi*term_v)
1417                    end if
1418 
1419                    if (optn==1) then
1420                      if(jdir(id)<=3) then
1421                        term_n = dn_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id))) + n_at
1422                      else
1423                        term_n = dn_at*gcart(eps1(jdir(id)))*gcart(eps2(jdir(id)))
1424                      end if
1425 
1426                      workn(re,ii,id) = workn(re,ii,id) - (sfr*term_n)
1427                      workn(im,ii,id) = workn(im,ii,id) - (sfi*term_n)
1428                    end if
1429 
1430                  end if
1431 !                End loop on ndir
1432 
1433                end do
1434 !              End skip G**2 outside cutoff
1435              end if
1436 !            End loop on n1, n2, n3
1437            end do
1438          end if ! this plane is selected
1439        end do
1440      end do
1441      ia1=ia2+1
1442    end do ! end loop itype
1443 
1444    ABI_FREE(phre_igia)
1445    ABI_FREE(phim_igia)
1446 
1447    if(ipert==natom+3.or.ipert==natom+4) then
1448 !    Set Vloc(G=0)=0:
1449      if (optv==1) then
1450        workv(re,1,:)=zero
1451        workv(im,1,:)=zero
1452      end if
1453    end if
1454 
1455 !  Identify unbalanced g-vectors
1456    if (qeq05) then  !This doesn't work in parallel
1457      ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2
1458      ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2
1459      ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2
1460      if (abs(abs(qphon(1))-half)<tol12) then
1461        if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max)
1462        if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min)
1463      end if
1464      if (abs(abs(qphon(2))-half)<tol12) then
1465        if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max)
1466        if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min)
1467      end if
1468      if (abs(abs(qphon(3))-half)<tol12) then
1469        if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max)
1470        if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min)
1471      end if
1472    end if
1473 
1474 
1475 !  Get 1st-order potential/density back to real space
1476 !  Non-symetrized non-zero elements have to be nullified
1477 !  Divide by unit cell volume
1478 
1479    if (optv==1.or.optn==1) then
1480 
1481      xnorm=one/ucvol
1482 !    Create fake mpi_enreg to wrap fourdp
1483      call initmpi_seq(mpi_enreg_fft)
1484      ABI_FREE(mpi_enreg_fft%distribfft)
1485      if (present(comm_fft)) then
1486        call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb)
1487        my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb
1488      else
1489        my_comm_fft=xmpi_comm_self;paral_kgb_fft=0;
1490        mpi_enreg_fft%distribfft => my_distribfft
1491      end if
1492 
1493      if (optv==1) then
1494        do id=1,ndir
1495 !        Eliminate unbalanced g-vectors
1496          if (qeq0) then       !q=0
1497            call zerosym(workv(:,:,id),2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft)
1498          else if (qeq05) then !q=1/2; this doesn't work in parallel
1499            call zerosym(workv(:,:,id),2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3)
1500          end if
1501          call fourdp(cplex,workv(:,:,id),atmvlocr1(:,id),1,mpi_enreg_fft,nfft,1,ngfft,0)
1502          atmvlocr1(:,id)=atmvlocr1(:,id)*xnorm
1503        end do
1504 
1505        !if (present(atmvlocg1)) atmvlocg1 = workv
1506        ABI_FREE(workv)
1507      end if
1508 
1509      if (optn==1) then
1510        do id=1,ndir
1511 !        Eliminate unbalanced g-vectors
1512          if (qeq0) then       !q=0
1513            call zerosym(workn(:,:,id),2,n1,n2,n3,comm_fft=my_comm_fft,distribfft=my_distribfft)
1514          else if (qeq05) then !q=1/2; this doesn't work in parallel
1515            call zerosym(workn(:,:,id),2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3)
1516          end if
1517          call fourdp(cplex,workn(:,:,id),atmrhor1(:,id),1,mpi_enreg_fft,nfft,1,ngfft,0)
1518          atmrhor1(:,id)=atmrhor1(:,id)*xnorm
1519        end do
1520        !if (present(atmrhog1)) atmrhog1 = workn
1521        ABI_FREE(workn)
1522      end if
1523 
1524 !    Destroy fake mpi_enreg
1525      call unset_mpi_enreg_fft(mpi_enreg_fft)
1526    end if
1527 
1528    if (.not.present(distribfft)) then
1529      call destroy_distribfft(my_distribfft)
1530      ABI_FREE(my_distribfft)
1531    end if
1532 
1533 !  End the condition of non-electric-field
1534  end if
1535 
1536  DBG_EXIT("COLL")
1537 
1538 end subroutine dfpt_atm2fft

ABINIT/m_atm2fft [ Modules ]

[ Top ] [ Modules ]

NAME

  m_atm2fft

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2022 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 .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_atm2fft
23 
24  use defs_basis
25  use m_abicore
26  use m_errors
27  use m_xmpi
28  use m_distribfft
29  use m_dtset
30 
31  use defs_abitypes, only : mpi_type
32  use m_time,        only : timab
33  use defs_datatypes,only : pseudopotential_type
34  use m_gtermcutoff, only : termcutoff
35  use m_pawtab,      only : pawtab_type
36  use m_fft,         only : zerosym, fourdp
37  use m_mpinfo,      only : set_mpi_enreg_fft, unset_mpi_enreg_fft, initmpi_seq
38 
39  implicit none
40 
41  private