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

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

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

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

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

ABINIT/m_atm2fft [ Modules ]

[ Top ] [ Modules ]

NAME

  m_atm2fft

FUNCTION

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 .

PARENTS

CHILDREN

SOURCE

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