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 (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
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
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
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
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*&
966
967 !              Vall(1) X Rhocore(1) term
968                eltfrxc_core(is1,is2_in)=eltfrxc_core(is1,is2_in)+0.25_dp*&
970                eltfrxc_core(is2_in,is1)=eltfrxc_core(is2_in,is1)+0.25_dp*&
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&
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 (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
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