gengclq.f90

Go to the documentation of this file.

! Copyright (C) 2002-2006 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

!BOP
! !ROUTINE: gengclq
! !INTERFACE:
subroutine gengclq
! !USES:
use modmain
use modtest
! !DESCRIPTION:
!   The Fock matrix elements
!   $$ V_{ij{\bf k}}\equiv\sum_{l{\bf k'}}\int
!    \frac{\Psi^{\dag}_{i{\bf k}}({\bf r})\cdot\Psi_{l{\bf k}'}({\bf r})
!    \Psi^{\dag}_{l{\bf k}'}({\bf r}')\cdot\Psi_{j{\bf k}}({\bf r}')}
!    {|{\bf r}-{\bf r'}|}\,d^3r\,d^3r' $$
!   contain a divergent term in the sum over ${\bf k}'$ which behaves as
!   $1/q^2$, where ${\bf q}\equiv{\bf k}-{\bf k}'$ is in the first Brillouin
!   zone. The resulting convergence with respect to the number of discrete
!   $q$-points, $N_q$, is very slow. This routine computes the regularised
!   Coulomb Green's function
!   \begin{align}
!    g({\bf q}_i)=\frac{4\pi}{V}\int_{V_i}\frac{1}{q^2}\,d^3q,
!   \end{align}
!   where the integral is over the small parallelepiped with volume
!   $V=\Omega_{\rm BZ}/N_q$ and centered on the discrete point ${\bf q}_i$.
!   This dramatically increases the rate of convergence of methods which involve
!   a summation over the $1/q^2$ part of the Coulomb interaction. The above
!   integral is evaluated numerically on increasingly finer grids and then
!   extrapolated to the continuum.
!
! !REVISION HISTORY:
!   Created August 2004 (JKD,SS)
!   Changed from genwiq2, July 2017 (JKD)
!EOP
!BOC
implicit none
! local variables
integer, parameter :: np=5,ns0=10,nss=20
integer ns,iq,i1,i2,i3,ip
real(8) d(3),sm,t1,t2
real(8) v1(3),v2(3),v3(3)
real(8) xa(np),ya(np)
! external functions
real(8), external :: polynm
! allocate global gclq array
if (allocated(gclq)) deallocate(gclq)
allocate(gclq(nqpt))
! begin loop over q-points, note that the vectors vqc are assumed to be in the
! first Brillouin zone
do iq=1,nqpt
! loop over different subdivisions
  ns=ns0
  do ip=1,np
! subdivision vectors in lattice coordinates
    d(1:3)=1.d0/dble(ngridq(1:3)*2*ns)
! compute the integral of 1/q²
    sm=0.d0
    do i1=-ns,ns-1
      t1=dble(i1)*d(1)
      v1(1:3)=vqc(1:3,iq)+t1*bvec(1:3,1)
      do i2=-ns,ns-1
        t1=dble(i2)*d(2)
        v2(1:3)=v1(1:3)+t1*bvec(1:3,2)
        do i3=-ns,ns-1
          t1=dble(i3)*d(3)
          v3(1:3)=v2(1:3)+t1*bvec(1:3,3)
          t2=v3(1)**2+v3(2)**2+v3(3)**2
          if (t2 > 1.d-14) sm=sm+1.d0/t2
        end do
      end do
    end do
    t1=1.d0/dble(2*ns)
    xa(ip)=t1
    ya(ip)=fourpi*sm*t1**3
! increment number of subdivisions
    ns=ns+nss
  end do
! extrapolate the volume element to zero with a polynomial
  gclq(iq)=polynm(0,np,xa,ya,0.d0)
end do
! zero the Green's function at q = 0 if required
if (t0gclq0) gclq(1)=0.d0
! write gclq to test file
call writetest(700,"regularised Coulomb Green''s function (gclq)",nv=nqpt, &
 tol=1.d-8,rva=gclq)
end subroutine
!EOC