TABLE OF CONTENTS
- ABINIT/m_vcoul
- m_vcoul/cmod_qpg
- m_vcoul/cutoff_cylinder
- m_vcoul/cutoff_sphere
- m_vcoul/cutoff_surface
- m_vcoul/vcoul_free
- m_vcoul/vcoul_init
- m_vcoul/vcoul_plot
- m_vcoul/vcoul_print
- m_vcoul/vcoul_t
ABINIT/m_vcoul [ Modules ]
NAME
m_vcoul
FUNCTION
This module contains the definition of the vcoul_t as well as procedures to calculate the Coulomb interaction in reciprocal space taking into account a possible cutoff in real space. Procedures to deal with the singularity for q-->0 are also provided.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (MG, FB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
PARENTS
CHILDREN
SOURCE
25 #if defined HAVE_CONFIG_H 26 #include "config.h" 27 #endif 28 29 #include "abi_common.h" 30 31 MODULE m_vcoul 32 33 use defs_basis 34 use m_abicore 35 use m_errors 36 use m_xmpi 37 use m_splines 38 use m_sort 39 40 use m_fstrings, only : sjoin, itoa 41 use m_special_funcs, only : abi_derf 42 use m_gwdefs, only : GW_TOLQ0 43 use m_io_tools, only : open_file 44 use m_numeric_tools, only : arth, geop, imin_loc, llsfit_svd, l2norm, OPERATOR(.x.), quadrature, isdiagmat 45 use m_bessel, only : CALJY0, CALJY1, CALCK0, CALCK1 46 use m_hide_lapack, only : matrginv 47 use m_geometry, only : normv, metric 48 use m_crystal, only : crystal_t 49 use m_bz_mesh, only : kmesh_t, get_BZ_item 50 use m_gsphere, only : gsphere_t 51 use m_paw_numeric, only : paw_jbessel 52 use m_dtfil, only : isfile 53 54 implicit none 55 56 private
m_vcoul/cmod_qpg [ Functions ]
[ Top ] [ m_vcoul ] [ Functions ]
NAME
cmod_qpg
FUNCTION
Set up table of lengths |q+G| for the Coulomb potential.
INPUTS
gvec(3,npwvec)=Reduced coordinates of the G vectors. gprimd(3,3)=Dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$) iq=Index specifying the q point. npwvec=Number of planewaves nq=Number of q points. q=Coordinates of q points.
OUTPUT
qplusg(npwvec)=Norm of q+G vector
PARENTS
m_ppmodel,m_vcoul
CHILDREN
calck0,paw_jbessel,quadrature
SOURCE
2133 subroutine cmod_qpg(nq,iq,q,npwvec,gvec,gprimd,qplusg) 2134 2135 2136 !This section has been created automatically by the script Abilint (TD). 2137 !Do not modify the following lines by hand. 2138 #undef ABI_FUNC 2139 #define ABI_FUNC 'cmod_qpg' 2140 !End of the abilint section 2141 2142 implicit none 2143 2144 !Arguments ------------------------------------ 2145 !scalars 2146 integer,intent(in) :: iq,npwvec,nq 2147 !arrays 2148 integer,intent(in) :: gvec(3,npwvec) 2149 real(dp),intent(in) :: gprimd(3,3),q(3,nq) 2150 real(dp),intent(out) :: qplusg(npwvec) 2151 2152 !Local variables ------------------------------ 2153 !scalars 2154 integer :: ig,ii 2155 !arrays 2156 real(dp) :: gmet(3,3),gpq(3) 2157 2158 !************************************************************************ 2159 2160 ! Compute reciprocal space metrics 2161 do ii=1,3 2162 gmet(ii,:)=gprimd(1,ii)*gprimd(1,:)+& 2163 & gprimd(2,ii)*gprimd(2,:)+& 2164 & gprimd(3,ii)*gprimd(3,:) 2165 end do 2166 2167 if (ALL(ABS(q(:,iq))<1.e-3)) then !FIXME avoid this, everything should be under the control of the programmer. 2168 ! * Treat q as it were zero except when G=0 2169 qplusg(1)=two_pi*SQRT(DOT_PRODUCT(q(:,iq),MATMUL(gmet,q(:,iq)))) 2170 do ig=2,npwvec 2171 gpq(:)=gvec(:,ig) 2172 qplusg(ig)=two_pi*SQRT(DOT_PRODUCT(gpq,MATMUL(gmet,gpq))) 2173 end do 2174 else 2175 do ig=1,npwvec 2176 gpq(:)=gvec(:,ig)+q(:,iq) 2177 qplusg(ig)=two_pi*SQRT(DOT_PRODUCT(gpq,MATMUL(gmet,gpq))) 2178 end do 2179 end if 2180 2181 end subroutine cmod_qpg 2182 2183 !---------------------------------------------------------------------- 2184 2185 function F2(xx) 2186 2187 2188 !This section has been created automatically by the script Abilint (TD). 2189 !Do not modify the following lines by hand. 2190 #undef ABI_FUNC 2191 #define ABI_FUNC 'F2' 2192 !End of the abilint section 2193 2194 real(dp),intent(in) :: xx 2195 real(dp) :: F2 2196 2197 !Local variables------------------------------- 2198 !scalars 2199 integer :: ierr 2200 real(dp) :: intr 2201 !************************************************************************ 2202 2203 zz_=xx 2204 call quadrature(F1,zero,rcut_,qopt_,intr,ierr,ntrial_,accuracy_,npts_) 2205 if (ierr/=0) then 2206 MSG_ERROR("Accuracy not reached") 2207 end if 2208 2209 F2=intr*COS(qpg_para_*xx) 2210 2211 end function F2
m_vcoul/cutoff_cylinder [ Functions ]
[ Top ] [ m_vcoul ] [ Functions ]
NAME
cutoff_cylinder
FUNCTION
Calculate the Fourier components of an effective Coulomb interaction
zeroed outside a finite cylindrical region.
Two methods are implemented:
method==1: The interaction in the (say) x-y plane is truncated outside the Wigner-Seitz
cell centered on the wire in the x-y plane. The interaction has infinite
extent along the z axis and the Fourier transform is singular only at the Gamma point.
Only orthorombic Bravais lattice are supported.
method==2: The interaction is truncated outside a cylinder of radius rcut. The cylinder has finite
extent along z. No singularity occurs.
INPUTS
boxcenter(3)= center of the wire in the x-y axis gvec(3,ng)=G vectors in reduced coordinates ng=number of G vectors qpt(3,nq)= q-points nq=number of q-points rprimd(3,3)=dimensional real space primitive translations (bohr) where: rprimd(i,j)=rprim(i,j)*acell(j) method=1 for Beigi approach (infinite cylinder with interaction truncated outside the W-S cell) 2 for Rozzi method (finite cylinder) comm=MPI communicator.
OUTPUT
vccut(ng,nq)= Fourier components of the effective Coulomb interaction
PARENTS
m_vcoul
CHILDREN
calck0,paw_jbessel,quadrature
SOURCE
1672 subroutine cutoff_cylinder(nq,qpt,ng,gvec,rcut,hcyl,pdir,boxcenter,rprimd,vccut,method,comm) 1673 1674 1675 !This section has been created automatically by the script Abilint (TD). 1676 !Do not modify the following lines by hand. 1677 #undef ABI_FUNC 1678 #define ABI_FUNC 'cutoff_cylinder' 1679 !End of the abilint section 1680 1681 implicit none 1682 1683 !Arguments ------------------------------------ 1684 !scalars 1685 integer,intent(in) :: ng,nq,method,comm 1686 real(dp),intent(in) :: rcut,hcyl 1687 !arrays 1688 integer,intent(in) :: gvec(3,ng),pdir(3) 1689 real(dp),intent(in) :: boxcenter(3),qpt(3,nq),rprimd(3,3) 1690 real(dp),intent(out) :: vccut(ng,nq) 1691 1692 !Local variables------------------------------- 1693 !scalars 1694 integer,parameter :: N0=1000 1695 integer :: icount,ig,igs,iq 1696 integer :: ntasks,ierr,my_start,my_stop 1697 real(dp) :: j0,j1,k0,k1,qpg2,qpg_xy,tmp 1698 real(dp) :: qpg_z,quad,rcut2,hcyl2,c1,c2,ucvol,SMALL 1699 logical :: q_is_zero 1700 character(len=500) :: msg 1701 !arrays 1702 real(dp) :: qpg(3),b1(3),b2(3),b3(3),gmet(3,3),rmet(3,3),gprimd(3,3),qc(3),gcart(3) 1703 real(dp),allocatable :: qcart(:,:) 1704 !************************************************************************ 1705 1706 ABI_UNUSED(pdir) 1707 ABI_UNUSED(boxcenter) 1708 ! 1709 ! =================================================== 1710 ! === Setup for the quadrature of matrix elements === 1711 ! =================================================== 1712 qopt_ =6 ! Quadrature method, see quadrature routine. 1713 ntrial_ =30 ! Max number of attempts. 1714 accuracy_=0.001 ! Fractional accuracy required. 1715 npts_ =6 ! Initial number of point (only for Gauss-Legendre method). 1716 SMALL =GW_TOLQ0 ! Below this value (q+G)_i is treated as zero. 1717 rcut_ =rcut ! Radial cutoff, used only if method==2 1718 hcyl_ =hcyl ! Lenght of cylinder along z, only if method==2 1719 1720 write(msg,'(3a,2(a,i5,a),a,f8.5)')ch10,& 1721 & ' cutoff_cylinder: Info on the quadrature method : ',ch10,& 1722 & ' Quadrature scheme = ',qopt_,ch10,& 1723 & ' Max number of attempts = ',ntrial_,ch10,& 1724 & ' Fractional accuracy = ',accuracy_ 1725 call wrtout(std_out,msg,'COLL') 1726 ! 1727 ! === From reduced to Cartesian coordinates === 1728 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1729 b1(:)=two_pi*gprimd(:,1) 1730 b2(:)=two_pi*gprimd(:,2) 1731 b3(:)=two_pi*gprimd(:,3) 1732 1733 ABI_MALLOC(qcart,(3,nq)) 1734 do iq=1,nq 1735 qcart(:,iq)=b1(:)*qpt(1,iq)+b2(:)*qpt(2,iq)+b3(:)*qpt(3,iq) 1736 end do 1737 1738 ntasks=nq*ng 1739 call xmpi_split_work(ntasks,comm,my_start,my_stop,msg,ierr) 1740 ! 1741 ! ================================================ 1742 ! === Different approaches according to method === 1743 ! ================================================ 1744 vccut(:,:)=zero 1745 1746 SELECT CASE (method) 1747 1748 CASE (1) 1749 ! === Infinite cylinder, interaction is zeroed outside the Wigner-Seitz cell === 1750 ! * Beigi"s expression holds only if the BZ is sampled only along z. 1751 write(msg,'(2(a,f8.4))')' cutoff_cylinder: Using Beigi''s Infinite cylinder ' 1752 call wrtout(std_out,msg,'COLL') 1753 if (ANY(qcart(1:2,:)>SMALL)) then 1754 write(std_out,*)' qcart = ',qcart(:,:) 1755 write(msg,'(5a)')& 1756 & ' found q-points with non zero components in the X-Y plane. ',ch10,& 1757 & ' This is not allowed, see Notes in cutoff_cylinder.F90. ',ch10,& 1758 & ' ACTION: Modify the q-point sampling. ' 1759 MSG_ERROR(msg) 1760 end if 1761 ! * Check if Bravais lattice is orthorombic and parallel to the Cartesian versors. 1762 ! In this case the intersection of the W-S cell with the x-y plane is a rectangle with -ha_<=x<=ha_ and -hb_<=y<=hb_ 1763 if ( (ANY(ABS(rprimd(2:3, 1))>tol6)).or.& 1764 & (ANY(ABS(rprimd(1:3:2,2))>tol6)).or.& 1765 & (ANY(ABS(rprimd(1:2, 3))>tol6)) & 1766 & ) then 1767 msg = ' Bravais lattice should be orthorombic and parallel to the cartesian versors ' 1768 MSG_ERROR(msg) 1769 end if 1770 1771 ha_=half*SQRT(DOT_PRODUCT(rprimd(:,1),rprimd(:,1))) 1772 hb_=half*SQRT(DOT_PRODUCT(rprimd(:,2),rprimd(:,2))) 1773 r0_=MIN(ha_,hb_)/N0 1774 ! 1775 ! === For each pair (q,G) evaluate the integral defining the cutoff Coulomb === 1776 ! * Here the code assumes that all q-vectors are non zero and q_xy/=0. 1777 do iq=1,nq 1778 igs=1 1779 ! * Skip singularity at Gamma, it will be treated "by hand" in csigme. 1780 q_is_zero = (normv(qpt(:,iq),gmet,'G')<GW_TOLQ0) 1781 !if (q_is_zero) igs=2 1782 qc(:)=qcart(:,iq) 1783 write(msg,'(2(a,i3))')' entering loop iq: ',iq,' with igs = ',igs 1784 call wrtout(std_out,msg,'COLL') 1785 1786 do ig=igs,ng 1787 icount=ig+(iq-1)*ng; if (icount<my_start.or.icount>my_stop) CYCLE 1788 1789 gcart(:)=b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig) 1790 qpg(:)=qc(:)+gcart(:) 1791 qpgx_=qpg(1) ; qpgy_=qpg(2) ; qpg_para_=ABS(qpg(3)) 1792 write(std_out,*)"qpgx_=",qpgx_ 1793 write(std_out,*)"qpgy_=",qpgy_ 1794 write(std_out,*)"qpg_para=",qpg_para_ 1795 1796 ! Avoid singularity in K_0{qpg_para_\rho) by using a small q along the periodic dimension. 1797 if (q_is_zero.and.qpg_para_<tol6) then 1798 qpg_para_ = tol6 1799 write(std_out,*)"setting qpg_para to=",qpg_para_ 1800 end if 1801 ! 1802 ! * Calculate $ 2\int_{WS} dxdy K_0{qpg_para_\rho) cos(x.qpg_x + y.qpg_y) $ 1803 ! where WS is the Wigner-Seitz cell. 1804 tmp=zero 1805 ! Difficult part, integrate on a small cirle of radius r0 using spherical coordinates 1806 !call quadrature(K0cos_dth_r0,zero,r0_,qopt_,quad,ierr,ntrial_,accuracy_,npts_) 1807 !if (ierr/=0) MSG_ERROR("Accuracy not reached") 1808 !write(std_out,'(i8,a,es14.6)')ig,' 1 ',quad 1809 !tmp=tmp+quad 1810 ! Add region with 0<=x<=r0 and y>=+-(SQRT(r0^2-x^2))since WS is rectangular 1811 !call quadrature(K0cos_dy_r0,zero,r0_,qopt_,quad,ierr,ntrial_,accuracy_,npts_) 1812 !if (ierr/=0) MSG_ERROR("Accuracy not reached") 1813 !write(std_out,'(i8,a,es14.6)')ig,' 2 ',quad 1814 !tmp=tmp+quad 1815 ! Get the in integral in the rectangle with x>=r0, should be the easiest but sometimes has problems to converge 1816 !call quadrature(K0cos_dy,r0_,ha_,qopt_,quad,ierr,ntrial_,accuracy_,npts_) 1817 !if (ierr/=0) MSG_ERROR("Accuracy not reached") 1818 !write(std_out,'(i8,a,es14.6)')ig,' 3 ',quad 1819 ! 1820 ! === More stable method: midpoint integration with Romberg extrapolation === 1821 call quadrature(K0cos_dy,zero,ha_,qopt_,quad,ierr,ntrial_,accuracy_,npts_) 1822 !write(std_out,'(i8,a,es14.6)')ig,' 3 ',quad 1823 if (ierr/=0) then 1824 MSG_ERROR("Accuracy not reached") 1825 end if 1826 ! === Store final result === 1827 ! * Factor two comes from the replacement WS -> (1,4) quadrant thanks to symmetries of the integrad. 1828 tmp=tmp+quad 1829 vccut(ig,iq)=two*(tmp*two) 1830 end do !ig 1831 end do !iq 1832 1833 CASE (2) 1834 ! === Finite cylinder of length hcyl, from Rozzi et al === 1835 ! TODO add check on hcyl value that should be smaller that 1/deltaq 1836 if (hcyl_<zero) then 1837 write(msg,'(a,f8.4)')' Negative value for cylinder length hcyl_=',hcyl_ 1838 MSG_BUG(msg) 1839 end if 1840 1841 if (ABS(hcyl_)>tol12) then 1842 write(msg,'(2(a,f8.4))')' cutoff_cylinder: using finite cylinder of length= ',hcyl_,' rcut= ',rcut_ 1843 call wrtout(std_out,msg,'COLL') 1844 hcyl2=hcyl_**2 1845 rcut2=rcut_**2 1846 1847 do iq=1,nq 1848 write(msg,'(a,i3)')' entering loop iq: ',iq 1849 call wrtout(std_out,msg,'COLL') 1850 qc(:)=qcart(:,iq) 1851 do ig=1,ng 1852 ! === No singularity occurs in finite cylinder, thus start from 1 === 1853 icount=ig+(iq-1)*ng ; if (icount<my_start.or.icount>my_stop) CYCLE 1854 gcart(:)=b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig) 1855 qpg(:)=qc(:)+gcart(:) 1856 qpg_para_=ABS(qpg(3)) ; qpg_perp_=SQRT(qpg(1)**2+qpg(2)**2) 1857 1858 if (qpg_perp_/=zero.and.qpg_para_/=zero) then 1859 ! $ 4\pi\int_0^{R_c} d\rho\rho j_o(qpg_perp_.\rho)\int_0^hcyl dz\cos(qpg_para_*z)/sqrt(\rho^2+z^2) $ 1860 call quadrature(F2,zero,rcut_,qopt_,quad,ierr,ntrial_,accuracy_,npts_) 1861 if (ierr/=0) then 1862 MSG_ERROR("Accuracy not reached") 1863 end if 1864 1865 vccut(ig,iq)=four_pi*quad 1866 1867 else if (qpg_perp_==zero.and.qpg_para_/=zero) then 1868 ! $ \int_0^h sin(qpg_para_.z)/\sqrt(rcut^2+z^2)dz $ 1869 call quadrature(F3,zero,hcyl_,qopt_,quad,ierr,ntrial_,accuracy_,npts_) 1870 if (ierr/=0) then 1871 MSG_ERROR("Accuracy not reached") 1872 end if 1873 1874 c1=one/qpg_para_**2-COS(qpg_para_*hcyl_)/qpg_para_**2-hcyl_*SIN(qpg_para_*hcyl_)/qpg_para_ 1875 c2=SIN(qpg_para_*hcyl_)*SQRT(hcyl2+rcut2) 1876 vccut(ig,iq)=four_pi*c1+four_pi*(c2-quad)/qpg_para_ 1877 1878 else if (qpg_perp_/=zero.and.qpg_para_==zero) then 1879 ! $ 4pi\int_0^rcut d\rho \rho J_o(qpg_perp_.\rho) ln((h+\sqrt(h^2+\rho^2))/\rho) $ 1880 call quadrature(F4,zero,rcut_,qopt_,quad,ierr,ntrial_,accuracy_,npts_) 1881 if (ierr/=0) then 1882 MSG_ERROR("Accuracy not reached") 1883 end if 1884 1885 vccut(ig,iq)=four_pi*quad 1886 1887 else if (qpg_perp_==zero.and.qpg_para_==zero) then 1888 ! Use lim q+G --> 0 1889 vccut(ig,iq)=two_pi*(-hcyl2+hcyl_*SQRT(hcyl2+rcut2)+rcut2*LOG((hcyl_+SQRT(hcyl_+SQRT(hcyl2+rcut2)))/rcut_)) 1890 1891 else 1892 MSG_BUG('You should not be here!') 1893 end if 1894 1895 end do !ig 1896 end do !iq 1897 1898 else 1899 ! === Infinite cylinder === 1900 call wrtout(std_out,' cutoff_cylinder: using Rozzi''s method with infinite cylinder ','COLL') 1901 do iq=1,nq 1902 write(msg,'(a,i3)')' entering loop iq: ',iq 1903 call wrtout(std_out,msg,'COLL') 1904 qc(:)=qcart(:,iq) 1905 do ig=1,ng 1906 icount=ig+(iq-1)*ng ; if (icount<my_start.or.icount>my_stop) CYCLE 1907 gcart(:)=b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig) 1908 qpg(:)=qc(:)+gcart(:) 1909 qpg2 =DOT_PRODUCT(qpg,qpg) 1910 qpg_z =ABS(qpg(3)) ; qpg_xy=SQRT(qpg(1)**2+qpg(2)**2) 1911 if (qpg_z>SMALL) then 1912 ! === Analytic expression === 1913 call CALCK0(qpg_z *rcut_,k0,1) 1914 call CALJY1(qpg_xy*rcut_,j1,0) 1915 call CALJY0(qpg_xy*rcut_,j0,0) 1916 call CALCK1(qpg_z *rcut_,k1,1) 1917 vccut(iq,ig)=(four_pi/qpg2)*(one+rcut_*qpg_xy*j1*k0-qpg_z*rcut_*j0*k1) 1918 else 1919 if (qpg_xy>SMALL) then 1920 ! === Integrate r*Jo(G_xy r)log(r) from 0 up to rcut_ === 1921 call quadrature(F5,zero,rcut_,qopt_,quad,ierr,ntrial_,accuracy_,npts_) 1922 if (ierr/=0) then 1923 MSG_ERROR("Accuracy not reached") 1924 end if 1925 vccut(ig,iq)=-four_pi*quad 1926 else 1927 ! === Analytic expression === 1928 vccut(ig,iq)=-pi*rcut_**2*(two*LOG(rcut_)-one) 1929 end if 1930 end if 1931 end do !ig 1932 end do !iq 1933 end if !finite/infinite 1934 1935 CASE DEFAULT 1936 MSG_BUG(sjoin('Wrong value for method:',itoa(method))) 1937 END SELECT 1938 ! 1939 ! === Collect vccut on each node === 1940 call xmpi_sum(vccut,comm,ierr) 1941 1942 ABI_FREE(qcart) 1943 1944 end subroutine cutoff_cylinder
m_vcoul/cutoff_sphere [ Functions ]
[ Top ] [ m_vcoul ] [ Functions ]
NAME
cutoff_sphere
FUNCTION
Calculate the Fourier transform of the Coulomb interaction with a spherical cutoff:
$ v_{cut}(G)= \frac{4\pi}{|q+G|^2} [ 1-cos(|q+G|*R_cut) ] $ (1)
INPUTS
gmet(3,3)=Metric in reciprocal space. gvec(3,ngvec)=G vectors in reduced coordinates. rcut=Cutoff radius of the sphere. ngvec=Number of G vectors nqpt=Number of q-points qpt(3,nqpt)=q-points where the cutoff Coulomb is required.
OUTPUT
vc_cut(ngvec,nqpt)=Fourier components of the effective Coulomb interaction.
NOTES
For |q|<small and G=0 we use 2pi.R_cut^2, namely we consider the limit q-->0 of Eq. (1)
PARENTS
m_vcoul
CHILDREN
calck0,paw_jbessel,quadrature
SOURCE
1569 subroutine cutoff_sphere(nqpt,qpt,ngvec,gvec,gmet,rcut,vc_cut) 1570 1571 1572 !This section has been created automatically by the script Abilint (TD). 1573 !Do not modify the following lines by hand. 1574 #undef ABI_FUNC 1575 #define ABI_FUNC 'cutoff_sphere' 1576 !End of the abilint section 1577 1578 implicit none 1579 1580 !Arguments ------------------------------------ 1581 !scalars 1582 integer,intent(in) :: ngvec,nqpt 1583 real(dp),intent(in) :: rcut 1584 !arrays 1585 integer,intent(in) :: gvec(3,ngvec) 1586 real(dp),intent(in) :: gmet(3,3),qpt(3,nqpt) 1587 real(dp),intent(out) :: vc_cut(ngvec,nqpt) 1588 1589 !Local variables------------------------------- 1590 !scalars 1591 integer :: ig,igs,iqpt 1592 real(dp) :: qpg 1593 logical :: ltest 1594 1595 !************************************************************************ 1596 1597 !------------------------------------------------------------ 1598 ! Code modification by Bruno Rousseau, Montreal, 06/11/2013 1599 !------------------------------------------------------------ 1600 ! 1601 ! In order for the code below to run in parallel using MPI 1602 ! within the gwls code (by Laflamme-Jansen, Cote and Rousseau) 1603 ! the ABI_CHECK below must be removed. 1604 ! 1605 ! The ABI_CHECK below will fail if G-vectors are distributed on 1606 ! many processors, as only the master process has G=(0,0,0). 1607 ! 1608 ! The test does not seem necessary; IF a process has G=(0,0,0) 1609 ! (namely, the master process), then it will be the first G vector. 1610 ! If a process does not have G=(0,0,0) (ie, all the other processes), 1611 ! then they don't need to worry about the G-> 0 limit. 1612 ! 1613 1614 ltest=ALL(gvec(:,1)==0) 1615 !ABI_CHECK(ltest,'The first G vector should be Gamma') 1616 1617 do iqpt=1,nqpt 1618 igs=1 1619 if (ltest .and. normv(qpt(:,iqpt),gmet,'G')<GW_TOLQ0) then ! For small q and G=0, use the limit q-->0. 1620 vc_cut(1,iqpt)=two_pi*rcut**2 1621 igs=2 1622 end if 1623 do ig=igs,ngvec 1624 qpg=normv(qpt(:,iqpt)+gvec(:,ig),gmet,'G') 1625 vc_cut(ig,iqpt)=four_pi*(one-COS(rcut*qpg))/qpg**2 1626 end do 1627 end do 1628 1629 end subroutine cutoff_sphere
m_vcoul/cutoff_surface [ Functions ]
[ Top ] [ m_vcoul ] [ Functions ]
NAME
cutoff_surface
FUNCTION
Calculate the Fourier components of an effective Coulomb interaction within a slab of thickness 2*rcut which is symmetric with respect to the xy plane. In this implementation rcut=L_z/2 where L_z is the periodicity along z
INPUTS
gprimd(3,3)=Dimensional primitive translations in reciprocal space ($\textrm{bohr}^{-1}$). gvec(3,ng)=G vectors in reduced coordinates. ng=Number of G vectors. qpt(3,nq)=q-points nq=Number of q-points gmet(3,3)=Metric in reciprocal space.
OUTPUT
vc_cut(ng,nq)=Fourier components of the effective Coulomb interaction.
NOTES
The Fourier expression for an interaction truncated along the z-direction
(i.e non-zero only if |z|<R) is :
vc(q.G) = 4pi/|q+G|^2 * [ 1 + e^{-((q+G)_xy)*R} * ( (q_z+G_z)/(q+G)_xy * sin((q_z+G_z)R) -
- cos((q_z+G_Z)R)) ] (1)
Equation (1) diverges when q_xy+G_xy --> 0 for any non zero q_z+G_z
However if we choose R=L/2, where L defines the periodicity along z,
and we limit ourselves to consider q-points such as q_z==0, then
sin((q_z+G_z)R)=sin(G_Z 2pi/L)=0 for every G.
Under these assumptions we obtain
v(q,G) = 4pi/|q+G|^2 [1-e^{-(q+G)_xy*L/2}\cos((q_z+G_z)R)]
which is always finite when G_z /=0 while it diverges as 4piR/(q+G)_xy as (q+G)_xy -->0
but only in the x-y plane.
PARENTS
m_vcoul
CHILDREN
calck0,paw_jbessel,quadrature
SOURCE
1995 subroutine cutoff_surface(nq,qpt,ng,gvec,gprimd,rcut,boxcenter,pdir,alpha,vc_cut,method) 1996 1997 1998 !This section has been created automatically by the script Abilint (TD). 1999 !Do not modify the following lines by hand. 2000 #undef ABI_FUNC 2001 #define ABI_FUNC 'cutoff_surface' 2002 !End of the abilint section 2003 2004 implicit none 2005 2006 !Arguments ------------------------------------ 2007 !scalars 2008 integer,intent(in) :: method,ng,nq 2009 real(dp),intent(in) :: rcut 2010 !arrays 2011 integer,intent(in) :: gvec(3,ng),pdir(3) 2012 real(dp),intent(in) :: alpha(3),boxcenter(3),gprimd(3,3),qpt(3,nq) 2013 real(dp),intent(out) :: vc_cut(ng,nq) 2014 2015 !Local variables------------------------------- 2016 !scalars 2017 integer :: ig,iq,igs 2018 real(dp),parameter :: SMALL=tol4 !@WC: was tol6 2019 real(dp) :: qpg2,qpg_para,qpg_perp 2020 character(len=500) :: msg 2021 !arrays 2022 real(dp) :: b1(3),b2(3),b3(3),gcart(3),qc(3),qpg(3) 2023 real(dp),allocatable :: qcart(:,:) 2024 2025 ! ************************************************************************* 2026 2027 ABI_UNUSED(pdir) 2028 ABI_UNUSED(boxcenter) 2029 2030 ! === From reduced to cartesian coordinates === 2031 b1(:)=two_pi*gprimd(:,1) 2032 b2(:)=two_pi*gprimd(:,2) 2033 b3(:)=two_pi*gprimd(:,3) 2034 ABI_MALLOC(qcart,(3,nq)) 2035 do iq=1,nq 2036 qcart(:,iq) = b1*qpt(1,iq) + b2*qpt(2,iq) + b3*qpt(3,iq) 2037 end do 2038 2039 ! === Different approaches according to method === 2040 vc_cut(:,:)=zero 2041 2042 SELECT CASE (method) 2043 2044 CASE (1) ! Beigi's expression. 2045 ! * q-points with non-zero component along the z-axis are not allowed if 2046 ! the simplified Eq.1 for the Coulomb interaction is used. 2047 if (ANY(ABS(qcart(3,:))>SMALL)) then 2048 write(std_out,*)qcart(:,:) 2049 write(msg,'(5a)')& 2050 & 'Found q-points with non-zero component along non-periodic direction ',ch10,& 2051 & 'This is not allowed, see Notes in cutoff_surface.F90 ',ch10,& 2052 & 'ACTION : Modify the q-point sampling ' 2053 MSG_ERROR(msg) 2054 end if 2055 ! 2056 ! === Calculate truncated Coulomb interaction for a infinite surface === 2057 ! * Here I suppose that all the input q-points are different from zero 2058 do iq=1,nq 2059 qc(:)=qcart(:,iq) 2060 igs=1; if (SQRT(DOT_PRODUCT(qc,qc))<tol16) igs=2 ! avoid (q=0, G=0) 2061 do ig=igs,ng 2062 gcart(:)=b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig) 2063 qpg(:)=qc(:)+gcart(:) 2064 qpg2 =DOT_PRODUCT(qpg(:),qpg(:)) 2065 qpg_para=SQRT(qpg(1)**2+qpg(2)**2) ; qpg_perp=qpg(3) 2066 ! if (abs(qpg_perp)<SMALL.and.qpg_para<SMALL) stop 'SMALL in cutoff_surface 2067 vc_cut(ig,iq)=four_pi/qpg2*(one-EXP(-qpg_para*rcut)*COS(qpg_perp*rcut)) 2068 end do 2069 end do 2070 2071 CASE (2) ! Rozzi"s method 2072 MSG_ERROR("Work in progress") 2073 ABI_UNUSED(alpha) ! just to keep alpha as an argument 2074 !alpha=?? ; ap1sqrt=SQRT(one+alpha**2) 2075 do iq=1,nq 2076 qc(:)=qcart(:,iq) 2077 do ig=1,ng 2078 gcart(:)=b1(:)*gvec(1,ig)+b2(:)*gvec(2,ig)+b3(:)*gvec(3,ig) 2079 qpg(:)=qc(:)+gcart(:) 2080 qpg2 =DOT_PRODUCT(qpg(:),qpg(:)) 2081 qpg_para=SQRT(qpg(1)**2+qpg(2)**2) ; qpg_perp =qpg(3) 2082 if (qpg_para>SMALL) then 2083 vc_cut(ig,iq)=four_pi/qpg2*(one+EXP(-qpg_para*rcut)*(qpg_perp/qpg_para*SIN(qpg_perp*rcut)-COS(qpg_perp*rcut))) 2084 else 2085 if (ABS(qpg_perp)>SMALL) then 2086 vc_cut(ig,iq)=four_pi/qpg_perp**2*(one-COS(qpg_perp*rcut)-qpg_perp*rcut*SIN(qpg_perp*rcut)) !& 2087 !& + 8*rcut*SIN(qpg_perp*rcut)/qpg_perp*LOG((alpha+ap1sqrt)*(one+ap1sqrt)/alpha) ! contribution due to finite surface 2088 else 2089 vc_cut(ig,iq)=-two_pi*rcut**2 2090 end if 2091 end if 2092 end do !ig 2093 end do !iq 2094 2095 CASE DEFAULT 2096 write(msg,'(a,i3)')' Wrong value of method: ',method 2097 MSG_BUG(msg) 2098 END SELECT 2099 2100 ABI_FREE(qcart) 2101 2102 end subroutine cutoff_surface
m_vcoul/vcoul_free [ Functions ]
[ Top ] [ m_vcoul ] [ Functions ]
NAME
vcoul_free
FUNCTION
Destroy a vcoul_t type
SIDE EFFECTS
Vcp<vcoul_t>=the datatype to be destroyed
PARENTS
bethe_salpeter,gwls_hamiltonian,mrgscr,screening,setup_sigma,sigma
CHILDREN
calck0,paw_jbessel,quadrature
SOURCE
1497 subroutine vcoul_free(Vcp) 1498 1499 1500 !This section has been created automatically by the script Abilint (TD). 1501 !Do not modify the following lines by hand. 1502 #undef ABI_FUNC 1503 #define ABI_FUNC 'vcoul_free' 1504 !End of the abilint section 1505 1506 implicit none 1507 1508 !Arguments ------------------------------------ 1509 !scalars 1510 type(vcoul_t),intent(inout) :: Vcp 1511 1512 ! ************************************************************************* 1513 1514 DBG_ENTER("COLL") 1515 1516 if (allocated(Vcp%qibz)) then 1517 ABI_FREE(Vcp%qibz) 1518 end if 1519 if (allocated(Vcp%qlwl)) then 1520 ABI_FREE(Vcp%qlwl) 1521 end if 1522 if (allocated(Vcp%vc_sqrt)) then 1523 ABI_FREE(Vcp%vc_sqrt) 1524 end if 1525 if (allocated(Vcp%vc_sqrt_resid)) then 1526 ABI_FREE(Vcp%vc_sqrt_resid) 1527 end if 1528 if (allocated(Vcp%vcqlwl_sqrt)) then 1529 ABI_FREE(Vcp%vcqlwl_sqrt) 1530 end if 1531 1532 DBG_EXIT("COLL") 1533 1534 end subroutine vcoul_free
m_vcoul/vcoul_init [ Functions ]
[ Top ] [ m_vcoul ] [ Functions ]
NAME
vcoul_init
FUNCTION
Perform general check and initialize the data type containing information on the cutoff technique Note Vcp%vc_sqrt_resid and Vcp%i_sz_resid are simply initialized at the same value as Vcp%vc_sqrt and Vcp%i_sz
INPUTS
Qmesh<kmesh_t>=Info on the q-point sampling. Kmesh<kmesh_t>=Info on the k-point sampling. Gsph<gsphere_t>=Info of the G sphere. %gmet(3,3)=Metric in reciprocal space. %gprimd(3,3)=Dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$) %gvec=G vectors rcut=Cutoff radius for the cylinder. icutcoul=Option of the cutoff technique. vcutgeo(3)= Info on the orientation and extension of the cutoff region. ng=Number of G-vectors to be used to describe the Coulomb interaction nqlwl=Number of point around Gamma for treatment of long-wavelength limit qlwl(3,nqlwl)= The nqlwl "small" q-points ngfft(18)=Information on the (fine) FFT grid used for the density. rprimd(3,3)=Direct lattice vectors. comm=MPI communicator.
OUTPUT
Vcp=<vcoul_t>=Datatype gathering information on the Coulomb interaction.
PARENTS
gwls_hamiltonian,mrgscr,setup_bse,setup_bse_interp,setup_screening setup_sigma
CHILDREN
calck0,paw_jbessel,quadrature
SOURCE
208 subroutine vcoul_init(Vcp,Gsph,Cryst,Qmesh,Kmesh,rcut,icutcoul,vcutgeo,ecut,ng,nqlwl,qlwl,ngfft,comm) 209 210 211 !This section has been created automatically by the script Abilint (TD). 212 !Do not modify the following lines by hand. 213 #undef ABI_FUNC 214 #define ABI_FUNC 'vcoul_init' 215 !End of the abilint section 216 217 implicit none 218 219 !Arguments ------------------------------------ 220 !scalars 221 integer,intent(in) :: ng,nqlwl,icutcoul,comm 222 real(dp),intent(in) :: rcut, ecut 223 type(kmesh_t),intent(in) :: Kmesh,Qmesh 224 type(gsphere_t),intent(in) :: Gsph 225 type(crystal_t),intent(in) :: Cryst 226 type(vcoul_t),intent(out) :: Vcp 227 !arrays 228 integer,intent(in) :: ngfft(18) 229 real(dp),intent(in) :: qlwl(3,nqlwl),vcutgeo(3) 230 231 !Local variables------------------------------- 232 !scalars 233 integer,parameter :: master=0,ncell=3 234 integer :: nmc_max=2500000 235 integer :: nmc,nseed 236 integer :: i1,i2,i3,ig,imc 237 integer :: ii,iqlwl,iq_bz,iq_ibz,npar,npt 238 integer :: opt_cylinder,opt_surface,test,rank,nprocs 239 integer, allocatable :: seed(:) 240 real(dp),parameter :: tolq0=1.d-3 241 real(dp) :: b1b1,b2b2,b3b3,b1b2,b2b3,b3b1 242 real(dp) :: bz_geometry_factor,bz_plane,check,dx,integ,q0_vol,q0_volsph 243 real(dp) :: qbz_norm,step,ucvol,intfauxgb, alfa 244 character(len=500) :: msg 245 !arrays 246 integer :: gamma_pt(3,1) 247 real(dp) :: a1(3),a2(3),a3(3),bb(3),b1(3),b2(3),b3(3),gmet(3,3),gprimd(3,3) 248 real(dp) :: qbz_cart(3),rmet(3,3) 249 real(dp),allocatable :: cov(:,:),par(:),qfit(:,:),sigma(:),var(:),qcart(:,:) 250 real(dp),allocatable :: vcfit(:,:),vcoul(:,:),vcoul_lwl(:,:),xx(:),yy(:) 251 real(dp),allocatable :: qran(:,:) 252 real(dp) :: lmin,vlength,qpg2,rcut2 253 real(dp) :: q0sph 254 real(dp) :: qtmp(3),qmin(3),qmin_cart(3),qpg(3) 255 real(dp) :: rprimd_sc(3,3),gprimd_sc(3,3),gmet_sc(3,3),rmet_sc(3,3),ucvol_sc 256 real(dp) :: qcart2red(3,3) 257 real(dp) :: qqgg(3) 258 259 ! ************************************************************************* 260 261 DBG_ENTER("COLL") 262 ! 263 ! === Test if the first q-point is zero === 264 ! FIXME this wont work if nqptdm/=0 265 !if (normv(Qmesh%ibz(:,1),gmet,'G')<GW_TOLQ0)) STOP 'vcoul_init, non zero first point ' 266 rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 267 268 call metric(gmet,gprimd,-1,rmet,Cryst%rprimd,ucvol) 269 ! 270 ! === Save dimension and other useful quantities in Vcp% === 271 Vcp%nfft = PRODUCT(ngfft(1:3)) ! Number of points in the FFT mesh. 272 Vcp%ng = ng ! Number of G-vectors in the Coulomb matrix elements. 273 Vcp%nqibz = Qmesh%nibz ! Number of irred q-point. 274 Vcp%nqlwl = nqlwl ! Number of small q-directions to deal with singularity and non Analytic behavior. 275 Vcp%rcut = rcut ! Cutoff radius for cylinder. 276 Vcp%hcyl = zero ! Length of finite cylinder (Rozzi"s method, default is Beigi). 277 Vcp%ucvol = ucvol ! Unit cell volume. 278 279 Vcp%rprimd = Cryst%rprimd(:,:) ! Dimensional direct lattice. 280 Vcp%boxcenter = zero ! boxcenter at the moment is supposed to be at the origin. 281 Vcp%vcutgeo = vcutgeo(:) ! Info on the orientation and extension of the cutoff region. 282 Vcp%ngfft(:) = ngfft(:) ! Info on the FFT mesh. 283 ! 284 ! === Define geometry and cutoff radius (if used) === 285 Vcp%mode='NONE' 286 if (icutcoul==14) Vcp%mode='MINIBZ-ERF' 287 if (icutcoul==15) Vcp%mode='MINIBZ-ERFC' 288 if (icutcoul==16) Vcp%mode='MINIBZ' 289 if (icutcoul==0) Vcp%mode='SPHERE' 290 if (icutcoul==1) Vcp%mode='CYLINDER' 291 if (icutcoul==2) Vcp%mode='SURFACE' 292 if (icutcoul==3) Vcp%mode='CRYSTAL' 293 if (icutcoul==4) Vcp%mode='ERF' 294 if (icutcoul==5) Vcp%mode='ERFC' 295 if (icutcoul==6) Vcp%mode='AUXILIARY_FUNCTION' 296 if (icutcoul==7) Vcp%mode='AUX_GB' 297 ! 298 ! === Calculate Fourier coefficient of Coulomb interaction === 299 ! * For the limit q-->0 we consider ng vectors due to a possible anisotropy in case of a cutoff interaction 300 ABI_MALLOC(Vcp%qibz,(3,Vcp%nqibz)) 301 Vcp%qibz = Qmesh%ibz(:,:) 302 ABI_MALLOC(Vcp%qlwl,(3,Vcp%nqlwl)) 303 Vcp%qlwl = qlwl(:,:) 304 305 ! =============================================== 306 ! == Calculation of the FT of the Coulomb term == 307 ! =============================================== 308 ABI_MALLOC(vcoul ,(ng,Vcp%nqibz)) 309 ABI_MALLOC(vcoul_lwl,(ng,Vcp%nqlwl)) 310 311 a1=Cryst%rprimd(:,1); b1=two_pi*gprimd(:,1) 312 a2=Cryst%rprimd(:,2); b2=two_pi*gprimd(:,2) 313 a3=Cryst%rprimd(:,3); b3=two_pi*gprimd(:,3) 314 b1b1=dot_product(b1,b1) ; b2b2=dot_product(b2,b2) ; b3b3=dot_product(b3,b3) 315 bb(1)=b1b1 ; bb(2)=b2b2 ; bb(3)=b3b3 316 b1b2=dot_product(b1,b2) ; b2b3=dot_product(b2,b3) ; b3b1=dot_product(b3,b1) 317 318 SELECT CASE (TRIM(Vcp%mode)) 319 320 CASE ('MINIBZ','MINIBZ-ERFC','MINIBZ-ERF') 321 322 ! Mimicking the BerkeleyGW technique 323 ! A Monte-Carlo sampling of each miniBZ surrounding each (q+G) point 324 ! However - extended to the multiple shifts 325 ! - with an adaptative number of MonteCarlo sampling points 326 327 ! 328 ! Supercell defined by the k-points 329 rprimd_sc(:,:) = MATMUL(Cryst%rprimd,Kmesh%kptrlatt) 330 call metric(gmet_sc,gprimd_sc,-1,rmet_sc,rprimd_sc,ucvol_sc) 331 332 qcart2red(:,:) = two_pi * Cryst%gprimd(:,:) 333 call matrginv(qcart2red,3,3) 334 335 ! 336 ! Find the largest sphere inside the miniBZ 337 ! in order to integrate the divergence analytically 338 q0sph = HUGE(1.0_dp) 339 do i1 = -ncell+1, ncell 340 qtmp(1) = dble(i1) * 0.5_dp 341 do i2 = -ncell+1, ncell 342 qtmp(2) = dble(i2) * 0.5_dp 343 do i3 = -ncell+1, ncell 344 qtmp(3) = dble(i3) * 0.5_dp 345 if( i1==0 .AND. i2==0 .AND. i3==0 ) cycle 346 vlength = normv(qtmp,gmet_sc,'G') 347 if (vlength < q0sph) then 348 q0sph = vlength 349 end if 350 enddo 351 enddo 352 enddo 353 354 355 ! 356 ! Setup the random vectors for the Monte Carlo sampling of the miniBZ 357 ! at q=0 358 ABI_MALLOC(qran,(3,nmc_max)) 359 call random_seed(size=nseed) 360 ABI_MALLOC(seed,(nseed)) 361 do i1=1,nseed 362 seed(i1) = NINT(SQRT(DBLE(i1)*103731)) 363 end do 364 call random_seed(put=seed) 365 call random_number(qran) 366 367 ! Overide the first "random vector" with 0 368 qran(:,1) = 0.0_dp 369 370 ! Fold qran into the Wignez-Seitz cell around q=0 371 do imc=2,nmc_max 372 lmin = HUGE(1.0_dp) 373 do i1 = -ncell+1, ncell 374 qtmp(1) = qran(1,imc) + dble(i1) 375 do i2 = -ncell+1, ncell 376 qtmp(2) = qran(2,imc) + dble(i2) 377 do i3 = -ncell+1, ncell 378 qtmp(3) = qran(3,imc) + dble(i3) 379 vlength = normv(qtmp,gmet_sc,'G') 380 if (vlength < lmin) then 381 lmin = vlength 382 ! Get the q-vector in cartersian coordinates 383 qmin_cart(:) = two_pi * MATMUL( gprimd_sc(:,:) , qtmp ) 384 ! Transform it back to the reciprocal space 385 qmin(:) = MATMUL( qcart2red , qmin_cart ) 386 end if 387 enddo 388 enddo 389 enddo 390 391 qran(1,imc) = qmin(1) 392 qran(2,imc) = qmin(2) 393 qran(3,imc) = qmin(3) 394 395 enddo 396 397 398 rcut2 = Vcp%rcut**2 399 vcoul(:,:) = zero 400 401 ! 402 ! Admittedly, it's a lot of duplication of code hereafter, 403 ! but I think it's the clearest and fastest way 404 ! 405 select case( TRIM(Vcp%mode) ) 406 case('MINIBZ') 407 408 do iq_ibz=1,Vcp%nqibz 409 do ig=1,ng 410 if( iq_ibz==1 .AND. ig==1 ) cycle 411 412 qpg(:) = Vcp%qibz(:,iq_ibz) + Gsph%gvec(:,ig) 413 qpg2 = normv(qpg,gmet,'G')**2 414 nmc = adapt_nmc(nmc_max,qpg2) 415 do imc=1,nmc 416 qpg(:) = Vcp%qibz(:,iq_ibz) + Gsph%gvec(:,ig) + qran(:,imc) 417 qpg2 = normv(qpg,gmet,'G')**2 418 vcoul(ig,iq_ibz) = vcoul(ig,iq_ibz) + four_pi / qpg2 / REAL(nmc,dp) 419 end do 420 end do 421 end do 422 423 ! Override q=1, ig=1 424 vcoul(1,1) = four_pi**2 * Kmesh%nbz * ucvol / ( 8.0_dp * pi**3 ) * q0sph 425 do imc=1,nmc_max 426 qpg(:) = Vcp%qibz(:,1) + Gsph%gvec(:,1) + qran(:,imc) 427 qpg2 = normv(qpg,gmet,'G')**2 428 if( qpg2 > q0sph**2 ) then 429 vcoul(1,1) = vcoul(1,1) + four_pi / qpg2 / REAL(nmc_max,dp) 430 end if 431 end do 432 433 vcoul_lwl(:,:) = zero 434 do iq_ibz=1,Vcp%nqlwl 435 do ig=1,ng 436 qpg(:) = Vcp%qlwl(:,iq_ibz) + Gsph%gvec(:,ig) 437 qpg2 = normv(qpg,gmet,'G')**2 438 nmc = adapt_nmc(nmc_max,qpg2) 439 do imc=1,nmc 440 qpg(:) = Vcp%qlwl(:,iq_ibz) + Gsph%gvec(:,ig) + qran(:,imc) 441 qpg2 = normv(qpg,gmet,'G')**2 442 vcoul_lwl(ig,iq_ibz) = vcoul_lwl(ig,iq_ibz) + four_pi / qpg2 / REAL(nmc,dp) 443 end do 444 end do 445 end do 446 447 448 case('MINIBZ-ERFC') 449 450 do iq_ibz=1,Vcp%nqibz 451 do ig=1,ng 452 if( iq_ibz==1 .AND. ig==1 ) cycle 453 454 qpg(:) = Vcp%qibz(:,iq_ibz) + Gsph%gvec(:,ig) 455 qpg2 = normv(qpg,gmet,'G')**2 456 nmc = adapt_nmc(nmc_max,qpg2) 457 do imc=1,nmc 458 qpg(:) = Vcp%qibz(:,iq_ibz) + Gsph%gvec(:,ig) + qran(:,imc) 459 qpg2 = normv(qpg,gmet,'G')**2 460 vcoul(ig,iq_ibz) = vcoul(ig,iq_ibz) + four_pi / qpg2 / REAL(nmc,dp) & 461 & * ( one - EXP( -0.25d0 * rcut2 * qpg2 ) ) 462 end do 463 end do 464 end do 465 466 ! Override q=1, ig=1 467 vcoul(1,1) = four_pi**2 * Kmesh%nbz * ucvol / ( 8.0_dp * pi**3 ) & 468 & * ( q0sph - SQRT(pi/rcut2) * abi_derf(0.5_dp*SQRT(rcut2)*q0sph) ) 469 do imc=1,nmc_max 470 qpg(:) = Vcp%qibz(:,1) + Gsph%gvec(:,1) + qran(:,imc) 471 qpg2 = normv(qpg,gmet,'G')**2 472 if( qpg2 > q0sph**2 ) then 473 vcoul(1,1) = vcoul(1,1) + four_pi / qpg2 / REAL(nmc_max,dp) & 474 & * ( one - EXP( -0.25d0 * rcut2 * qpg2 ) ) 475 end if 476 end do 477 478 vcoul_lwl(:,:) = zero 479 do iq_ibz=1,Vcp%nqlwl 480 do ig=1,ng 481 qpg(:) = Vcp%qlwl(:,iq_ibz) + Gsph%gvec(:,ig) 482 qpg2 = normv(qpg,gmet,'G')**2 483 nmc = adapt_nmc(nmc_max,qpg2) 484 do imc=1,nmc 485 qpg(:) = Vcp%qlwl(:,iq_ibz) + Gsph%gvec(:,ig) + qran(:,imc) 486 qpg2 = normv(qpg,gmet,'G')**2 487 vcoul_lwl(ig,iq_ibz) = vcoul_lwl(ig,iq_ibz) + four_pi / qpg2 / REAL(nmc,dp) & 488 & * ( one - EXP( -0.25d0 * rcut2 * qpg2 ) ) 489 end do 490 end do 491 end do 492 493 494 case('MINIBZ-ERF') 495 496 do iq_ibz=1,Vcp%nqibz 497 do ig=1,ng 498 if( iq_ibz==1 .AND. ig==1 ) cycle 499 500 qpg(:) = Vcp%qibz(:,iq_ibz) + Gsph%gvec(:,ig) 501 qpg2 = normv(qpg,gmet,'G')**2 502 nmc = adapt_nmc(nmc_max,qpg2) 503 do imc=1,nmc 504 qpg(:) = Vcp%qibz(:,iq_ibz) + Gsph%gvec(:,ig) + qran(:,imc) 505 qpg2 = normv(qpg,gmet,'G')**2 506 vcoul(ig,iq_ibz) = vcoul(ig,iq_ibz) + four_pi / qpg2 / REAL(nmc,dp) & 507 & * EXP( -0.25d0 * rcut2 * qpg2 ) 508 end do 509 end do 510 end do 511 512 ! Override q=1, ig=1 513 vcoul(1,1) = four_pi**2 * Kmesh%nbz * ucvol / ( 8.0_dp * pi**3 ) & 514 & * SQRT(pi/rcut2) * abi_derf(0.5_dp*SQRT(rcut2)*q0sph) 515 516 do imc=1,nmc_max 517 qpg(:) = Vcp%qibz(:,1) + Gsph%gvec(:,1) + qran(:,imc) 518 qpg2 = normv(qpg,gmet,'G')**2 519 if( qpg2 > q0sph**2 ) then 520 vcoul(1,1) = vcoul(1,1) + four_pi / qpg2 / REAL(nmc_max,dp) * EXP( -0.25d0 * rcut2 * qpg2 ) 521 end if 522 end do 523 524 vcoul_lwl(:,:) = zero 525 do iq_ibz=1,Vcp%nqlwl 526 do ig=1,ng 527 qpg(:) = Vcp%qlwl(:,iq_ibz) + Gsph%gvec(:,ig) 528 qpg2 = normv(qpg,gmet,'G')**2 529 nmc = adapt_nmc(nmc_max,qpg2) 530 do imc=1,nmc 531 qpg(:) = Vcp%qlwl(:,iq_ibz) + Gsph%gvec(:,ig) + qran(:,imc) 532 qpg2 = normv(qpg,gmet,'G')**2 533 vcoul_lwl(ig,iq_ibz) = vcoul_lwl(ig,iq_ibz) + four_pi / qpg2 / REAL(nmc,dp) & 534 & * EXP( -0.25d0 * rcut2 * qpg2 ) 535 end do 536 end do 537 end do 538 539 end select 540 541 Vcp%i_sz = vcoul(1,1) 542 543 ABI_FREE(qran) 544 ABI_FREE(seed) 545 546 CASE ('SPHERE') 547 ! TODO check that L-d > R_c > d 548 ! * A non-positive value of rcut imposes the recipe of Spencer & Alavi, PRB 77, 193110 (2008) [[cite:Spencer2008]]. 549 if (Vcp%rcut<tol12) then 550 Vcp%rcut = (ucvol*Kmesh%nbz*3.d0/four_pi)**third 551 write(msg,'(2a,2x,f8.4,a)')ch10,& 552 & ' Using a calculated value for rcut = ',Vcp%rcut, ' to have the same volume as the BvK crystal ' 553 call wrtout(std_out,msg,'COLL') 554 end if 555 556 Vcp%vcutgeo=zero 557 call cutoff_sphere(Qmesh%nibz,Qmesh%ibz,ng,Gsph%gvec,gmet,Vcp%rcut,vcoul) 558 559 ! q-points for optical limit. 560 call cutoff_sphere(Vcp%nqlwl,Vcp%qlwl,ng,Gsph%gvec,gmet,Vcp%rcut,vcoul_lwl) 561 ! 562 ! === Treat the limit q--> 0 === 563 ! * The small cube is approximated by a sphere, while vc(q=0)=2piR**2. 564 ! * if a single q-point is used, the expression for the volume is exact. 565 Vcp%i_sz=two_pi*Vcp%rcut**2 566 call vcoul_print(Vcp,unit=ab_out) 567 568 CASE ('CYLINDER') 569 test=COUNT(ABS(Vcp%vcutgeo)>tol6) 570 ABI_CHECK(test==1,'Wrong cutgeo for cylinder') 571 572 ! === Beigi method is the default one, i.e infinite cylinder of radius rcut === 573 ! * Negative values to use Rozzi method with finite cylinder of extent hcyl. 574 opt_cylinder=1; Vcp%hcyl=zero; Vcp%pdir(:)=0 575 do ii=1,3 576 check=Vcp%vcutgeo(ii) 577 if (ABS(check)>tol6) then 578 Vcp%pdir(ii)=1 579 if (check<zero) then ! use Rozzi's method. 580 Vcp%hcyl=ABS(check)*SQRT(SUM(Cryst%rprimd(:,ii)**2)) 581 opt_cylinder=2 582 end if 583 end if 584 end do 585 586 test=COUNT(Vcp%pdir==1) 587 ABI_CHECK((test==1),'Wrong pdir for cylinder') 588 if (Vcp%pdir(3)/=1) then 589 MSG_ERROR("The cylinder must be along the z-axis") 590 end if 591 592 call cutoff_cylinder(Qmesh%nibz,Qmesh%ibz,ng,Gsph%gvec,Vcp%rcut,Vcp%hcyl,Vcp%pdir,& 593 & Vcp%boxcenter,Cryst%rprimd,vcoul,opt_cylinder,comm) 594 595 ! q-points for optical limit. 596 call cutoff_cylinder(Vcp%nqlwl,Vcp%qlwl,ng,Gsph%gvec,Vcp%rcut,Vcp%hcyl,Vcp%pdir,& 597 & Vcp%boxcenter,Cryst%rprimd,vcoul_lwl,opt_cylinder,comm) 598 ! 599 ! === If Beigi, treat the limit q--> 0 === 600 if (opt_cylinder==1) then 601 npar=8; npt=100 ; gamma_pt=RESHAPE((/0,0,0/),(/3,1/)) 602 ABI_MALLOC(qfit,(3,npt)) 603 ABI_MALLOC(vcfit,(1,npt)) 604 if (Qmesh%nibz==1) then 605 MSG_ERROR("nqibz == 1 not supported when Beigi's method is used") 606 endif 607 qfit(:,:)=zero 608 step=half/(npt*(Qmesh%nibz-1)) ; qfit(3,:)=arth(tol6,step,npt) 609 !step=(half/(Qmesh%nibz-1)/tol6)**(one/npt) ; qfit(3,:)=geop(tol6,step,npt) 610 call cutoff_cylinder(npt,qfit,1,gamma_pt,Vcp%rcut,Vcp%hcyl,Vcp%pdir,Vcp%boxcenter,Cryst%rprimd,vcfit,opt_cylinder,comm) 611 612 ABI_MALLOC(xx,(npt)) 613 ABI_MALLOC(yy,(npt)) 614 ABI_MALLOC(sigma,(npt)) 615 ABI_MALLOC(par,(npar)) 616 ABI_MALLOC(var,(npar)) 617 ABI_MALLOC(cov,(npar,npar)) 618 do ii=1,npt 619 xx(ii)=normv(qfit(:,ii),gmet,'G') 620 end do 621 ABI_FREE(qfit) 622 sigma=one ; yy(:)=vcfit(1,:) 623 ABI_FREE(vcfit) 624 !call llsfit_svd(xx,yy,sigma,npar,K0fit,chisq,par,var,cov,info) 625 !do ii=1,npt 626 ! write(99,*)xx(ii),yy(ii),DOT_PRODUCT(par,K0fit(xx(ii),npar)) 627 !end do 628 bz_plane=l2norm(b1.x.b2) 629 !integ=K0fit_int(xx(npt),par,npar) 630 !write(std_out,*)' SVD fit : chi-square',chisq 631 !write(std_out,*)' fit-parameters : ',par 632 !write(std_out,*)' variance ',var 633 !write(std_out,*)' bz_plane ',bz_plane 634 !write(std_out,*)' SCD integ ',integ 635 ! Here Im assuming homogeneous mesh 636 dx=(xx(2)-xx(1)) 637 integ=yy(2)*dx*3.0/2.0 638 do ii=3,npt-2 639 integ=integ+yy(ii)*dx 640 end do 641 integ=integ+yy(npt-1)*dx*3.0/2.0 642 write(std_out,*)' simple integral',integ 643 q0_volsph=(two_pi)**3/(Kmesh%nbz*ucvol) 644 q0_vol=bz_plane*two*xx(npt) 645 write(std_out,*)' q0 sphere : ',q0_volsph,' q0_vol cyl ',q0_vol 646 Vcp%i_sz=bz_plane*two*integ/q0_vol 647 write(std_out,*)' spherical approximation ',four_pi*7.44*q0_volsph**(-two_thirds) 648 write(std_out,*)' Cylindrical cutoff value ',Vcp%i_sz 649 !Vcp%i_sz=four_pi*7.44*q0_vol**(-two_thirds) 650 ABI_FREE(xx) 651 ABI_FREE(yy) 652 ABI_FREE(sigma) 653 ABI_FREE(par) 654 ABI_FREE(var) 655 ABI_FREE(cov) 656 else 657 ! In Rozzi"s method the lim q+G --> 0 is finite. 658 Vcp%i_sz=vcoul(1,1) 659 end if 660 661 call vcoul_print(Vcp,unit=ab_out ) 662 663 CASE ('SURFACE') 664 665 test=COUNT(Vcp%vcutgeo/=zero) 666 ABI_CHECK(test==2,"Wrong vcutgeo") 667 ! 668 ! === Default is Beigi"s method === 669 opt_surface=1; Vcp%alpha(:)=zero 670 if (ANY(Vcp%vcutgeo<zero)) opt_surface=2 671 Vcp%pdir(:)=zero 672 do ii=1,3 673 check=Vcp%vcutgeo(ii) 674 if (ABS(check)>zero) then ! Use Rozzi"s method with a finite surface along x-y 675 Vcp%pdir(ii)=1 676 if (check<zero) Vcp%alpha(ii)=normv(check*Cryst%rprimd(:,ii),rmet,'R') 677 end if 678 end do 679 680 ! Beigi"s method: the surface must be along x-y and R must be L_Z/2. 681 if (opt_surface==1) then 682 ABI_CHECK(ALL(Vcp%pdir == (/1,1,0/)),"Surface must be in the x-y plane") 683 Vcp%rcut = half*SQRT(DOT_PRODUCT(a3,a3)) 684 end if 685 686 call cutoff_surface(Qmesh%nibz,Qmesh%ibz,ng,Gsph%gvec,gprimd,Vcp%rcut,& 687 & Vcp%boxcenter,Vcp%pdir,Vcp%alpha,vcoul,opt_surface) 688 689 ! q-points for optical limit. 690 call cutoff_surface(Vcp%nqlwl,Vcp%qlwl,ng,Gsph%gvec,gprimd,Vcp%rcut,& 691 & Vcp%boxcenter,Vcp%pdir,Vcp%alpha,vcoul_lwl,opt_surface) 692 ! 693 ! === If Beigi, treat the limit q--> 0 === 694 if (opt_surface==1) then 695 ! Integrate numerically in the plane close to 0 696 npt=100 ! Number of points in 1D 697 gamma_pt=RESHAPE((/0,0,0/),(/3,1/)) ! Gamma point 698 ABI_MALLOC(qfit,(3,npt)) 699 ABI_MALLOC(qcart,(3,npt)) 700 ABI_MALLOC(vcfit,(1,npt)) 701 if (Qmesh%nibz==1) then 702 MSG_ERROR("nqibz == 1 not supported when Beigi's method is used") 703 endif 704 qfit(:,:)=zero 705 qcart(:,:)=zero 706 ! Size of the third vector 707 bz_plane=l2norm(b3) 708 q0_volsph=(two_pi)**3/(Kmesh%nbz*ucvol) 709 ! radius that gives the same volume as q0_volsph 710 ! Let's assume that c is perpendicular to the plane 711 ! We also assume isotropic BZ around gamma 712 step=sqrt((q0_volsph/bz_plane)/pi)/npt 713 714 !step=half/(npt*(Qmesh%nibz-1)) 715 ! Let's take qpoints along 1 line, the vcut does depend only on the norm 716 qcart(1,:)=arth(tol6,step,npt) 717 718 do ii = 1,npt 719 qfit(:,ii) = MATMUL(TRANSPOSE(Cryst%rprimd),qcart(:,ii))/(2*pi) 720 end do 721 722 call cutoff_surface(npt,qfit,1,gamma_pt,gprimd,Vcp%rcut,& 723 & Vcp%boxcenter,Vcp%pdir,Vcp%alpha,vcfit,opt_surface) 724 725 ABI_MALLOC(xx,(npt)) 726 ABI_MALLOC(yy,(npt)) 727 ABI_MALLOC(sigma,(npt)) 728 do ii=1,npt 729 !xx(ii)=qfit(1,:) 730 xx(ii)=normv(qfit(:,ii),gmet,'G') 731 end do 732 ABI_FREE(qfit) 733 sigma=one 734 yy(:)=vcfit(1,:) 735 !yy(:)=one 736 ABI_FREE(vcfit) 737 dx=(xx(2)-xx(1)) 738 ! integ = \int dr r f(r) 739 integ=xx(2)*yy(2)*dx*3.0/2.0 740 do ii=3,npt-2 741 integ=integ+xx(ii)*yy(ii)*dx 742 end do 743 integ=integ+xx(npt-1)*yy(npt-1)*dx*3.0/2.0 744 write(std_out,*)' simple integral',integ 745 q0_vol=bz_plane*pi*xx(npt)**2 746 write(std_out,*)' q0 sphere : ',q0_volsph,' q0_vol cyl ',q0_vol 747 Vcp%i_sz=bz_plane*2*pi*integ/q0_vol 748 write(std_out,*)' spherical approximation ',four_pi*7.44*q0_volsph**(-two_thirds) 749 write(std_out,*)' Cylindrical cutoff value ',Vcp%i_sz 750 !Vcp%i_sz=four_pi*7.44*q0_vol**(-two_thirds) 751 ABI_FREE(xx) 752 ABI_FREE(yy) 753 else 754 ! In Rozzi"s method the lim q+G --> 0 is finite. 755 Vcp%i_sz=vcoul(1,1) 756 end if 757 758 CASE ('AUXILIARY_FUNCTION') 759 ! 760 ! Numerical integration of the exact-exchange divergence through the 761 ! auxiliary function of Carrier et al. PRB 75, 205126 (2007) [[cite:Carrier2007]]. 762 ! 763 do iq_ibz=1,Vcp%nqibz 764 call cmod_qpg(Vcp%nqibz,iq_ibz,Vcp%qibz,ng,Gsph%gvec,gprimd,vcoul(:,iq_ibz)) 765 766 if (iq_ibz==1) then ! The singularity is treated using vcoul_lwl. 767 vcoul(1, iq_ibz) = zero 768 vcoul(2:,iq_ibz) = four_pi/vcoul(2:,iq_ibz)**2 769 else 770 vcoul(:,iq_ibz) = four_pi/vcoul(:,iq_ibz)**2 771 end if 772 end do 773 774 ! q-points for optical limit. 775 do iqlwl=1,Vcp%nqlwl 776 call cmod_qpg(Vcp%nqlwl,iqlwl,Vcp%qlwl,ng,Gsph%gvec,gprimd,vcoul_lwl(:,iqlwl)) 777 end do 778 vcoul_lwl = four_pi/vcoul_lwl**2 779 ! 780 ! === Integration of 1/q^2 singularity === 781 ! * We use the auxiliary function from PRB 75, 205126 (2007) [[cite:Carrier2007]] 782 q0_vol=(two_pi)**3/(Kmesh%nbz*ucvol) ; bz_geometry_factor=zero 783 784 ! It might be useful to perform the integration using the routine quadrature 785 ! so that we can control the accuracy of the integral and improve the 786 ! numerical stability of the GW results. 787 do iq_bz=1,Qmesh%nbz 788 qbz_cart(:)=Qmesh%bz(1,iq_bz)*b1(:)+Qmesh%bz(2,iq_bz)*b2(:)+Qmesh%bz(3,iq_bz)*b3(:) 789 qbz_norm=SQRT(SUM(qbz_cart(:)**2)) 790 if (qbz_norm>tolq0) bz_geometry_factor=bz_geometry_factor-faux(Qmesh%bz(:,iq_bz)) 791 end do 792 793 bz_geometry_factor = bz_geometry_factor + integratefaux()*Qmesh%nbz 794 795 write(msg,'(2a,2x,f8.4)')ch10,& 796 & ' integrate q->0 : numerical BZ geometry factor = ',bz_geometry_factor*q0_vol**(2./3.) 797 call wrtout(std_out,msg,'COLL') 798 799 Vcp%i_sz=four_pi*bz_geometry_factor ! Final result stored here 800 801 CASE ('AUX_GB') 802 ! 803 do iq_ibz=1,Vcp%nqibz 804 call cmod_qpg(Vcp%nqibz,iq_ibz,Vcp%qibz,ng,Gsph%gvec,gprimd,vcoul(:,iq_ibz)) 805 806 if (iq_ibz==1) then ! The singularity is treated using vcoul_lwl. 807 vcoul(1, iq_ibz) = zero 808 vcoul(2:,iq_ibz) = four_pi/vcoul(2:,iq_ibz)**2 809 else 810 vcoul(:,iq_ibz) = four_pi/vcoul(:,iq_ibz)**2 811 end if 812 end do 813 814 ! q-points for optical limit. 815 do iqlwl=1,Vcp%nqlwl 816 call cmod_qpg(Vcp%nqlwl,iqlwl,Vcp%qlwl,ng,Gsph%gvec,gprimd,vcoul_lwl(:,iqlwl)) 817 end do 818 vcoul_lwl = four_pi/vcoul_lwl**2 819 ! 820 ! === Integration of 1/q^2 singularity === 821 ! * We use the auxiliary function of a Gygi-Baldereschi variant [[cite:Gigy1986]] 822 q0_vol=(two_pi)**3/(Kmesh%nbz*ucvol) ; bz_geometry_factor=zero 823 ! * the choice of alfa (the width of gaussian) is somehow empirical 824 alfa = 150.0/ecut 825 !write(msg,'(2a,2x,f12.4,2x,f12.4)')ch10, ' alfa, ecutsigx = ', alfa, ecut 826 !call wrtout(std_out,msg,'COLL') 827 828 do iq_bz=1,Qmesh%nbz 829 do ig = 1,ng 830 qqgg(:) = Qmesh%bz(:,iq_bz) + Gsph%gvec(:,ig) 831 qqgg(:) = qqgg(1)*b1(:) + qqgg(2)*b2(:) + qqgg(3)*b3(:) 832 qbz_norm=SQRT(SUM(qqgg(:)**2)) 833 if (qbz_norm>tolq0*0.01) then 834 bz_geometry_factor = bz_geometry_factor - EXP(-alfa*qbz_norm**2)/qbz_norm**2 835 end if 836 end do 837 end do 838 839 intfauxgb = ucvol/four_pi/SQRT(0.5*two_pi*alfa) 840 !write(msg,'(2a,2x,f12.4,2x,f12.4)')ch10, ' bz_geometry_factor, f(q) integral = ', & 841 !& bz_geometry_factor, intfauxgb 842 !call wrtout(std_out,msg,'COLL') 843 844 bz_geometry_factor = bz_geometry_factor + intfauxgb*Qmesh%nbz 845 846 write(msg,'(2a,2x,f8.4)')ch10,& 847 & ' integrate q->0 : numerical BZ geometry factor = ',bz_geometry_factor*q0_vol**(2./3.) 848 call wrtout(std_out,msg,'COLL') 849 850 Vcp%i_sz=four_pi*bz_geometry_factor 851 852 CASE ('CRYSTAL') 853 854 do iq_ibz=1,Vcp%nqibz 855 call cmod_qpg(Vcp%nqibz,iq_ibz,Vcp%qibz,ng,Gsph%gvec,gprimd,vcoul(:,iq_ibz)) 856 857 if (iq_ibz==1) then ! The singularity is treated using vcoul_lwl. 858 vcoul(1, iq_ibz) = zero 859 vcoul(2:,iq_ibz) = four_pi/vcoul(2:,iq_ibz)**2 860 else 861 vcoul(:,iq_ibz) = four_pi/vcoul(:,iq_ibz)**2 862 end if 863 end do 864 865 ! q-points for optical limit. 866 do iqlwl=1,Vcp%nqlwl 867 call cmod_qpg(Vcp%nqlwl,iqlwl,Vcp%qlwl,ng,Gsph%gvec,gprimd,vcoul_lwl(:,iqlwl)) 868 end do 869 vcoul_lwl = four_pi/vcoul_lwl**2 870 ! 871 ! === Integration of 1/q^2 singularity === 872 ! * We use the auxiliary function from PRB 75, 205126 (2007) [[cite:Carrier2007]] 873 q0_vol=(two_pi)**3/(Kmesh%nbz*ucvol) ; bz_geometry_factor=zero 874 875 ! $$ MG: this is to restore the previous implementation, it will facilitate the merge 876 ! Analytic integration of 4pi/q^2 over the volume element: 877 ! $4pi/V \int_V d^3q 1/q^2 =4pi bz_geometric_factor V^(-2/3)$ 878 ! i_sz=4*pi*bz_geometry_factor*q0_vol**(-two_thirds) where q0_vol= V_BZ/N_k 879 ! bz_geometry_factor: sphere=7.79, fcc=7.44, sc=6.188, bcc=6.946, wz=5.255 880 ! (see gwa.pdf, appendix A.4) 881 Vcp%i_sz=four_pi*7.44*q0_vol**(-two_thirds) 882 883 CASE ('ERF') 884 ! Modified long-range only Coulomb interaction thanks to the error function: 885 ! * Vc = erf(r/rcut)/r 886 ! * The singularity is treated using vcoul_lwl. 887 do iq_ibz=1,Qmesh%nibz 888 call cmod_qpg(Qmesh%nibz,iq_ibz,Qmesh%ibz,ng,Gsph%gvec,gprimd,vcoul(:,iq_ibz)) 889 890 ! * The Fourier transform of the error function reads 891 if (iq_ibz==1) then 892 vcoul(1, iq_ibz) = zero 893 vcoul(2:,iq_ibz) = four_pi/(vcoul(2:,iq_ibz)**2) * EXP( -0.25d0 * (Vcp%rcut*vcoul(2:,iq_ibz))**2 ) 894 else 895 vcoul(:,iq_ibz) = four_pi/(vcoul(:, iq_ibz)**2) * EXP( -0.25d0 * (Vcp%rcut*vcoul(: ,iq_ibz))**2 ) 896 end if 897 end do 898 899 ! q-points for optical limit. 900 do iqlwl=1,Vcp%nqlwl 901 call cmod_qpg(Vcp%nqlwl,iqlwl,Vcp%qlwl,ng,Gsph%gvec,gprimd,vcoul_lwl(:,iqlwl)) 902 end do 903 vcoul_lwl = four_pi/(vcoul_lwl**2) * EXP( -0.25d0 * (Vcp%rcut*vcoul_lwl)**2 ) 904 905 ! 906 ! === Treat 1/q^2 singularity === 907 ! * We use the auxiliary function from PRB 75, 205126 (2007) [[cite:Carrier2007]] 908 q0_vol=(two_pi)**3/(Kmesh%nbz*ucvol) ; bz_geometry_factor=zero 909 do iq_bz=1,Qmesh%nbz 910 qbz_cart(:)=Qmesh%bz(1,iq_bz)*b1(:)+Qmesh%bz(2,iq_bz)*b2(:)+Qmesh%bz(3,iq_bz)*b3(:) 911 qbz_norm=SQRT(SUM(qbz_cart(:)**2)) 912 if (qbz_norm>tolq0) bz_geometry_factor=bz_geometry_factor-faux(Qmesh%bz(:,iq_bz)) 913 end do 914 915 bz_geometry_factor = bz_geometry_factor + integratefaux()*Qmesh%nbz 916 917 write(msg,'(2a,2x,f12.4)')ch10,& 918 & ' integrate q->0 : numerical BZ geometry factor = ',bz_geometry_factor*q0_vol**(2./3.) 919 call wrtout(std_out,msg,'COLL') 920 Vcp%i_sz=four_pi*bz_geometry_factor ! Final result stored here 921 922 CASE ('ERFC') 923 ! * Use a modified short-range only Coulomb interaction thanks to the complementar error function: 924 ! $ V_c = [1-erf(r/r_{cut})]/r $ 925 ! * The Fourier transform of the error function reads 926 ! vcoul=four_pi/(vcoul**2) * ( 1.d0 - exp( -0.25d0 * (Vcp%rcut*vcoul)**2 ) ) 927 do iq_ibz=1,Qmesh%nibz 928 call cmod_qpg(Qmesh%nibz,iq_ibz,Qmesh%ibz,ng,Gsph%gvec,gprimd,vcoul(:,iq_ibz)) 929 930 if (iq_ibz==1) then 931 vcoul(1 ,iq_ibz) = zero 932 vcoul(2:,iq_ibz) = four_pi/(vcoul(2:,iq_ibz)**2) * ( one - EXP( -0.25d0 * (Vcp%rcut*vcoul(2:,iq_ibz))**2 ) ) 933 else 934 vcoul(:, iq_ibz) = four_pi/(vcoul(:, iq_ibz)**2) * ( one - EXP( -0.25d0 * (Vcp%rcut*vcoul(:, iq_ibz))**2 ) ) 935 end if 936 end do 937 ! 938 ! q-points for optical limit. 939 do iqlwl=1,Vcp%nqlwl 940 call cmod_qpg(Vcp%nqlwl,iqlwl,Vcp%qlwl,ng,Gsph%gvec,gprimd,vcoul_lwl(:,iqlwl)) 941 end do 942 vcoul_lwl = four_pi/(vcoul_lwl**2) * ( one - EXP( -0.25d0 * (Vcp%rcut*vcoul_lwl)**2 ) ) 943 ! 944 ! === Treat 1/q^2 singularity === 945 ! * There is NO singularity in this case. 946 Vcp%i_sz=pi*Vcp%rcut**2 ! Final result stored here 947 948 CASE DEFAULT 949 MSG_BUG(sjoin('Unknown cutoff mode:',Vcp%mode)) 950 END SELECT 951 ! 952 ! === Store final results in complex array === 953 ! * Rozzi"s cutoff can give real negative values 954 955 ABI_MALLOC(Vcp%vc_sqrt,(ng,Vcp%nqibz)) 956 ABI_MALLOC(Vcp%vc_sqrt_resid,(ng,Vcp%nqibz)) 957 Vcp%vc_sqrt=CMPLX(vcoul,zero) 958 Vcp%vc_sqrt=SQRT(Vcp%vc_sqrt) 959 Vcp%vc_sqrt_resid=Vcp%vc_sqrt 960 Vcp%i_sz_resid=Vcp%i_sz 961 call vcoul_plot(Vcp,Qmesh,Gsph,ng,vcoul,comm) 962 ABI_FREE(vcoul) 963 964 ABI_MALLOC(Vcp%vcqlwl_sqrt,(ng,Vcp%nqlwl)) 965 Vcp%vcqlwl_sqrt=CMPLX(vcoul_lwl,zero) 966 Vcp%vcqlwl_sqrt=SQRT(Vcp%vcqlwl_sqrt) 967 ABI_FREE(vcoul_lwl) 968 969 call vcoul_print(Vcp,unit=std_out) 970 971 DBG_EXIT("COLL") 972 973 contains !=============================================================== 974 975 real(dp) function integratefaux() 976 977 978 !This section has been created automatically by the script Abilint (TD). 979 !Do not modify the following lines by hand. 980 #undef ABI_FUNC 981 #define ABI_FUNC 'integratefaux' 982 !End of the abilint section 983 984 implicit none 985 986 !Arguments ------------------------------------ 987 !MG here there is a bug in abilint since integratefaux is chaged to integrate, dont know why! 988 !real(dp) :: integrate 989 990 !Local variables------------------------------- 991 integer,parameter :: nref=3 992 integer,parameter :: nq=50 993 integer :: ierr,iq,iqx1,iqy1,iqz1,iqx2,iqy2,iqz2,miniqy1,maxiqy1,nqhalf 994 real(dp) :: invnq,invnq3,qq,weightq,weightxy,weightxyz 995 real(dp) :: qq1(3),qq2(3),bb4sinpiqq_2(3,nq),sin2piqq(nq),bb4sinpiqq2_2(3,0:nq),sin2piqq2(3,0:nq) 996 997 ! nq is the number of sampling points along each axis for the numerical integration 998 ! nref is the area where the mesh is refined 999 1000 integratefaux=zero 1001 invnq=one/DBLE(nq) 1002 invnq3=invnq**3 1003 nqhalf=nq/2 1004 1005 !In order to speed up the calculation, precompute the sines 1006 do iq=1,nq 1007 qq=DBLE(iq)*invnq-half 1008 bb4sinpiqq_2(:,iq)=bb(:)*four*SIN(pi*qq)**2 ; sin2piqq(iq)=SIN(two_pi*qq) 1009 end do 1010 1011 do iqx1=1,nq 1012 if(modulo(iqx1,nprocs)/=rank) cycle 1013 qq1(1)=DBLE(iqx1)*invnq-half 1014 !Here take advantage of the q <=> -q symmetry : arrange the sampling of qx, qy space to avoid duplicating calculations. 1015 !Need weights to do this ... 1016 !do iqy1=1,nq 1017 miniqy1=nqhalf+1 ; maxiqy1=nq 1018 if(iqx1>=nqhalf)miniqy1=nqhalf 1019 if(iqx1>nqhalf .and. iqx1<nq)maxiqy1=nq-1 1020 do iqy1=miniqy1,maxiqy1 1021 qq1(2)=DBLE(iqy1)*invnq-half 1022 !By default, the factor of two is for the q <=> -q symmetry 1023 weightq=invnq3*two 1024 !But not all qx qy lines have a symmetric one ... 1025 if( (iqx1==nqhalf .or. iqx1==nq) .and. (iqy1==nqhalf .or. iqy1==nq))weightq=weightq*half 1026 1027 do iqz1=1,nq 1028 qq1(3)=DBLE(iqz1)*invnq-half 1029 ! 1030 ! Refine the mesh for the point close to the origin 1031 if( abs(iqx1-nqhalf)<=nref .and. abs(iqy1-nqhalf)<=nref .and. abs(iqz1-nqhalf)<=nref ) then 1032 !Note that the set of point is symmetric around the central point, while weights are taken into account 1033 do iq=0,nq 1034 qq2(:)=qq1(:)+ (DBLE(iq)*invnq-half)*invnq 1035 bb4sinpiqq2_2(:,iq)=bb(:)*four*SIN(pi*qq2(:))**2 ; sin2piqq2(:,iq)=SIN(two_pi*qq2(:)) 1036 enddo 1037 do iqx2=0,nq 1038 qq2(1)=qq1(1) + (DBLE(iqx2)*invnq-half ) *invnq 1039 do iqy2=0,nq 1040 qq2(2)=qq1(2) + (DBLE(iqy2)*invnq-half ) *invnq 1041 weightxy=invnq3*weightq 1042 if(iqx2==0 .or. iqx2==nq)weightxy=weightxy*half 1043 if(iqy2==0 .or. iqy2==nq)weightxy=weightxy*half 1044 do iqz2=0,nq 1045 qq2(3)=qq1(3) + (DBLE(iqz2)*invnq-half ) *invnq 1046 weightxyz=weightxy 1047 if(iqz2==0 .or. iqz2==nq)weightxyz=weightxy*half 1048 ! 1049 ! Treat the remaining divergence in the origin as if it would be a spherical 1050 ! integration of 1/q^2 1051 if( iqx1/=nqhalf .or. iqy1/=nqhalf .or. iqz1/=nqhalf .or. iqx2/=nqhalf .or. iqy2/=nqhalf .or. iqz2/=nqhalf ) then 1052 ! integratefaux=integratefaux+ faux(qq2) *invnq**6 1053 integratefaux=integratefaux+ faux_fast(qq2,bb4sinpiqq2_2(1,iqx2),bb4sinpiqq2_2(2,iqy2),bb4sinpiqq2_2(3,iqz2),& 1054 & sin2piqq2(1,iqx2),sin2piqq2(2,iqy2),sin2piqq2(3,iqz2)) * weightxyz 1055 else 1056 integratefaux=integratefaux + 7.7955* ( (two_pi)**3/ucvol*invnq3*invnq3 )**(-2./3.) *invnq3*invnq3 1057 end if 1058 end do 1059 end do 1060 end do 1061 else 1062 ! integratefaux=integratefaux+faux(qq1)*invnq**3 1063 integratefaux=integratefaux+ faux_fast(qq1,bb4sinpiqq_2(1,iqx1),bb4sinpiqq_2(2,iqy1),bb4sinpiqq_2(3,iqz1),& 1064 & sin2piqq(iqx1),sin2piqq(iqy1),sin2piqq(iqz1)) * weightq 1065 end if 1066 end do 1067 end do 1068 end do 1069 1070 call xmpi_sum(integratefaux,comm,ierr) 1071 1072 end function integratefaux 1073 1074 function faux(qq) 1075 1076 1077 !This section has been created automatically by the script Abilint (TD). 1078 !Do not modify the following lines by hand. 1079 #undef ABI_FUNC 1080 #define ABI_FUNC 'faux' 1081 !End of the abilint section 1082 1083 real(dp),intent(in) :: qq(3) 1084 real(dp) :: faux 1085 real(dp) :: bb4sinpiqq1_2,bb4sinpiqq2_2,bb4sinpiqq3_2,sin2piqq1,sin2piqq2,sin2piqq3 1086 1087 bb4sinpiqq1_2=b1b1*four*SIN(pi*qq(1))**2 1088 bb4sinpiqq2_2=b2b2*four*SIN(pi*qq(2))**2 1089 bb4sinpiqq3_2=b3b3*four*SIN(pi*qq(3))**2 1090 sin2piqq1=SIN(two_pi*qq(1)) 1091 sin2piqq2=SIN(two_pi*qq(2)) 1092 sin2piqq3=SIN(two_pi*qq(3)) 1093 1094 faux=faux_fast(qq,bb4sinpiqq1_2,bb4sinpiqq2_2,bb4sinpiqq3_2,sin2piqq1,sin2piqq2,sin2piqq3) 1095 1096 end function faux 1097 1098 function faux_fast(qq,bb4sinpiqq1_2,bb4sinpiqq2_2,bb4sinpiqq3_2,sin2piqq1,sin2piqq2,sin2piqq3) 1099 1100 1101 !This section has been created automatically by the script Abilint (TD). 1102 !Do not modify the following lines by hand. 1103 #undef ABI_FUNC 1104 #define ABI_FUNC 'faux_fast' 1105 !End of the abilint section 1106 1107 real(dp),intent(in) :: qq(3) 1108 real(dp) :: faux_fast 1109 real(dp) :: bb4sinpiqq1_2,bb4sinpiqq2_2,bb4sinpiqq3_2,sin2piqq1,sin2piqq2,sin2piqq3 1110 1111 faux_fast= bb4sinpiqq1_2 + bb4sinpiqq2_2 + bb4sinpiqq3_2 & 1112 & +two*( b1b2 * sin2piqq1*sin2piqq2 & 1113 & +b2b3 * sin2piqq2*sin2piqq3 & 1114 & +b3b1 * sin2piqq3*sin2piqq1 & 1115 & ) 1116 if(Vcp%rcut>tol6)then 1117 faux_fast=two_pi*two_pi/faux_fast * exp( -0.25d0*Vcp%rcut**2* sum( ( qq(1)*b1(:)+qq(2)*b2(:)+qq(3)*b3(:) )**2 ) ) 1118 else 1119 faux_fast=two_pi*two_pi/faux_fast 1120 endif 1121 1122 end function faux_fast 1123 1124 1125 function adapt_nmc(nmc_max,qpg2) result(nmc) 1126 1127 1128 !This section has been created automatically by the script Abilint (TD). 1129 !Do not modify the following lines by hand. 1130 #undef ABI_FUNC 1131 #define ABI_FUNC 'adapt_nmc' 1132 !End of the abilint section 1133 1134 real(dp),intent(in) :: qpg2 1135 integer,intent(in) :: nmc_max 1136 integer :: nmc 1137 1138 ! Empirical law to decrease the Monte Carlo sampling 1139 ! for large q+G, for which the accuracy is not an issue 1140 nmc = NINT( nmc_max / ( 1.0_dp + 1.0_dp * qpg2**6 ) ) 1141 nmc = MIN(nmc_max,nmc) 1142 nmc = MAX(1,nmc) 1143 1144 end function adapt_nmc 1145 1146 1147 end subroutine vcoul_init
m_vcoul/vcoul_plot [ Functions ]
[ Top ] [ m_vcoul ] [ Functions ]
NAME
vcoul_plot
FUNCTION
Plot vccut(q,G) as a function of |q+G|. Calculate also vc in real space.
INPUTS
OUTPUT
PARENTS
m_vcoul
CHILDREN
calck0,paw_jbessel,quadrature
SOURCE
1171 subroutine vcoul_plot(Vcp,Qmesh,Gsph,ng,vc,comm) 1172 1173 1174 !This section has been created automatically by the script Abilint (TD). 1175 !Do not modify the following lines by hand. 1176 #undef ABI_FUNC 1177 #define ABI_FUNC 'vcoul_plot' 1178 !End of the abilint section 1179 1180 implicit none 1181 1182 !Arguments ------------------------------------ 1183 !scalars 1184 integer,intent(in) :: ng,comm 1185 type(kmesh_t),intent(in) :: Qmesh 1186 type(vcoul_t),intent(in) :: Vcp 1187 type(gsphere_t),intent(in) :: Gsph 1188 !arrays 1189 real(dp),intent(in) :: vc(ng,Qmesh%nibz) 1190 1191 !Local variables------------------------------- 1192 !scalars 1193 integer,parameter :: master=0 1194 integer :: icount,idx_Sm1G,ierr,ig,igs,ii,iq_bz,iq_ibz,iqg,ir,isym,itim 1195 integer :: my_start,my_stop,nqbz,nqibz,nr,ntasks,rank,unt 1196 real(dp) :: arg,fact,l1,l2,l3,lmax,step,tmp,vcft,vc_bare 1197 character(len=500) :: msg 1198 character(len=fnlen) :: filnam 1199 !arrays 1200 integer,allocatable :: insort(:) 1201 real(dp) :: b1(3),b2(3),b3(3),gmet(3,3),gprimd(3,3),qbz(3),qpgc(3) 1202 real(dp),allocatable :: qpg_mod(:),rr(:,:,:),vcr(:,:),vcr_cut(:,:) 1203 1204 !************************************************************************ 1205 1206 if (TRIM(Vcp%mode)/='CYLINDER') RETURN 1207 1208 rank = xmpi_comm_rank(comm) 1209 1210 nqibz=Qmesh%nibz; nqbz=Qmesh%nbz 1211 gmet=Gsph%gmet; gprimd=Gsph%gprimd 1212 1213 b1(:)=two_pi*gprimd(:,1) 1214 b2(:)=two_pi*gprimd(:,2) 1215 b3(:)=two_pi*gprimd(:,3) 1216 1217 ! === Compare in Fourier space the true Coulomb with the cutted one === 1218 if (rank==master) then 1219 ABI_MALLOC(insort,(nqibz*ng)) 1220 ABI_MALLOC(qpg_mod,(nqibz*ng)) 1221 iqg=1 1222 do iq_ibz=1,nqibz 1223 do ig=1,ng 1224 qpg_mod(iqg)=normv(Qmesh%ibz(:,iq_ibz)+Gsph%gvec(:,ig),gmet,'g') 1225 insort(iqg)=iqg; iqg=iqg+1 1226 end do 1227 end do 1228 call sort_dp(nqibz*ng,qpg_mod,insort,tol14) 1229 1230 filnam='_VCoulFT_' 1231 call isfile(filnam,'new') 1232 if (open_file(filnam,msg,newunit=unt,status='new',form='formatted') /=0) then 1233 MSG_ERROR(msg) 1234 end if 1235 write(unt,'(a,i3,a,i6,a)')& 1236 & '# |q+G| q-point (Tot no.',nqibz,') Gvec (',ng,') vc_bare(q,G) vc_cutoff(q,G) ' 1237 1238 do iqg=1,nqibz*ng 1239 iq_ibz=(insort(iqg)-1)/ng +1 1240 ig=(insort(iqg))-(iq_ibz-1)*ng 1241 vc_bare=zero 1242 if (qpg_mod(iqg)>tol16) vc_bare=four_pi/qpg_mod(iqg)**2 1243 write(unt,'(f12.6,2x,3f8.4,2x,3i6,2x,2es14.6)')& 1244 & qpg_mod(iqg),Qmesh%ibz(:,iq_ibz),Gsph%gvec(:,ig),vc_bare,vc(ig,iq_ibz) 1245 end do 1246 1247 close(unt) 1248 ABI_FREE(insort) 1249 ABI_FREE(qpg_mod) 1250 end if ! rank==master 1251 1252 ! === Fourier transform back to real space just to check cutoff implementation === 1253 ntasks=nqbz*ng 1254 call xmpi_split_work(ntasks,comm,my_start,my_stop,msg,ierr) 1255 if (ierr/=0) then 1256 MSG_WARNING(msg) 1257 end if 1258 1259 l1=SQRT(SUM(Vcp%rprimd(:,1)**2)) 1260 l2=SQRT(SUM(Vcp%rprimd(:,2)**2)) 1261 l3=SQRT(SUM(Vcp%rprimd(:,3)**2)) 1262 1263 nr=50 1264 lmax=MAX(l1,l2,l3) ; step=lmax/(nr-1) 1265 fact=one/(Vcp%ucvol*nqbz) 1266 ! 1267 ! numb coding 1268 ABI_MALLOC(rr,(3,nr,3)) 1269 rr=zero 1270 do ii=1,3 1271 do ir=1,nr 1272 rr(ii,ir,ii)=(ir-1)*step 1273 end do 1274 end do 1275 1276 ABI_MALLOC(vcr,(nr,3)) 1277 ABI_MALLOC(vcr_cut,(nr,3)) 1278 vcr=zero; vcr_cut=zero 1279 1280 do iq_bz=1,nqbz 1281 call get_BZ_item(Qmesh,iq_bz,qbz,iq_ibz,isym,itim) 1282 if (ABS(qbz(1))<0.01) qbz(1)=zero 1283 if (ABS(qbz(2))<0.01) qbz(2)=zero 1284 if (ABS(qbz(3))<0.01) qbz(3)=zero 1285 igs=1 ; if (ALL(qbz(:)==zero)) igs=2 1286 do ig=igs,ng 1287 icount=ig+(iq_bz-1)*ng 1288 if (icount<my_start.or.icount>my_stop) CYCLE 1289 idx_Sm1G=Gsph%rottbm1(ig,itim,isym) ! IS{^-1}G 1290 vcft=vc(idx_Sm1G,iq_ibz) 1291 qpgc(:)=qbz(:)+Gsph%gvec(:,ig) ; qpgc(:)=b1(:)*qpgc(1)+b2(:)*qpgc(2)+b3(:)*qpgc(3) 1292 tmp=SQRT(DOT_PRODUCT(qpgc,qpgc)) ; tmp=tmp**2 1293 do ii=1,3 1294 do ir=1,nr 1295 arg=DOT_PRODUCT(rr(:,ir,ii),qpgc) 1296 vcr_cut(ir,ii)=vcr_cut(ir,ii) + vcft*COS(arg) 1297 vcr (ir,ii)=vcr (ir,ii) + four_pi/tmp*COS(arg) 1298 end do 1299 end do 1300 end do !ig 1301 end do !iq_ibz 1302 1303 call xmpi_sum_master(vcr_cut,master,comm,ierr) 1304 call xmpi_sum_master(vcr ,master,comm,ierr) 1305 1306 if (rank==master) then 1307 filnam='_VCoulR_' 1308 call isfile(filnam,'new') 1309 if (open_file(filnam,msg,newunit=unt,status='new',form='formatted') /= 0) then 1310 MSG_ERROR(msg) 1311 end if 1312 write(unt,'(a)')'# length vc_bare(r) vc_cut(r) ' 1313 do ir=1,nr 1314 write(unt,'(7es18.6)')(ir-1)*step,(fact*vcr(ir,ii),fact*vcr_cut(ir,ii),ii=1,3) 1315 end do 1316 close(unt) 1317 end if 1318 1319 ABI_FREE(rr) 1320 ABI_FREE(vcr) 1321 ABI_FREE(vcr_cut) 1322 1323 end subroutine vcoul_plot
m_vcoul/vcoul_print [ Functions ]
[ Top ] [ m_vcoul ] [ Functions ]
NAME
vcoul_print
FUNCTION
Print the content of a Coulomb datatype.
INPUTS
Vcp<vcoul_t>=The datatype whose content has to be printed. [unit]=The unit number for output [prtvol]=Verbosity level [mode_paral]=Either "COLL" or "PERS".
PARENTS
m_vcoul
CHILDREN
calck0,paw_jbessel,quadrature
SOURCE
1349 subroutine vcoul_print(Vcp,unit,prtvol,mode_paral) 1350 1351 1352 !This section has been created automatically by the script Abilint (TD). 1353 !Do not modify the following lines by hand. 1354 #undef ABI_FUNC 1355 #define ABI_FUNC 'vcoul_print' 1356 !End of the abilint section 1357 1358 implicit none 1359 1360 !Arguments ------------------------------------ 1361 !scalars 1362 integer,intent(in),optional :: prtvol,unit 1363 character(len=4),intent(in),optional :: mode_paral 1364 type(vcoul_t),intent(in) :: Vcp 1365 1366 !Local variables------------------------------- 1367 !scalars 1368 integer :: ii,my_unt,my_prtvol,iqlwl 1369 character(len=4) :: my_mode 1370 character(len=500) :: msg 1371 1372 ! ************************************************************************* 1373 1374 my_unt =std_out; if (PRESENT(unit )) my_unt =unit 1375 my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral 1376 my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol 1377 1378 SELECT CASE (Vcp%mode) 1379 1380 CASE ('MINIBZ') 1381 write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',TRIM(Vcp%mode) 1382 call wrtout(my_unt,msg,my_mode) 1383 1384 CASE ('MINIBZ-ERF') 1385 write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',TRIM(Vcp%mode) 1386 call wrtout(my_unt,msg,my_mode) 1387 write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,& 1388 & ' === Error function cutoff === ',ch10,ch10,& 1389 & ' Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10 1390 call wrtout(my_unt,msg,my_mode) 1391 1392 CASE ('MINIBZ-ERFC') 1393 write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',TRIM(Vcp%mode) 1394 call wrtout(my_unt,msg,my_mode) 1395 write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,& 1396 & ' === Complement Error function cutoff === ',ch10,ch10,& 1397 & ' Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10 1398 call wrtout(my_unt,msg,my_mode) 1399 1400 CASE ('SPHERE') 1401 write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,& 1402 & ' === Spherical cutoff === ',ch10,ch10,& 1403 & ' Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10,& 1404 & ' Volume of the sphere .. ',four_pi/three*Vcp%rcut**3,' [Bohr^3] ' 1405 !FB: This has no meaning here! & ' Sphere centered at .... ',Vcp%boxcenter,' (r.l.u) ',ch10 1406 !MG It might be useful if the system is not centered on the origin because in this case the 1407 ! matrix elements of the Coulomb have to be multiplied by a phase depending on boxcenter. 1408 ! I still have to decide if it is useful to code this possibility and which variable use to 1409 ! define the center (boxcenter is used in the tddft part). 1410 call wrtout(my_unt,msg,my_mode) 1411 1412 CASE ('CYLINDER') 1413 ii=imin_loc(ABS(Vcp%pdir-1)) 1414 write(msg,'(5a,f10.4,3a,i2,2a,3f10.2,a)')ch10,& 1415 & ' === Cylindrical cutoff === ',ch10,ch10,& 1416 & ' Cutoff radius ............... ',Vcp%rcut,' [Bohr] ',ch10,& 1417 & ' Axis parallel to direction... ',ii,ch10,& 1418 & ' Passing through point ....... ',Vcp%boxcenter,' (r.l.u) ' 1419 call wrtout(my_unt,msg,my_mode) 1420 1421 write(msg,'(2a)')' Infinite length ....... ',ch10 1422 if (Vcp%hcyl/=zero) write(msg,'(a,f8.5,2a)')' Finite length of ....... ',Vcp%hcyl,' [Bohr] ',ch10 1423 call wrtout(my_unt,msg,my_mode) 1424 1425 CASE ('SURFACE') 1426 write(msg,'(5a,f10.4,3a,3f10.2,2a)')ch10,& 1427 & ' === Surface cutoff === ',ch10,ch10,& 1428 & ' Cutoff radius .................... ',Vcp%rcut,' [Bohr] ',ch10,& 1429 & ' Central plane passing through .... ',Vcp%boxcenter,' (r.l.u) ',ch10 1430 call wrtout(my_unt,msg,my_mode) 1431 !write(msg,'(a)')' Infinite length .......' 1432 !if (Vcp%hcyl/=zero) write(msg,'(a,f8.5,a)')' Finite length of .......',Vcp%hcyl,' [Bohr] ' 1433 !call wrtout(my_unt,msg,my_mode) 1434 1435 CASE ('AUXILIARY_FUNCTION') 1436 write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',TRIM(Vcp%mode) 1437 call wrtout(my_unt,msg,my_mode) 1438 1439 CASE ('AUX_GB') 1440 write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',TRIM(Vcp%mode) 1441 call wrtout(my_unt,msg,my_mode) 1442 1443 CASE ('CRYSTAL') 1444 write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',TRIM(Vcp%mode) 1445 call wrtout(my_unt,msg,my_mode) 1446 1447 CASE ('ERF') 1448 write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,& 1449 & ' === Error function cutoff === ',ch10,ch10,& 1450 & ' Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10 1451 call wrtout(my_unt,msg,my_mode) 1452 1453 CASE ('ERFC') 1454 write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,& 1455 & ' === Complement Error function cutoff === ',ch10,ch10,& 1456 & ' Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10 1457 call wrtout(my_unt,msg,my_mode) 1458 1459 CASE DEFAULT 1460 MSG_BUG(sjoin('Unknown cutoff mode: ', Vcp%mode)) 1461 END SELECT 1462 1463 if (Vcp%nqlwl>0) then 1464 write(msg,'(a,i3)')" q-points for optical limit: ",Vcp%nqlwl 1465 call wrtout(my_unt,msg,my_mode) 1466 do iqlwl=1,Vcp%nqlwl 1467 write(msg,'(1x,i5,a,2x,3f12.6)') iqlwl,')',Vcp%qlwl(:,iqlwl) 1468 call wrtout(my_unt,msg,my_mode) 1469 end do 1470 end if 1471 1472 !TODO add additional information 1473 1474 end subroutine vcoul_print
m_vcoul/vcoul_t [ Types ]
NAME
vcoul_t
FUNCTION
This data type contains the square root of the Fourier components of the Coulomb interaction calculated taking into account a possible cutoff. Info on the particular geometry used for the cutoff as well as quantities required to deal with the Coulomb divergence.
SOURCE
70 type,public :: vcoul_t 71 72 ! TODO: Remove it 73 integer :: nfft 74 ! Number of points in FFT grid 75 76 integer :: ng 77 ! Number of G-vectors 78 79 integer :: nqibz 80 ! Number of irreducible q-points 81 82 integer :: nqlwl 83 ! Number of small q-points around Gamma 84 85 real(dp) :: alpha(3) 86 ! Lenght of the finite surface 87 88 real(dp) :: rcut 89 ! Cutoff radius 90 91 real(dp) :: i_sz 92 ! Value of the integration of the Coulomb singularity 4\pi/V_BZ \int_BZ d^3q 1/q^2 93 94 real(dp) :: i_sz_resid 95 ! Residual difference between the i_sz in the sigma self-energy for exchange, 96 ! and the i_sz already present in the generalized Kohn-Sham eigenenergies 97 ! Initialized to the same value as i_sz 98 99 real(dp) :: hcyl 100 ! Length of the finite cylinder along the periodic dimension 101 102 real(dp) :: ucvol 103 ! Volume of the unit cell 104 105 character(len=50) :: mode 106 ! String defining the cutoff mode, possible values are: sphere,cylinder,surface,crystal 107 108 integer :: pdir(3) 109 ! 1 if the system is periodic along this direction 110 111 ! TODO: Remove it 112 integer :: ngfft(18) 113 ! Information on the FFT grid 114 115 real(dp) :: boxcenter(3) 116 ! 1 if the point in inside the cutoff region 0 otherwise 117 ! Reduced coordinates of the center of the box (input variable) 118 119 real(dp) :: vcutgeo(3) 120 ! For each reduced direction gives the length of the finite system 121 ! 0 if the system is infinite along that particular direction 122 ! negative value to indicate that a finite size has to be used 123 124 real(dp) :: rprimd(3,3) 125 ! Lattice vectors in real space. 126 127 real(dp),allocatable :: qibz(:,:) 128 ! qibz(3,nqibz) 129 ! q-points in the IBZ. 130 131 real(dp),allocatable :: qlwl(:,:) 132 ! qibz(3,nqlwl) 133 ! q-points for the treatment of the Coulomb singularity. 134 135 complex(gwpc),allocatable :: vc_sqrt(:,:) 136 ! vc_sqrt(ng,nqibz) 137 ! Square root of the Coulomb interaction in reciprocal space. 138 ! A cut might be applied. 139 140 complex(gwpc),allocatable :: vcqlwl_sqrt(:,:) 141 ! vcqs_sqrt(ng,nqlwl) 142 ! Square root of the Coulomb term calculated for small q-points 143 144 complex(gwpc),allocatable :: vc_sqrt_resid(:,:) 145 ! vc_sqrt_resid(ng,nqibz) 146 ! Square root of the residual difference between the Coulomb interaction in the sigma self-energy for exchange, 147 ! and the Coulomb interaction already present in the generalized Kohn-Sham eigenenergies (when they come from an hybrid) 148 ! Given in reciprocal space. At the call to vcoul_init, it is simply initialized at the value of vc_sqrt(:,:), 149 ! and only later modified. 150 ! A cut might be applied. 151 152 end type vcoul_t 153 154 public :: vcoul_init ! Main creation method. 155 public :: vcoul_plot ! Plot vc in real and reciprocal space. 156 public :: vcoul_print ! Report info on the object. 157 public :: vcoul_free ! Destruction method. 158 public :: cmod_qpg ! FT of the long ranged Coulomb interaction.