TABLE OF CONTENTS


ABINIT/d2kindstr2 [ Functions ]

[ Top ] [ Functions ]

NAME

 d2kindstr2

FUNCTION

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

INPUTS

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

OUTPUT

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

NOTES

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

SOURCE

1559 subroutine d2kindstr2(cwavef,ecut,ecutsm,effmass_free,ekinout,gmet,gprimd,&
1560 &            istwfk,kg_k,kpt,npw,nspinor)
1561 
1562 !Arguments ------------------------------------
1563 !scalars
1564  integer,intent(in) :: istwfk,npw,nspinor
1565  real(dp),intent(in) :: ecut,ecutsm,effmass_free
1566 !arrays
1567  integer,intent(in) :: kg_k(3,npw)
1568  real(dp),intent(in) :: cwavef(2,npw*nspinor),gmet(3,3),gprimd(3,3),kpt(3)
1569  real(dp),intent(inout) :: ekinout(36) !vz_i
1570 
1571 !Local variables-------------------------------
1572 !scalars
1573  integer,parameter :: im=2,re=1
1574  integer :: ig,igs,ii,ispinor,istr1,istr2,ka,kb,kd,kg
1575  real(dp) :: d2fkin,d2fsm,d2kinacc,d2kpg2,dfkin,dfsm,dkpg21,dkpg22,ecutsm_inv
1576  real(dp) :: fsm,gpk1,gpk2,gpk3,htpisq,kpg2,term,xx
1577 !arrays
1578  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1579  real(dp) :: d2gm(3,3),dgm01(3,3),dgm10(3,3)
1580 
1581 ! *************************************************************************
1582 !
1583 !htpisq is (1/2) (2 Pi) **2:
1584    htpisq=0.5_dp*(two_pi)**2
1585 
1586    ecutsm_inv=0.0_dp
1587    if(ecutsm>1.0d-20)ecutsm_inv=1/ecutsm
1588 
1589 !Loop over 2nd strain index
1590    do istr2=1,6
1591 !  Loop over 1st strain index, upper triangle only
1592      do istr1=1,istr2
1593 
1594        ka=idx(2*istr1-1);kb=idx(2*istr1);kg=idx(2*istr2-1);kd=idx(2*istr2)
1595 
1596        do ii = 1,3
1597          dgm01(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
1598          dgm10(:,ii)=-(gprimd(kg,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kg,ii))
1599        end do
1600 
1601        d2gm(:,:)=0._dp
1602        do ii = 1,3
1603          if(ka==kg) d2gm(:,ii)=d2gm(:,ii)&
1604 &         +gprimd(kb,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kb,ii)
1605          if(ka==kd) d2gm(:,ii)=d2gm(:,ii)&
1606 &         +gprimd(kb,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(kb,ii)
1607          if(kb==kg) d2gm(:,ii)=d2gm(:,ii)&
1608 &         +gprimd(ka,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(ka,ii)
1609          if(kb==kd) d2gm(:,ii)=d2gm(:,ii)&
1610 &         +gprimd(ka,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(ka,ii)
1611        end do
1612        d2gm(:,:)=0.5_dp*d2gm(:,:)
1613 
1614        d2kinacc=0._dp
1615 
1616 !    loop on spinor index
1617        do ispinor=1,nspinor
1618          igs=(ispinor-1)*npw
1619 !      loop on plane waves
1620          do ig=1,npw
1621            gpk1=dble(kg_k(1,ig))+kpt(1)
1622            gpk2=dble(kg_k(2,ig))+kpt(2)
1623            gpk3=dble(kg_k(3,ig))+kpt(3)
1624            kpg2=htpisq*&
1625 &           ( gmet(1,1)*gpk1**2+         &
1626 &           gmet(2,2)*gpk2**2+         &
1627 &           gmet(3,3)*gpk3**2          &
1628 &           +2.0_dp*(gpk1*gmet(1,2)*gpk2+  &
1629 &           gpk1*gmet(1,3)*gpk3+  &
1630 &           gpk2*gmet(2,3)*gpk3 )  )
1631            dkpg21=htpisq*&
1632 &           ( dgm01(1,1)*gpk1**2+         &
1633 &           dgm01(2,2)*gpk2**2+         &
1634 &           dgm01(3,3)*gpk3**2          &
1635 &           +2.0_dp*(gpk1*dgm01(1,2)*gpk2+  &
1636 &           gpk1*dgm01(1,3)*gpk3+  &
1637 &           gpk2*dgm01(2,3)*gpk3 )  )
1638            dkpg22=htpisq*&
1639 &           ( dgm10(1,1)*gpk1**2+         &
1640 &           dgm10(2,2)*gpk2**2+         &
1641 &           dgm10(3,3)*gpk3**2          &
1642 &           +2.0_dp*(gpk1*dgm10(1,2)*gpk2+  &
1643 &           gpk1*dgm10(1,3)*gpk3+  &
1644 &           gpk2*dgm10(2,3)*gpk3 )  )
1645            d2kpg2=htpisq*&
1646 &           ( d2gm(1,1)*gpk1**2+         &
1647 &           d2gm(2,2)*gpk2**2+         &
1648 &           d2gm(3,3)*gpk3**2          &
1649 &           +2.0_dp*(gpk1*d2gm(1,2)*gpk2+  &
1650 &           gpk1*d2gm(1,3)*gpk3+  &
1651 &           gpk2*d2gm(2,3)*gpk3 )  )
1652 
1653            if(kpg2>ecut-tol12)then
1654              dfkin=0._dp
1655              d2fkin=0._dp
1656            elseif(kpg2>ecut-ecutsm)then
1657 !          This kinetic cutoff smoothing function and its xx derivatives
1658 !          were produced with Mathematica and the fortran code has been
1659 !          numerically checked against Mathematica.
1660              xx=(ecut-kpg2)*ecutsm_inv
1661              fsm=1.0_dp/(xx**2*(3+xx*(1+xx*(-6+3*xx))))
1662              dfsm=-3.0_dp*(-1+xx)**2*xx*(2+5*xx)*fsm**2
1663              d2fsm=6.0_dp*xx**2*(9+xx*(8+xx*(-52+xx*(-3+xx*(137+xx*&
1664 &             (-144+45*xx))))))*fsm**3
1665              dfkin=fsm-ecutsm_inv*kpg2*dfsm
1666              d2fkin=ecutsm_inv*(-2.0_dp*dfsm+ecutsm_inv*kpg2*d2fsm)
1667            else
1668              dfkin=1._dp
1669              d2fkin=0._dp
1670            end if
1671 
1672 !        accumulate kinetic energy 2nd derivative with wavefunction components
1673            term=d2fkin*dkpg21*dkpg22 + dfkin*d2kpg2
1674            if(istwfk==2 .and. ig/=1)term=2.0_dp*term
1675            if(istwfk>2)term=2.0_dp*term
1676            d2kinacc=d2kinacc + term*(cwavef(re,ig+igs)**2 + cwavef(im,ig+igs)**2)
1677 
1678          end do  !ig
1679        end do !ispinor
1680 
1681        ekinout(istr1+6*(istr2-1))=d2kinacc/effmass_free
1682 
1683      end do !istr1
1684    end do !istr2
1685 
1686   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=information 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

SOURCE

1718 subroutine dfpt_eltfrhar(eltfrhar,rprimd,gsqcut,mpi_enreg,nfft,ngfft,rhog)
1719 
1720 !Arguments ------------------------------------
1721 !scalars
1722  integer,intent(in) :: nfft
1723  real(dp),intent(in) :: gsqcut
1724  type(MPI_type),intent(in) :: mpi_enreg
1725 !arrays
1726  integer,intent(in) :: ngfft(18)
1727  real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3)
1728  real(dp),intent(out) :: eltfrhar(6,6)
1729 
1730 !Local variables-------------------------------
1731 !scalars
1732  integer,parameter :: im=2,re=1
1733  integer :: i1,i2,i23,i3,id2,id3,ierr,ig,ig2,ig3,ii,ii1,ing,istr1,istr2,jj
1734  integer :: ka,kb,kd,kg,me_fft,n1,n2,n3,nproc_fft
1735  real(dp),parameter :: tolfix=1.000000001_dp
1736  real(dp) :: cutoff,d2eacc,d2etot,d2gs,deacc01,deacc10,dgs01,dgs10,eacc,fact,gs
1737  real(dp) :: term,ucvol
1738 !arrays
1739  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1740  integer :: id(3)
1741  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
1742  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1743  real(dp) :: d2gm(3,3),dgm01(3,3),dgm10(3,3),gmet(3,3),gprimd(3,3),gqr(3)
1744  real(dp) :: rmet(3,3),tsec(2)
1745  real(dp),allocatable :: gq(:,:)
1746 
1747 ! *************************************************************************
1748 
1749 !Compute gmet, gprimd and ucvol from rprimd
1750  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1751 
1752  eltfrhar(:,:)=0.0_dp
1753 
1754  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1755  me_fft=ngfft(11)
1756  nproc_fft=ngfft(10)
1757 
1758 !Get the distrib associated with this fft_grid
1759  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1760 
1761 !Initialize a few quantities
1762  fact=0.5_dp*ucvol/pi
1763  cutoff=gsqcut*tolfix
1764 
1765 !In order to speed the routine, precompute the components of g+q
1766 !Also check if the booked space was large enough...
1767  ABI_MALLOC(gq,(3,max(n1,n2,n3)))
1768  do ii=1,3
1769    id(ii)=ngfft(ii)/2+2
1770    do ing=1,ngfft(ii)
1771      ig=ing-(ing/id(ii))*ngfft(ii)-1
1772      gq(ii,ing)=ig
1773    end do
1774  end do
1775 
1776 !Loop over 2nd strain index
1777  do istr2=1,6
1778 !  Loop over 1st strain index, upper triangle only
1779    do istr1=1,istr2
1780 
1781      ka=idx(2*istr1-1);kb=idx(2*istr1);kg=idx(2*istr2-1);kd=idx(2*istr2)
1782 
1783      do ii = 1,3
1784        dgm01(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii))
1785        dgm10(:,ii)=-(gprimd(kg,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kg,ii))
1786      end do
1787 
1788      d2gm(:,:)=0._dp
1789      do ii = 1,3
1790        if(ka==kg) d2gm(:,ii)=d2gm(:,ii)&
1791 &       +gprimd(kb,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kb,ii)
1792        if(ka==kd) d2gm(:,ii)=d2gm(:,ii)&
1793 &       +gprimd(kb,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(kb,ii)
1794        if(kb==kg) d2gm(:,ii)=d2gm(:,ii)&
1795 &       +gprimd(ka,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(ka,ii)
1796        if(kb==kd) d2gm(:,ii)=d2gm(:,ii)&
1797 &       +gprimd(ka,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(ka,ii)
1798      end do
1799      d2gm(:,:)=0.5_dp*d2gm(:,:)
1800 
1801 !    initialize energy accumulator
1802      eacc=0._dp
1803      deacc01=0._dp
1804      deacc10=0._dp
1805      d2eacc=0._dp
1806 
1807      id2=n2/2+2
1808      id3=n3/2+2
1809 !    Triple loop on each dimension
1810      do i3=1,n3
1811        ig3=i3-(i3/id3)*n3-1
1812        gqr(3)=gq(3,i3)
1813        do i2=1,n2
1814          if (fftn2_distrib(i2)==me_fft) then
1815            gqr(2)=gq(2,i2)
1816            ig2=i2-(i2/id2)*n2-1
1817            i23=n1*(ffti2_local(i2)-1 +(n2/nproc_fft)*(i3-1))
1818 !          Do the test that eliminates the Gamma point outside
1819 !          of the inner loop
1820            ii1=1
1821            if(i23==0 .and. ig2==0 .and. ig3==0)then
1822              ii1=2
1823            end if
1824 
1825 !          Final inner loop on the first dimension
1826 !          (note the lower limit)
1827            do i1=ii1,n1
1828              gqr(1)=gq(1,i1)
1829              gs=(gmet(1,1)*gqr(1)*gqr(1)+gmet(2,2)*gqr(2)*gqr(2)+&
1830 &             gmet(3,3)*gqr(3)*gqr(3)+2._dp*&
1831 &             (gmet(1,2)*gqr(1)*gqr(2) + gmet(1,3)*gqr(1)*gqr(3)+&
1832 &             gmet(2,3)*gqr(2)*gqr(3)) )
1833              ii=i1+i23
1834              if(gs<=cutoff)then
1835                dgs01=(dgm01(1,1)*gqr(1)*gqr(1)+dgm01(2,2)*gqr(2)*gqr(2)+&
1836 &               dgm01(3,3)*gqr(3)*gqr(3)+2._dp*&
1837 &               (dgm01(1,2)*gqr(1)*gqr(2) + dgm01(1,3)*gqr(1)*gqr(3)+&
1838 &               dgm01(2,3)*gqr(2)*gqr(3)) )
1839                dgs10=(dgm10(1,1)*gqr(1)*gqr(1)+dgm10(2,2)*gqr(2)*gqr(2)+&
1840 &               dgm10(3,3)*gqr(3)*gqr(3)+2._dp*&
1841 &               (dgm10(1,2)*gqr(1)*gqr(2) + dgm10(1,3)*gqr(1)*gqr(3)+&
1842 &               dgm10(2,3)*gqr(2)*gqr(3)) )
1843                d2gs =(d2gm(1,1)*gqr(1)*gqr(1)+d2gm(2,2)*gqr(2)*gqr(2)+&
1844 &               d2gm(3,3)*gqr(3)*gqr(3)+2._dp*&
1845 &               (d2gm(1,2)*gqr(1)*gqr(2) + d2gm(1,3)*gqr(1)*gqr(3)+&
1846 &               d2gm(2,3)*gqr(2)*gqr(3)) )
1847 
1848                term=(rhog(re,ii)**2+rhog(im,ii)**2)/gs
1849                eacc=eacc+term
1850                deacc01=deacc01+dgs01*term/gs
1851                deacc10=deacc10+dgs10*term/gs
1852                d2eacc=d2eacc+(-d2gs+2._dp*dgs01*dgs10/gs)*term/gs
1853              end if
1854 
1855 !            End loop on i1
1856            end do
1857          end if
1858 !        End loop on i2
1859        end do
1860 !      End loop on i3
1861      end do
1862 
1863 !    Add contributions taking account diagonal strain terms (from ucvol
1864 !    derivatives)
1865      d2etot=d2eacc
1866      if(istr1<=3) d2etot=d2etot+deacc10
1867      if(istr2<=3) d2etot=d2etot+deacc01
1868      if(istr1<=3 .and. istr2<=3) d2etot=d2etot+eacc
1869 
1870      eltfrhar(istr1,istr2)=fact*d2etot
1871 
1872 !    End loop on istr1
1873    end do
1874 !  End loop in istr2
1875  end do
1876 
1877  ABI_FREE(gq)
1878 
1879 !Init mpi_comm
1880  call timab(48,1,tsec)
1881  call xmpi_sum(eltfrhar,mpi_enreg%comm_fft,ierr)
1882  call timab(48,2,tsec)
1883 
1884 !Fill in lower triangle
1885  do jj=2,6
1886    do ii=1,jj-1
1887      eltfrhar(jj,ii)=eltfrhar(ii,jj)
1888    end do
1889  end do
1890 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_mem*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
  mband_mem=maximum number of bands in memory
  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

SOURCE

1372 subroutine dfpt_eltfrkin(cg,eltfrkin,ecut,ecutsm,effmass_free,&
1373 &  istwfk,kg,kptns,mband,mband_mem,mgfft,mkmem,mpi_enreg,&
1374 &  mpw,nband,nkpt,ngfft,npwarr,nspinor,nsppol,occ,rprimd,wtk)
1375 
1376 !Arguments ------------------------------------
1377 !scalars
1378  integer,intent(in) :: mband,mband_mem,mgfft,mkmem,mpw,nkpt,nspinor,nsppol
1379  real(dp),intent(in) :: ecut,ecutsm,effmass_free
1380  type(MPI_type),intent(in) :: mpi_enreg
1381 !arrays
1382  integer,intent(in) :: istwfk(nkpt),kg(3,mpw*mkmem),nband(nkpt*nsppol)
1383  integer,intent(in) :: ngfft(18),npwarr(nkpt)
1384  real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem*nsppol),kptns(3,nkpt)
1385  real(dp),intent(in) :: occ(mband*nkpt*nsppol),rprimd(3,3),wtk(nkpt)
1386  real(dp),intent(out) :: eltfrkin(6,6)
1387 
1388 !Local variables-------------------------------
1389 !scalars
1390  integer :: bdtot_index,iband,icg,ierr,ii,ikg
1391  integer :: ikpt,index,ipw,isppol,istwf_k,jj,master,me,n1,n2
1392  integer :: n3,nband_k,nkinout,npw_k,spaceComm
1393  integer :: nband_me, iband_me
1394  real(dp) :: ucvol
1395 !arrays
1396  integer,allocatable :: gbound(:,:),kg_k(:,:)
1397  real(dp) :: gmet(3,3),gprimd(3,3),kpoint(3),rmet(3,3),tsec(2)
1398  real(dp),allocatable :: cwavef(:,:),ekinout(:)
1399  real(dp),allocatable :: eltfrkink(:,:)
1400 
1401 ! *************************************************************************
1402 
1403  DBG_ENTER("COLL")
1404 
1405 !Default for sequential use
1406  master=0
1407 !Init mpi_comm
1408  spaceComm=mpi_enreg%comm_cell
1409  me=mpi_enreg%me_kpt
1410 
1411 !Compute gmet, gprimd and ucvol from rprimd
1412  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1413 
1414  eltfrkin(:,:)=0.0_dp
1415  bdtot_index=0
1416  icg=0
1417 
1418  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1419  ABI_MALLOC(kg_k,(3,mpw))
1420  ABI_MALLOC(cwavef,(2,mpw*nspinor))
1421  ABI_MALLOC(eltfrkink,(6,6))
1422 
1423 !Define k-points distribution
1424 
1425 !LOOP OVER SPINS
1426  do isppol=1,nsppol
1427    ikg=0
1428 
1429 !  Loop over k points
1430    do ikpt=1,nkpt
1431 
1432      nband_k=nband(ikpt+(isppol-1)*nkpt)
1433      istwf_k=istwfk(ikpt)
1434      npw_k=npwarr(ikpt)
1435 
1436 !    Skip this k-point if not the proper processor
1437      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
1438        bdtot_index=bdtot_index+nband_k
1439        cycle
1440      end if
1441 
1442 ! find number of bands I will actually treat
1443      nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me)
1444 
1445      ABI_MALLOC(gbound,(2*mgfft+8,2))
1446      kpoint(:)=kptns(:,ikpt)
1447 
1448      kg_k(:,:) = 0
1449 
1450 
1451 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(ikg,kg,kg_k,npw_k)
1452      do ipw=1,npw_k
1453        kg_k(1,ipw)=kg(1,ipw+ikg)
1454        kg_k(2,ipw)=kg(2,ipw+ikg)
1455        kg_k(3,ipw)=kg(3,ipw+ikg)
1456      end do
1457 
1458      call sphereboundary(gbound,istwf_k,kg_k,mgfft,npw_k)
1459 
1460      index=1+icg
1461 
1462      eltfrkink(:,:)=0.0_dp
1463 
1464      nkinout=6*6
1465      ABI_MALLOC(ekinout,(nkinout))
1466      ekinout(:)=zero
1467 
1468      iband_me = 0
1469      do iband=1,nband_k
1470 
1471        if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= me) cycle
1472        iband_me = iband_me + 1
1473 
1474        cwavef(:,1:npw_k*nspinor)=cg(:,1+(iband_me-1)*npw_k*nspinor+icg:iband_me*npw_k*nspinor+icg)
1475 
1476        call d2kindstr2(cwavef,ecut,ecutsm,effmass_free,ekinout,gmet,gprimd,&
1477 &       istwf_k,kg_k,kpoint,npw_k,nspinor)
1478 
1479        eltfrkink(:,:)=eltfrkink(:,:)+ occ(iband+bdtot_index)* reshape(ekinout(:), (/6,6/) )
1480 
1481      end do !iband
1482 
1483      ABI_FREE(ekinout)
1484 
1485      eltfrkin(:,:)=eltfrkin(:,:)+wtk(ikpt)*eltfrkink(:,:)
1486 
1487      ABI_FREE(gbound)
1488 
1489      bdtot_index=bdtot_index+nband_k
1490 
1491      if (mkmem/=0) then
1492 !      Handle case in which kg, cg, are kept in core
1493        icg=icg+npw_k*nspinor*nband_me
1494        ikg=ikg+npw_k
1495      end if
1496 
1497    end do
1498  end do  ! End loops on isppol and ikpt
1499 
1500 !Fill in lower triangle
1501  do jj=2,6
1502    do ii=1,jj-1
1503      eltfrkin(jj,ii)=eltfrkin(ii,jj)
1504    end do
1505  end do
1506 
1507 !Accumulate eltfrkin on all proc.
1508  call timab(48,1,tsec)
1509  call xmpi_sum(eltfrkin,spaceComm,ierr)
1510  call timab(48,2,tsec)
1511 
1512  ABI_FREE(cwavef)
1513  ABI_FREE(eltfrkink)
1514  ABI_FREE(kg_k)
1515 
1516  DBG_EXIT("COLL")
1517 
1518   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=information 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.

SOURCE

1032 subroutine dfpt_eltfrloc(atindx,eltfrloc,gmet,gprimd,gsqcut,mgfft,&
1033 &  mpi_enreg,mqgrid,natom,nattyp,nfft,ngfft,ntypat,ph1d,qgrid,rhog,vlspl)
1034 
1035 !Arguments ------------------------------------
1036 !scalars
1037  integer,intent(in) :: mgfft,mqgrid,natom,nfft,ntypat
1038  real(dp),intent(in) :: gsqcut
1039  type(MPI_type),intent(in) :: mpi_enreg
1040 !arrays
1041  integer,intent(in) :: atindx(natom),nattyp(ntypat),ngfft(18)
1042  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),ph1d(2,3*(2*mgfft+1)*natom)
1043  real(dp),intent(in) :: qgrid(mqgrid),rhog(2,nfft),vlspl(mqgrid,2,ntypat)
1044  real(dp),intent(out) :: eltfrloc(6+3*natom,6)
1045 
1046 !Local variables-------------------------------
1047 !scalars
1048  integer,parameter :: im=2,re=1
1049  integer :: i1,i2,i3,ia,ia1,ia2,id1,id2,id3,ielt,ieltx,ierr,ig1,ig2,ig3,ii
1050  integer :: is1,is2,itypat,jj,ka,kb,kd,kg,me_fft,n1,n2,n3,nproc_fft
1051  real(dp),parameter :: tolfix=1.0000001_dp
1052  real(dp) :: aa,bb,cc,cutoff,d2g
1053  real(dp) :: dd,dg1,dg2,diff,dq
1054  real(dp) :: dq2div6,dqdiv6,dqm1,ee,ff,gmag,gsquar
1055  real(dp) :: sfi,sfr,term,term1
1056 !real(dp) :: ph1_elt,ph2_elt,ph3_elt,phi_elt,phr_elt
1057  real(dp) :: term2,term3,term4,term5,vion1,vion2,vion3
1058 !arrays
1059  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1060  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
1061  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
1062  real(dp) :: dgm(3,3,6),tsec(2)
1063  real(dp),allocatable :: d2gm(:,:,:,:),elt_work(:,:)
1064 
1065 ! *************************************************************************
1066 
1067  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1068  me_fft=ngfft(11) ; nproc_fft=ngfft(10)
1069 
1070 !Get the distrib associated with this fft_grid
1071  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
1072 
1073 !-----
1074 !Compute 1st and 2nd derivatives of metric tensor wrt all strain components
1075 !and store for use in inner loop below.
1076  ABI_MALLOC(d2gm,(3,3,6,6))
1077 
1078 !Loop over 2nd strain index
1079  do is2=1,6
1080    kg=idx(2*is2-1);kd=idx(2*is2)
1081    do jj = 1,3
1082      dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj))
1083    end do
1084 !  Loop over 1st strain index, upper triangle only
1085    do is1=1,is2
1086      ka=idx(2*is1-1);kb=idx(2*is1)
1087      d2gm(:,:,is1,is2)=0._dp
1088      do jj = 1,3
1089        if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1090 &       +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj)
1091        if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1092 &       +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj)
1093        if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1094 &       +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj)
1095        if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1096 &       +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj)
1097      end do
1098    end do !is1
1099  end do !is2
1100 
1101 !Zero out array to permit accumulation over atom types below:
1102  eltfrloc(:,:)=0.0_dp
1103 
1104  dq=(qgrid(mqgrid)-qgrid(1))/dble(mqgrid-1)
1105  dqm1=1.0_dp/dq
1106  dqdiv6=dq/6.0_dp
1107  dq2div6=dq**2/6.0_dp
1108  cutoff=gsqcut*tolfix
1109  id1=n1/2+2
1110  id2=n2/2+2
1111  id3=n3/2+2
1112 
1113  ia1=1
1114  do itypat=1,ntypat
1115 !  ia1,ia2 sets range of loop over atoms:
1116    ia2=ia1+nattyp(itypat)-1
1117    ii=0
1118    do i3=1,n3
1119      ig3=i3-(i3/id3)*n3-1
1120      do i2=1,n2
1121        if (fftn2_distrib(i2)==me_fft) then
1122          ig2=i2-(i2/id2)*n2-1
1123          do i1=1,n1
1124            ig1=i1-(i1/id1)*n1-1
1125 
1126            ii=ii+1
1127 !          Skip G=0:
1128            if (ig1==0 .and. ig2==0 .and. ig3==0) cycle
1129 
1130 !          Skip G**2 outside cutoff:
1131            gsquar=gsq_elt(ig1,ig2,ig3)
1132            if (gsquar<=cutoff) then
1133              gmag=sqrt(gsquar)
1134 
1135 !            Compute vion(G) for given type of atom
1136              jj=1+int(gmag*dqm1)
1137              diff=gmag-qgrid(jj)
1138 
1139 !            Evaluate spline fit from q^2 V(q) to get V(q):
1140 !            (p. 86 Numerical Recipes, Press et al; NOTE error in book for sign
1141 !             of "aa" term in derivative; also see splfit routine).
1142              bb = diff*dqm1
1143              aa = 1.0_dp-bb
1144              cc = aa*(aa**2-1.0_dp)*dq2div6
1145              dd = bb*(bb**2-1.0_dp)*dq2div6
1146              term1 = (aa*vlspl(jj,1,itypat)+bb*vlspl(jj+1,1,itypat) +&
1147 &             cc*vlspl(jj,2,itypat)+dd*vlspl(jj+1,2,itypat))
1148              vion1=term1 / gsquar
1149 
1150 !            Also get dV(q)/dq:
1151 !            (note correction of Numerical Recipes sign error
1152 !             before (3._dp*aa**2-1._dp)
1153              ee= vlspl(jj+1,1,itypat)-vlspl(jj,1,itypat)
1154              ff=  (3._dp*bb**2-1._dp)*vlspl(jj+1,2,itypat) &
1155 &             - (3._dp*aa**2-1._dp)*vlspl(jj,2,itypat)
1156              term2 = ee*dqm1 + ff*dqdiv6
1157              vion2 = term2/gsquar - 2._dp*term1/(gsquar*gmag)
1158 
1159 !            Also get V''(q)
1160              term3=aa*vlspl(jj,2,itypat)+bb*vlspl(jj+1,2,itypat)
1161              vion3 = (term3 - 4.0_dp*term2/gmag + 6._dp*term1/gsquar)/gsquar
1162 
1163 !            Assemble structure factor over all atoms of given type:
1164              sfr=zero;sfi=zero
1165              do ia=ia1,ia2
1166                sfr=sfr+phre_elt(ig1,ig2,ig3,ia)
1167                sfi=sfi-phimag_elt(ig1,ig2,ig3,ia)
1168              end do
1169              term=(rhog(re,ii)*sfr+rhog(im,ii)*sfi)
1170 
1171 !            Loop over 2nd strain index
1172              do is2=1,6
1173                dg2=0.5_dp*dgsqds_elt(ig1,ig2,ig3,is2)/gmag
1174 !              Loop over 1st strain index, upper triangle only
1175                do is1=1,is2
1176                  dg1=0.5_dp*dgsqds_elt(ig1,ig2,ig3,is1)/gmag
1177                  d2g=(0.25_dp*d2gsqds_elt(ig1,ig2,ig3,is1,is2)-dg1*dg2)/gmag
1178                  eltfrloc(is1,is2)=eltfrloc(is1,is2)+&
1179 &                 term*(vion3*dg1*dg2+vion2*d2g)
1180                  if(is2<=3)&
1181 &                 eltfrloc(is1,is2)=eltfrloc(is1,is2)-term*vion2*dg1
1182                  if(is1<=3)&
1183 &                 eltfrloc(is1,is2)=eltfrloc(is1,is2)-term*vion2*dg2
1184                  if(is1<=3 .and. is2<=3)&
1185 &                 eltfrloc(is1,is2)=eltfrloc(is1,is2)+term*vion1
1186                end do !is1
1187 
1188 !              Internal strain section - loop over current atoms
1189                do ia=ia1,ia2
1190                  if(is2 <=3) then
1191                    term4=vion2*dg2-vion1
1192                  else
1193                    term4=vion2*dg2
1194                  end if
1195                  term5=-two_pi*(rhog(re,ii)*phimag_elt(ig1,ig2,ig3,ia)&
1196 &                 +rhog(im,ii)*phre_elt(ig1,ig2,ig3,ia))*term4
1197                  eltfrloc(7+3*(ia-1),is2)=eltfrloc(7+3*(ia-1),is2)+term5*dble(ig1)
1198                  eltfrloc(8+3*(ia-1),is2)=eltfrloc(8+3*(ia-1),is2)+term5*dble(ig2)
1199                  eltfrloc(9+3*(ia-1),is2)=eltfrloc(9+3*(ia-1),is2)+term5*dble(ig3)
1200                end do
1201 
1202              end do !is2
1203 
1204 !            End skip G**2 outside cutoff:
1205            end if
1206 
1207 !          End loop on n1, n2, n3. There is a "cycle" inside the loop
1208          end do
1209        end if
1210      end do
1211    end do
1212 
1213 !  End loop on type of atoms
1214    ia1=ia2+1
1215  end do
1216 !Init mpi_comm
1217  call timab(48,1,tsec)
1218  call xmpi_sum(eltfrloc,mpi_enreg%comm_fft,ierr)
1219  call timab(48,2,tsec)
1220 
1221 !Fill in lower triangle
1222  do is2=2,6
1223    do is1=1,is2-1
1224      eltfrloc(is2,is1)=eltfrloc(is1,is2)
1225    end do
1226  end do
1227 
1228 !The indexing array atindx is used to reestablish the correct
1229 !order of atoms
1230  ABI_MALLOC(elt_work,(6+3*natom,6))
1231  elt_work(1:6,1:6)=eltfrloc(1:6,1:6)
1232  do ia=1,natom
1233    ielt=7+3*(ia-1)
1234    ieltx=7+3*(atindx(ia)-1)
1235    elt_work(ielt:ielt+2,1:6)=eltfrloc(ieltx:ieltx+2,1:6)
1236  end do
1237  eltfrloc(:,:)=elt_work(:,:)
1238 
1239  ABI_FREE(d2gm)
1240  ABI_FREE(elt_work)
1241 
1242  contains
1243 
1244 !Real and imaginary parts of phase.
1245    function phr_elt(x1,y1,x2,y2,x3,y3)
1246 
1247    real(dp) :: phr_elt,x1,x2,x3,y1,y2,y3
1248    phr_elt=(x1*x2-y1*y2)*x3-(y1*x2+x1*y2)*y3
1249  end function phr_elt
1250 
1251    function phi_elt(x1,y1,x2,y2,x3,y3)
1252 
1253    real(dp):: phi_elt,x1,x2,x3,y1,y2,y3
1254    phi_elt=(x1*x2-y1*y2)*y3+(y1*x2+x1*y2)*x3
1255  end function phi_elt
1256 
1257    function ph1_elt(nri,ig1,ia)
1258 
1259    real(dp):: ph1_elt
1260    integer :: nri,ig1,ia
1261    ph1_elt=ph1d(nri,ig1+1+n1+(ia-1)*(2*n1+1))
1262  end function ph1_elt
1263 
1264    function ph2_elt(nri,ig2,ia)
1265 
1266    real(dp):: ph2_elt
1267    integer :: nri,ig2,ia
1268    ph2_elt=ph1d(nri,ig2+1+n2+(ia-1)*(2*n2+1)+natom*(2*n1+1))
1269  end function ph2_elt
1270 
1271    function ph3_elt(nri,ig3,ia)
1272 
1273    real(dp):: ph3_elt
1274    integer :: nri,ig3,ia
1275    ph3_elt=ph1d(nri,ig3+1+n3+(ia-1)*(2*n3+1)+natom*(2*n1+1+2*n2+1))
1276  end function ph3_elt
1277 
1278    function phre_elt(ig1,ig2,ig3,ia)
1279 
1280    real(dp):: phre_elt
1281    integer :: ig1,ig2,ig3,ia
1282    phre_elt=phr_elt(ph1_elt(re,ig1,ia),ph1_elt(im,ig1,ia),&
1283 &   ph2_elt(re,ig2,ia),ph2_elt(im,ig2,ia),ph3_elt(re,ig3,ia),ph3_elt(im,ig3,ia))
1284  end function phre_elt
1285 
1286    function phimag_elt(ig1,ig2,ig3,ia)
1287 
1288    real(dp) :: phimag_elt
1289    integer :: ig1,ig2,ig3,ia
1290    phimag_elt=phi_elt(ph1_elt(re,ig1,ia),ph1_elt(im,ig1,ia),&
1291 &   ph2_elt(re,ig2,ia),ph2_elt(im,ig2,ia),ph3_elt(re,ig3,ia),ph3_elt(im,ig3,ia))
1292  end function phimag_elt
1293 
1294    function gsq_elt(i1,i2,i3)
1295 
1296    real(dp) :: gsq_elt
1297    integer :: i1,i2,i3
1298 !Define G^2 based on G space metric gmet.
1299    gsq_elt=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+&
1300 &   dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+&
1301 &   dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
1302  end function gsq_elt
1303 
1304    function dgsqds_elt(i1,i2,i3,is)
1305 
1306    real(dp) :: dgsqds_elt
1307    integer :: i1,i2,i3,is
1308 !Define dG^2/ds based on G space metric derivative
1309    dgsqds_elt=dble(i1*i1)*dgm(1,1,is)+dble(i2*i2)*dgm(2,2,is)+&
1310 &   dble(i3*i3)*dgm(3,3,is)+&
1311 &   dble(i1*i2)*(dgm(1,2,is)+dgm(2,1,is))+&
1312 &   dble(i1*i3)*(dgm(1,3,is)+dgm(3,1,is))+&
1313 &   dble(i2*i3)*(dgm(2,3,is)+dgm(3,2,is))
1314  end function dgsqds_elt
1315 
1316    function d2gsqds_elt(i1,i2,i3,is1,is2)
1317 
1318    real(dp) :: d2gsqds_elt
1319    integer :: i1,i2,i3,is1,is2
1320 !Define 2dG^2/ds1ds2  based on G space metric derivative
1321    d2gsqds_elt=dble(i1*i1)*d2gm(1,1,is1,is2)+&
1322 &   dble(i2*i2)*d2gm(2,2,is1,is2)+dble(i3*i3)*d2gm(3,3,is1,is2)+&
1323 &   dble(i1*i2)*(d2gm(1,2,is1,is2)+d2gm(2,1,is1,is2))+&
1324 &   dble(i1*i3)*(d2gm(1,3,is1,is2)+d2gm(3,1,is1,is2))+&
1325 &   dble(i2*i3)*(d2gm(2,3,is1,is2)+d2gm(3,2,is1,is2))
1326  end function d2gsqds_elt
1327 
1328 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.

SOURCE

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

SOURCE

2337 subroutine dfpt_ewald(dyew,gmet,my_natom,natom,qphon,rmet,sumg0,typat,ucvol,xred,zion, &
2338 &                 mpi_atmtab,comm_atom ) ! optional arguments (parallelism))
2339 
2340 !Arguments -------------------------------
2341 !scalars
2342  integer,intent(in) :: my_natom,natom,sumg0
2343  real(dp),intent(in) :: ucvol
2344 !arrays
2345  integer,intent(in) :: typat(natom)
2346  integer,optional,intent(in) :: comm_atom
2347  integer,optional,target,intent(in) :: mpi_atmtab(:)
2348  real(dp),intent(in) :: gmet(3,3),qphon(3),rmet(3,3),xred(3,natom),zion(*)
2349  real(dp),intent(out) :: dyew(2,3,natom,3,natom)
2350 
2351 !Local variables -------------------------
2352 !nr, ng affect convergence of sums (nr=3,ng=5 is not good enough):
2353 !scalars
2354  integer,parameter :: im=2,ng=10,nr=6,re=1
2355  integer :: ia,ia0,ib,ierr,ig1,ig2,ig3,ii,ir1,ir2,ir3,mu,my_comm_atom,nu
2356  logical :: my_atmtab_allocated,paral_atom
2357  real(dp) :: arg,arga,argb,c1i,c1r,da1,da2,da3,derfc_arg
2358  real(dp) :: direct,dot1,dot2,dot3,dotr1,dotr2,dotr3
2359  real(dp) :: eta,fac,gdot12,gdot13,gdot23,gsq,gsum,norm1
2360  real(dp) :: r1,r2,r3,rdot12,rdot13,rdot23,recip,reta
2361  real(dp) :: reta3m,rmagn,rsq,term,term1,term2
2362  real(dp) :: term3
2363  character(len=500) :: message
2364 !arrays
2365  real(dp) :: tsec(2)
2366  integer,pointer :: my_atmtab(:)
2367  real(dp) :: gpq(3),rq(3)
2368 
2369 ! *************************************************************************
2370 
2371 !Set up parallelism over atoms
2372  paral_atom=(present(comm_atom).and.(my_natom/=natom))
2373  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2374  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2375  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
2376 
2377 !Compute eta for approximately optimized summations:
2378  direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
2379 & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
2380  recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
2381 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
2382  eta=pi*(dble(ng)/dble(nr))*sqrt(1.69_dp*recip/direct)
2383 
2384 !Test Ewald s summation
2385 !eta=1.2_dp*eta
2386 
2387 !Sum terms over g space:
2388  fac=pi**2/eta
2389  gsum=zero
2390  da1=zero
2391  da2=zero
2392  da3=zero
2393  dyew(:,:,:,:,:)=zero
2394  do ig3=-ng,ng
2395    do ig2=-ng,ng
2396      do ig1=-ng,ng
2397        gpq(1)=dble(ig1)+qphon(1)
2398        gpq(2)=dble(ig2)+qphon(2)
2399        gpq(3)=dble(ig3)+qphon(3)
2400        gdot12=gmet(2,1)*gpq(1)*gpq(2)
2401        gdot13=gmet(3,1)*gpq(1)*gpq(3)
2402        gdot23=gmet(3,2)*gpq(2)*gpq(3)
2403        dot1=gmet(1,1)*gpq(1)**2+gdot12+gdot13
2404        dot2=gmet(2,2)*gpq(2)**2+gdot12+gdot23
2405        dot3=gmet(3,3)*gpq(3)**2+gdot13+gdot23
2406        gsq=dot1+dot2+dot3
2407 !      Skip q=0:
2408        if (gsq<1.0d-20) then
2409          if (sumg0==1) then
2410            write(message,'(5a)')&
2411 &           'The phonon wavelength should not be zero : ',ch10,&
2412 &           'there are non-analytical terms that the code cannot handle.',ch10,&
2413 &           'Action : subtract this wavelength from the input.'
2414            ABI_ERROR(message)
2415          end if
2416        else
2417          arg=fac*gsq
2418 !        Larger arg gives 0 contribution:
2419          if (arg <= 80._dp) then
2420            term=exp(-arg)/gsq
2421            do ia0=1,my_natom
2422              ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2423              arga=two_pi*(gpq(1)*xred(1,ia)+gpq(2)*xred(2,ia)+gpq(3)*xred(3,ia))
2424              do ib=1,ia
2425                argb=two_pi*(gpq(1)*xred(1,ib)+gpq(2)*xred(2,ib)+gpq(3)*xred(3,ib))
2426                arg=arga-argb
2427                c1r=cos(arg)*term
2428                c1i=sin(arg)*term
2429 
2430                do mu=1,3
2431                  do nu=1,mu
2432                    dyew(re,mu,ia,nu,ib)=dyew(re,mu,ia,nu,ib)+gpq(mu)*gpq(nu)*c1r
2433                    dyew(im,mu,ia,nu,ib)=dyew(im,mu,ia,nu,ib)+gpq(mu)*gpq(nu)*c1i
2434                  end do
2435                end do
2436 
2437              end do
2438            end do
2439          end if
2440 !        Endif g/=0 :
2441        end if
2442 !      End triple loop over G s:
2443      end do
2444    end do
2445  end do
2446 
2447 !End G summation by accounting for some common factors.
2448 !(for the charges:see end of routine)
2449  norm1=4.0_dp*pi/ucvol
2450  do ia0=1,my_natom
2451    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2452    do ib=1,ia
2453      do mu=1,3
2454        do nu=1,mu
2455          dyew(:,mu,ia,nu,ib)=dyew(:,mu,ia,nu,ib)*norm1
2456        end do
2457      end do
2458    end do
2459  end do
2460 
2461 !Do sums over real space:
2462  reta=sqrt(eta)
2463  reta3m=-eta*reta
2464  fac=4._dp/3.0_dp/sqrt(pi)
2465  do ir3=-nr,nr
2466    do ir2=-nr,nr
2467      do ir1=-nr,nr
2468        arg=-two_pi*(qphon(1)*ir1+qphon(2)*ir2+qphon(3)*ir3)
2469        c1r=cos(arg)*reta3m
2470        c1i=sin(arg)*reta3m
2471        do ia0=1,my_natom
2472          ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2473          do ib=1,ia
2474            r1=dble(ir1)+xred(1,ia)-xred(1,ib)
2475            r2=dble(ir2)+xred(2,ia)-xred(2,ib)
2476            r3=dble(ir3)+xred(3,ia)-xred(3,ib)
2477            rdot12=rmet(2,1)*r1*r2
2478            rdot13=rmet(3,1)*r1*r3
2479            rdot23=rmet(3,2)*r2*r3
2480            dotr1=rmet(1,1)*r1**2+rdot12+rdot13
2481            dotr2=rmet(2,2)*r2**2+rdot12+rdot23
2482            dotr3=rmet(3,3)*r3**2+rdot13+rdot23
2483            rsq=dotr1+dotr2+dotr3
2484            rmagn=sqrt(rsq)
2485 !          Avoid zero denominators in term :
2486            if (rmagn>=1.0d-12) then
2487              arg=reta*rmagn
2488              term=zero
2489              if (arg<8.0_dp) then
2490 !              Note: erfc(8) is about 1.1e-29,
2491 !              so don t bother with larger arg.
2492 !              Also: exp(-64) is about 1.6e-28,
2493 !              so don t bother with larger arg**2 in exp.
2494                derfc_arg = abi_derfc(arg)
2495                term=derfc_arg/arg**3
2496                term1=2.0_dp/sqrt(pi)*exp(-arg**2)/arg**2
2497                term2=-(term+term1)
2498                term3=(3*term+term1*(3.0_dp+2.0_dp*arg**2))/rsq
2499                rq(1)=rmet(1,1)*r1+rmet(1,2)*r2+rmet(1,3)*r3
2500                rq(2)=rmet(2,1)*r1+rmet(2,2)*r2+rmet(2,3)*r3
2501                rq(3)=rmet(3,1)*r1+rmet(3,2)*r2+rmet(3,3)*r3
2502                do mu=1,3
2503                  do nu=1,mu
2504                    dyew(re,mu,ia,nu,ib)=dyew(re,mu,ia,nu,ib)+&
2505 &                   c1r*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
2506                    dyew(im,mu,ia,nu,ib)=dyew(im,mu,ia,nu,ib)+&
2507 &                   c1i*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
2508                  end do
2509                end do
2510              end if
2511            else
2512              if (ia/=ib)then
2513                write(message,'(a,a,a,a,a,i5,a,i5,a)')&
2514 &               'The distance between two atoms vanishes.',ch10,&
2515 &               'This is not allowed.',ch10,&
2516 &               'Action: check the input for the atoms number',ia,' and',ib,'.'
2517                ABI_ERROR(message)
2518              else
2519                do mu=1,3
2520                  do nu=1,mu
2521                    dyew(re,mu,ia,nu,ib)=dyew(re,mu,ia,nu,ib)+&
2522 &                   fac*reta3m*rmet(mu,nu)
2523                  end do
2524                end do
2525              end if
2526            end if
2527 
2528          end do ! End loop over ib:
2529        end do ! End loop over ia:
2530      end do ! End triple loop over real space points:
2531    end do
2532  end do
2533 
2534 !Take account of the charges
2535 !write(std_out,*)' '
2536  do ia0=1,my_natom
2537    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2538    do ib=1,ia
2539      do mu=1,3
2540        do nu=1,mu
2541          do ii=1,2
2542 !          write(std_out,*)dyew(ii,mu,ia,nu,ib)
2543            dyew(ii,mu,ia,nu,ib)=dyew(ii,mu,ia,nu,ib)*&
2544 &           zion(typat(ia))*zion(typat(ib))
2545          end do
2546        end do
2547      end do
2548    end do
2549  end do
2550 
2551 !Symmetrize with respect to the directions
2552  do ia0=1,my_natom
2553    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2554    do ib=1,ia
2555      do mu=1,3
2556        do nu=1,mu
2557          dyew(re,nu,ia,mu,ib)=dyew(re,mu,ia,nu,ib)
2558          dyew(im,nu,ia,mu,ib)=dyew(im,mu,ia,nu,ib)
2559        end do
2560      end do
2561    end do
2562  end do
2563 
2564 !In case of parallelism over atoms: communicate
2565  if (paral_atom) then
2566    call timab(48,1,tsec)
2567    call xmpi_sum(dyew,my_comm_atom,ierr)
2568    call timab(48,2,tsec)
2569  end if
2570 
2571 !Fill the upper part of the matrix, with the hermitian conjugate
2572  do ia=1,natom
2573    do ib=1,ia
2574      do nu=1,3
2575        do mu=1,3
2576          dyew(re,mu,ib,nu,ia)=dyew(re,mu,ia,nu,ib)
2577          dyew(im,mu,ib,nu,ia)=-dyew(im,mu,ia,nu,ib)
2578        end do
2579      end do
2580    end do
2581  end do
2582 
2583 !Destroy atom table used for parallelism
2584  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2585 
2586 end subroutine dfpt_ewald

ABINIT/dfpt_ewalddq [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_ewalddq

FUNCTION

 Compute the first q-gradient of Ewald contribution to the dynamical matrix, at a given q wavevector.
 If q=0 is asked, sumg0 should be put to 0. Otherwise, it should be put to 1.

COPYRIGHT

 Copyright (C) 1998-2022 ABINIT group (MR, MS)
  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

 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

 dyewdq(2,3,natom,3,natom,3)= First q-gradient of Ewald part of the dynamical matrix,
    second energy derivative wrt xred(3,natom), Hartrees.

SOURCE

2626 subroutine dfpt_ewalddq(dyewdq,gmet,my_natom,natom,qphon,rmet,sumg0,typat,ucvol,xred,zion, &
2627 &                 mpi_atmtab,comm_atom ) ! optional arguments (parallelism))
2628 
2629 !Arguments -------------------------------
2630 !scalars
2631  integer,intent(in) :: my_natom,natom,sumg0
2632  real(dp),intent(in) :: ucvol
2633 !arrays
2634  integer,intent(in) :: typat(natom)
2635  integer,optional,intent(in) :: comm_atom
2636  integer,optional,target,intent(in) :: mpi_atmtab(:)
2637  real(dp),intent(in) :: gmet(3,3),qphon(3),rmet(3,3),xred(3,natom),zion(*)
2638  real(dp),intent(out) :: dyewdq(2,3,natom,3,natom,3)
2639 
2640 !Local variables -------------------------
2641 !nr, ng affect convergence of sums (nr=3,ng=5 is not good enough):
2642 !scalars
2643  integer,parameter :: im=2,nng=10,nnr=6,re=1
2644  integer ::ia,ia0,ib,ierr,ig1,ig2,ig3,ii,iq,ir1,ir2,ir3,mu,my_comm_atom,newg,newr,ng,nr,nu
2645  logical :: my_atmtab_allocated,paral_atom
2646  real(dp) :: arg,arga,argb,c1i,c1r,delag,delbg,derfc_arg
2647  real(dp) :: direct,dot1,dot2,dot3,dotr1,dotr2,dotr3
2648  real(dp) :: eta,fac,fac2,gdot12,gdot13,gdot23,gsq,gpqdq,gterms,norm1
2649  real(dp) :: r1,r2,r3,rdot12,rdot13,rdot23,recip,reta
2650  real(dp) :: reta3m,rmagn,rsq,term,term1,term2,term3
2651  character(len=500) :: message
2652 !arrays
2653  real(dp) :: tsec(2)
2654  integer,pointer :: my_atmtab(:)
2655  real(dp) :: dakk(3),gpq(3),rq(3)
2656 
2657 ! *************************************************************************
2658 
2659 !Set up parallelism over atoms
2660  paral_atom=(present(comm_atom).and.(my_natom/=natom))
2661  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2662  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2663  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
2664 
2665 !Compute eta for approximately optimized summations:
2666  direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
2667 & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
2668  recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
2669 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
2670  eta=pi*(dble(nng)/dble(nnr))*sqrt(1.69_dp*recip/direct)
2671 
2672 !Test Ewald s summation
2673 !eta=1.2_dp*eta
2674 
2675 !Sum over G space, done shell after shell until all
2676 !contributions are too small.
2677  fac=pi**2.d0/eta
2678  fac2=2.d0*fac
2679  dyewdq(:,:,:,:,:,:)=zero
2680  ng=0
2681  do
2682    ng=ng+1
2683    newg=0
2684 
2685    do ig3=-ng,ng
2686      do ig2=-ng,ng
2687        do ig1=-ng,ng
2688 
2689 !        Exclude shells previously summed over
2690          if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng&
2691 &         .or. ng==1 ) then
2692 
2693            gpq(1)=dble(ig1)+qphon(1)
2694            gpq(2)=dble(ig2)+qphon(2)
2695            gpq(3)=dble(ig3)+qphon(3)
2696            gdot12=gmet(2,1)*gpq(1)*gpq(2)
2697            gdot13=gmet(3,1)*gpq(1)*gpq(3)
2698            gdot23=gmet(3,2)*gpq(2)*gpq(3)
2699            dot1=gmet(1,1)*gpq(1)**2+gdot12+gdot13
2700            dot2=gmet(2,2)*gpq(2)**2+gdot12+gdot23
2701            dot3=gmet(3,3)*gpq(3)**2+gdot13+gdot23
2702            gsq=dot1+dot2+dot3
2703   !        Skip q=0:
2704            if (gsq<1.0d-20) then
2705              if (sumg0==1) then
2706                write(message,'(3a)')&
2707   &             'The G=0 term has no contributions at first order in q: ',ch10,&
2708   &             'Action : sumg0=0 '
2709                ABI_ERROR(message)
2710              end if
2711            else
2712              arg=fac*gsq
2713   !          Larger arg gives 0 contribution:
2714              if (arg <= 80._dp) then
2715 !              When any term contributes then include next shell
2716                newg=1
2717                term=exp(-arg)/gsq
2718                do ia0=1,my_natom
2719                  ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2720                  arga=two_pi*(gpq(1)*xred(1,ia)+gpq(2)*xred(2,ia)+gpq(3)*xred(3,ia))
2721                  do ib=1,ia
2722                    argb=two_pi*(gpq(1)*xred(1,ib)+gpq(2)*xred(2,ib)+gpq(3)*xred(3,ib))
2723                    arg=arga-argb
2724                    c1r=cos(arg)*term
2725                    c1i=sin(arg)*term
2726 
2727                    do iq=1,3
2728                      gpqdq=gmet(iq,1)*gpq(1)+gmet(iq,2)*gpq(2)+gmet(iq,3)*gpq(3)
2729                      do mu=1,3
2730                        delag=zero; if(iq==mu) delag=one
2731                        do nu=1,mu
2732                          delbg=zero; if(iq==nu) delbg=one
2733                          term1=delag*gpq(nu)+delbg*gpq(mu)
2734                          term2=gpq(mu)*gpq(nu)*gpqdq
2735                          term3=fac2*term2
2736                          term2=two*term2/gsq
2737                          gterms=term1-term2-term3
2738                          dyewdq(re,mu,ia,nu,ib,iq)=dyewdq(re,mu,ia,nu,ib,iq)+gterms*c1r
2739                          dyewdq(im,mu,ia,nu,ib,iq)=dyewdq(im,mu,ia,nu,ib,iq)+gterms*c1i
2740                        end do
2741                      end do
2742                    end do
2743 
2744                  end do
2745                end do
2746              end if
2747   !          Endif g/=0 :
2748            end if
2749          end if
2750   !        End triple loop over G s:
2751        end do
2752      end do
2753    end do
2754 
2755 !  Check if new shell must be calculated
2756    if (newg==0) exit
2757  end do !  End the loop on ng (new shells). Note that there is one exit from this loop.
2758 
2759 !End G summation by accounting for some common factors.
2760 !(for the charges:see end of routine)
2761  norm1=4.0_dp*pi/ucvol
2762  do ia0=1,my_natom
2763    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2764    do ib=1,ia
2765      do iq=1,3
2766        do mu=1,3
2767          do nu=1,mu
2768            dyewdq(:,mu,ia,nu,ib,iq)=dyewdq(:,mu,ia,nu,ib,iq)*norm1
2769          end do
2770        end do
2771      end do
2772    end do
2773  end do
2774 
2775 !Do sums over real space:
2776  reta=sqrt(eta)
2777  reta3m=-eta*reta
2778  fac=4._dp/3.0_dp/sqrt(pi)
2779  nr=0
2780  do
2781    nr=nr+1
2782    newr=0
2783 
2784    do ir3=-nr,nr
2785      do ir2=-nr,nr
2786        do ir1=-nr,nr
2787          if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
2788 
2789            arg=two_pi*(qphon(1)*ir1+qphon(2)*ir2+qphon(3)*ir3)
2790            c1r=cos(arg)*reta3m
2791            c1i=sin(arg)*reta3m
2792            do ia0=1,my_natom
2793              ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2794              do ib=1,ia
2795                r1=dble(ir1)+xred(1,ib)-xred(1,ia)
2796                r2=dble(ir2)+xred(2,ib)-xred(2,ia)
2797                r3=dble(ir3)+xred(3,ib)-xred(3,ia)
2798                dakk(:)=two_pi*(/r1,r2,r3/)
2799                rdot12=rmet(2,1)*r1*r2
2800                rdot13=rmet(3,1)*r1*r3
2801                rdot23=rmet(3,2)*r2*r3
2802                dotr1=rmet(1,1)*r1**2+rdot12+rdot13
2803                dotr2=rmet(2,2)*r2**2+rdot12+rdot23
2804                dotr3=rmet(3,3)*r3**2+rdot13+rdot23
2805                rsq=dotr1+dotr2+dotr3
2806                rmagn=sqrt(rsq)
2807 !              Avoid zero denominators in term :
2808                if (rmagn>=1.0d-12) then
2809                  arg=reta*rmagn
2810                  term=zero
2811                  if (arg<8.0_dp) then
2812 !                  Note: erfc(8) is about 1.1e-29,
2813 !                  so don t bother with larger arg.
2814 !                  Also: exp(-64) is about 1.6e-28,
2815 !                  so don t bother with larger arg**2 in exp.
2816                    newr=1
2817                    derfc_arg = abi_derfc(arg)
2818                    term=derfc_arg/arg**3
2819                    term1=2.0_dp/sqrt(pi)*exp(-arg**2)/arg**2
2820                    term2=-(term+term1)
2821                    term3=(3*term+term1*(3.0_dp+2.0_dp*arg**2))/rsq
2822                    rq(1)=rmet(1,1)*r1+rmet(1,2)*r2+rmet(1,3)*r3
2823                    rq(2)=rmet(2,1)*r1+rmet(2,2)*r2+rmet(2,3)*r3
2824                    rq(3)=rmet(3,1)*r1+rmet(3,2)*r2+rmet(3,3)*r3
2825                    do iq=1,3
2826                      do mu=1,3
2827 !                       do nu=1,3
2828                        do nu=1,mu
2829                          dyewdq(re,mu,ia,nu,ib,iq)=dyewdq(re,mu,ia,nu,ib,iq)-&
2830 &                         c1i*dakk(iq)*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
2831                          dyewdq(im,mu,ia,nu,ib,iq)=dyewdq(im,mu,ia,nu,ib,iq)+&
2832 &                         c1r*dakk(iq)*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
2833                        end do
2834                      end do
2835                    end do
2836                  end if
2837                else
2838                  if (ia/=ib)then
2839                    write(message,'(a,a,a,a,a,i5,a,i5,a)')&
2840 &                   'The distance between two atoms vanishes.',ch10,&
2841 &                   'This is not allowed.',ch10,&
2842 &                   'Action: check the input for the atoms number',ia,' and',ib,'.'
2843                    ABI_ERROR(message)
2844                  end if
2845                end if
2846 
2847              end do ! End loop over ib:
2848            end do ! End loop over ia:
2849          end if
2850        end do ! End triple loop over real space points:
2851      end do
2852    end do
2853 
2854 !  Check if new shell must be calculated
2855    if(newr==0) exit
2856  end do !  End loop on nr (new shells). Note that there is an exit within the loop
2857 
2858 !Take account of the charges
2859 !write(std_out,*)' '
2860  do ia0=1,my_natom
2861    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2862    do ib=1,ia
2863      do iq=1,3
2864        do mu=1,3
2865          do nu=1,mu
2866            do ii=1,2
2867              dyewdq(ii,mu,ia,nu,ib,iq)=dyewdq(ii,mu,ia,nu,ib,iq)*&
2868 &             zion(typat(ia))*zion(typat(ib))
2869            end do
2870          end do
2871        end do
2872      end do
2873    end do
2874  end do
2875 
2876 !Symmetrize with respect to the directions
2877  do ia0=1,my_natom
2878    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2879    do ib=1,ia
2880      do iq=1,3
2881        do mu=1,3
2882          do nu=1,mu
2883            dyewdq(re,nu,ia,mu,ib,iq)=dyewdq(re,mu,ia,nu,ib,iq)
2884            dyewdq(im,nu,ia,mu,ib,iq)=dyewdq(im,mu,ia,nu,ib,iq)
2885          end do
2886        end do
2887      end do
2888    end do
2889  end do
2890 
2891 !In case of parallelism over atoms: communicate
2892  if (paral_atom) then
2893    call timab(48,1,tsec)
2894    call xmpi_sum(dyewdq,my_comm_atom,ierr)
2895    call timab(48,2,tsec)
2896  end if
2897 
2898 !Fill the upper part of the matrix, with the hermitian conjugate
2899  do ia=1,natom
2900    do ib=1,ia
2901      do iq=1,3
2902        do nu=1,3
2903          do mu=1,3
2904            dyewdq(re,mu,ib,nu,ia,iq)=dyewdq(re,mu,ia,nu,ib,iq)
2905            dyewdq(im,mu,ib,nu,ia,iq)=-dyewdq(im,mu,ia,nu,ib,iq)
2906          end do
2907        end do
2908      end do
2909    end do
2910  end do
2911 
2912 !Destroy atom table used for parallelism
2913  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2914 
2915 end subroutine dfpt_ewalddq

ABINIT/dfpt_ewalddqdq [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_ewalddqdq

FUNCTION

 Compute the second q-gradient of Ewald contribution to the dynamical matrix, at a given q wavevector,
 sumed over the second atomic sublattice.
 If q=0 is asked, sumg0 should be put to 0. Otherwise, it should be put to 1.

COPYRIGHT

 Copyright (C) 1998-2022 ABINIT group (MR, MS)
  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

 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

 dyewdqdq(2,3,natom,3,3,3)= First q-gradient of Ewald part of the dynamical matrix,
    sumed over second atomic sublattice.

SOURCE

2956 subroutine dfpt_ewalddqdq(dyewdqdq,gmet,my_natom,natom,qphon,rmet,sumg0,typat,ucvol,xred,zion, &
2957 &                 mpi_atmtab,comm_atom ) ! optional arguments (parallelism))
2958 
2959 !Arguments -------------------------------
2960 !scalars
2961  integer,intent(in) :: my_natom,natom,sumg0
2962  real(dp),intent(in) :: ucvol
2963 !arrays
2964  integer,intent(in) :: typat(natom)
2965  integer,optional,intent(in) :: comm_atom
2966  integer,optional,target,intent(in) :: mpi_atmtab(:)
2967  real(dp),intent(in) :: gmet(3,3),qphon(3),rmet(3,3),xred(3,natom),zion(*)
2968  real(dp),intent(out) :: dyewdqdq(2,3,natom,3,3,3)
2969 
2970 !Local variables -------------------------
2971 !nr, ng affect convergence of sums (nr=3,ng=5 is not good enough):
2972 !scalars
2973  integer,parameter :: im=2,nng=10,nnr=6,re=1
2974  integer :: ia,ia0,ib,ierr,ig1,ig2,ig3,ii,iq1,iq2,ir1,ir2,ir3,mu,my_comm_atom,newg,newr,ng,nr,nu
2975  logical :: my_atmtab_allocated,paral_atom
2976  real(dp) :: arg,arga,argb,c1i,c1r,delad,delag,delbd,delbg,derfc_arg
2977  real(dp) :: direct,dot1,dot2,dot3,dotr1,dotr2,dotr3
2978  real(dp) :: eta,fac,fac2,fac8,fac2sqr,gdot12,gdot13,gdot23,gsq,gsqsq,gpqdq1,gpqdq2,gterms,g0term,norm1
2979  real(dp) :: r1,r2,r3,rdot12,rdot13,rdot23,recip,reta
2980  real(dp) :: reta3m,rmagn,rsq,term,term1,term2,term3
2981  character(len=500) :: message
2982 !arrays
2983  integer,pointer :: my_atmtab(:)
2984  real(dp) :: dakk(3),gpq(3),rq(3)
2985  real(dp) :: tsec(2)
2986  real(dp),allocatable :: work(:,:,:,:,:,:,:)
2987 
2988 ! *************************************************************************
2989 
2990 !Set up parallelism over atoms
2991  paral_atom=(present(comm_atom).and.(my_natom/=natom))
2992  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2993  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
2994  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
2995 
2996 !Compute eta for approximately optimized summations:
2997  direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
2998 & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
2999  recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
3000 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
3001  eta=pi*(dble(nng)/dble(nnr))*sqrt(1.69_dp*recip/direct)
3002 ! eta=1.0_dp
3003 
3004 !Test Ewald s summation
3005 !eta=1.2_dp*eta
3006 
3007 !Sum terms over g space:
3008  fac=pi**2.0_dp/eta
3009  fac2=2.0_dp*fac
3010  fac8=4.0_dp*fac2
3011  fac2sqr=fac2*fac2
3012  ABI_MALLOC(work,(2,3,natom,3,natom,3,3))
3013  work(:,:,:,:,:,:,:)=zero
3014  ng=0
3015  do
3016    ng=ng+1
3017    newg=0
3018 
3019    do ig3=-ng,ng
3020      do ig2=-ng,ng
3021        do ig1=-ng,ng
3022 
3023 !        Exclude shells previously summed over
3024          if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng&
3025 &         .or. ng==1 ) then
3026 
3027            gpq(1)=dble(ig1)+qphon(1)
3028            gpq(2)=dble(ig2)+qphon(2)
3029            gpq(3)=dble(ig3)+qphon(3)
3030            gdot12=gmet(2,1)*gpq(1)*gpq(2)
3031            gdot13=gmet(3,1)*gpq(1)*gpq(3)
3032            gdot23=gmet(3,2)*gpq(2)*gpq(3)
3033            dot1=gmet(1,1)*gpq(1)**2+gdot12+gdot13
3034            dot2=gmet(2,2)*gpq(2)**2+gdot12+gdot23
3035            dot3=gmet(3,3)*gpq(3)**2+gdot13+gdot23
3036            gsq=dot1+dot2+dot3
3037            gsqsq=gsq*gsq
3038 !          Skip q=0:
3039            if (gsq<1.0d-20) then
3040 
3041 !            At second order in q there is a nonvanishing G=0 contribution in the longwave limit
3042              if (sumg0==1) then
3043                do ia0=1,my_natom
3044                  ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
3045                  do ib=1,ia
3046                    do iq2=1,3
3047                      do iq1=1,3
3048                        do mu=1,3
3049                          delag=zero; if(iq1==mu) delag=one
3050                          delad=zero; if(iq2==mu) delad=one
3051                          do nu=1,mu
3052                            delbg=zero; if(iq1==nu) delbg=one
3053                            delbd=zero; if(iq2==nu) delbd=one
3054                            g0term=-fac*(delad*delbg+delbd*delag)
3055                            work(re,mu,ia,nu,ib,iq1,iq2)=work(re,mu,ia,nu,ib,iq1,iq2)+g0term
3056                          end do
3057                        end do
3058                      end do
3059                    end do
3060                  end do
3061                end do
3062              end if
3063 
3064            else
3065              arg=fac*gsq
3066 !            Larger arg gives 0 contribution:
3067              if (arg <= 80._dp) then
3068                newg=1
3069                term=exp(-arg)/gsq
3070                do ia0=1,my_natom
3071                  ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
3072                  arga=two_pi*(gpq(1)*xred(1,ia)+gpq(2)*xred(2,ia)+gpq(3)*xred(3,ia))
3073                  do ib=1,ia
3074                    argb=two_pi*(gpq(1)*xred(1,ib)+gpq(2)*xred(2,ib)+gpq(3)*xred(3,ib))
3075                    arg=arga-argb
3076                    c1r=cos(arg)*term
3077                    c1i=sin(arg)*term
3078 
3079                    do iq2=1,3
3080                      gpqdq2=gmet(iq2,1)*gpq(1)+gmet(iq2,2)*gpq(2)+gmet(iq2,3)*gpq(3)
3081                      do iq1=1,3
3082                        gpqdq1=gmet(iq1,1)*gpq(1)+gmet(iq1,2)*gpq(2)+gmet(iq1,3)*gpq(3)
3083                        do mu=1,3
3084                          delag=zero; if(iq1==mu) delag=one
3085                          delad=zero; if(iq2==mu) delad=one
3086                          do nu=1,mu
3087                            delbg=zero; if(iq1==nu) delbg=one
3088                            delbd=zero; if(iq2==nu) delbd=one
3089 
3090                            term1=gpqdq2*(delag*gpq(nu)+delbg*gpq(mu))
3091                            term1=term1+gpqdq1*(delad*gpq(nu)+delbd*gpq(mu))
3092                            term1=term1+gpq(mu)*gpq(nu)*gmet(iq1,iq2)
3093                            term1=-term1*(fac2+2.0_dp/gsq)
3094 
3095                            term2=delag*delbd + delbg*delad
3096 
3097                            term3=gpqdq1*gpqdq2*gpq(mu)*gpq(nu)
3098                            term3=term3*(fac8/gsq + fac2sqr + 8.0_dp/gsqsq)
3099 
3100                            gterms=term1+term2+term3
3101                            work(re,mu,ia,nu,ib,iq1,iq2)=work(re,mu,ia,nu,ib,iq1,iq2)+gterms*c1r
3102                            work(im,mu,ia,nu,ib,iq1,iq2)=work(im,mu,ia,nu,ib,iq1,iq2)+gterms*c1i
3103                          end do
3104                        end do
3105                      end do
3106                    end do
3107                  end do
3108                end do
3109              end if
3110 !            Endif g/=0 :
3111            end if
3112          end if
3113 !        End triple loop over G s:
3114        end do
3115      end do
3116    end do
3117 
3118 !  Check if new shell must be calculated
3119    if (newg==0) exit
3120  end do !  End the loop on ng (new shells). Note that there is one exit from this loop.
3121 
3122 !End G summation by accounting for some common factors.
3123 !(for the charges:see end of routine)
3124  norm1=4.0_dp*pi/ucvol
3125  do ia0=1,my_natom
3126    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
3127    do ib=1,ia
3128      do iq2=1,3
3129        do iq1=1,3
3130          do mu=1,3
3131            do nu=1,mu
3132              work(:,mu,ia,nu,ib,iq1,iq2)=work(:,mu,ia,nu,ib,iq1,iq2)*norm1
3133            end do
3134          end do
3135        end do
3136      end do
3137    end do
3138  end do
3139 
3140 !Do sums over real space:
3141  reta=sqrt(eta)
3142  reta3m=eta*reta
3143  fac=4._dp/3.0_dp/sqrt(pi)
3144  nr=0
3145  do
3146    nr=nr+1
3147    newr=0
3148 
3149    do ir3=-nr,nr
3150      do ir2=-nr,nr
3151        do ir1=-nr,nr
3152          if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
3153 
3154            arg=two_pi*(qphon(1)*ir1+qphon(2)*ir2+qphon(3)*ir3)
3155            c1r=cos(arg)*reta3m
3156            c1i=sin(arg)*reta3m
3157            do ia0=1,my_natom
3158              ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
3159              do ib=1,ia
3160                r1=dble(ir1)+xred(1,ib)-xred(1,ia)
3161                r2=dble(ir2)+xred(2,ib)-xred(2,ia)
3162                r3=dble(ir3)+xred(3,ib)-xred(3,ia)
3163                dakk(:)=two_pi*(/r1,r2,r3/)
3164                rdot12=rmet(2,1)*r1*r2
3165                rdot13=rmet(3,1)*r1*r3
3166                rdot23=rmet(3,2)*r2*r3
3167                dotr1=rmet(1,1)*r1**2+rdot12+rdot13
3168                dotr2=rmet(2,2)*r2**2+rdot12+rdot23
3169                dotr3=rmet(3,3)*r3**2+rdot13+rdot23
3170                rsq=dotr1+dotr2+dotr3
3171                rmagn=sqrt(rsq)
3172 !              Avoid zero denominators in term :
3173                if (rmagn>=1.0d-12) then
3174                  arg=reta*rmagn
3175                  term=zero
3176                  if (arg<8.0_dp) then
3177 !                  Note: erfc(8) is about 1.1e-29,
3178 !                  so don t bother with larger arg.
3179 !                  Also: exp(-64) is about 1.6e-28,
3180 !                  so don t bother with larger arg**2 in exp.
3181                    newr=1
3182                    derfc_arg = abi_derfc(arg)
3183                    term=derfc_arg/arg**3
3184                    term1=2.0_dp/sqrt(pi)*exp(-arg**2)/arg**2
3185                    term2=-(term+term1)
3186                    term3=(3*term+term1*(3.0_dp+2.0_dp*arg**2))/rsq
3187                    rq(1)=rmet(1,1)*r1+rmet(1,2)*r2+rmet(1,3)*r3
3188                    rq(2)=rmet(2,1)*r1+rmet(2,2)*r2+rmet(2,3)*r3
3189                    rq(3)=rmet(3,1)*r1+rmet(3,2)*r2+rmet(3,3)*r3
3190                    do iq2=1,3
3191                      do iq1=1,3
3192                        do mu=1,3
3193 !                         do nu=1,3
3194                          do nu=1,mu
3195                            work(re,mu,ia,nu,ib,iq1,iq2)=work(re,mu,ia,nu,ib,iq1,iq2)+&
3196 &                           c1r*dakk(iq1)*dakk(iq2)*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
3197                            work(im,mu,ia,nu,ib,iq1,iq2)=work(im,mu,ia,nu,ib,iq1,iq2)+&
3198 &                           c1i*dakk(iq1)*dakk(iq2)*(rq(mu)*rq(nu)*term3+rmet(mu,nu)*term2)
3199                          end do
3200                        end do
3201                      end do
3202                    end do
3203                  end if
3204                else
3205                  if (ia/=ib)then
3206                    write(message,'(a,a,a,a,a,i5,a,i5,a)')&
3207 &                   'The distance between two atoms vanishes.',ch10,&
3208 &                   'This is not allowed.',ch10,&
3209 &                   'Action: check the input for the atoms number',ia,' and',ib,'.'
3210                    ABI_ERROR(message)
3211                  end if
3212                end if
3213 
3214              end do ! End loop over ib:
3215            end do ! End loop over ia:
3216          end if
3217        end do ! End triple loop over real space points:
3218      end do
3219    end do
3220 
3221 !  Check if new shell must be calculated
3222    if(newr==0) exit
3223  end do !  End loop on nr (new shells). Note that there is an exit within the loop
3224 
3225 !Take account of the charges
3226 !write(std_out,*)' '
3227  do ia0=1,my_natom
3228    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
3229    do ib=1,ia
3230      do iq2=1,3
3231        do iq1=1,3
3232          do mu=1,3
3233            do nu=1,mu
3234              do ii=1,2
3235                work(ii,mu,ia,nu,ib,iq1,iq2)=work(ii,mu,ia,nu,ib,iq1,iq2)*&
3236 &               zion(typat(ia))*zion(typat(ib))
3237              end do
3238            end do
3239          end do
3240        end do
3241      end do
3242    end do
3243  end do
3244 
3245 !Symmetrize with respect to the directions
3246  do ia0=1,my_natom
3247    ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
3248    do ib=1,ia
3249      do iq2=1,3
3250        do iq1=1,3
3251          do mu=1,3
3252            do nu=1,mu
3253              work(re,nu,ia,mu,ib,iq1,iq2)=work(re,mu,ia,nu,ib,iq1,iq2)
3254              work(im,nu,ia,mu,ib,iq1,iq2)=work(im,mu,ia,nu,ib,iq1,iq2)
3255            end do
3256          end do
3257        end do
3258      end do
3259    end do
3260  end do
3261 
3262 !In case of parallelism over atoms: communicate
3263  if (paral_atom) then
3264    call timab(48,1,tsec)
3265    call xmpi_sum(work,my_comm_atom,ierr)
3266    call timab(48,2,tsec)
3267  end if
3268 
3269 !Fill the upper part of the matrix, with the hermitian conjugate
3270  do ia=1,natom
3271    do ib=1,ia
3272      do iq2=1,3
3273        do iq1=1,3
3274          do nu=1,3
3275            do mu=1,3
3276              work(re,mu,ib,nu,ia,iq1,iq2)=work(re,mu,ia,nu,ib,iq1,iq2)
3277              work(im,mu,ib,nu,ia,iq1,iq2)=-work(im,mu,ia,nu,ib,iq1,iq2)
3278            end do
3279          end do
3280        end do
3281      end do
3282    end do
3283  end do
3284 
3285 !Perform the summation over the second atomic sublattice
3286  dyewdqdq(:,:,:,:,:,:)=zero
3287  do ia=1,natom
3288    do iq2=1,3
3289      do iq1=1,3
3290        do nu=1,3
3291          do mu=1,3
3292            do ib=1,natom
3293              dyewdqdq(re,mu,ia,nu,iq1,iq2)=dyewdqdq(re,mu,ia,nu,iq1,iq2) + &
3294            & work(re,mu,ia,nu,ib,iq1,iq2)
3295              dyewdqdq(im,mu,ia,nu,iq1,iq2)=dyewdqdq(im,mu,ia,nu,iq1,iq2) + &
3296            & work(im,mu,ia,nu,ib,iq1,iq2)
3297            end do
3298          end do
3299        end do
3300      end do
3301    end do
3302  end do
3303  ABI_FREE(work)
3304 
3305 !Destroy atom table used for parallelism
3306  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
3307 
3308 end subroutine dfpt_ewalddqdq

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

SOURCE

1922 subroutine elt_ewald(elteew,gmet,gprimd,my_natom,natom,ntypat,rmet,rprimd,&
1923 &                 typat,ucvol,xred,zion, &
1924 &                 mpi_atmtab,comm_atom) ! optional arguments (parallelism)
1925 
1926 !Arguments ------------------------------------
1927 !scalars
1928  integer,intent(in) :: my_natom,natom,ntypat
1929  real(dp),intent(in) :: ucvol
1930 !arrays
1931  integer,intent(in) :: typat(natom)
1932  integer,optional,intent(in) :: comm_atom
1933  integer,optional,target,intent(in) :: mpi_atmtab(:)
1934  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
1935  real(dp),intent(in) :: xred(3,natom),zion(ntypat)
1936  real(dp),intent(out) :: elteew(6+3*natom,6)
1937 
1938 !Local variables-------------------------------
1939 !scalars
1940  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
1941  logical :: my_atmtab_allocated,paral_atom
1942  real(dp) :: arg,ch,chsq,cos_term,d2derfc,d2gss,d2r,d2rs,dderfc,derfc_arg
1943  real(dp) :: dgss1,dgss2,direct,dr1,dr2,drs1,drs2,eew,eta,fac,fraca1,fraca2
1944  real(dp) :: fraca3,fracb1,fracb2,fracb3,gsq,gsum,r1,r2,r3,recip,reta
1945  real(dp) :: rmagn,rsq,sin_term,sumg,summi,summr,sumr,t1,term
1946  character(len=500) :: message
1947 !arrays
1948  integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/)
1949  integer,pointer :: my_atmtab(:)
1950  real(dp) :: d2gm(3,3,6,6),d2ris(3),d2rm(3,3,6,6),dgm(3,3,6),dris(3),drm(3,3,6)
1951  real(dp) :: t2(3),ts2(3),tsec(2),tt(3)
1952  real(dp),allocatable :: d2sumg(:,:),d2sumr(:,:),drhoisi(:,:),drhoisr(:,:)
1953  real(dp),allocatable :: mpibuf(:)
1954 
1955 ! *************************************************************************
1956 
1957 !DEBUG
1958 !write(std_out,*)' elt_ewald : enter '
1959 !stop
1960 !ENDDEBUG
1961 
1962 !Compute 1st and 2nd derivatives of metric tensor wrt all strain components
1963 !and store for use in inner loop below.
1964 
1965 !Set up parallelism over atoms
1966  paral_atom=(present(comm_atom).and.(my_natom/=natom))
1967  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
1968  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1969  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
1970 
1971 !Loop over 2nd strain index
1972  do is2=1,6
1973    kg=idx(2*is2-1);kd=idx(2*is2)
1974    do jj = 1,3
1975      drm(:,jj,is2)=rprimd(kg,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kg,jj)
1976      dgm(:,jj,is2)=-(gprimd(kg,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kg,jj))
1977    end do
1978 
1979 !  Loop over 1st strain index, upper triangle only
1980    do is1=1,is2
1981 
1982      ka=idx(2*is1-1);kb=idx(2*is1)
1983      d2rm(:,:,is1,is2)=zero
1984      d2gm(:,:,is1,is2)=zero
1985      do jj = 1,3
1986        if(ka==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
1987 &       +rprimd(kb,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(kb,jj)
1988        if(ka==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
1989 &       +rprimd(kb,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(kb,jj)
1990        if(kb==kg) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
1991 &       +rprimd(ka,:)*rprimd(kd,jj)+rprimd(kd,:)*rprimd(ka,jj)
1992        if(kb==kd) d2rm(:,jj,is1,is2)=d2rm(:,jj,is1,is2)&
1993 &       +rprimd(ka,:)*rprimd(kg,jj)+rprimd(kg,:)*rprimd(ka,jj)
1994 
1995        if(ka==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1996 &       +gprimd(kb,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(kb,jj)
1997        if(ka==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
1998 &       +gprimd(kb,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(kb,jj)
1999        if(kb==kg) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
2000 &       +gprimd(ka,:)*gprimd(kd,jj)+gprimd(kd,:)*gprimd(ka,jj)
2001        if(kb==kd) d2gm(:,jj,is1,is2)=d2gm(:,jj,is1,is2)&
2002 &       +gprimd(ka,:)*gprimd(kg,jj)+gprimd(kg,:)*gprimd(ka,jj)
2003      end do
2004    end do !is1
2005  end do !is2
2006 
2007 !Add up total charge and sum of $charge^2$ in cell
2008  chsq=zero
2009  ch=zero
2010  do ia=1,natom
2011    ch=ch+zion(typat(ia))
2012    chsq=chsq+zion(typat(ia))**2
2013  end do
2014 
2015 !Compute eta, the Ewald summation convergence parameter,
2016 !for approximately optimized summations:
2017  direct=rmet(1,1)+rmet(1,2)+rmet(1,3)+rmet(2,1)+&
2018 & rmet(2,2)+rmet(2,3)+rmet(3,1)+rmet(3,2)+rmet(3,3)
2019  recip=gmet(1,1)+gmet(1,2)+gmet(1,3)+gmet(2,1)+&
2020 & gmet(2,2)+gmet(2,3)+gmet(3,1)+gmet(3,2)+gmet(3,3)
2021 !Here, a bias is introduced, because G-space summation scales
2022 !better than r space summation ! Note : debugging is the most
2023 !easier at fixed eta.
2024 ! eta=pi*200._dp/33.0_dp*sqrt(1.69_dp*recip/direct)
2025 eta=1.0_dp
2026 
2027 !Conduct reciprocal space summations
2028  fac=pi**2/eta ; gsum=zero
2029  ABI_MALLOC(d2sumg,(6+3*natom,6))
2030  ABI_MALLOC(drhoisr,(3,natom))
2031  ABI_MALLOC(drhoisi,(3,natom))
2032  d2sumg(:,:)=zero
2033 
2034 !Sum over G space, done shell after shell until all
2035 !contributions are too small.
2036  ng=0
2037  do
2038    ng=ng+1
2039    newg=0
2040 
2041    do ig3=-ng,ng
2042      do ig2=-ng,ng
2043        do ig1=-ng,ng
2044 
2045 !        Exclude shells previously summed over
2046          if(abs(ig1)==ng .or. abs(ig2)==ng .or. abs(ig3)==ng&
2047 &         .or. ng==1 ) then
2048 
2049 !          gsq is G dot G = |G|^2
2050            gsq=gmet(1,1)*dble(ig1*ig1)+gmet(2,2)*dble(ig2*ig2)+&
2051 &           gmet(3,3)*dble(ig3*ig3)+2._dp*(gmet(2,1)*dble(ig1*ig2)+&
2052 &           gmet(3,1)*dble(ig1*ig3)+gmet(3,2)*dble(ig3*ig2))
2053 
2054 !          Skip g=0:
2055            if (gsq>1.0d-20) then
2056              arg=fac*gsq
2057 
2058 !            Larger arg gives 0 contribution because of exp(-arg)
2059              if (arg <= 80._dp) then
2060 !              When any term contributes then include next shell
2061                newg=1
2062                term=exp(-arg)/gsq
2063                summr = zero
2064                summi = zero
2065 !              Note that if reduced atomic coordinates xred drift outside
2066 !              of unit cell (outside [0,1)) it is irrelevant in the following
2067 !              term, which only computes a phase.
2068                do ia=1,natom
2069                  arg=two_pi*(ig1*xred(1,ia)+ig2*xred(2,ia)+ig3*xred(3,ia))
2070 !                Sum real and imaginary parts (avoid complex variables)
2071                  cos_term=cos(arg)
2072                  sin_term=sin(arg)
2073                  summr=summr+zion(typat(ia))*cos_term
2074                  summi=summi+zion(typat(ia))*sin_term
2075                  drhoisr(1,ia)=-two_pi*zion(typat(ia))*sin_term*dble(ig1)
2076                  drhoisi(1,ia)= two_pi*zion(typat(ia))*cos_term*dble(ig1)
2077                  drhoisr(2,ia)=-two_pi*zion(typat(ia))*sin_term*dble(ig2)
2078                  drhoisi(2,ia)= two_pi*zion(typat(ia))*cos_term*dble(ig2)
2079                  drhoisr(3,ia)=-two_pi*zion(typat(ia))*sin_term*dble(ig3)
2080                  drhoisi(3,ia)= two_pi*zion(typat(ia))*cos_term*dble(ig3)
2081                end do
2082 
2083 !              The following two checks avoid an annoying
2084 !              underflow error message
2085                if (abs(summr)<1.d-16) summr=zero
2086                if (abs(summi)<1.d-16) summi=zero
2087 
2088 !              The product of term and summr**2 or summi**2 below
2089 !              can underflow if not for checks above
2090                t1=term*(summr*summr+summi*summi)
2091                gsum=gsum+t1
2092 !              Loop over 2nd strain index
2093                do is2=1,6
2094                  dgss2=dgm(1,1,is2)*dble(ig1*ig1)+dgm(2,2,is2)*dble(ig2*ig2)+&
2095 &                 dgm(3,3,is2)*dble(ig3*ig3)+2._dp*(dgm(2,1,is2)*dble(ig1*ig2)+&
2096 &                 dgm(3,1,is2)*dble(ig1*ig3)+dgm(3,2,is2)*dble(ig3*ig2))
2097 !                Loop over 1st strain index, upper triangle only
2098                  do is1=1,is2
2099                    dgss1=dgm(1,1,is1)*dble(ig1*ig1)+dgm(2,2,is1)*dble(ig2*ig2)+&
2100 &                   dgm(3,3,is1)*dble(ig3*ig3)+2._dp*(dgm(2,1,is1)*dble(ig1*ig2)+&
2101 &                   dgm(3,1,is1)*dble(ig1*ig3)+dgm(3,2,is1)*dble(ig3*ig2))
2102 
2103                    d2gss=d2gm(1,1,is1,is2)*dble(ig1*ig1)+&
2104 &                   d2gm(2,2,is1,is2)*dble(ig2*ig2)+&
2105 &                   d2gm(3,3,is1,is2)*dble(ig3*ig3)+2._dp*&
2106 &                   (d2gm(2,1,is1,is2)*dble(ig1*ig2)+&
2107 &                   d2gm(3,1,is1,is2)*dble(ig1*ig3)+&
2108 &                   d2gm(3,2,is1,is2)*dble(ig3*ig2))
2109 
2110                    d2sumg(is1,is2)=d2sumg(is1,is2)+&
2111 &                   t1*((fac**2 + 2.0_dp*fac/gsq + 2.0_dp/(gsq**2))*dgss1*dgss2 -&
2112 &                   0.5_dp*(fac + 1.0_dp/gsq)*d2gss)
2113                    if(is1<=3) d2sumg(is1,is2)=d2sumg(is1,is2)+&
2114 &                   t1*(fac + 1.0_dp/gsq)*dgss2
2115                    if(is2<=3) d2sumg(is1,is2)=d2sumg(is1,is2)+&
2116 &                   t1*(fac + 1.0_dp/gsq)*dgss1
2117                    if(is1<=3 .and. is2<=3) d2sumg(is1,is2)=d2sumg(is1,is2)+t1
2118 
2119                  end do !is1
2120 
2121 !                Internal strain contributions
2122                  do ia=1,natom
2123                    js=7+3*(ia-1)
2124                    t2(:)=2.0_dp*term*(summr*drhoisr(:,ia)+summi*drhoisi(:,ia))
2125                    d2sumg(js:js+2,is2)=d2sumg(js:js+2,is2)-&
2126 &                   (fac + 1.0_dp/gsq)*dgss2*t2(:)
2127                    if(is2<=3) d2sumg(js:js+2,is2)=d2sumg(js:js+2,is2)-t2(:)
2128                  end do
2129                end do !is2
2130 
2131              end if ! End condition of not larger than 80.0
2132            end if ! End skip g=0
2133          end if ! End triple loop over G s and associated new shell condition
2134        end do
2135      end do
2136    end do
2137 
2138 !  Check if new shell must be calculated
2139    if (newg==0) exit
2140  end do !  End the loop on ng (new shells). Note that there is one exit from this loop.
2141 
2142  sumg=gsum/(two_pi*ucvol)
2143  d2sumg(:,:)=d2sumg(:,:)/(two_pi*ucvol)
2144 
2145  ABI_FREE(drhoisr)
2146  ABI_FREE(drhoisi)
2147 !Stress tensor is now computed elsewhere (ewald2) hence do not need
2148 !length scale gradients (used to compute them here).
2149 
2150 !Conduct real space summations
2151  reta=sqrt(eta)
2152  fac=2._dp*sqrt(eta/pi)
2153  ABI_MALLOC(d2sumr,(6+3*natom,6))
2154  sumr=zero;d2sumr(:,:)=zero
2155 
2156 !In the following a summation is being conducted over all
2157 !unit cells (ir1, ir2, ir3) so it is appropriate to map all
2158 !reduced coordinates xred back into [0,1).
2159 !
2160 !Loop on shells in r-space as was done in g-space
2161  nr=0
2162  do
2163    nr=nr+1
2164    newr=0
2165 !
2166    do ir3=-nr,nr
2167      do ir2=-nr,nr
2168        do ir1=-nr,nr
2169          if( abs(ir3)==nr .or. abs(ir2)==nr .or. abs(ir1)==nr .or. nr==1 )then
2170 
2171            do ia0=1,my_natom
2172              ia=ia0;if(paral_atom)ia=my_atmtab(ia0)
2173              js=7+3*(ia-1)
2174 !            Map reduced coordinate xred(mu,ia) into [0,1)
2175              fraca1=xred(1,ia)-aint(xred(1,ia))+0.5_dp-sign(0.5_dp,xred(1,ia))
2176              fraca2=xred(2,ia)-aint(xred(2,ia))+0.5_dp-sign(0.5_dp,xred(2,ia))
2177              fraca3=xred(3,ia)-aint(xred(3,ia))+0.5_dp-sign(0.5_dp,xred(3,ia))
2178              do ib=1,natom
2179                fracb1=xred(1,ib)-aint(xred(1,ib))+0.5_dp-sign(0.5_dp,xred(1,ib))
2180                fracb2=xred(2,ib)-aint(xred(2,ib))+0.5_dp-sign(0.5_dp,xred(2,ib))
2181                fracb3=xred(3,ib)-aint(xred(3,ib))+0.5_dp-sign(0.5_dp,xred(3,ib))
2182                r1=dble(ir1)+fracb1-fraca1
2183                r2=dble(ir2)+fracb2-fraca2
2184                r3=dble(ir3)+fracb3-fraca3
2185                rsq=rmet(1,1)*r1*r1+rmet(2,2)*r2*r2+rmet(3,3)*r3*r3+&
2186 &               2.0_dp*(rmet(2,1)*r2*r1+rmet(3,2)*r3*r2+rmet(3,1)*r1*r3)
2187 
2188 !              Avoid zero denominators in 'term':
2189                if (rsq>=1.0d-24) then
2190 
2191 !                Note: erfc(8) is about 1.1e-29,
2192 !                so do not bother with larger arg.
2193 !                Also: exp(-64) is about 1.6e-28,
2194 !                so do not bother with larger arg**2 in exp.
2195                  term=zero
2196                  if (eta*rsq<64.0_dp) then
2197                    newr=1
2198                    rmagn=sqrt(rsq)
2199                    arg=reta*rmagn
2200 !                  derfc computes the complementary error function
2201 !                  dderfc is the derivative of the complementary error function
2202 !                  d2derfc is the 2nd derivative of the complementary error function
2203                    dderfc=-fac*exp(-eta*rsq)
2204                    d2derfc=-2._dp*eta*rmagn*dderfc
2205                    derfc_arg = abi_derfc(arg)
2206                    term=derfc_arg/rmagn
2207                    sumr=sumr+zion(typat(ia))*zion(typat(ib))*term
2208                    tt(:)=rmet(:,1)*r1+rmet(:,2)*r2+rmet(:,3)*r3
2209                    dris(:)=tt(:)/rmagn
2210 !                  Loop over 2nd strain index
2211                    do is2=1,6
2212                      drs2=drm(1,1,is2)*r1*r1+drm(2,2,is2)*r2*r2+&
2213 &                     drm(3,3,is2)*r3*r3+&
2214 &                     2.0_dp*(drm(2,1,is2)*r2*r1+drm(3,2,is2)*r3*r2+&
2215 &                     drm(3,1,is2)*r1*r3)
2216                      dr2=0.5_dp*drs2/rmagn
2217 !                    Loop over 1st strain index, upper triangle only
2218                      do is1=1,is2
2219                        drs1=drm(1,1,is1)*r1*r1+drm(2,2,is1)*r2*r2+&
2220 &                       drm(3,3,is1)*r3*r3+&
2221 &                       2.0_dp*(drm(2,1,is1)*r2*r1+drm(3,2,is1)*r3*r2+&
2222 &                       drm(3,1,is1)*r1*r3)
2223                        dr1=0.5_dp*drs1/rmagn
2224                        d2rs=d2rm(1,1,is1,is2)*r1*r1+d2rm(2,2,is1,is2)*r2*r2+&
2225 &                       d2rm(3,3,is1,is2)*r3*r3+&
2226 &                       2.0_dp*(d2rm(2,1,is1,is2)*r2*r1+d2rm(3,2,is1,is2)*r3*r2+&
2227 &                       d2rm(3,1,is1,is2)*r1*r3)
2228                        d2r=(0.25_dp*d2rs-dr1*dr2)/rmagn
2229                        d2sumr(is1,is2)=d2sumr(is1,is2)+&
2230 &                       zion(typat(ia))*zion(typat(ib))*&
2231 &                       ((d2derfc-2.0_dp*dderfc/rmagn+2.0_dp*derfc_arg/rsq)*dr1*dr2+&
2232 &                       (dderfc-derfc_arg/rmagn)*d2r)/rmagn
2233                      end do !is1
2234 !                    Internal strain contribution
2235                      ts2(:)=drm(:,1,is2)*r1+drm(:,2,is2)*r2+drm(:,3,is2)*r3
2236                      d2ris(:)=ts2(:)/rmagn-0.5_dp*drs2*tt(:)/(rsq*rmagn)
2237 
2238                      d2sumr(js:js+2,is2)=d2sumr(js:js+2,is2)-&
2239 &                     2.0_dp*zion(typat(ia))*zion(typat(ib))*&
2240 &                     ((d2derfc-2.0_dp*dderfc/rmagn+2.0_dp*derfc_arg/rsq)*dr2*dris(:)+&
2241 &                     (dderfc-derfc_arg/rmagn)*d2ris(:))/rmagn
2242                    end do !is2
2243                  end if
2244                end if ! End avoid zero denominators in'term'
2245 
2246              end do ! end loop over ib:
2247            end do ! end loop over ia:
2248          end if ! end triple loop over real space points and associated condition of new shell
2249        end do
2250      end do
2251    end do
2252 
2253 !  Check if new shell must be calculated
2254    if(newr==0) exit
2255  end do !  End loop on nr (new shells). Note that there is an exit within the loop
2256 
2257 !In case of parallelism over atoms: communicate
2258  if (paral_atom) then
2259    call timab(48,1,tsec)
2260    ABI_MALLOC(mpibuf,((6+3*natom)*6+1))
2261    mpibuf(1:(6+3*natom)*6)=reshape(d2sumr(:,:),shape=(/((6+3*natom)*6)/))
2262    mpibuf((6+3*natom)*6+1)=sumr
2263    call xmpi_sum(mpibuf,my_comm_atom,ierr)
2264    sumr=mpibuf((6+3*natom)*6+1)
2265    d2sumr(:,:)=reshape(mpibuf(1:(6+3*natom)*6),shape=(/(6+3*natom),6/))
2266    ABI_FREE(mpibuf)
2267    call timab(48,2,tsec)
2268  end if
2269 
2270  sumr=0.5_dp*sumr
2271  d2sumr(:,:)=0.5_dp*d2sumr(:,:)
2272  fac=pi*ch**2/(2.0_dp*eta*ucvol)
2273 
2274 !Finally assemble Ewald energy, eew
2275  eew=sumg+sumr-chsq*reta/sqrt(pi)-fac
2276 
2277  elteew(:,:)=d2sumg(:,:)+d2sumr(:,:)
2278 
2279 !Additional term for all strains diagonal (from "fac" term in eew)
2280  elteew(1:3,1:3)=elteew(1:3,1:3)-fac
2281 
2282 !Fill in lower triangle
2283  do is2=2,6
2284    do is1=1,is2-1
2285      elteew(is2,is1)=elteew(is1,is2)
2286    end do
2287  end do
2288 
2289  ABI_FREE(d2sumg)
2290  ABI_FREE(d2sumr)
2291 
2292 !Output the final values of ng and nr
2293  write(message, '(a,i4,a,i4)' )' elt_ewald : nr and ng are ',nr,' and ',ng
2294  call wrtout(std_out,message,'COLL')
2295 
2296 !Destroy atom table used for parallelism
2297  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2298 
2299 end subroutine elt_ewald

ABINIT/eltxccore [ Functions ]

[ Top ] [ Functions ]

NAME

 eltxccore

FUNCTION

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

COPYRIGHT

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

SOURCE

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

ABINIT/m_dfpt_elt [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dfpt_elt

FUNCTION

COPYRIGHT

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

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_dfpt_elt
23 
24  use defs_basis
25  use m_abicore
26  use m_errors
27  use m_xmpi
28  use m_dtset
29 
30  use defs_datatypes, only : pseudopotential_type
31  use defs_abitypes, only : MPI_type
32  use m_time,        only : timab
33  use m_special_funcs,  only : abi_derfc
34  use m_geometry,    only : metric
35  use m_cgtools,     only : dotprod_vn
36  use m_pawtab,      only : pawtab_type,pawtab_free,pawtab_nullify
37  use m_pawrad,      only : pawrad_type,pawrad_init,pawrad_free
38  use m_pawpsp,      only : pawpsp_cg
39  use m_paw_numeric, only : paw_spline
40  use m_spacepar,    only : redgr
41  use m_atm2fft,     only : atm2fft, dfpt_atm2fft
42  use m_mkcore,      only : dfpt_mkcore
43  use m_dfpt_mkvxcstr, only : dfpt_mkvxcstr
44  use m_paral_atom,   only : get_my_atmtab, free_my_atmtab
45  use m_mpinfo,       only : ptabs_fourdp, proc_distrb_cycle, proc_distrb_nband
46  use m_fftcore,      only : sphereboundary
47  use m_fft,          only : fourdp
48 
49  implicit none
50 
51  private