TABLE OF CONTENTS


ABINIT/d2kindstr2 [ Functions ]

[ Top ] [ Functions ]

NAME

 d2kindstr2

FUNCTION

 compute expectation value of the second derivatives of the kinetic energy
 wrt strain for one band and kpoint

INPUTS

  cwavef(2,npw*nspinor)=wavefunction for current band
  ecut=cut-off energy for plane wave basis sphere (Ha)
  ecutsm=smearing energy for plane wave kinetic energy (Ha)
  effmass_free=effective mass for electrons (1. in common case)
  gmet(3,3)=reciprocal lattice metric tensor ($\textrm{Bohr}^{-2}$)
  gprimd(3,3)=primitive vectors in reciprocal space
  istwfk=information about wavefunction storage
  kg_k(3,npw)=integer coordinates of planewaves in basis sphere.
  kpt(3)=reduced coordinates of k point
  npw=number of plane waves at kpt.
  nspinor=number of spinorial components of the wavefunction

OUTPUT

  ekinout(36)=expectation values of the second strain derivatives
   of the (modified) kinetic energy

NOTES

 Usually, the kinetic energy expression is $(1/2) (2 \pi)^2 (k+G)^2 $
 However, the present implementation allows for a modification
 of this kinetic energy, in order to obtain smooth total energy
 curves with respect to the cut-off energy or the cell size and shape.
 Thus the usual expression is kept if it is lower then ecut-ecutsm,
 zero is returned beyond ecut, and in between, the kinetic
 energy is DIVIDED by a smearing factor (to make it infinite at the
 cut-off energy). The smearing factor is $x^2 (3-2x)$, where
 x = (ecut- unmodified energy)/ecutsm.

PARENTS

      dfpt_eltfrkin

CHILDREN

SOURCE

1694 subroutine d2kindstr2(cwavef,ecut,ecutsm,effmass_free,ekinout,gmet,gprimd,&
1695 &            istwfk,kg_k,kpt,npw,nspinor)
1696 
1697 
1698 !This section has been created automatically by the script Abilint (TD).
1699 !Do not modify the following lines by hand.
1700 #undef ABI_FUNC
1701 #define ABI_FUNC 'd2kindstr2'
1702 !End of the abilint section
1703 
1704  implicit none
1705 
1706 !Arguments ------------------------------------
1707 !scalars
1708  integer,intent(in) :: istwfk,npw,nspinor
1709  real(dp),intent(in) :: ecut,ecutsm,effmass_free
1710 !arrays
1711  integer,intent(in) :: kg_k(3,npw)
1712  real(dp),intent(in) :: cwavef(2,npw*nspinor),gmet(3,3),gprimd(3,3),kpt(3)
1713  real(dp),intent(inout) :: ekinout(36) !vz_i
1714 
1715 !Local variables-------------------------------
1716 !scalars
1717  integer,parameter :: im=2,re=1
1718  integer :: ig,igs,ii,ispinor,istr1,istr2,ka,kb,kd,kg
1719  real(dp) :: d2fkin,d2fsm,d2kinacc,d2kpg2,dfkin,dfsm,dkpg21,dkpg22,ecutsm_inv
1720  real(dp) :: fsm,gpk1,gpk2,gpk3,htpisq,kpg2,term,xx
1721 !arrays
1722  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1723  real(dp) :: d2gm(3,3),dgm01(3,3),dgm10(3,3)
1724 
1725 ! *************************************************************************
1726 !
1727 !htpisq is (1/2) (2 Pi) **2:
1728    htpisq=0.5_dp*(two_pi)**2
1729 
1730    ecutsm_inv=0.0_dp
1731    if(ecutsm>1.0d-20)ecutsm_inv=1/ecutsm
1732 
1733 !Loop over 2nd strain index
1734    do istr2=1,6
1735 !  Loop over 1st strain index, upper triangle only
1736      do istr1=1,istr2
1737 
1738        ka=idx(2*istr1-1);kb=idx(2*istr1);kg=idx(2*istr2-1);kd=idx(2*istr2)
1739 
1740        do ii = 1,3
1741          dgm01(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
1742          dgm10(:,ii)=-(gprimd(kg,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kg,ii))
1743        end do
1744 
1745        d2gm(:,:)=0._dp
1746        do ii = 1,3
1747          if(ka==kg) d2gm(:,ii)=d2gm(:,ii)&
1748 &         +gprimd(kb,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kb,ii)
1749          if(ka==kd) d2gm(:,ii)=d2gm(:,ii)&
1750 &         +gprimd(kb,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(kb,ii)
1751          if(kb==kg) d2gm(:,ii)=d2gm(:,ii)&
1752 &         +gprimd(ka,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(ka,ii)
1753          if(kb==kd) d2gm(:,ii)=d2gm(:,ii)&
1754 &         +gprimd(ka,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(ka,ii)
1755        end do
1756        d2gm(:,:)=0.5_dp*d2gm(:,:)
1757 
1758        d2kinacc=0._dp
1759 
1760 !    loop on spinor index
1761        do ispinor=1,nspinor
1762          igs=(ispinor-1)*npw
1763 !      loop on plane waves
1764          do ig=1,npw
1765            gpk1=dble(kg_k(1,ig))+kpt(1)
1766            gpk2=dble(kg_k(2,ig))+kpt(2)
1767            gpk3=dble(kg_k(3,ig))+kpt(3)
1768            kpg2=htpisq*&
1769 &           ( gmet(1,1)*gpk1**2+         &
1770 &           gmet(2,2)*gpk2**2+         &
1771 &           gmet(3,3)*gpk3**2          &
1772 &           +2.0_dp*(gpk1*gmet(1,2)*gpk2+  &
1773 &           gpk1*gmet(1,3)*gpk3+  &
1774 &           gpk2*gmet(2,3)*gpk3 )  )
1775            dkpg21=htpisq*&
1776 &           ( dgm01(1,1)*gpk1**2+         &
1777 &           dgm01(2,2)*gpk2**2+         &
1778 &           dgm01(3,3)*gpk3**2          &
1779 &           +2.0_dp*(gpk1*dgm01(1,2)*gpk2+  &
1780 &           gpk1*dgm01(1,3)*gpk3+  &
1781 &           gpk2*dgm01(2,3)*gpk3 )  )
1782            dkpg22=htpisq*&
1783 &           ( dgm10(1,1)*gpk1**2+         &
1784 &           dgm10(2,2)*gpk2**2+         &
1785 &           dgm10(3,3)*gpk3**2          &
1786 &           +2.0_dp*(gpk1*dgm10(1,2)*gpk2+  &
1787 &           gpk1*dgm10(1,3)*gpk3+  &
1788 &           gpk2*dgm10(2,3)*gpk3 )  )
1789            d2kpg2=htpisq*&
1790 &           ( d2gm(1,1)*gpk1**2+         &
1791 &           d2gm(2,2)*gpk2**2+         &
1792 &           d2gm(3,3)*gpk3**2          &
1793 &           +2.0_dp*(gpk1*d2gm(1,2)*gpk2+  &
1794 &           gpk1*d2gm(1,3)*gpk3+  &
1795 &           gpk2*d2gm(2,3)*gpk3 )  )
1796 
1797            if(kpg2>ecut-tol12)then
1798              dfkin=0._dp
1799              d2fkin=0._dp
1800            elseif(kpg2>ecut-ecutsm)then
1801 !          This kinetic cutoff smoothing function and its xx derivatives
1802 !          were produced with Mathematica and the fortran code has been
1803 !          numerically checked against Mathematica.
1804              xx=(ecut-kpg2)*ecutsm_inv
1805              fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
1806              dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
1807              d2fsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*&
1808 &             (-144+45*xx))))))*fsm**3
1809              dfkin=fsm-ecutsm_inv*kpg2*dfsm
1810              d2fkin=ecutsm_inv*(-2.0_dp*dfsm+ecutsm_inv*kpg2*d2fsm)
1811            else
1812              dfkin=1._dp
1813              d2fkin=0._dp
1814            end if
1815 
1816 !        accumulate kinetic energy 2nd derivative with wavefunction components
1817            term=d2fkin*dkpg21*dkpg22 + dfkin*d2kpg2
1818            if(istwfk==2 .and. ig/=1)term=2.0_dp*term
1819            if(istwfk>2)term=2.0_dp*term
1820            d2kinacc=d2kinacc + term*(cwavef(re,ig+igs)**2 + cwavef(im,ig+igs)**2)
1821 
1822          end do  !ig
1823        end do !ispinor
1824 
1825        ekinout(istr1+6*(istr2-1))=d2kinacc/effmass_free
1826 
1827      end do !istr1
1828    end do !istr2
1829 
1830   end subroutine d2kindstr2

ABINIT/dfpt_eltfrhar [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_eltfrhar

FUNCTION

 Compute the frozen-wavefunction hartree enegy contribution to the elastic tensor

INPUTS

  rprimd(3,3)=dimensional primitive translation vectors (bohr)
  gsqcut =Fourier cutoff on G^2 for "large sphere" of radius double
   that of the basis sphere--appropriate for charge density rho(G),
   Hartree potential, and pseudopotentials
  mpi_enreg=informations about MPI parallelization
  nfft =(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT,
    see ~abinit/doc/variables/vargs.htm#ngfft
  rhog(2,nfft)=total electron density in G space

OUTPUT

  eltfrhar(6,6)=non-symmetrized kinetic energy contribution to the
                    elastic tensor

NOTES

 *based largely on hartre.f

PARENTS

      respfn

CHILDREN

      metric,ptabs_fourdp,timab,xmpi_sum

SOURCE

1868 subroutine dfpt_eltfrhar(eltfrhar,rprimd,gsqcut,mpi_enreg,nfft,ngfft,rhog)
1869 
1870 
1871 !This section has been created automatically by the script Abilint (TD).
1872 !Do not modify the following lines by hand.
1873 #undef ABI_FUNC
1874 #define ABI_FUNC 'dfpt_eltfrhar'
1875 !End of the abilint section
1876 
1877  implicit none
1878 
1879 !Arguments ------------------------------------
1880 !scalars
1881  integer,intent(in) :: nfft
1882  real(dp),intent(in) :: gsqcut
1883  type(MPI_type),intent(in) :: mpi_enreg
1884 !arrays
1885  integer,intent(in) :: ngfft(18)
1886  real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3)
1887  real(dp),intent(out) :: eltfrhar(6,6)
1888 
1889 !Local variables-------------------------------
1890 !scalars
1891  integer,parameter :: im=2,re=1
1892  integer :: i1,i2,i23,i3,id2,id3,ierr,ig,ig2,ig3,ii,ii1,ing,istr1,istr2,jj
1893  integer :: ka,kb,kd,kg,me_fft,n1,n2,n3,nproc_fft
1894  real(dp),parameter :: tolfix=1.000000001_dp
1895  real(dp) :: cutoff,d2eacc,d2etot,d2gs,deacc01,deacc10,dgs01,dgs10,eacc,fact,gs
1896  real(dp) :: term,ucvol
1897 !arrays
1898  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1899  integer :: id(3)
1900  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
1901  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1902  real(dp) :: d2gm(3,3),dgm01(3,3),dgm10(3,3),gmet(3,3),gprimd(3,3),gqr(3)
1903  real(dp) :: rmet(3,3),tsec(2)
1904  real(dp),allocatable :: gq(:,:)
1905 
1906 ! *************************************************************************
1907 
1908 !Compute gmet, gprimd and ucvol from rprimd
1909  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1910 
1911  eltfrhar(:,:)=0.0_dp
1912 
1913  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1914  me_fft=ngfft(11)
1915  nproc_fft=ngfft(10)
1916 
1917 !Get the distrib associated with this fft_grid
1918  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1919 
1920 !Initialize a few quantities
1921  fact=0.5_dp*ucvol/pi
1922  cutoff=gsqcut*tolfix
1923 
1924 !In order to speed the routine, precompute the components of g+q
1925 !Also check if the booked space was large enough...
1926  ABI_ALLOCATE(gq,(3,max(n1,n2,n3)))
1927  do ii=1,3
1928    id(ii)=ngfft(ii)/2+2
1929    do ing=1,ngfft(ii)
1930      ig=ing-(ing/id(ii))*ngfft(ii)-1
1931      gq(ii,ing)=ig
1932    end do
1933  end do
1934 
1935 !Loop over 2nd strain index
1936  do istr2=1,6
1937 !  Loop over 1st strain index, upper triangle only
1938    do istr1=1,istr2
1939 
1940      ka=idx(2*istr1-1);kb=idx(2*istr1);kg=idx(2*istr2-1);kd=idx(2*istr2)
1941 
1942      do ii = 1,3
1943        dgm01(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
1944        dgm10(:,ii)=-(gprimd(kg,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kg,ii))
1945      end do
1946 
1947      d2gm(:,:)=0._dp
1948      do ii = 1,3
1949        if(ka==kg) d2gm(:,ii)=d2gm(:,ii)&
1950 &       +gprimd(kb,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kb,ii)
1951        if(ka==kd) d2gm(:,ii)=d2gm(:,ii)&
1952 &       +gprimd(kb,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(kb,ii)
1953        if(kb==kg) d2gm(:,ii)=d2gm(:,ii)&
1954 &       +gprimd(ka,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(ka,ii)
1955        if(kb==kd) d2gm(:,ii)=d2gm(:,ii)&
1956 &       +gprimd(ka,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(ka,ii)
1957      end do
1958      d2gm(:,:)=0.5_dp*d2gm(:,:)
1959 
1960 !    initialize energy accumulator
1961      eacc=0._dp
1962      deacc01=0._dp
1963      deacc10=0._dp
1964      d2eacc=0._dp
1965 
1966      id2=n2/2+2
1967      id3=n3/2+2
1968 !    Triple loop on each dimension
1969      do i3=1,n3
1970        ig3=i3-(i3/id3)*n3-1
1971        gqr(3)=gq(3,i3)
1972        do i2=1,n2
1973          if (fftn2_distrib(i2)==me_fft) then
1974            gqr(2)=gq(2,i2)
1975            ig2=i2-(i2/id2)*n2-1
1976            i23=n1*(ffti2_local(i2)-1 +(n2/nproc_fft)*(i3-1))
1977 !          Do the test that eliminates the Gamma point outside
1978 !          of the inner loop
1979            ii1=1
1980            if(i23==0 .and. ig2==0 .and. ig3==0)then
1981              ii1=2
1982            end if
1983 
1984 !          Final inner loop on the first dimension
1985 !          (note the lower limit)
1986            do i1=ii1,n1
1987              gqr(1)=gq(1,i1)
1988              gs=(gmet(1,1)*gqr(1)*gqr(1)+gmet(2,2)*gqr(2)*gqr(2)+&
1989 &             gmet(3,3)*gqr(3)*gqr(3)+2._dp*&
1990 &             (gmet(1,2)*gqr(1)*gqr(2) + gmet(1,3)*gqr(1)*gqr(3)+&
1991 &             gmet(2,3)*gqr(2)*gqr(3)) )
1992              ii=i1+i23
1993              if(gs<=cutoff)then
1994                dgs01=(dgm01(1,1)*gqr(1)*gqr(1)+dgm01(2,2)*gqr(2)*gqr(2)+&
1995 &               dgm01(3,3)*gqr(3)*gqr(3)+2._dp*&
1996 &               (dgm01(1,2)*gqr(1)*gqr(2) + dgm01(1,3)*gqr(1)*gqr(3)+&
1997 &               dgm01(2,3)*gqr(2)*gqr(3)) )
1998                dgs10=(dgm10(1,1)*gqr(1)*gqr(1)+dgm10(2,2)*gqr(2)*gqr(2)+&
1999 &               dgm10(3,3)*gqr(3)*gqr(3)+2._dp*&
2000 &               (dgm10(1,2)*gqr(1)*gqr(2) + dgm10(1,3)*gqr(1)*gqr(3)+&
2001 &               dgm10(2,3)*gqr(2)*gqr(3)) )
2002                d2gs =(d2gm(1,1)*gqr(1)*gqr(1)+d2gm(2,2)*gqr(2)*gqr(2)+&
2003 &               d2gm(3,3)*gqr(3)*gqr(3)+2._dp*&
2004 &               (d2gm(1,2)*gqr(1)*gqr(2) + d2gm(1,3)*gqr(1)*gqr(3)+&
2005 &               d2gm(2,3)*gqr(2)*gqr(3)) )
2006 
2007                term=(rhog(re,ii)**2+rhog(im,ii)**2)/gs
2008                eacc=eacc+term
2009                deacc01=deacc01+dgs01*term/gs
2010                deacc10=deacc10+dgs10*term/gs
2011                d2eacc=d2eacc+(-d2gs+2._dp*dgs01*dgs10/gs)*term/gs
2012              end if
2013 
2014 !            End loop on i1
2015            end do
2016          end if
2017 !        End loop on i2
2018        end do
2019 !      End loop on i3
2020      end do
2021 
2022 !    Add contributions taking account diagonal strain terms (from ucvol
2023 !    derivatives)
2024      d2etot=d2eacc
2025      if(istr1<=3) d2etot=d2etot+deacc10
2026      if(istr2<=3) d2etot=d2etot+deacc01
2027      if(istr1<=3 .and. istr2<=3) d2etot=d2etot+eacc
2028 
2029      eltfrhar(istr1,istr2)=fact*d2etot
2030 
2031 !    End loop on istr1
2032    end do
2033 !  End loop in istr2
2034  end do
2035 
2036  ABI_DEALLOCATE(gq)
2037 
2038 !Init mpi_comm
2039  call timab(48,1,tsec)
2040  call xmpi_sum(eltfrhar,mpi_enreg%comm_fft,ierr)
2041  call timab(48,2,tsec)
2042 
2043 !Fill in lower triangle
2044  do jj=2,6
2045    do ii=1,jj-1
2046      eltfrhar(jj,ii)=eltfrhar(ii,jj)
2047    end do
2048  end do
2049 end subroutine dfpt_eltfrhar

ABINIT/dfpt_eltfrkin [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_eltfrkin

FUNCTION

 Compute the frozen-wavefunction kinetic enegy contribution to the
 elastic tensor

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol)=<G|Cnk>=Fourier coefficients of wavefunction
  ecut=cut-off energy for plane wave basis sphere (Ha)
  ecutsm=smearing energy for plane wave kinetic energy (Ha) (NOT NEEDED !)
  effmass_free=effective mass for electrons (1. in common case)
  istwfk(nkpt)=input option parameter that describes the storage of wfs
  kg(3,mpw*mkmem)=work array for coordinates of G vectors in basis
  kptns(3,nkpt)=coordinates of k points in terms of reciprocal space
   primitive translations
  mband=maximum number of bands
  mgfft=maximum size of 1D FFTs
  mkmem=number of k points treated by this node.
  mpi_enreg=information about MPI parallelization
  mpw=maximum dimension for number of planewaves
  nband(nkpt*nsppol)=number of bands being considered per k point
  nkpt=number of k points
  ngfft(18)=contain all needed information about 3D FFT, i
    see ~abinit/doc/variables/vargs.htm#ngfft
  npwarr(nkpt)=number of planewaves at each k point, and boundary
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for polarized
  occ(mband*nkpt*nsppol)=occupation numbers of bands (usually 2)
    at each k point
  rprimd(3,3)=dimensional real space primitive translations (bohr)
  wtk(nkpt)=k point weights

OUTPUT

  eltfrkin(6,6)=non-symmetrized kinetic energy contribution to the
                    elastic tensor

PARENTS

      respfn

CHILDREN

      d2kindstr2,metric,sphereboundary,timab,xmpi_sum

SOURCE

1500 subroutine dfpt_eltfrkin(cg,eltfrkin,ecut,ecutsm,effmass_free,&
1501 &  istwfk,kg,kptns,mband,mgfft,mkmem,mpi_enreg,&
1502 &  mpw,nband,nkpt,ngfft,npwarr,nspinor,nsppol,occ,rprimd,wtk)
1503 
1504 
1505 !This section has been created automatically by the script Abilint (TD).
1506 !Do not modify the following lines by hand.
1507 #undef ABI_FUNC
1508 #define ABI_FUNC 'dfpt_eltfrkin'
1509 !End of the abilint section
1510 
1511  implicit none
1512 
1513 !Arguments ------------------------------------
1514 !scalars
1515  integer,intent(in) :: mband,mgfft,mkmem,mpw,nkpt,nspinor,nsppol
1516  real(dp),intent(in) :: ecut,ecutsm,effmass_free
1517  type(MPI_type),intent(in) :: mpi_enreg
1518 !arrays
1519  integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),nband(nkpt*nsppol)
1520  integer,intent(in) :: ngfft(18),npwarr(nkpt)
1521  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),kptns(3,nkpt)
1522  real(dp),intent(in) :: occ(mband*nkpt*nsppol),rprimd(3,3),wtk(nkpt)
1523  real(dp),intent(out) :: eltfrkin(6,6)
1524 
1525 !Local variables-------------------------------
1526 !scalars
1527  integer :: bdtot_index,iband,icg,ierr,ii,ikg
1528  integer :: ikpt,index,ipw,isppol,istwf_k,jj,master,me,n1,n2
1529  integer :: n3,nband_k,nkinout,npw_k,spaceComm
1530  real(dp) :: ucvol
1531 !arrays
1532  integer,allocatable :: gbound(:,:),kg_k(:,:)
1533  real(dp) :: gmet(3,3),gprimd(3,3),kpoint(3),rmet(3,3),tsec(2)
1534  real(dp),allocatable :: cwavef(:,:),ekinout(:)
1535  real(dp),allocatable :: eltfrkink(:,:)
1536 
1537 ! *************************************************************************
1538 
1539  DBG_ENTER("COLL")
1540 
1541 !Default for sequential use
1542  master=0
1543 !Init mpi_comm
1544  spaceComm=mpi_enreg%comm_cell
1545  me=mpi_enreg%me_kpt
1546 
1547 !Compute gmet, gprimd and ucvol from rprimd
1548  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1549 
1550  eltfrkin(:,:)=0.0_dp
1551  bdtot_index=0
1552  icg=0
1553 
1554  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1555  ABI_ALLOCATE(kg_k,(3,mpw))
1556  ABI_ALLOCATE(cwavef,(2,mpw*nspinor))
1557  ABI_ALLOCATE(eltfrkink,(6,6))
1558 
1559 !Define k-points distribution
1560 
1561 !LOOP OVER SPINS
1562  do isppol=1,nsppol
1563    ikg=0
1564 
1565 !  Loop over k points
1566    do ikpt=1,nkpt
1567 
1568      nband_k=nband(ikpt+(isppol-1)*nkpt)
1569      istwf_k=istwfk(ikpt)
1570      npw_k=npwarr(ikpt)
1571 
1572 !    Skip this k-point if not the proper processor
1573      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
1574        bdtot_index=bdtot_index+nband_k
1575        cycle
1576      end if
1577 
1578      ABI_ALLOCATE(gbound,(2*mgfft+8,2))
1579      kpoint(:)=kptns(:,ikpt)
1580 
1581      kg_k(:,:) = 0
1582 
1583 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(ikg,kg,kg_k,npw_k)
1584      do ipw=1,npw_k
1585        kg_k(1,ipw)=kg(1,ipw+ikg)
1586        kg_k(2,ipw)=kg(2,ipw+ikg)
1587        kg_k(3,ipw)=kg(3,ipw+ikg)
1588      end do
1589 
1590      call sphereboundary(gbound,istwf_k,kg_k,mgfft,npw_k)
1591 
1592      index=1+icg
1593 
1594      eltfrkink(:,:)=0.0_dp
1595 
1596      nkinout=6*6
1597      ABI_ALLOCATE(ekinout,(nkinout))
1598      ekinout(:)=zero
1599 
1600      do iband=1,nband_k
1601 
1602        if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= me) cycle
1603 
1604        cwavef(:,1:npw_k*nspinor)=cg(:,1+(iband-1)*npw_k*nspinor+icg:iband*npw_k*nspinor+icg)
1605 
1606        call d2kindstr2(cwavef,ecut,ecutsm,effmass_free,ekinout,gmet,gprimd,&
1607 &       istwf_k,kg_k,kpoint,npw_k,nspinor)
1608 
1609        eltfrkink(:,:)=eltfrkink(:,:)+ occ(iband+bdtot_index)* reshape(ekinout(:), (/6,6/) )
1610 
1611      end do !iband
1612 
1613      ABI_DEALLOCATE(ekinout)
1614 
1615      eltfrkin(:,:)=eltfrkin(:,:)+wtk(ikpt)*eltfrkink(:,:)
1616 
1617      ABI_DEALLOCATE(gbound)
1618 
1619      bdtot_index=bdtot_index+nband_k
1620 
1621      if (mkmem/=0) then
1622 !      Handle case in which kg, cg, are kept in core
1623        icg=icg+npw_k*nspinor*nband_k
1624        ikg=ikg+npw_k
1625      end if
1626 
1627    end do
1628  end do  ! End loops on isppol and ikpt
1629 
1630 !Fill in lower triangle
1631  do jj=2,6
1632    do ii=1,jj-1
1633      eltfrkin(jj,ii)=eltfrkin(ii,jj)
1634    end do
1635  end do
1636 
1637 !Accumulate eltfrkin on all proc.
1638  call timab(48,1,tsec)
1639  call xmpi_sum(eltfrkin,spaceComm,ierr)
1640  call timab(48,2,tsec)
1641 
1642  ABI_DEALLOCATE(cwavef)
1643  ABI_DEALLOCATE(eltfrkink)
1644  ABI_DEALLOCATE(kg_k)
1645 
1646  DBG_EXIT("COLL")
1647 
1648   contains

ABINIT/dfpt_eltfrloc [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_eltfrloc

FUNCTION

 Compute the frozen-wavefunction local pseudopotential contribution
 to the elastic tensor and the internal strain (derivative wrt one
 cartesian strain component and one reduced-coordinate atomic displacement).

INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
  gmet(3,3)=reciprocal space metric (bohr^-2)
  gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr**-1)
  gsqcut=cutoff on G^2 based on ecut
  mgfft=maximum size of 1D FFTs
  mpi_enreg=informations about MPI parallelization
  mqgrid=dimensioned number of q grid points for local psp spline
  natom=number of atoms in unit cell
  nattyp(ntypat)=number of atoms of each type
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT,
    see ~abinit/doc/variables/vargs.htm#ngfft
  ntypat=number of types of atoms
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  qgrid(mqgrid)=q point array for local psp spline fits
  rhog(2,nfft)=electron density in G space
  vlspl(mqgrid,2,ntypat)=q^2 v(q) spline for each type of atom.

OUTPUT

  eltfrloc(6+3*natom,6)=non-symmetrized local pseudopotenial contribution
   to the elastic tensor and internal strain.

PARENTS

      respfn

CHILDREN

      ptabs_fourdp,timab,xmpi_sum

SOURCE

1076 subroutine dfpt_eltfrloc(atindx,eltfrloc,gmet,gprimd,gsqcut,mgfft,&
1077 &  mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,ntypat,ph1d,qgrid,rhog,vlspl)
1078 
1079 
1080 !This section has been created automatically by the script Abilint (TD).
1081 !Do not modify the following lines by hand.
1082 #undef ABI_FUNC
1083 #define ABI_FUNC 'dfpt_eltfrloc'
1084 !End of the abilint section
1085 
1086  implicit none
1087 
1088 !Arguments ------------------------------------
1089 !scalars
1090  integer,intent(in) :: mgfft,mqgrid,natom,nfft,ntypat
1091  real(dp),intent(in) :: gsqcut
1092  type(MPI_type),intent(in) :: mpi_enreg
1093 !arrays
1094  integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18)
1095  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom)
1096  real(dp),intent(in) :: qgrid(mqgrid),rhog(2,nfft),vlspl(mqgrid,2,ntypat)
1097  real(dp),intent(out) :: eltfrloc(6+3*natom,6)
1098 
1099 !Local variables-------------------------------
1100 !scalars
1101  integer,parameter :: im=2,re=1
1102  integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ielt,ieltx,ierr,ig1,ig2,ig3,ii
1103  integer :: is1,is2,itypat,jj,ka,kb,kd,kg,me_fft,n1,n2,n3,nproc_fft
1104  real(dp),parameter :: tolfix=1.0000001_dp
1105  real(dp) :: aa,bb,cc,cutoff,d2g
1106  real(dp) :: dd,dg1,dg2,diff,dq
1107  real(dp) :: dq2div6,dqdiv6,dqm1,ee,ff,gmag,gsquar
1108  real(dp) :: sfi,sfr,term,term1
1109 !real(dp) :: ph1_elt,ph2_elt,ph3_elt,phi_elt,phr_elt
1110  real(dp) :: term2,term3,term4,term5,vion1,vion2,vion3
1111 !arrays
1112  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1113  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
1114  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1115  real(dp) :: dgm(3,3,6),tsec(2)
1116  real(dp),allocatable :: d2gm(:,:,:,:),elt_work(:,:)
1117 
1118 ! *************************************************************************
1119 
1120  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1121  me_fft=ngfft(11) ; nproc_fft=ngfft(10)
1122 
1123 !Get the distrib associated with this fft_grid
1124  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1125 
1126 !-----
1127 !Compute 1st and 2nd derivatives of metric tensor wrt all strain components
1128 !and store for use in inner loop below.
1129  ABI_ALLOCATE(d2gm,(3,3,6,6))
1130 
1131 !Loop over 2nd strain index
1132  do is2=1,6
1133    kg=idx(2*is2-1);kd=idx(2*is2)
1134    do jj = 1,3
1135      dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj))
1136    end do
1137 !  Loop over 1st strain index, upper triangle only
1138    do is1=1,is2
1139      ka=idx(2*is1-1);kb=idx(2*is1)
1140      d2gm(:,:,is1,is2)=0._dp
1141      do jj = 1,3
1142        if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1143 &       +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj)
1144        if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1145 &       +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj)
1146        if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1147 &       +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj)
1148        if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1149 &       +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj)
1150      end do
1151    end do !is1
1152  end do !is2
1153 
1154 !Zero out array to permit accumulation over atom types below:
1155  eltfrloc(:,:)=0.0_dp
1156 
1157  dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
1158  dqm1=1.0_dp/dq
1159  dqdiv6=dq/6.0_dp
1160  dq2div6=dq**2/6.0_dp
1161  cutoff=gsqcut*tolfix
1162  id1=n1/2+2
1163  id2=n2/2+2
1164  id3=n3/2+2
1165 
1166  ia1=1
1167  do itypat=1,ntypat
1168 !  ia1,ia2 sets range of loop over atoms:
1169    ia2=ia1+nattyp(itypat)-1
1170    ii=0
1171    do i3=1,n3
1172      ig3=i3-(i3/id3)*n3-1
1173      do i2=1,n2
1174        if (fftn2_distrib(i2)==me_fft) then
1175          ig2=i2-(i2/id2)*n2-1
1176          do i1=1,n1
1177            ig1=i1-(i1/id1)*n1-1
1178 
1179            ii=ii+1
1180 !          Skip G=0:
1181            if (ig1==0 .and. ig2==0 .and. ig3==0) cycle
1182 
1183 !          Skip G**2 outside cutoff:
1184            gsquar=gsq_elt(ig1,ig2,ig3)
1185            if (gsquar<=cutoff) then
1186              gmag=sqrt(gsquar)
1187 
1188 !            Compute vion(G) for given type of atom
1189              jj=1+int(gmag*dqm1)
1190              diff=gmag-qgrid(jj)
1191 
1192 !            Evaluate spline fit from q^2 V(q) to get V(q):
1193 !            (p. 86 Numerical Recipes, Press et al; NOTE error in book for sign
1194 !             of "aa" term in derivative; also see splfit routine).
1195              bb = diff*dqm1
1196              aa = 1.0_dp-bb
1197              cc = aa*(aa**2-1.0_dp)*dq2div6
1198              dd = bb*(bb**2-1.0_dp)*dq2div6
1199              term1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) +&
1200 &             cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat))
1201              vion1=term1 / gsquar
1202 
1203 !            Also get dV(q)/dq:
1204 !            (note correction of Numerical Recipes sign error
1205 !             before (3._dp*aa**2-1._dp)
1206              ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat)
1207              ff=  (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) &
1208 &             - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat)
1209              term2 = ee*dqm1 + ff*dqdiv6
1210              vion2 = term2/gsquar - 2._dp*term1/(gsquar*gmag)
1211 
1212 !            Also get V''(q)
1213              term3=aa*vlspl(jj,2,itypat)+bb*vlspl(jj+1,2,itypat)
1214              vion3 = (term3 - 4.0_dp*term2/gmag + 6._dp*term1/gsquar)/gsquar
1215 
1216 !            Assemble structure factor over all atoms of given type:
1217              sfr=zero;sfi=zero
1218              do ia=ia1,ia2
1219                sfr=sfr+phre_elt(ig1,ig2,ig3,ia)
1220                sfi=sfi-phimag_elt(ig1,ig2,ig3,ia)
1221              end do
1222              term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi)
1223 
1224 !            Loop over 2nd strain index
1225              do is2=1,6
1226                dg2=0.5_dp*dgsqds_elt(ig1,ig2,ig3,is2)/gmag
1227 !              Loop over 1st strain index, upper triangle only
1228                do is1=1,is2
1229                  dg1=0.5_dp*dgsqds_elt(ig1,ig2,ig3,is1)/gmag
1230                  d2g=(0.25_dp*d2gsqds_elt(ig1,ig2,ig3,is1,is2)-dg1*dg2)/gmag
1231                  eltfrloc(is1,is2)=eltfrloc(is1,is2)+&
1232 &                 term*(vion3*dg1*dg2+vion2*d2g)
1233                  if(is2<=3)&
1234 &                 eltfrloc(is1,is2)=eltfrloc(is1,is2)-term*vion2*dg1
1235                  if(is1<=3)&
1236 &                 eltfrloc(is1,is2)=eltfrloc(is1,is2)-term*vion2*dg2
1237                  if(is1<=3 .and. is2<=3)&
1238 &                 eltfrloc(is1,is2)=eltfrloc(is1,is2)+term*vion1
1239                end do !is1
1240 
1241 !              Internal strain section - loop over current atoms
1242                do ia=ia1,ia2
1243                  if(is2 <=3) then
1244                    term4=vion2*dg2-vion1
1245                  else
1246                    term4=vion2*dg2
1247                  end if
1248                  term5=-two_pi*(rhog(re,ii)*phimag_elt(ig1,ig2,ig3,ia)&
1249 &                 +rhog(im,ii)*phre_elt(ig1,ig2,ig3,ia))*term4
1250                  eltfrloc(7+3*(ia-1),is2)=eltfrloc(7+3*(ia-1),is2)+term5*dble(ig1)
1251                  eltfrloc(8+3*(ia-1),is2)=eltfrloc(8+3*(ia-1),is2)+term5*dble(ig2)
1252                  eltfrloc(9+3*(ia-1),is2)=eltfrloc(9+3*(ia-1),is2)+term5*dble(ig3)
1253                end do
1254 
1255              end do !is2
1256 
1257 !            End skip G**2 outside cutoff:
1258            end if
1259 
1260 !          End loop on n1, n2, n3. There is a "cycle" inside the loop
1261          end do
1262        end if
1263      end do
1264    end do
1265 
1266 !  End loop on type of atoms
1267    ia1=ia2+1
1268  end do
1269 !Init mpi_comm
1270  call timab(48,1,tsec)
1271  call xmpi_sum(eltfrloc,mpi_enreg%comm_fft,ierr)
1272  call timab(48,2,tsec)
1273 
1274 !Fill in lower triangle
1275  do is2=2,6
1276    do is1=1,is2-1
1277      eltfrloc(is2,is1)=eltfrloc(is1,is2)
1278    end do
1279  end do
1280 
1281 !The indexing array atindx is used to reestablish the correct
1282 !order of atoms
1283  ABI_ALLOCATE(elt_work,(6+3*natom,6))
1284  elt_work(1:6,1:6)=eltfrloc(1:6,1:6)
1285  do ia=1,natom
1286    ielt=7+3*(ia-1)
1287    ieltx=7+3*(atindx(ia)-1)
1288    elt_work(ielt:ielt+2,1:6)=eltfrloc(ieltx:ieltx+2,1:6)
1289  end do
1290  eltfrloc(:,:)=elt_work(:,:)
1291 
1292  ABI_DEALLOCATE(d2gm)
1293  ABI_DEALLOCATE(elt_work)
1294 
1295  contains
1296 
1297 !Real and imaginary parts of phase.
1298    function phr_elt(x1,y1,x2,y2,x3,y3)
1299 
1300 
1301 !This section has been created automatically by the script Abilint (TD).
1302 !Do not modify the following lines by hand.
1303 #undef ABI_FUNC
1304 #define ABI_FUNC 'phr_elt'
1305 !End of the abilint section
1306 
1307    real(dp) :: phr_elt,x1,x2,x3,y1,y2,y3
1308    phr_elt=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
1309  end function phr_elt
1310 
1311    function phi_elt(x1,y1,x2,y2,x3,y3)
1312 
1313 
1314 !This section has been created automatically by the script Abilint (TD).
1315 !Do not modify the following lines by hand.
1316 #undef ABI_FUNC
1317 #define ABI_FUNC 'phi_elt'
1318 !End of the abilint section
1319 
1320    real(dp):: phi_elt,x1,x2,x3,y1,y2,y3
1321    phi_elt=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
1322  end function phi_elt
1323 
1324    function ph1_elt(nri,ig1,ia)
1325 
1326 
1327 !This section has been created automatically by the script Abilint (TD).
1328 !Do not modify the following lines by hand.
1329 #undef ABI_FUNC
1330 #define ABI_FUNC 'ph1_elt'
1331 !End of the abilint section
1332 
1333    real(dp):: ph1_elt
1334    integer :: nri,ig1,ia
1335    ph1_elt=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1))
1336  end function ph1_elt
1337 
1338    function ph2_elt(nri,ig2,ia)
1339 
1340 
1341 !This section has been created automatically by the script Abilint (TD).
1342 !Do not modify the following lines by hand.
1343 #undef ABI_FUNC
1344 #define ABI_FUNC 'ph2_elt'
1345 !End of the abilint section
1346 
1347    real(dp):: ph2_elt
1348    integer :: nri,ig2,ia
1349    ph2_elt=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1))
1350  end function ph2_elt
1351 
1352    function ph3_elt(nri,ig3,ia)
1353 
1354 
1355 !This section has been created automatically by the script Abilint (TD).
1356 !Do not modify the following lines by hand.
1357 #undef ABI_FUNC
1358 #define ABI_FUNC 'ph3_elt'
1359 !End of the abilint section
1360 
1361    real(dp):: ph3_elt
1362    integer :: nri,ig3,ia
1363    ph3_elt=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
1364  end function ph3_elt
1365 
1366    function phre_elt(ig1,ig2,ig3,ia)
1367 
1368 
1369 !This section has been created automatically by the script Abilint (TD).
1370 !Do not modify the following lines by hand.
1371 #undef ABI_FUNC
1372 #define ABI_FUNC 'phre_elt'
1373 !End of the abilint section
1374 
1375    real(dp):: phre_elt
1376    integer :: ig1,ig2,ig3,ia
1377    phre_elt=phr_elt(ph1_elt(re,ig1,ia),ph1_elt(im,ig1,ia),&
1378 &   ph2_elt(re,ig2,ia),ph2_elt(im,ig2,ia),ph3_elt(re,ig3,ia),ph3_elt(im,ig3,ia))
1379  end function phre_elt
1380 
1381    function phimag_elt(ig1,ig2,ig3,ia)
1382 
1383 
1384 !This section has been created automatically by the script Abilint (TD).
1385 !Do not modify the following lines by hand.
1386 #undef ABI_FUNC
1387 #define ABI_FUNC 'phimag_elt'
1388 !End of the abilint section
1389 
1390    real(dp) :: phimag_elt
1391    integer :: ig1,ig2,ig3,ia
1392    phimag_elt=phi_elt(ph1_elt(re,ig1,ia),ph1_elt(im,ig1,ia),&
1393 &   ph2_elt(re,ig2,ia),ph2_elt(im,ig2,ia),ph3_elt(re,ig3,ia),ph3_elt(im,ig3,ia))
1394  end function phimag_elt
1395 
1396    function gsq_elt(i1,i2,i3)
1397 
1398 
1399 !This section has been created automatically by the script Abilint (TD).
1400 !Do not modify the following lines by hand.
1401 #undef ABI_FUNC
1402 #define ABI_FUNC 'gsq_elt'
1403 !End of the abilint section
1404 
1405    real(dp) :: gsq_elt
1406    integer :: i1,i2,i3
1407 !Define G^2 based on G space metric gmet.
1408    gsq_elt=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+&
1409 &   dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+&
1410 &   dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
1411  end function gsq_elt
1412 
1413    function dgsqds_elt(i1,i2,i3,is)
1414 
1415 
1416 !This section has been created automatically by the script Abilint (TD).
1417 !Do not modify the following lines by hand.
1418 #undef ABI_FUNC
1419 #define ABI_FUNC 'dgsqds_elt'
1420 !End of the abilint section
1421 
1422    real(dp) :: dgsqds_elt
1423    integer :: i1,i2,i3,is
1424 !Define dG^2/ds based on G space metric derivative
1425    dgsqds_elt=dble(i1*i1)*dgm(1,1,is)+dble(i2*i2)*dgm(2,2,is)+&
1426 &   dble(i3*i3)*dgm(3,3,is)+&
1427 &   dble(i1*i2)*(dgm(1,2,is)+dgm(2,1,is))+&
1428 &   dble(i1*i3)*(dgm(1,3,is)+dgm(3,1,is))+&
1429 &   dble(i2*i3)*(dgm(2,3,is)+dgm(3,2,is))
1430  end function dgsqds_elt
1431 
1432    function d2gsqds_elt(i1,i2,i3,is1,is2)
1433 
1434 
1435 !This section has been created automatically by the script Abilint (TD).
1436 !Do not modify the following lines by hand.
1437 #undef ABI_FUNC
1438 #define ABI_FUNC 'd2gsqds_elt'
1439 !End of the abilint section
1440 
1441    real(dp) :: d2gsqds_elt
1442    integer :: i1,i2,i3,is1,is2
1443 !Define 2dG^2/ds1ds2  based on G space metric derivative
1444    d2gsqds_elt=dble(i1*i1)*d2gm(1,1,is1,is2)+&
1445 &   dble(i2*i2)*d2gm(2,2,is1,is2)+dble(i3*i3)*d2gm(3,3,is1,is2)+&
1446 &   dble(i1*i2)*(d2gm(1,2,is1,is2)+d2gm(2,1,is1,is2))+&
1447 &   dble(i1*i3)*(d2gm(1,3,is1,is2)+d2gm(3,1,is1,is2))+&
1448 &   dble(i2*i3)*(d2gm(2,3,is1,is2)+d2gm(3,2,is1,is2))
1449  end function d2gsqds_elt
1450 
1451 end subroutine dfpt_eltfrloc

ABINIT/dfpt_eltfrxc [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_eltfrxc

FUNCTION

 Compute the 2nd derivatives of exchange-correlation energy
 with respect to all pairs of strain and strain-atomic displacement
 for the frozen wavefunction contribution to the elastic
 and internal strain tensors

INPUTS

  atindx(natom)=index table for atoms ordered by type
  dtset <type(dataset_type)>=all input variables for this dataset
   | natom=number of atoms in unit cell
   | nfft=(effective) number of FFT grid points (for this processor)
   | nspden=number of spin-density components
   | ntypat=number of types of atoms in cell.
   | typat(natom)=integer type for each atom in cell
  enxc=exchange and correlation energy (hartree)
  gsqcut=Fourier cutoff on G^2 for "large sphere" of radius double that of the basis sphere
  kxc(nfft,nkxc)=exchange and correlation kernel
  mgfft=maximum size of 1D FFTs
  mpi_enreg=information about MPI parallelization
  ngfft(18)=contain all needed information about 3D FFT,
     see ~abinit/doc/variables/vargs.htm#ngfft
  ngfftf(18)= -PAW ONLY- contain all needed information about 3D FFT for the fine grid
              (ngfftf=ngfft for norm-conserving potential runs)
  nkxc=2nd dimension of kxc
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  n3xccc=dimension of xccc3d (0 if no core charge, nfft otherwise)
  nhat(nfft,nspden*nhatdim)= -PAW only- compensation density
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  ph1d(2,3*(2*mgfft+1)*natom)=1-dim phase (structure factor) information
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rhor(nfft,nspden)=electron density in r space
   (if spin polarized, array contains total density in first half and
    spin-up density in second half)
   (for non-collinear magnetism, first element: total density,
    3 next ones: mx,my,mz)
  rprimd(3,3)=dimensional primitive translation vectors (bohr)
  usexcnhat= -PAW only- 1 if nhat density has to be taken into account in Vxc
  vxc(nfft,nspden)=xc potential (spin up in first half and spin down in
   second half if nspden=2)
  xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
  xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives,
   for each type of atom, from psp
  xccc3d(n3xccc)=3D core electron density for XC core correction, bohr^-3
  xred(3,natom)=reduced coordinates for atoms in unit cell

OUTPUT

  eltfrxc(6+3*natom,6) = xc frozen wavefunction contribution to the
   elastic tensor

SIDE EFFECTS

NOTES

      Much of the code in versions of this routine prior to 4.4.5
      has been transfered to its child eltxccore.

PARENTS

      respfn

CHILDREN

      atm2fft,dfpt_atm2fft,dfpt_mkcore,dfpt_mkvxcstr,dotprod_vn,eltxccore
      fourdp,metric,paw_spline,pawpsp_cg,pawrad_free,pawrad_init,pawtab_free
      pawtab_nullify,redgr,timab,xmpi_sum

SOURCE

138 subroutine dfpt_eltfrxc(atindx,dtset,eltfrxc,enxc,gsqcut,kxc,mpi_enreg,mgfft,&
139 & nattyp,nfft,ngfft,ngfftf,nhat,nkxc,n3xccc,pawtab,ph1d,psps,rhor,rprimd,&
140 & usexcnhat,vxc,xccc3d,xred)
141 
142 
143 !This section has been created automatically by the script Abilint (TD).
144 !Do not modify the following lines by hand.
145 #undef ABI_FUNC
146 #define ABI_FUNC 'dfpt_eltfrxc'
147 !End of the abilint section
148 
149  implicit none
150 
151 !Arguments ------------------------------------
152 !type
153 !scalars
154  integer,intent(in) :: mgfft,n3xccc,nfft,nkxc,usexcnhat
155  real(dp),intent(in) :: enxc,gsqcut
156  type(MPI_type),intent(in) :: mpi_enreg
157  type(dataset_type),intent(in) :: dtset
158  type(pseudopotential_type),intent(inout) :: psps
159 !arrays
160  integer,intent(in) :: atindx(dtset%natom),nattyp(dtset%ntypat),ngfft(18)
161  integer,intent(in) :: ngfftf(18)
162  real(dp),intent(in) :: nhat(nfft,dtset%nspden*psps%usepaw)
163  real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom),rprimd(3,3)
164  real(dp),intent(in) :: vxc(nfft,dtset%nspden),xccc3d(n3xccc)
165  real(dp),intent(in) :: xred(3,dtset%natom)
166  real(dp),intent(in),target :: rhor(nfft,dtset%nspden)
167  real(dp),intent(inout) :: kxc(nfft,nkxc)
168  real(dp),intent(out) :: eltfrxc(6+3*dtset%natom,6)
169  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
170 
171 !Local variables-------------------------------
172 !scalars
173  integer,parameter :: mshift=401
174  integer :: cplex,fgga,ia,idir,ielt,ieltx,ierr,ifft,ii,ipert,is1,is2,ispden,ispden_c,jj,ka,kb
175  integer :: kd,kg,n1,n1xccc,n2,n3,n3xccc_loc,optatm,optdyfr,opteltfr,optgr
176  integer :: option,optn,optn2,optstr,optv
177  real(dp) :: d2eacc,d2ecdgs2,d2exdgs2,d2gsds1ds2,d2gstds1ds2,decdgs,dexdgs
178  real(dp) :: dgsds10,dgsds20,dgstds10,dgstds20,rstep,spnorm,tmp0,tmp0t
179  real(dp) :: ucvol,valuei,yp1,ypn
180  type(pawrad_type) :: core_mesh
181 !arrays
182  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
183  real(dp) :: corstr(6),dummy6(0),dummy_in(0,0)
184  real(dp) :: dummy_out1(0),dummy_out2(0),dummy_out3(0),dummy_out4(0),dummy_out5(0),dummy_out6(0)
185  real(dp) :: eltfrxc_test1(6+3*dtset%natom,6),eltfrxc_test2(6+3*dtset%natom,6)
186  real(dp) :: gmet(3,3),gprimd(3,3),qphon(3),rmet(3,3),tsec(2)
187  real(dp) :: strn_dummy6(0), strv_dummy6(0)
188  real(dp),allocatable :: d2gm(:,:,:,:),dgm(:,:,:),eltfrxc_tmp(:,:)
189  real(dp),allocatable :: eltfrxc_tmp2(:,:),elt_work(:,:),rho0_redgr(:,:,:)
190  real(dp),allocatable :: vxc10(:,:),vxc10_core(:),vxc10_coreg(:,:)
191  real(dp),allocatable :: vxc1is_core(:),vxc1is_coreg(:,:),vxc_core(:)
192  real(dp),allocatable :: vxc_coreg(:,:),work(:),workgr(:,:),xccc1d(:,:,:)
193  real(dp),allocatable :: xccc3d1(:),xccc3d1_temp(:,:),xcccrc(:)
194  real(dp),pointer :: rhor_(:,:)
195  type(pawtab_type),allocatable :: pawtab_test(:)
196 
197 ! *************************************************************************
198 
199 !Initialize variables
200  cplex=1
201  qphon(:)=zero
202  n1=ngfft(1)
203  n2=ngfft(2)
204  n3=ngfft(3)
205 
206  n1xccc = psps%n1xccc
207  if(psps%usepaw==0)then
208    ABI_ALLOCATE(xcccrc,(dtset%ntypat))
209    ABI_ALLOCATE(xccc1d,(n1xccc,6,dtset%ntypat))
210    xcccrc = psps%xcccrc
211    xccc1d = psps%xccc1d
212  end if
213 
214  if (usexcnhat==0.and.dtset%usepaw==1) then
215    ABI_ALLOCATE(rhor_,(nfft,dtset%nspden))
216    rhor_(:,:) = rhor(:,:)-nhat(:,:)
217  else
218    rhor_ => rhor
219  end if
220 
221 !HACK - should be fixed globally
222  if(n1xccc==0) then
223    n3xccc_loc=0
224  else
225    n3xccc_loc=n3xccc
226  end if
227 
228  fgga=0 ; if(nkxc==7.or.nkxc==19) fgga=1
229 
230  ABI_ALLOCATE(eltfrxc_tmp,(6+3*dtset%natom,6))
231  ABI_ALLOCATE(eltfrxc_tmp2,(6+3*dtset%natom,6))
232  ABI_ALLOCATE(vxc10,(nfft,dtset%nspden))
233  ABI_ALLOCATE(xccc3d1,(cplex*nfft))
234 
235  if(n1xccc/=0) then
236    ABI_ALLOCATE(vxc_core,(nfft))
237    ABI_ALLOCATE(vxc10_core,(nfft))
238    ABI_ALLOCATE(vxc1is_core,(nfft))
239 
240    if(dtset%nspden==1) then
241      vxc_core(:)=vxc(:,1)
242    else
243      vxc_core(:)=0.5_dp*(vxc(:,1)+vxc(:,2))
244    end if
245  end if
246 
247 !Compute gmet, gprimd and ucvol from rprimd
248  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
249 
250 !For GGA case, prepare quantities needed to evaluate contributions
251 !arising from the strain dependence of the gradient operator itself
252 
253  if(fgga==1) then
254    ABI_ALLOCATE(rho0_redgr,(3,nfft,dtset%nspden))
255    ABI_ALLOCATE(work,(nfft))
256    ABI_ALLOCATE(workgr,(nfft,3))
257 
258 !  Set up metric tensor derivatives
259    ABI_ALLOCATE(dgm,(3,3,6))
260    ABI_ALLOCATE(d2gm,(3,3,6,6))
261 !  Loop over 2nd strain index
262    do is2=1,6
263      kg=idx(2*is2-1);kd=idx(2*is2)
264      do jj = 1,3
265        dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj))
266      end do
267 
268 !    Loop over 1st strain index
269      do is1=1,6
270        ka=idx(2*is1-1);kb=idx(2*is1)
271        d2gm(:,:,is1,is2)=0._dp
272        do jj = 1,3
273          if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
274 &         +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj)
275          if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
276 &         +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj)
277          if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
278 &         +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj)
279          if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
280 &         +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj)
281        end do
282        d2gm(:,:,is1,is2)=0.5_dp*d2gm(:,:,is1,is2)
283      end do
284    end do
285 
286 !  Compute the reduced gradients of the zero-order charge density.
287 !  Note that in the spin-polarized case, we are computing the reduced
288 !  gradients of 2 X the spin-up or spin-down charge.  This simplifies
289 !  subsequent code for the non-spin-polarized case.
290    if(dtset%nspden==1) then
291      work(:)=rhor_(:,1)
292    else
293      work(:)=2.0_dp*rhor_(:,2)
294    end if
295    if(n1xccc/=0) then
296      work(:)=work(:)+xccc3d(:)
297    end if
298    call redgr (work,workgr,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb)
299    do ifft=1,nfft
300      rho0_redgr(:,ifft,1)=workgr(ifft,:)
301    end do
302    if(dtset%nspden==2) then
303      work(:)=2.0_dp*(rhor_(:,1)-rhor_(:,2))
304      if(n1xccc/=0) then
305        work(:)=work(:)+xccc3d(:)
306      end if
307      call redgr(work,workgr,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb)
308      do ifft=1,nfft
309        rho0_redgr(:,ifft,2)=workgr(ifft,:)
310      end do
311    end if
312    ABI_DEALLOCATE(work)
313    ABI_DEALLOCATE(workgr)
314  end if !GGA
315 
316 
317 !Null the elastic tensor accumulator
318  eltfrxc(:,:)=zero;eltfrxc_tmp(:,:)=zero;eltfrxc_tmp2(:,:) = zero
319 
320 !Normalization factor
321  if(dtset%nspden==1) then
322    spnorm=one
323  else
324    spnorm=half
325  end if
326 
327 !Big loop over 2nd strain index
328  do is2=1,6
329 
330 !  Translate strain index as needed by dfpt_mkcore below.
331    if(is2<=3) then
332      ipert=dtset%natom+3
333      idir=is2
334    else
335      ipert=dtset%natom+4
336      idir=is2-3
337    end if
338 
339 !  Generate first-order core charge for is2 strain if core charges are present.
340    if(n1xccc/=0)then
341 
342      if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
343 !      Calculation in Reciprocal space for paw or NC with nc_xccc_gspace
344        ABI_ALLOCATE(xccc3d1_temp,(cplex*nfft,1))
345        xccc3d1_temp = zero
346 
347        call dfpt_atm2fft(atindx,cplex,gmet,gprimd,gsqcut,is2,ipert,&
348 &       mgfft,psps%mqgrid_vl,dtset%natom,1,nfft,ngfftf,dtset%ntypat,&
349 &       ph1d,psps%qgrid_vl,qphon,dtset%typat,ucvol,psps%usepaw,xred,psps,pawtab,&
350 &       atmrhor1=xccc3d1_temp,optn2_in=1,&
351 &       comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
352 &       paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
353        xccc3d1(:) = xccc3d1_temp(:,1)
354        ABI_DEALLOCATE(xccc3d1_temp)
355 
356      else
357 !      Calculation in direct space for norm conserving:
358        call dfpt_mkcore(cplex,idir,ipert,dtset%natom,dtset%ntypat,n1,n1xccc,&
359 &       n2,n3,qphon,rprimd,dtset%typat,ucvol,&
360 &       xcccrc,xccc1d,xccc3d1,xred)
361      end if
362    else
363      xccc3d1(:)=zero
364    end if
365 
366 !  Compute the first-order potentials.
367 !  Standard first-order potential for LDA and GGA with core charge
368    if(fgga==0 .or. (fgga==1 .and. n1xccc/=0)) then
369      option=0
370      call dfpt_mkvxcstr(cplex,idir,ipert,kxc,mpi_enreg,dtset%natom,nfft,ngfft,nhat,&
371 &     dummy_in,nkxc,dtset%nspden,n3xccc_loc,option,mpi_enreg%paral_kgb,qphon,rhor,rhor,&
372 &     rprimd,dtset%usepaw,usexcnhat,vxc10,xccc3d1)
373      if(n1xccc/=0)then
374        if(dtset%nspden==1) then
375          vxc10_core(:)=vxc10(:,1)
376          vxc1is_core(:)=vxc10(:,1)
377        else
378          vxc10_core(:)=0.5_dp*(vxc10(:,1)+vxc10(:,2))
379          vxc1is_core(:)=0.5_dp*(vxc10(:,1)+vxc10(:,2))
380        end if
381      end if
382    end if
383 
384 !  For GGA, first-order potential with doubled gradient operator strain
385 !  derivative terms needed for elastic tensor but not internal strain.
386    if(fgga==1) then
387      option=2
388      call dfpt_mkvxcstr(cplex,idir,ipert,kxc,mpi_enreg,dtset%natom,nfft,ngfft,nhat,&
389 &     dummy_in,nkxc,dtset%nspden,n3xccc_loc,option,mpi_enreg%paral_kgb,qphon,rhor,rhor,&
390 &     rprimd,dtset%usepaw,usexcnhat,vxc10,xccc3d1)
391      if(n1xccc/=0)then
392        if(dtset%nspden==1) then
393          vxc10_core(:)=vxc10(:,1)
394        else
395          vxc10_core(:)=0.5_dp*(vxc10(:,1)+vxc10(:,2))
396        end if
397      end if
398    end if
399 
400 
401 !  Additional term for diagonal strains.
402    if(is2<=3) then
403      vxc10(:,:)=vxc10(:,:)+vxc(:,:)
404      if(n1xccc/=0) then
405        vxc10_core(:)=vxc10_core(:)+2.0_dp*vxc_core(:)
406        vxc1is_core(:)=vxc1is_core(:)+vxc_core(:)
407      end if
408    end if
409 
410 !  For GGA, compute the contributions from the strain derivatives acting
411 !  on the gradient operators.
412    if(fgga==1) then
413 
414      if (dtset%nspden==1) then
415        do ifft=1,nfft
416 !        Collect the needed derivatives of Exc.  The factors introduced
417 !        deal with the difference between density as used here and
418 !        spin density as used with these kxc terms in other contexts.
419          dexdgs  =half   *kxc(ifft,2)
420          d2exdgs2=quarter*kxc(ifft,4)
421 !        Loop over 1st strain index
422          do is1=1,6
423 !          The notation here is .gs... for the derivatives of the squared-
424 !          gradient of (2X) each spin density, and .gst... for the total density.
425            dgsds10=zero;dgsds20=zero;d2gsds1ds2=zero
426            do jj=1,3
427              do ii=1,3
428                tmp0=rho0_redgr(ii,ifft,1)*rho0_redgr(jj,ifft,1)
429                dgsds10=dgsds10+dgm(ii,jj,is1)*tmp0
430                dgsds20=dgsds20+dgm(ii,jj,is2)*tmp0
431                d2gsds1ds2=d2gsds1ds2+d2gm(ii,jj,is1,is2)*tmp0
432              end do
433            end do
434 !          Volume derivative terms added
435            if(is1<=3) d2gsds1ds2=d2gsds1ds2+dgsds20
436            if(is2<=3) d2gsds1ds2=d2gsds1ds2+dgsds10
437 !          Add the gradient derivative terms to eltfrxc.
438            eltfrxc(is1,is2)=eltfrxc(is1,is2)+d2exdgs2*dgsds10*dgsds20+dexdgs*d2gsds1ds2
439          end do !is1
440        end do !ifft
441 
442      else ! nspden==2
443 
444        do ispden=1,dtset%nspden
445          ispden_c=dtset%nspden-ispden+1
446 
447          do ifft=1,nfft
448 
449 !          Collect the needed derivatives of Exc.  The factors introduced
450 !          deal with the difference between density as used here and
451 !          spin density as used with these kxc terms in other contexts.
452            dexdgs  =quarter       *kxc(ifft,3+ispden)
453            d2exdgs2=quarter*eighth*kxc(ifft,7+ispden)
454            decdgs  =eighth        *kxc(ifft,10)
455            d2ecdgs2=eighth*eighth *kxc(ifft,13)
456 
457 !          Loop over 1st strain index
458            do is1=1,6
459 
460 !            The notation here is .gs... for the derivatives of the squared-
461 !            gradient of (2X) each spin density, and .gst... for the total
462 !            density.  Note the hack that the the total density is given
463 !            by the same expression for either the non-polarized or spin-
464 !            polarized case, implemented with the "complementary" index ispden_c
465 !            in the expression for tmp0t below.
466              dgsds10=zero;dgsds20=zero;d2gsds1ds2=zero
467              dgstds10=zero;dgstds20=zero;d2gstds1ds2=zero
468              do jj=1,3
469                do ii=1,3
470                  tmp0=rho0_redgr(ii,ifft,ispden)*rho0_redgr(jj,ifft,ispden)
471                  tmp0t=(rho0_redgr(ii,ifft,ispden)+rho0_redgr(ii,ifft,ispden_c))&
472 &                 *(rho0_redgr(jj,ifft,ispden)+rho0_redgr(jj,ifft,ispden_c))
473                  dgsds10=dgsds10+dgm(ii,jj,is1)*tmp0
474                  dgsds20=dgsds20+dgm(ii,jj,is2)*tmp0
475                  dgstds10=dgstds10+dgm(ii,jj,is1)*tmp0t
476                  dgstds20=dgstds20+dgm(ii,jj,is2)*tmp0t
477                  d2gsds1ds2=d2gsds1ds2+d2gm(ii,jj,is1,is2)*tmp0
478                  d2gstds1ds2=d2gstds1ds2+d2gm(ii,jj,is1,is2)*tmp0t
479                end do
480              end do
481 !            Volume derivative terms added
482              if(is1<=3) then
483                d2gsds1ds2=d2gsds1ds2+dgsds20
484                d2gstds1ds2=d2gstds1ds2+dgstds20
485              end if
486              if(is2<=3) then
487                d2gsds1ds2=d2gsds1ds2+dgsds10
488                d2gstds1ds2=d2gstds1ds2+dgstds10
489              end if
490 
491 !            Add the gradient derivative terms to eltfrxc.
492              eltfrxc(is1,is2)=eltfrxc(is1,is2)+spnorm*&
493 &             (d2exdgs2*(dgsds10 *dgsds20) + dexdgs*d2gsds1ds2&
494 &             +d2ecdgs2*(dgstds10*dgstds20)+ decdgs*d2gstds1ds2)
495 
496            end do !is1
497          end do !ifft
498        end do !ispden
499 
500      end if ! nspden
501 
502    end if !GGA
503 
504 !  Compute valence electron 1st-order charge contributions.  Recall that
505 !  the diagonal strain derivatives of the valence charge are minus the
506 !  zero-order density.  The explicit symmetrization avoids the need
507 !  to store vxc10 for strain indices other than is2.
508 
509    call dotprod_vn(1,rhor_,d2eacc,valuei,nfft,nfft,dtset%nspden,1,&
510 &   vxc10,ucvol)
511    do is1=1,3
512      eltfrxc_tmp(is1,is2)=eltfrxc_tmp(is1,is2)-0.5_dp*d2eacc
513      eltfrxc_tmp(is2,is1)=eltfrxc_tmp(is2,is1)-0.5_dp*d2eacc
514    end do
515 
516 !  Compute additional core contributions from is1 perturbation
517 !  Internal strain terms calculated here.
518    if(n1xccc/=0) then
519 
520      if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then
521 !      Calculation in Reciprocal space for paw or NC with nc_xccc_gspace
522        optatm=0;optdyfr=0;optgr=0;optstr=0;optv=0;optn=n3xccc/nfft;optn2=1;opteltfr=1
523        ABI_ALLOCATE(vxc10_coreg,(2,nfft))
524        ABI_ALLOCATE(vxc_coreg,(2,nfft))
525        ABI_ALLOCATE(vxc1is_coreg,(2,nfft))
526 
527        vxc10_coreg(:,:)=zero;vxc10_coreg(:,:)=zero;vxc1is_coreg(:,:)=zero;
528 
529 !      Fourier transform of Vxc_core/vxc10_core to use in atm2fft (reciprocal space calculation)
530        call fourdp(1,vxc10_coreg,vxc10_core,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0)
531        call fourdp(1,vxc_coreg,vxc_core,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0)
532        call fourdp(1,vxc1is_coreg,vxc1is_core,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0)
533 
534        call atm2fft(atindx,dummy_out1,dummy_out2,dummy_out3,dummy_out4,eltfrxc_tmp2,dummy_in,gmet,gprimd,&
535 &       dummy_out5,dummy_out6,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%ntypat,&
536 &       optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1d,psps%qgrid_vl,dtset%qprtrb,&
537 &       dummy_in,strn_dummy6,strv_dummy6,ucvol,psps%usepaw,vxc_coreg,vxc10_coreg,vxc1is_coreg,dtset%vprtrb,psps%vlspl,is2_in=is2,&
538 &       comm_fft=mpi_enreg%comm_fft,me_g0=mpi_enreg%me_g0,&
539 &       paral_kgb=mpi_enreg%paral_kgb,distribfft=mpi_enreg%distribfft)
540 
541 !       The indexing array atindx is used to reestablish the correct order of atoms
542        ABI_ALLOCATE(elt_work,(6+3*dtset%natom,6))
543        elt_work(1:6,1:6)=eltfrxc_tmp2(1:6,1:6)
544        do ia=1,dtset%natom
545          ielt=7+3*(ia-1)
546          ieltx=7+3*(atindx(ia)-1)
547          elt_work(ielt:ielt+2,1:6)=eltfrxc_tmp2(ieltx:ieltx+2,1:6)
548        end do
549        eltfrxc_tmp2(:,:)=elt_work(:,:)
550        ABI_DEALLOCATE(elt_work)
551 
552 
553        ABI_DEALLOCATE(vxc10_coreg)
554        ABI_DEALLOCATE(vxc_coreg)
555        ABI_DEALLOCATE(vxc1is_coreg)
556        eltfrxc(:,:)= eltfrxc(:,:) + eltfrxc_tmp2(:,:)
557 
558      else
559 
560        call eltxccore(eltfrxc,is2,mpi_enreg%my_natom,dtset%natom,nfft,dtset%ntypat,&
561 &       n1,n1xccc,n2,n3,rprimd,dtset%typat,ucvol,vxc_core,vxc10_core,vxc1is_core,&
562 &       xcccrc,xccc1d,xred,mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
563 
564 !DEBUG
565 !      TEST ZONE (DO NOT REMOVE) USE TO RECIPROCAL SPACE IN NC CASE
566        if (dtset%userid==567) then
567          eltfrxc_test1(:,is2)=zero;eltfrxc_test2(:,is2)=zero
568          call eltxccore(eltfrxc_test1,is2,mpi_enreg%my_natom,dtset%natom,nfft,dtset%ntypat,&
569 &         n1,n1xccc,n2,n3,rprimd,dtset%typat,ucvol,vxc_core,vxc10_core,vxc1is_core,&
570 &         xcccrc,xccc1d,xred,mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
571 !        if (is2==1) print*,"elt-frxc from eltxccore",is2,eltfrxc_test1(1,1)*ucvol/dble(nfft)
572          ABI_DATATYPE_ALLOCATE(pawtab_test,(dtset%ntypat))
573          call pawtab_nullify(pawtab_test)
574          do jj=1,dtset%ntypat
575            pawtab_test(jj)%mqgrid=psps%mqgrid_vl
576            ABI_ALLOCATE(pawtab_test(jj)%tcorespl,(pawtab_test(jj)%mqgrid,2))
577            rstep=xcccrc(jj)/dble(n1xccc-1)
578            call pawrad_init(mesh=core_mesh,mesh_size=n1xccc,mesh_type=1,rstep=rstep)
579            call pawpsp_cg(pawtab_test(jj)%dncdq0,pawtab_test(jj)%d2ncdq0,psps%mqgrid_vl,psps%qgrid_vl,&
580 &           pawtab_test(jj)%tcorespl(:,1),core_mesh,xccc1d(:,1,jj),yp1,ypn)
581            call paw_spline(psps%qgrid_vl,pawtab_test(jj)%tcorespl(:,1),psps%mqgrid_vl,yp1,ypn,pawtab_test(jj)%tcorespl(:,2))
582 !          if (is2==1) then
583 !            do ii=1,n1xccc;write(100+jj,*) (ii-1)*rstep,xccc1d(ii,1,jj);enddo
584 !            do ii=1,psps%mqgrid_vl;write(200+jj,*) psps%qgrid_vl(ii),pawtab_test(jj)%tcorespl(ii,1);enddo
585 !          end if
586          end do
587          ABI_ALLOCATE(vxc10_coreg,(2,nfft))
588          ABI_ALLOCATE(vxc_coreg,(2,nfft))
589          ABI_ALLOCATE(vxc1is_coreg,(2,nfft))
590          vxc10_coreg(:,:)=zero;vxc10_coreg(:,:)=zero;vxc1is_coreg(:,:)=zero;
591          call fourdp(1,vxc10_coreg,vxc10_core,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0)
592          call fourdp(1,vxc_coreg,vxc_core,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0)
593          call fourdp(1,vxc1is_coreg,vxc1is_core,-1,mpi_enreg,nfft,ngfft,mpi_enreg%paral_kgb,0)
594          optatm=0;optdyfr=0;optgr=0;optstr=0;optv=0;optn=1;optn2=1;opteltfr=1;corstr=zero
595          call atm2fft(atindx,dummy_out1,dummy_out2,dummy_out3,dummy_out4,eltfrxc_test2,dummy_in,gmet,gprimd,&
596 &         dummy_out5,dummy_out6,gsqcut,mgfft,psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%ntypat,&
597 &         optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab_test,ph1d,psps%qgrid_vl,dtset%qprtrb,&
598 &         dummy_in,corstr,dummy6,ucvol,psps%usepaw,vxc_coreg,vxc10_coreg,vxc1is_coreg,dtset%vprtrb,psps%vlspl,is2_in=is2)
599          ABI_DEALLOCATE(vxc10_coreg)
600          ABI_DEALLOCATE(vxc_coreg)
601          ABI_DEALLOCATE(vxc1is_coreg)
602          call pawrad_free(core_mesh)
603          call pawtab_free(pawtab_test)
604          ABI_DATATYPE_DEALLOCATE(pawtab_test)
605          eltfrxc(:,:)= eltfrxc(:,:)+eltfrxc_test2(:,:)
606 !        if (is2==1) print*,"cor-str from atm2fft",is2,corstr*ucvol
607 !        if (is2==1) print*,"elt-frxc from atm2fft  ",is2,eltfrxc_test2(1,1)
608        end if
609 !DEBUG
610 
611      end if
612    end if
613 
614 !  Additional term for diagonal strains
615    if(is2<=3) then
616      do is1=1,3
617        eltfrxc_tmp(is1,is2)=eltfrxc_tmp(is1,is2)+enxc
618      end do
619    end if
620  end do !is2 outermost strain loop
621 
622 !Accumulate eltfrxc accross processors
623  call timab(48,1,tsec)
624  call xmpi_sum(eltfrxc,mpi_enreg%comm_fft,ierr)
625  call timab(48,2,tsec)
626 
627  !Normalize accumulated 2nd derivatives in NC case
628  if(psps%usepaw==1)then
629    eltfrxc(:,:)=eltfrxc_tmp(:,:)+eltfrxc
630  else
631    eltfrxc(:,:)=eltfrxc_tmp(:,:)+eltfrxc*ucvol/dble(nfft)
632  end if
633 
634  ABI_DEALLOCATE(eltfrxc_tmp)
635  ABI_DEALLOCATE(eltfrxc_tmp2)
636  ABI_DEALLOCATE(vxc10)
637  ABI_DEALLOCATE(xccc3d1)
638  if(psps%usepaw==0)then
639    ABI_DEALLOCATE(xccc1d)
640    ABI_DEALLOCATE(xcccrc)
641  end if
642  if (usexcnhat==0.and.dtset%usepaw==1) then
643    ABI_DEALLOCATE(rhor_)
644  end if
645 
646  if(n1xccc/=0) then
647    ABI_DEALLOCATE(vxc_core)
648    ABI_DEALLOCATE(vxc10_core)
649    ABI_DEALLOCATE(vxc1is_core)
650  end if
651 
652  if(fgga==1) then
653    ABI_DEALLOCATE(rho0_redgr)
654    ABI_DEALLOCATE(dgm)
655    ABI_DEALLOCATE(d2gm)
656  end if
657 
658 end subroutine dfpt_eltfrxc

ABINIT/dfpt_ewald [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_ewald

FUNCTION

 Compute ewald contribution to the dynamical matrix, at a given q wavevector.
 Note: the q=0 part should be subtracted, by another call to
 the present routine, with q=0. The present routine correspond
 to the quantity C_bar defined in Eq.(24) or (27) in Phys. Rev. B 55, 10355 (1997) [[cite:Gonze1997a]].
 The two calls correspond to Eq.(23) of the same paper.
 If q=0 is asked, sumg0 should be put to 0. Otherwise, it should be put to 1.

INPUTS

 gmet(3,3)=metric tensor in reciprocal space (length units **-2)
 mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
 comm_atom=--optional-- MPI communicator over atoms
 my_natom=number of atoms treated by current processor
 natom=number of atoms in unit cell
 qphon(3)=phonon wavevector (same system of coordinates as the
          reciprocal lattice vectors)
 rmet(3,3)=metric tensor in real space (length units squared)
 sumg0: if=1, the sum in reciprocal space must include g=0,
   if=0, this contribution must be skipped (q=0 singularity)
 typat(natom)=integer label of each type of atom (1,2,...)
 ucvol=unit cell volume in (whatever length scale units)**3
 xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
 zion(ntypat)=charge on each type of atom (real number)

OUTPUT

 dyew(2,3,natom,3,natom)= Ewald part of the dynamical matrix,
    second energy derivative wrt xred(3,natom), Hartrees.

PARENTS

      respfn

CHILDREN

      free_my_atmtab,get_my_atmtab,timab,xmpi_sum

SOURCE

2516 subroutine dfpt_ewald(dyew,gmet,my_natom,natom,qphon,rmet,sumg0,typat,ucvol,xred,zion, &
2517 &                 mpi_atmtab,comm_atom ) ! optional arguments (parallelism))
2518 
2519 
2520 !This section has been created automatically by the script Abilint (TD).
2521 !Do not modify the following lines by hand.
2522 #undef ABI_FUNC
2523 #define ABI_FUNC 'dfpt_ewald'
2524 !End of the abilint section
2525 
2526  implicit none
2527 
2528 !Arguments -------------------------------
2529 !scalars
2530  integer,intent(in) :: my_natom,natom,sumg0
2531  real(dp),intent(in) :: ucvol
2532 !arrays
2533  integer,intent(in) :: typat(natom)
2534  integer,optional,intent(in) :: comm_atom
2535  integer,optional,target,intent(in) :: mpi_atmtab(:)
2536  real(dp),intent(in) :: gmet(3,3),qphon(3),rmet(3,3),xred(3,natom),zion(*)
2537  real(dp),intent(out) :: dyew(2,3,natom,3,natom)
2538 
2539 !Local variables -------------------------
2540 !nr, ng affect convergence of sums (nr=3,ng=5 is not good enough):
2541 !scalars
2542  integer,parameter :: im=2,ng=10,nr=6,re=1
2543  integer :: ia,ia0,ib,ierr,ig1,ig2,ig3,ii,ir1,ir2,ir3,mu,my_comm_atom,nu
2544  logical :: my_atmtab_allocated,paral_atom
2545  real(dp) :: arg,arga,argb,c1i,c1r,da1,da2,da3,derfc_arg
2546  real(dp) :: direct,dot1,dot2,dot3,dotr1,dotr2,dotr3
2547  real(dp) :: eta,fac,gdot12,gdot13,gdot23,gsq,gsum,norm1
2548  real(dp) :: r1,r2,r3,rdot12,rdot13,rdot23,recip,reta
2549  real(dp) :: reta3m,rmagn,rsq,term,term1,term2
2550  real(dp) :: term3
2551  character(len=500) :: message
2552 !arrays
2553  real(dp) :: tsec(2)
2554  integer,pointer :: my_atmtab(:)
2555  real(dp) :: gpq(3),rq(3)
2556 
2557 ! *************************************************************************
2558 
2559 !Set up parallelism over atoms
2560  paral_atom=(present(comm_atom).and.(my_natom/=natom))
2561  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2562  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2563  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
2564 
2565 !Compute eta for approximately optimized summations:
2566  direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
2567 & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
2568  recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
2569 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
2570  eta=pi*(dble(ng)/dble(nr))*sqrt(1.69_dp*recip/direct)
2571 
2572 !Test Ewald s summation
2573 !eta=1.2_dp*eta
2574 
2575 !Sum terms over g space:
2576  fac=pi**2/eta
2577  gsum=zero
2578  da1=zero
2579  da2=zero
2580  da3=zero
2581  dyew(:,:,:,:,:)=zero
2582  do ig3=-ng,ng
2583    do ig2=-ng,ng
2584      do ig1=-ng,ng
2585        gpq(1)=dble(ig1)+qphon(1)
2586        gpq(2)=dble(ig2)+qphon(2)
2587        gpq(3)=dble(ig3)+qphon(3)
2588        gdot12=gmet(2,1)*gpq(1)*gpq(2)
2589        gdot13=gmet(3,1)*gpq(1)*gpq(3)
2590        gdot23=gmet(3,2)*gpq(2)*gpq(3)
2591        dot1=gmet(1,1)*gpq(1)**2+gdot12+gdot13
2592        dot2=gmet(2,2)*gpq(2)**2+gdot12+gdot23
2593        dot3=gmet(3,3)*gpq(3)**2+gdot13+gdot23
2594        gsq=dot1+dot2+dot3
2595 !      Skip q=0:
2596        if (gsq<1.0d-20) then
2597          if (sumg0==1) then
2598            write(message,'(5a)')&
2599 &           'The phonon wavelength should not be zero : ',ch10,&
2600 &           'there are non-analytical terms that the code cannot handle.',ch10,&
2601 &           'Action : subtract this wavelength from the input.'
2602            MSG_ERROR(message)
2603          end if
2604        else
2605          arg=fac*gsq
2606 !        Larger arg gives 0 contribution:
2607          if (arg <= 80._dp) then
2608            term=exp(-arg)/gsq
2609            do ia0=1,my_natom
2610              ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2611              arga=two_pi*(gpq(1)*xred(1,ia)+gpq(2)*xred(2,ia)+gpq(3)*xred(3,ia))
2612              do ib=1,ia
2613                argb=two_pi*(gpq(1)*xred(1,ib)+gpq(2)*xred(2,ib)+gpq(3)*xred(3,ib))
2614                arg=arga-argb
2615                c1r=cos(arg)*term
2616                c1i=sin(arg)*term
2617 
2618                do mu=1,3
2619                  do nu=1,mu
2620                    dyew(re,mu,ia,nu,ib)=dyew(re,mu,ia,nu,ib)+gpq(mu)*gpq(nu)*c1r
2621                    dyew(im,mu,ia,nu,ib)=dyew(im,mu,ia,nu,ib)+gpq(mu)*gpq(nu)*c1i
2622                  end do
2623                end do
2624 
2625              end do
2626            end do
2627          end if
2628 !        Endif g/=0 :
2629        end if
2630 !      End triple loop over G s:
2631      end do
2632    end do
2633  end do
2634 
2635 !End G summation by accounting for some common factors.
2636 !(for the charges:see end of routine)
2637  norm1=4.0_dp*pi/ucvol
2638  do ia0=1,my_natom
2639    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2640    do ib=1,ia
2641      do mu=1,3
2642        do nu=1,mu
2643          dyew(:,mu,ia,nu,ib)=dyew(:,mu,ia,nu,ib)*norm1
2644        end do
2645      end do
2646    end do
2647  end do
2648 
2649 
2650 !Do sums over real space:
2651  reta=sqrt(eta)
2652  reta3m=-eta*reta
2653  fac=4._dp/3.0_dp/sqrt(pi)
2654  do ir3=-nr,nr
2655    do ir2=-nr,nr
2656      do ir1=-nr,nr
2657        arg=-two_pi*(qphon(1)*ir1+qphon(2)*ir2+qphon(3)*ir3)
2658        c1r=cos(arg)*reta3m
2659        c1i=sin(arg)*reta3m
2660        do ia0=1,my_natom
2661          ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2662          do ib=1,ia
2663            r1=dble(ir1)+xred(1,ia)-xred(1,ib)
2664            r2=dble(ir2)+xred(2,ia)-xred(2,ib)
2665            r3=dble(ir3)+xred(3,ia)-xred(3,ib)
2666            rdot12=rmet(2,1)*r1*r2
2667            rdot13=rmet(3,1)*r1*r3
2668            rdot23=rmet(3,2)*r2*r3
2669            dotr1=rmet(1,1)*r1**2+rdot12+rdot13
2670            dotr2=rmet(2,2)*r2**2+rdot12+rdot23
2671            dotr3=rmet(3,3)*r3**2+rdot13+rdot23
2672            rsq=dotr1+dotr2+dotr3
2673            rmagn=sqrt(rsq)
2674 !          Avoid zero denominators in term :
2675            if (rmagn>=1.0d-12) then
2676              arg=reta*rmagn
2677              term=zero
2678              if (arg<8.0_dp) then
2679 !              Note: erfc(8) is about 1.1e-29,
2680 !              so don t bother with larger arg.
2681 !              Also: exp(-64) is about 1.6e-28,
2682 !              so don t bother with larger arg**2 in exp.
2683                derfc_arg = abi_derfc(arg)
2684                term=derfc_arg/arg**3
2685                term1=2.0_dp/sqrt(pi)*exp(-arg**2)/arg**2
2686                term2=-(term+term1)
2687                term3=(3*term+term1*(3.0_dp+2.0_dp*arg**2))/rsq
2688                rq(1)=rmet(1,1)*r1+rmet(1,2)*r2+rmet(1,3)*r3
2689                rq(2)=rmet(2,1)*r1+rmet(2,2)*r2+rmet(2,3)*r3
2690                rq(3)=rmet(3,1)*r1+rmet(3,2)*r2+rmet(3,3)*r3
2691                do mu=1,3
2692                  do nu=1,mu
2693                    dyew(re,mu,ia,nu,ib)=dyew(re,mu,ia,nu,ib)+&
2694 &                   c1r*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
2695                    dyew(im,mu,ia,nu,ib)=dyew(im,mu,ia,nu,ib)+&
2696 &                   c1i*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
2697                  end do
2698                end do
2699              end if
2700            else
2701              if (ia/=ib)then
2702                write(message,'(a,a,a,a,a,i5,a,i5,a)')&
2703 &               'The distance between two atoms vanishes.',ch10,&
2704 &               'This is not allowed.',ch10,&
2705 &               'Action: check the input for the atoms number',ia,' and',ib,'.'
2706                MSG_ERROR(message)
2707              else
2708                do mu=1,3
2709                  do nu=1,mu
2710                    dyew(re,mu,ia,nu,ib)=dyew(re,mu,ia,nu,ib)+&
2711 &                   fac*reta3m*rmet(mu,nu)
2712                  end do
2713                end do
2714              end if
2715            end if
2716 
2717          end do ! End loop over ib:
2718        end do ! End loop over ia:
2719      end do ! End triple loop over real space points:
2720    end do
2721  end do
2722 
2723 !Take account of the charges
2724 !write(std_out,*)' '
2725  do ia0=1,my_natom
2726    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2727    do ib=1,ia
2728      do mu=1,3
2729        do nu=1,mu
2730          do ii=1,2
2731 !          write(std_out,*)dyew(ii,mu,ia,nu,ib)
2732            dyew(ii,mu,ia,nu,ib)=dyew(ii,mu,ia,nu,ib)*&
2733 &           zion(typat(ia))*zion(typat(ib))
2734          end do
2735        end do
2736      end do
2737    end do
2738  end do
2739 
2740 !Symmetrize with respect to the directions
2741  do ia0=1,my_natom
2742    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2743    do ib=1,ia
2744      do mu=1,3
2745        do nu=1,mu
2746          dyew(re,nu,ia,mu,ib)=dyew(re,mu,ia,nu,ib)
2747          dyew(im,nu,ia,mu,ib)=dyew(im,mu,ia,nu,ib)
2748        end do
2749      end do
2750    end do
2751  end do
2752 
2753 !In case of parallelism over atoms: communicate
2754  if (paral_atom) then
2755    call timab(48,1,tsec)
2756    call xmpi_sum(dyew,my_comm_atom,ierr)
2757    call timab(48,2,tsec)
2758  end if
2759 
2760 !Fill the upper part of the matrix, with the hermitian conjugate
2761  do ia=1,natom
2762    do ib=1,ia
2763      do nu=1,3
2764        do mu=1,3
2765          dyew(re,mu,ib,nu,ia)=dyew(re,mu,ia,nu,ib)
2766          dyew(im,mu,ib,nu,ia)=-dyew(im,mu,ia,nu,ib)
2767        end do
2768      end do
2769    end do
2770  end do
2771 
2772 !Destroy atom table used for parallelism
2773  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2774 
2775 end subroutine dfpt_ewald

ABINIT/elt_ewald [ Functions ]

[ Top ] [ Functions ]

NAME

 elt_ewald

FUNCTION

 Compute 2nd derivatives of Ewald energy wrt strain for frozen wavefunction
 contributions to elastic tensor

INPUTS

 gmet(3,3)=metric tensor in reciprocal space (bohr^-2)
 gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr^-1)
 mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
 comm_atom=--optional-- MPI communicator over atoms
 my_natom=number of atoms treated by current processor
 natom=number of atoms in unit cell
 ntypat=numbe of type of atoms
 rmet(3,3)=metric tensor in real space (bohr^2)
 rprimd(3,3)=dimensional primitive translation vectors (bohr)
 typat(natom)=integer label of each type of atom (1,2,...)
 ucvol=unit cell volume (bohr^3)
 xred(3,natom)=relative coords of atoms in unit cell (dimensionless)
 zion(ntypat)=charge on each type of atom (real number)

OUTPUT

 elteew(6+3*natom,6)=2nd derivatives of Ewald energy wrt strain

PARENTS

      respfn

CHILDREN

      free_my_atmtab,get_my_atmtab,timab,wrtout,xmpi_sum

SOURCE

2087 subroutine elt_ewald(elteew,gmet,gprimd,my_natom,natom,ntypat,rmet,rprimd,&
2088 &                 typat,ucvol,xred,zion, &
2089 &                 mpi_atmtab,comm_atom) ! optional arguments (parallelism)
2090 
2091 
2092 !This section has been created automatically by the script Abilint (TD).
2093 !Do not modify the following lines by hand.
2094 #undef ABI_FUNC
2095 #define ABI_FUNC 'elt_ewald'
2096 !End of the abilint section
2097 
2098  implicit none
2099 
2100 !Arguments ------------------------------------
2101 !scalars
2102  integer,intent(in) :: my_natom,natom,ntypat
2103  real(dp),intent(in) :: ucvol
2104 !arrays
2105  integer,intent(in) :: typat(natom)
2106  integer,optional,intent(in) :: comm_atom
2107  integer,optional,target,intent(in) :: mpi_atmtab(:)
2108  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
2109  real(dp),intent(in) :: xred(3,natom),zion(ntypat)
2110  real(dp),intent(out) :: elteew(6+3*natom,6)
2111 
2112 !Local variables-------------------------------
2113 !scalars
2114  integer :: ia,ia0,ib,ierr,ig1,ig2,ig3,ir1,ir2,ir3,is1,is2,jj,js,ka,kb,kd,kg,my_comm_atom,newg,newr,ng,nr
2115  logical :: my_atmtab_allocated,paral_atom
2116  real(dp) :: arg,ch,chsq,cos_term,d2derfc,d2gss,d2r,d2rs,dderfc,derfc_arg
2117  real(dp) :: dgss1,dgss2,direct,dr1,dr2,drs1,drs2,eew,eta,fac,fraca1,fraca2
2118  real(dp) :: fraca3,fracb1,fracb2,fracb3,gsq,gsum,r1,r2,r3,recip,reta
2119  real(dp) :: rmagn,rsq,sin_term,sumg,summi,summr,sumr,t1,term
2120  character(len=500) :: message
2121 !arrays
2122  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
2123  integer,pointer :: my_atmtab(:)
2124  real(dp) :: d2gm(3,3,6,6),d2ris(3),d2rm(3,3,6,6),dgm(3,3,6),dris(3),drm(3,3,6)
2125  real(dp) :: t2(3),ts2(3),tsec(2),tt(3)
2126  real(dp),allocatable :: d2sumg(:,:),d2sumr(:,:),drhoisi(:,:),drhoisr(:,:)
2127  real(dp),allocatable :: mpibuf(:)
2128 
2129 ! *************************************************************************
2130 
2131 !DEBUG
2132 !write(std_out,*)' elt_ewald : enter '
2133 !stop
2134 !ENDDEBUG
2135 
2136 !Compute 1st and 2nd derivatives of metric tensor wrt all strain components
2137 !and store for use in inner loop below.
2138 
2139 !Set up parallelism over atoms
2140  paral_atom=(present(comm_atom).and.(my_natom/=natom))
2141  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2142  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2143  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
2144 
2145 !Loop over 2nd strain index
2146  do is2=1,6
2147    kg=idx(2*is2-1);kd=idx(2*is2)
2148    do jj = 1,3
2149      drm(:,jj,is2)=rprimd(kg,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kg,jj)
2150      dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj))
2151    end do
2152 
2153 !  Loop over 1st strain index, upper triangle only
2154    do is1=1,is2
2155 
2156      ka=idx(2*is1-1);kb=idx(2*is1)
2157      d2rm(:,:,is1,is2)=zero
2158      d2gm(:,:,is1,is2)=zero
2159      do jj = 1,3
2160        if(ka==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
2161 &       +rprimd(kb,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kb,jj)
2162        if(ka==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
2163 &       +rprimd(kb,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(kb,jj)
2164        if(kb==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
2165 &       +rprimd(ka,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(ka,jj)
2166        if(kb==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
2167 &       +rprimd(ka,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(ka,jj)
2168 
2169        if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
2170 &       +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj)
2171        if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
2172 &       +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj)
2173        if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
2174 &       +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj)
2175        if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
2176 &       +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj)
2177      end do
2178    end do !is1
2179  end do !is2
2180 
2181 !Add up total charge and sum of $charge^2$ in cell
2182  chsq=zero
2183  ch=zero
2184  do ia=1,natom
2185    ch=ch+zion(typat(ia))
2186    chsq=chsq+zion(typat(ia))**2
2187  end do
2188 
2189 !Compute eta, the Ewald summation convergence parameter,
2190 !for approximately optimized summations:
2191  direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
2192 & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
2193  recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
2194 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
2195 !Here, a bias is introduced, because G-space summation scales
2196 !better than r space summation ! Note : debugging is the most
2197 !easier at fixed eta.
2198  eta=pi*200._dp/33.0_dp*sqrt(1.69_dp*recip/direct)
2199 
2200 !Conduct reciprocal space summations
2201  fac=pi**2/eta ; gsum=zero
2202  ABI_ALLOCATE(d2sumg,(6+3*natom,6))
2203  ABI_ALLOCATE(drhoisr,(3,natom))
2204  ABI_ALLOCATE(drhoisi,(3,natom))
2205  d2sumg(:,:)=zero
2206 
2207 !Sum over G space, done shell after shell until all
2208 !contributions are too small.
2209  ng=0
2210  do
2211    ng=ng+1
2212    newg=0
2213 
2214    do ig3=-ng,ng
2215      do ig2=-ng,ng
2216        do ig1=-ng,ng
2217 
2218 !        Exclude shells previously summed over
2219          if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng&
2220 &         .or. ng==1 ) then
2221 
2222 !          gsq is G dot G = |G|^2
2223            gsq=gmet(1,1)*dble(ig1*ig1)+gmet(2,2)*dble(ig2*ig2)+&
2224 &           gmet(3,3)*dble(ig3*ig3)+2._dp*(gmet(2,1)*dble(ig1*ig2)+&
2225 &           gmet(3,1)*dble(ig1*ig3)+gmet(3,2)*dble(ig3*ig2))
2226 
2227 !          Skip g=0:
2228            if (gsq>1.0d-20) then
2229              arg=fac*gsq
2230 
2231 !            Larger arg gives 0 contribution because of exp(-arg)
2232              if (arg <= 80._dp) then
2233 !              When any term contributes then include next shell
2234                newg=1
2235                term=exp(-arg)/gsq
2236                summr = zero
2237                summi = zero
2238 !              Note that if reduced atomic coordinates xred drift outside
2239 !              of unit cell (outside [0,1)) it is irrelevant in the following
2240 !              term, which only computes a phase.
2241                do ia=1,natom
2242                  arg=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia))
2243 !                Sum real and imaginary parts (avoid complex variables)
2244                  cos_term=cos(arg)
2245                  sin_term=sin(arg)
2246                  summr=summr+zion(typat(ia))*cos_term
2247                  summi=summi+zion(typat(ia))*sin_term
2248                  drhoisr(1,ia)=-two_pi*zion(typat(ia))*sin_term*dble(ig1)
2249                  drhoisi(1,ia)= two_pi*zion(typat(ia))*cos_term*dble(ig1)
2250                  drhoisr(2,ia)=-two_pi*zion(typat(ia))*sin_term*dble(ig2)
2251                  drhoisi(2,ia)= two_pi*zion(typat(ia))*cos_term*dble(ig2)
2252                  drhoisr(3,ia)=-two_pi*zion(typat(ia))*sin_term*dble(ig3)
2253                  drhoisi(3,ia)= two_pi*zion(typat(ia))*cos_term*dble(ig3)
2254                end do
2255 
2256 !              The following two checks avoid an annoying
2257 !              underflow error message
2258                if (abs(summr)<1.d-16) summr=zero
2259                if (abs(summi)<1.d-16) summi=zero
2260 
2261 !              The product of term and summr**2 or summi**2 below
2262 !              can underflow if not for checks above
2263                t1=term*(summr*summr+summi*summi)
2264                gsum=gsum+t1
2265 !              Loop over 2nd strain index
2266                do is2=1,6
2267                  dgss2=dgm(1,1,is2)*dble(ig1*ig1)+dgm(2,2,is2)*dble(ig2*ig2)+&
2268 &                 dgm(3,3,is2)*dble(ig3*ig3)+2._dp*(dgm(2,1,is2)*dble(ig1*ig2)+&
2269 &                 dgm(3,1,is2)*dble(ig1*ig3)+dgm(3,2,is2)*dble(ig3*ig2))
2270 !                Loop over 1st strain index, upper triangle only
2271                  do is1=1,is2
2272                    dgss1=dgm(1,1,is1)*dble(ig1*ig1)+dgm(2,2,is1)*dble(ig2*ig2)+&
2273 &                   dgm(3,3,is1)*dble(ig3*ig3)+2._dp*(dgm(2,1,is1)*dble(ig1*ig2)+&
2274 &                   dgm(3,1,is1)*dble(ig1*ig3)+dgm(3,2,is1)*dble(ig3*ig2))
2275 
2276                    d2gss=d2gm(1,1,is1,is2)*dble(ig1*ig1)+&
2277 &                   d2gm(2,2,is1,is2)*dble(ig2*ig2)+&
2278 &                   d2gm(3,3,is1,is2)*dble(ig3*ig3)+2._dp*&
2279 &                   (d2gm(2,1,is1,is2)*dble(ig1*ig2)+&
2280 &                   d2gm(3,1,is1,is2)*dble(ig1*ig3)+&
2281 &                   d2gm(3,2,is1,is2)*dble(ig3*ig2))
2282 
2283                    d2sumg(is1,is2)=d2sumg(is1,is2)+&
2284 &                   t1*((fac**2 + 2.0_dp*fac/gsq + 2.0_dp/(gsq**2))*dgss1*dgss2 -&
2285 &                   0.5_dp*(fac + 1.0_dp/gsq)*d2gss)
2286                    if(is1<=3) d2sumg(is1,is2)=d2sumg(is1,is2)+&
2287 &                   t1*(fac + 1.0_dp/gsq)*dgss2
2288                    if(is2<=3) d2sumg(is1,is2)=d2sumg(is1,is2)+&
2289 &                   t1*(fac + 1.0_dp/gsq)*dgss1
2290                    if(is1<=3 .and. is2<=3) d2sumg(is1,is2)=d2sumg(is1,is2)+t1
2291 
2292                  end do !is1
2293 
2294 !                Internal strain contributions
2295                  do ia=1,natom
2296                    js=7+3*(ia-1)
2297                    t2(:)=2.0_dp*term*(summr*drhoisr(:,ia)+summi*drhoisi(:,ia))
2298                    d2sumg(js:js+2,is2)=d2sumg(js:js+2,is2)-&
2299 &                   (fac + 1.0_dp/gsq)*dgss2*t2(:)
2300                    if(is2<=3) d2sumg(js:js+2,is2)=d2sumg(js:js+2,is2)-t2(:)
2301                  end do
2302                end do !is2
2303 
2304              end if ! End condition of not larger than 80.0
2305            end if ! End skip g=0
2306          end if ! End triple loop over G s and associated new shell condition
2307        end do
2308      end do
2309    end do
2310 
2311 !  Check if new shell must be calculated
2312    if (newg==0) exit
2313  end do !  End the loop on ng (new shells). Note that there is one exit from this loop.
2314 
2315  sumg=gsum/(two_pi*ucvol)
2316  d2sumg(:,:)=d2sumg(:,:)/(two_pi*ucvol)
2317 
2318  ABI_DEALLOCATE(drhoisr)
2319  ABI_DEALLOCATE(drhoisi)
2320 !Stress tensor is now computed elsewhere (ewald2) hence do not need
2321 !length scale gradients (used to compute them here).
2322 
2323 !Conduct real space summations
2324  reta=sqrt(eta)
2325  fac=2._dp*sqrt(eta/pi)
2326  ABI_ALLOCATE(d2sumr,(6+3*natom,6))
2327  sumr=zero;d2sumr(:,:)=zero
2328 
2329 !In the following a summation is being conducted over all
2330 !unit cells (ir1, ir2, ir3) so it is appropriate to map all
2331 !reduced coordinates xred back into [0,1).
2332 !
2333 !Loop on shells in r-space as was done in g-space
2334  nr=0
2335  do
2336    nr=nr+1
2337    newr=0
2338 !
2339    do ir3=-nr,nr
2340      do ir2=-nr,nr
2341        do ir1=-nr,nr
2342          if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
2343 
2344            do ia0=1,my_natom
2345              ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2346              js=7+3*(ia-1)
2347 !            Map reduced coordinate xred(mu,ia) into [0,1)
2348              fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia))
2349              fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia))
2350              fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia))
2351              do ib=1,natom
2352                fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib))
2353                fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib))
2354                fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib))
2355                r1=dble(ir1)+fracb1-fraca1
2356                r2=dble(ir2)+fracb2-fraca2
2357                r3=dble(ir3)+fracb3-fraca3
2358                rsq=rmet(1,1)*r1*r1+rmet(2,2)*r2*r2+rmet(3,3)*r3*r3+&
2359 &               2.0_dp*(rmet(2,1)*r2*r1+rmet(3,2)*r3*r2+rmet(3,1)*r1*r3)
2360 
2361 !              Avoid zero denominators in 'term':
2362                if (rsq>=1.0d-24) then
2363 
2364 !                Note: erfc(8) is about 1.1e-29,
2365 !                so do not bother with larger arg.
2366 !                Also: exp(-64) is about 1.6e-28,
2367 !                so do not bother with larger arg**2 in exp.
2368                  term=zero
2369                  if (eta*rsq<64.0_dp) then
2370                    newr=1
2371                    rmagn=sqrt(rsq)
2372                    arg=reta*rmagn
2373 !                  derfc computes the complementary error function
2374 !                  dderfc is the derivative of the complementary error function
2375 !                  d2derfc is the 2nd derivative of the complementary error function
2376                    dderfc=-fac*exp(-eta*rsq)
2377                    d2derfc=-2._dp*eta*rmagn*dderfc
2378                    derfc_arg = abi_derfc(arg)
2379                    term=derfc_arg/rmagn
2380                    sumr=sumr+zion(typat(ia))*zion(typat(ib))*term
2381                    tt(:)=rmet(:,1)*r1+rmet(:,2)*r2+rmet(:,3)*r3
2382                    dris(:)=tt(:)/rmagn
2383 !                  Loop over 2nd strain index
2384                    do is2=1,6
2385                      drs2=drm(1,1,is2)*r1*r1+drm(2,2,is2)*r2*r2+&
2386 &                     drm(3,3,is2)*r3*r3+&
2387 &                     2.0_dp*(drm(2,1,is2)*r2*r1+drm(3,2,is2)*r3*r2+&
2388 &                     drm(3,1,is2)*r1*r3)
2389                      dr2=0.5_dp*drs2/rmagn
2390 !                    Loop over 1st strain index, upper triangle only
2391                      do is1=1,is2
2392                        drs1=drm(1,1,is1)*r1*r1+drm(2,2,is1)*r2*r2+&
2393 &                       drm(3,3,is1)*r3*r3+&
2394 &                       2.0_dp*(drm(2,1,is1)*r2*r1+drm(3,2,is1)*r3*r2+&
2395 &                       drm(3,1,is1)*r1*r3)
2396                        dr1=0.5_dp*drs1/rmagn
2397                        d2rs=d2rm(1,1,is1,is2)*r1*r1+d2rm(2,2,is1,is2)*r2*r2+&
2398 &                       d2rm(3,3,is1,is2)*r3*r3+&
2399 &                       2.0_dp*(d2rm(2,1,is1,is2)*r2*r1+d2rm(3,2,is1,is2)*r3*r2+&
2400 &                       d2rm(3,1,is1,is2)*r1*r3)
2401                        d2r=(0.25_dp*d2rs-dr1*dr2)/rmagn
2402                        d2sumr(is1,is2)=d2sumr(is1,is2)+&
2403 &                       zion(typat(ia))*zion(typat(ib))*&
2404 &                       ((d2derfc-2.0_dp*dderfc/rmagn+2.0_dp*derfc_arg/rsq)*dr1*dr2+&
2405 &                       (dderfc-derfc_arg/rmagn)*d2r)/rmagn
2406                      end do !is1
2407 !                    Internal strain contribution
2408                      ts2(:)=drm(:,1,is2)*r1+drm(:,2,is2)*r2+drm(:,3,is2)*r3
2409                      d2ris(:)=ts2(:)/rmagn-0.5_dp*drs2*tt(:)/(rsq*rmagn)
2410 
2411                      d2sumr(js:js+2,is2)=d2sumr(js:js+2,is2)-&
2412 &                     2.0_dp*zion(typat(ia))*zion(typat(ib))*&
2413 &                     ((d2derfc-2.0_dp*dderfc/rmagn+2.0_dp*derfc_arg/rsq)*dr2*dris(:)+&
2414 &                     (dderfc-derfc_arg/rmagn)*d2ris(:))/rmagn
2415                    end do !is2
2416                  end if
2417                end if ! End avoid zero denominators in'term'
2418 
2419              end do ! end loop over ib:
2420            end do ! end loop over ia:
2421          end if ! end triple loop over real space points and associated condition of new shell
2422        end do
2423      end do
2424    end do
2425 
2426 !  Check if new shell must be calculated
2427    if(newr==0) exit
2428  end do !  End loop on nr (new shells). Note that there is an exit within the loop
2429 
2430 !In case of parallelism over atoms: communicate
2431  if (paral_atom) then
2432    call timab(48,1,tsec)
2433    ABI_ALLOCATE(mpibuf,((6+3*natom)*6+1))
2434    mpibuf(1:(6+3*natom)*6)=reshape(d2sumr(:,:),shape=(/((6+3*natom)*6)/))
2435    mpibuf((6+3*natom)*6+1)=sumr
2436    call xmpi_sum(mpibuf,my_comm_atom,ierr)
2437    sumr=mpibuf((6+3*natom)*6+1)
2438    d2sumr(:,:)=reshape(mpibuf(1:(6+3*natom)*6),shape=(/(6+3*natom),6/))
2439    ABI_DEALLOCATE(mpibuf)
2440    call timab(48,2,tsec)
2441  end if
2442 
2443  sumr=0.5_dp*sumr
2444  d2sumr(:,:)=0.5_dp*d2sumr(:,:)
2445  fac=pi*ch**2/(2.0_dp*eta*ucvol)
2446 
2447 !Finally assemble Ewald energy, eew
2448  eew=sumg+sumr-chsq*reta/sqrt(pi)-fac
2449 
2450  elteew(:,:)=d2sumg(:,:)+d2sumr(:,:)
2451 
2452 !Additional term for all strains diagonal (from "fac" term in eew)
2453  elteew(1:3,1:3)=elteew(1:3,1:3)-fac
2454 
2455 !Fill in lower triangle
2456  do is2=2,6
2457    do is1=1,is2-1
2458      elteew(is2,is1)=elteew(is1,is2)
2459    end do
2460  end do
2461 
2462  ABI_DEALLOCATE(d2sumg)
2463  ABI_DEALLOCATE(d2sumr)
2464 
2465 !Output the final values of ng and nr
2466  write(message, '(a,i4,a,i4)' )' elt_ewald : nr and ng are ',nr,' and ',ng
2467  call wrtout(std_out,message,'COLL')
2468 
2469 !Destroy atom table used for parallelism
2470  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2471 
2472 end subroutine elt_ewald

ABINIT/eltxccore [ Functions ]

[ Top ] [ Functions ]

NAME

 eltxccore

FUNCTION

 Compute the core charge contributions to the 2nd derivatives
 of the exchange-correlation energy with respect to all pairs of
 strain or strain and atomic displacement for the frozen wavefunction
 contribution to the elastic tensor. 1st-order potentials representing
 the perturbation by one strain are supplied, and the routine loops
 over the second strain and over all atomic displacements.

COPYRIGHT

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

INPUTS

  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell.
  nfft=number of fft grid points
  ntypat=number of types of atoms in cell.
  n1,n2,n3=fft grid dimensions.
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  rprimd(3,3)=dimensional primitive translation vectors (bohr)
  typat(natom)=integer type for each atom in cell
  ucvol=unit cell volume (bohr**3).
  vxc_core(nfft)=spin-averaged xc potential
  vxc10_core(nfft)=spin-averaged 1st-order xc potential for elastic tensor
  vxc1is_core(nfft)=spin-averaged 1st-order xc potential for internal strain
  xcccrc(ntypat)=XC core correction cutoff radius (bohr) for each atom type
  xccc1d(n1xccc,6,ntypat)=1D core charge function and five derivatives,
   for each type of atom, from psp
  xred(3,natom)=reduced coordinates for atoms in unit cell

OUTPUT

  eltfrxc(6+3*natom,6) = xc frozen wavefunction contribution to the
   elastic tensor

SIDE EFFECTS

  eltfrxc(6+3*natom,6) = xc frozen wavefunction contribution to the
   elastic and internal-strain tensor.  One column is incremented
   by the core contribution.

NOTES

 Note that this routine is related to the mkcore.f routine

PARENTS

      dfpt_eltfrxc

CHILDREN

      free_my_atmtab,get_my_atmtab,timab,xmpi_sum

SOURCE

 719 subroutine eltxccore(eltfrxc,is2_in,my_natom,natom,nfft,ntypat,&
 720 & n1,n1xccc,n2,n3,rprimd,typat,ucvol,vxc_core,vxc10_core,vxc1is_core,&
 721 & xcccrc,xccc1d,xred, &
 722 & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
 723 
 724 
 725 !This section has been created automatically by the script Abilint (TD).
 726 !Do not modify the following lines by hand.
 727 #undef ABI_FUNC
 728 #define ABI_FUNC 'eltxccore'
 729 !End of the abilint section
 730 
 731  implicit none
 732 
 733 !Arguments ------------------------------------
 734 !scalars
 735  integer,intent(in) :: is2_in,n1,n1xccc,n2,n3,my_natom,natom,nfft,ntypat
 736  integer,optional,intent(in) :: comm_atom
 737  real(dp),intent(in) :: ucvol
 738 !arrays
 739  integer,intent(in) :: typat(natom)
 740  integer,optional,target,intent(in) :: mpi_atmtab(:)
 741  real(dp),intent(in) :: rprimd(3,3),vxc10_core(nfft),vxc1is_core(nfft)
 742  real(dp),intent(in) :: vxc_core(nfft),xccc1d(n1xccc,6,ntypat)
 743  real(dp),intent(in) :: xcccrc(ntypat),xred(3,natom)
 744  real(dp),intent(inout) :: eltfrxc(6+3*natom,6)
 745 
 746 !Local variables-------------------------------
 747 !scalars
 748  integer,parameter :: mshift=401
 749  integer :: i1,i2,i3,iat,iatom,ierr,ifft,is1,is2,ishift,ishift1,ishift2
 750  integer :: ishift3,ixp,jj,js,ka,kb,kd,kg,mu,my_comm_atom,nu
 751  logical :: my_atmtab_allocated,paral_atom
 752  real(dp) :: aa,bb,cc,d2rss,dd,delta,delta2div6,deltam1,diff
 753  real(dp) :: difmag,difmag2,difmag2_fact,difmag2_part,drss1,drss2,func1
 754  real(dp) :: func2,range,range2,rangem1,rdiff1,rdiff2,rdiff3
 755  real(dp) :: term1,term2,yy
 756  character(len=500) :: message
 757 !arrays
 758  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
 759  integer :: igrid(3),ii(mshift,3),irange(3),ngfft(3)
 760  integer,pointer :: my_atmtab(:)
 761  real(dp) :: drm(3,3,6),eltfrxc_core(6+3*natom,6),lencp(3),rmet(3,3),rrdiff(mshift,3)
 762  real(dp) :: scale(3),tau(3),ts2(3),tsec(2),tt(3)
 763  real(dp),allocatable :: d2rm(:,:,:,:)
 764 
 765 ! *************************************************************************
 766 
 767 !Compute lengths of cross products for pairs of primitive
 768 !translation vectors (used in setting index search range below)
 769  lencp(1)=cross_elt(rprimd(1,2),rprimd(2,2),rprimd(3,2),&
 770 & rprimd(1,3),rprimd(2,3),rprimd(3,3))
 771  lencp(2)=cross_elt(rprimd(1,3),rprimd(2,3),rprimd(3,3),&
 772 & rprimd(1,1),rprimd(2,1),rprimd(3,1))
 773  lencp(3)=cross_elt(rprimd(1,1),rprimd(2,1),rprimd(3,1),&
 774 & rprimd(1,2),rprimd(2,2),rprimd(3,2))
 775 
 776 !Set up parallelism over atoms
 777  paral_atom=(present(comm_atom).and.(my_natom/=natom))
 778  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
 779  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
 780  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
 781 
 782 !Compute factor R1.(R2xR3)/|R2xR3| etc for 1, 2, 3
 783 !(recall ucvol=R1.(R2xR3))
 784  scale(:)=ucvol/lencp(:)
 785 
 786 !Compute metric tensor in real space rmet
 787  do nu=1,3
 788    rmet(:,nu)=rprimd(1,:)*rprimd(1,nu)+rprimd(2,:)*rprimd(2,nu)+&
 789 &   rprimd(3,:)*rprimd(3,nu)
 790  end do
 791 
 792 !Compute 1st and 2nd derivatives of metric tensor wrt all strain components
 793 !and store for use in inner loop below.
 794 
 795  ABI_ALLOCATE(d2rm,(3,3,6,6))
 796 
 797 !Loop over 2nd strain index
 798  do is2=1,6
 799    kg=idx(2*is2-1);kd=idx(2*is2)
 800    do jj = 1,3
 801      drm(:,jj,is2)=rprimd(kg,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kg,jj)
 802    end do
 803 
 804 !  Loop over 1st strain index
 805    do is1=1,6
 806 
 807      ka=idx(2*is1-1);kb=idx(2*is1)
 808      d2rm(:,:,is1,is2)=0._dp
 809      do jj = 1,3
 810        if(ka==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
 811 &       +rprimd(kb,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kb,jj)
 812        if(ka==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
 813 &       +rprimd(kb,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(kb,jj)
 814        if(kb==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
 815 &       +rprimd(ka,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(ka,jj)
 816        if(kb==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
 817 &       +rprimd(ka,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(ka,jj)
 818      end do
 819    end do !is1
 820  end do !is2
 821 
 822  ngfft(1)=n1
 823  ngfft(2)=n2
 824  ngfft(3)=n3
 825  delta=1.0_dp/(n1xccc-1)
 826  deltam1=n1xccc-1
 827  delta2div6=delta**2/6.0_dp
 828 
 829 !Loop over atoms in unit cell
 830  eltfrxc_core(:,:)=zero
 831 
 832  do iat=1,my_natom
 833    iatom=iat;if (paral_atom) iatom=my_atmtab(iat)
 834    js=7+3*(iatom-1)
 835 !  Set search range (density cuts off perfectly beyond range)
 836 !  Cycle if no range.
 837    range=0.0_dp
 838    range=xcccrc(typat(iatom))
 839    if(range<1.d-16) cycle
 840 
 841    range2=range**2
 842    rangem1=1.0_dp/range
 843 
 844 !  Consider each component in turn
 845    do mu=1,3
 846      tau(mu)=mod(xred(mu,iatom)+1._dp-aint(xred(mu,iatom)),1._dp)
 847 
 848 !    Use tau to find nearest grid point along R(mu)
 849 !    (igrid=0 is the origin; shift by 1 to agree with usual index)
 850      igrid(mu)=nint(tau(mu)*dble(ngfft(mu)))
 851 
 852 !    Use range to compute an index range along R(mu)
 853 !    (add 1 to make sure it covers full range)
 854      irange(mu)=1+nint((range/scale(mu))*dble(ngfft(mu)))
 855 
 856 !    Check that the largest range is smallest than the maximum
 857 !    allowed one
 858      if(2*irange(mu)+1 > mshift)then
 859        write(message, '(a,i0,a)' )' The range around atom',iatom,' is too large.'
 860        MSG_BUG(message)
 861      end if
 862 
 863 !    Set up a counter that explore the relevant range
 864 !    of points around the atom
 865      ishift=0
 866      do ixp=igrid(mu)-irange(mu),igrid(mu)+irange(mu)
 867        ishift=ishift+1
 868        ii(ishift,mu)=1+mod(ngfft(mu)+mod(ixp,ngfft(mu)),ngfft(mu))
 869        rrdiff(ishift,mu)=dble(ixp)/dble(ngfft(mu))-tau(mu)
 870      end do
 871 
 872 !    End loop on mu
 873    end do
 874 
 875 !  Conduct triple loop over restricted range of grid points for iatom
 876 
 877    do ishift3=1,1+2*irange(3)
 878 !    map back to [1,ngfft(3)] for usual fortran index in unit cell
 879      i3=ii(ishift3,3)
 880 !    find vector from atom location to grid point (reduced)
 881      rdiff3=rrdiff(ishift3,3)
 882 
 883      do ishift2=1,1+2*irange(2)
 884        i2=ii(ishift2,2)
 885        rdiff2=rrdiff(ishift2,2)
 886 !      Prepare the computation of difmag2
 887        difmag2_part=rmet(3,3)*rdiff3**2+rmet(2,2)*rdiff2**2&
 888 &       +2.0_dp*rmet(3,2)*rdiff3*rdiff2
 889        difmag2_fact=2.0_dp*(rmet(3,1)*rdiff3+rmet(2,1)*rdiff2)
 890 
 891        do ishift1=1,1+2*irange(1)
 892          rdiff1=rrdiff(ishift1,1)
 893 
 894 !        Compute (rgrid-tau-Rprim)**2
 895          difmag2= difmag2_part+rdiff1*(difmag2_fact+rmet(1,1)*rdiff1)
 896 
 897 !        Only accept contribution inside defined range
 898          if (difmag2<range2) then
 899 
 900 !          Prepare computation of core charge function and derivatives,
 901 !          using splines
 902            difmag=sqrt(difmag2)
 903            if (difmag>=1.0d-10) then
 904              i1=ii(ishift1,1)
 905              yy=difmag*rangem1
 906 
 907 !            Compute index of yy over 1 to n1xccc scale
 908              jj=1+int(yy*(n1xccc-1))
 909              diff=yy-(jj-1)*delta
 910 
 911 !            Will evaluate spline fit (p. 86 Numerical Recipes, Press et al;
 912 !            NOTE error in book for sign of "aa" term in derivative;
 913 !            also see splfit routine).
 914              bb = diff*deltam1
 915              aa = 1.0_dp-bb
 916              cc = aa*(aa**2-1.0_dp)*delta2div6
 917              dd = bb*(bb**2-1.0_dp)*delta2div6
 918 
 919 !            Evaluate spline fit of 1st der of core charge density
 920 !            from xccc1d(:,2,:) and (:,4,:)
 921              func1=aa*xccc1d(jj,2,typat(iatom))+bb*xccc1d(jj+1,2,typat(iatom)) +&
 922 &             cc*xccc1d(jj,4,typat(iatom))+dd*xccc1d(jj+1,4,typat(iatom))
 923              term1=func1*rangem1
 924 !            Evaluate spline fit of 2nd der of core charge density
 925 !            from xccc1d(:,3,:) and (:,5,:)
 926              func2=aa*xccc1d(jj,3,typat(iatom))+bb*xccc1d(jj+1,3,typat(iatom)) +&
 927 &             cc*xccc1d(jj,5,typat(iatom))+dd*xccc1d(jj+1,5,typat(iatom))
 928              term2=func2*rangem1**2
 929 
 930              ifft=i1+n1*(i2-1+n2*(i3-1))
 931              tt(:)=rmet(:,1)*rdiff1+rmet(:,2)*rdiff2+rmet(:,3)*rdiff3
 932 
 933 !            Add contributions to 2nd derivative tensor
 934              drss2=&
 935 &             (rdiff1*(drm(1,1,is2_in)*rdiff1+drm(1,2,is2_in)*rdiff2&
 936 &             +drm(1,3,is2_in)*rdiff3)&
 937 &             +rdiff2*(drm(2,1,is2_in)*rdiff1+drm(2,2,is2_in)*rdiff2&
 938 &             +drm(2,3,is2_in)*rdiff3)&
 939 &             +rdiff3*(drm(3,1,is2_in)*rdiff1+drm(3,2,is2_in)*rdiff2&
 940 &             +drm(3,3,is2_in)*rdiff3))
 941 
 942 !            Loop over 1st strain index
 943              do is1=1,6
 944 
 945                drss1=&
 946 &               (rdiff1*(drm(1,1,is1)*rdiff1+drm(1,2,is1)*rdiff2&
 947 &               +drm(1,3,is1)*rdiff3)&
 948 &               +rdiff2*(drm(2,1,is1)*rdiff1+drm(2,2,is1)*rdiff2&
 949 &               +drm(2,3,is1)*rdiff3)&
 950 &               +rdiff3*(drm(3,1,is1)*rdiff1+drm(3,2,is1)*rdiff2&
 951 &               +drm(3,3,is1)*rdiff3))
 952 
 953                d2rss=&
 954 &               (rdiff1*(d2rm(1,1,is1,is2_in)*rdiff1+d2rm(1,2,is1,is2_in)*rdiff2&
 955 &               +d2rm(1,3,is1,is2_in)*rdiff3)&
 956 &               +rdiff2*(d2rm(2,1,is1,is2_in)*rdiff1+d2rm(2,2,is1,is2_in)*rdiff2&
 957 &               +d2rm(2,3,is1,is2_in)*rdiff3)&
 958 &               +rdiff3*(d2rm(3,1,is1,is2_in)*rdiff1+d2rm(3,2,is1,is2_in)*rdiff2&
 959 &               +d2rm(3,3,is1,is2_in)*rdiff3))
 960 
 961 !              Vall(0) X Rhocore(2) term
 962                eltfrxc_core(is1,is2_in)=eltfrxc_core(is1,is2_in)+0.25_dp*&
 963 &               (vxc_core(ifft)*(term1*(d2rss/difmag&
 964 &               -drss1*drss2/difmag**3)&
 965 &               +term2*drss1*drss2/difmag**2))
 966 
 967 !              Vall(1) X Rhocore(1) term
 968                eltfrxc_core(is1,is2_in)=eltfrxc_core(is1,is2_in)+0.25_dp*&
 969 &               vxc10_core(ifft)*drss1*term1/difmag
 970                eltfrxc_core(is2_in,is1)=eltfrxc_core(is2_in,is1)+0.25_dp*&
 971 &               vxc10_core(ifft)*drss1*term1/difmag
 972 
 973 !              End loop in is1
 974              end do
 975 !            Internal strain contributions
 976              ts2(:)=drm(:,1,is2_in)*rdiff1+drm(:,2,is2_in)*rdiff2&
 977 &             +drm(:,3,is2_in)*rdiff3
 978 
 979              eltfrxc_core(js:js+2,is2_in)=eltfrxc_core(js:js+2,is2_in)&
 980 &             -(vxc1is_core(ifft)*term1/difmag&
 981 &             +0.5_dp*vxc_core(ifft)*(term2-term1/difmag)*drss2/difmag**2)*tt(:)&
 982 &             -(vxc_core(ifft)*term1/difmag)*ts2(:)
 983 
 984 !            End of the condition for the distance not to vanish
 985            end if
 986 
 987 !          End of condition to be inside the range
 988          end if
 989 
 990 !        End loop on ishift1
 991        end do
 992 
 993 !      End loop on ishift2
 994      end do
 995 
 996 !    End loop on ishift3
 997    end do
 998 
 999 !  End loop on atoms
1000  end do
1001 
1002 !In case of parallelism over atoms: communicate
1003  if (paral_atom) then
1004    call timab(48,1,tsec)
1005    call xmpi_sum(eltfrxc_core,my_comm_atom,ierr)
1006    call timab(48,2,tsec)
1007  end if
1008 
1009 !Add core contribution to XC elastic tensor
1010  eltfrxc(:,:)=eltfrxc(:,:)+eltfrxc_core(:,:)
1011 
1012 !Destroy atom table used for parallelism
1013  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1014 
1015  ABI_DEALLOCATE(d2rm)
1016 
1017  contains
1018 
1019    function cross_elt(xx,yy,zz,aa,bb,cc)
1020 !Define magnitude of cross product of two vectors
1021 
1022 !This section has been created automatically by the script Abilint (TD).
1023 !Do not modify the following lines by hand.
1024 #undef ABI_FUNC
1025 #define ABI_FUNC 'cross_elt'
1026 !End of the abilint section
1027 
1028    real(dp) :: cross_elt
1029    real(dp),intent(in) :: xx,yy,zz,aa,bb,cc
1030    cross_elt=sqrt((yy*cc-zz*bb)**2+(zz*aa-xx*cc)**2+(xx*bb-yy*aa)**2)
1031  end function cross_elt
1032 
1033 end subroutine eltxccore

ABINIT/m_dfpt_elt [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dfpt_elt

FUNCTION

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DRH, DCA, XG, GM, AR, MB)
  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_dfpt_elt
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use m_abicore
33  use m_errors
34  use m_xmpi
35 
36  use m_time,        only : timab
37  use m_special_funcs,  only : abi_derfc
38  use m_geometry,    only : metric
39  use m_cgtools,     only : dotprod_vn
40  use m_pawtab,      only : pawtab_type,pawtab_free,pawtab_nullify
41  use m_pawrad,      only : pawrad_type,pawrad_init,pawrad_free
42  use m_pawpsp,      only : pawpsp_cg
43  use m_paw_numeric, only : paw_spline
44  use m_spacepar,    only : redgr
45  use m_atm2fft,     only : atm2fft, dfpt_atm2fft
46  use m_mkcore,      only : dfpt_mkcore
47  use m_dfpt_mkvxcstr, only : dfpt_mkvxcstr
48  use m_paral_atom, only : get_my_atmtab, free_my_atmtab
49  use m_mpinfo,  only : ptabs_fourdp, proc_distrb_cycle
50  use m_fftcore,      only : sphereboundary
51  use m_fft,             only : fourdp
52 
53  implicit none
54 
55  private