TABLE OF CONTENTS
- ABINIT/d2kindstr2
- ABINIT/dfpt_eltfrhar
- ABINIT/dfpt_eltfrkin
- ABINIT/dfpt_eltfrloc
- ABINIT/dfpt_eltfrxc
- ABINIT/dfpt_ewald
- ABINIT/dfpt_ewalddq
- ABINIT/dfpt_ewalddqdq
- ABINIT/elt_ewald
- ABINIT/eltxccore
- ABINIT/m_dfpt_elt
ABINIT/d2kindstr2 [ 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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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