singularity.F90

Go to the documentation of this file.

!! Copyright (C) 2019-2022 N. Tancogne-Dejean
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module singularity_oct_m
  use comm_oct_m
  use debug_oct_m
  use distributed_oct_m
  use electron_space_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lattice_vectors_oct_m
  use messages_oct_m
  use mpi_oct_m
  use parser_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
  public ::                          &
    singularity_t,             &
    singularity_init,          &
    singularity_end
 
  integer, public, parameter ::        &
    SINGULARITY_NONE                  = 0, &
    singularity_general               = 1, &
    singularity_gygi                  = 2, &
    singularity_sphere                = 3
 
 
  type singularity_t
    !For the treatment of the singularity in solids
    integer :: coulomb_singularity = 0
    real(real64), allocatable :: Fk(:)
    real(real64) :: FF
    real(real64) :: energy
  end type singularity_t
 
contains
 
  subroutine singularity_init(this, namespace, space, st, kpoints)
    type(singularity_t),       intent(inout) :: this
    type(namespace_t),         intent(in)    :: namespace
    class(space_t),            intent(in)    :: space
    type(states_elec_t),       intent(in)    :: st
    type(kpoints_t),           intent(in)    :: kpoints
 
    integer :: default
 
    push_sub(singularity_init)
 
    this%energy = m_zero
 
    if (.not. allocated(this%Fk)) then
      safe_allocate(this%Fk(1:st%nik))
      this%Fk(1:st%nik) = m_zero
      this%FF = m_zero
 
      if (space%is_periodic()) then
 
        !%Variable HFSingularity
        !%Type integer
        !%Default general
        !%Section Hamiltonian::XC
        !%Description
        !% (Experimental) This variable selects the method used for the treatment of the
        !% singularity of the Coulomb potential in Hatree-Fock and hybrid-functional DFT calculations.
        !% This shoulbe be only applied for periodic systems and is only
        !% used for FFT kernels of the Poisson solvers.
        !%Option none 0
        !% The singularity is replaced by zero.
        !%Option general 1
        !% The general treatment of the singularity, as described in Carrier et al, PRB 75 205126 (2007).
        !% This is the default option
        !%Option fcc 2
        !% The treatment of the singulariy as described in Gygi and Baldereschi, PRB 34, 4405 (1986).
        !% This is formally valid for cubic systems only.
        !%Option spherical_bz 3
        !% The divergence in q=0 is treated analytically assuming a spherical Brillouin zone
        !%End
 
        default = singularity_general
        if (space%dim == 2 .or. space%dim > 3) default = singularity_none
 
        call parse_variable(namespace, 'HFSingularity', default, this%coulomb_singularity)
        call messages_print_var_option('HFSingularity', this%coulomb_singularity, namespace=namespace)
 
        if(this%coulomb_singularity /= singularity_none) then
          if(this%coulomb_singularity /= singularity_general .and. space%dim == 1) then
            call messages_not_implemented("HFSingularity /= general for 1D")
          end if
 
          if(space%dim == 2) then
            call messages_not_implemented("HFSingularity /= none for 2D")
          end if
        end if
 
        if (this%coulomb_singularity /= singularity_none) then
          call singularity_correction(this, namespace, space, st, kpoints)
        end if
      end if
    end if
 
    pop_sub(singularity_init)
  end subroutine singularity_init
 
  subroutine singularity_end(this)
    type(singularity_t), intent(inout) :: this
 
    push_sub(singularity_end)
 
    safe_deallocate_a(this%Fk)
    this%coulomb_singularity = -1
 
    pop_sub(singularity_end)
  end subroutine singularity_end
 
  !This routine implements the general tratment of the singularity for periodic solids,
  !as described in Carrier et al. PRB 75, 205126 (2007)
  subroutine singularity_correction(this, namespace, space, st, kpoints)
    type(singularity_t),       intent(inout) :: this
    type(namespace_t),         intent(in)    :: namespace
    class(space_t),            intent(in)    :: space
    type(states_elec_t),       intent(in)    :: st
    type(kpoints_t),           intent(in)    :: kpoints
 
    integer :: ik, ik2, ikpoint, nk, nsteps
    integer :: ikx, iky, ikz, istep, kpt_start, kpt_end
    real(real64) :: length
    real(real64) :: kpoint(space%dim), qpoint(space%dim)
    real(real64) :: kvol_element
    type(distributed_t) :: dist_kpt
    integer :: default_nk, default_step
    real(real64), parameter :: singul_cnst = 7.7955541794415_real64 !The constant is 4*pi*(3/(4*pi))^1/3
    real(real64), parameter :: log2_minus_gamma = 0.11593151565841244881_real64 !The constant is log(2) - gamma, where gamma is the Euler-Mascheroni constant
    push_sub(singularity_correction)
 
    call profiling_in("COULOMB_SINGULARITY")
 
    !At the moment this is only implemented in 3D and in 1D.
    assert(space%dim == 3 .or. this%coulomb_singularity == singularity_general)
 
    kpt_start = st%d%kpt%start
    kpt_end = st%d%kpt%end
 
    if (.not. st%d%kpt%parallel) then
      call distributed_init(dist_kpt, st%nik, mpi_world%comm, "singularity")
      kpt_start = dist_kpt%start
      kpt_end = dist_kpt%end
    end if
 
    do ik = kpt_start, kpt_end
      ikpoint = st%d%get_kpoint_index(ik)
      kpoint = kpoints%get_point(ikpoint, absolute_coordinates = .false.)
 
      this%Fk(ik) = m_zero
 
      do ik2 = 1, kpoints%full%npoints
        qpoint = kpoint - kpoints%full%red_point(:, ik2)
 
        !We remove potential umklapp
        qpoint = qpoint - anint(qpoint + 1e-5_real64)
 
        if (all(abs(qpoint) < 1e-6_real64)) cycle
 
        this%Fk(ik) = this%Fk(ik) + aux_funct(qpoint) * kpoints%full%weight(ik2)
      end do
      select case(space%dim)
      case(1)
        this%Fk(ik) = this%Fk(ik)/kpoints%latt%rcell_volume
      case(3)
        this%Fk(ik) = this%Fk(ik)*m_four*m_pi/kpoints%latt%rcell_volume
      end select
    end do
 
    if (st%d%kpt%parallel) then
      call comm_allreduce(st%d%kpt%mpi_grp, this%Fk)
    else
      call comm_allreduce(dist_kpt%mpi_grp, this%Fk)
      call distributed_end(dist_kpt)
    end if
 
    if (this%coulomb_singularity == singularity_general) then
      !%Variable HFSingularityNk
      !%Type integer
      !%Default 60 in 3D, 1200 in 1D
      !%Section Hamiltonian::XC
      !%Description
      !% Number of k-point used (total number of k-points) is (2*Nk+1)^3) in the numerical integration
      !% of the auxiliary function f(q). See PRB 75, 205126 (2007) for more details.
      !% Only for HFSingularity=general.
      !% Also used in 1D.
      !%End
      default_nk = 60
      if(space%dim == 1) default_nk = 1200
      call parse_variable(namespace, 'HFSingularityNk', default_nk, nk)
      if (abs(nk/m_three-nint(nk/m_three)) > m_epsilon) then
        message(1) = 'HFSingularity_Nk must be a multiple of 3.'
        call messages_fatal(1, namespace=namespace)
      end if
 
      !%Variable HFSingularityNsteps
      !%Type integer
      !%Default 7 in 3D, 15 in 1D
      !%Section Hamiltonian::XC
      !%Description
      !% Number of grid refinement steps in the numerical integration of the auxiliary function f(q).
      !% See PRB 75, 205126 (2007) for more details. Only for HFSingularity=general.
      !% Also used in 1D.
      !%End
      default_step = 7
      if(space%dim == 1) default_step = 15
      call parse_variable(namespace, 'HFSingularityNsteps', default_step, nsteps)
 
      select case(space%dim)
      case(1)
        this%FF = m_zero
        length = m_one
        kvol_element = (m_one/(m_two*nk+m_one))*((m_two*m_pi))/kpoints%latt%rcell_volume
        qpoint = m_zero
        do istep = 1, nsteps
 
          do ikx = 0, nk
            qpoint(1) = ikx/(m_two*nk)*length
 
            if(abs(ikx)<=nk/3) cycle
            this%FF = this%FF + aux_funct(qpoint)*kvol_element
          end do
          if(istep<nsteps) then
            length = length/m_three
            kvol_element = kvol_element/m_three
          end if
        end do
 
        !We have a factor two because we used the fact that f(q)=f(-q)
        !We multiply by 1/((2*pi)^1)
        this%FF = this%FF*m_two/((m_two*m_pi))
 
        !We add the remaining part
        !The constant is log(2) - gamma, where gamma is the Euler-Mascheroni constant
        this%FF = this%FF + m_two * (m_pi)/kpoints%latt%rcell_volume * length &
          * (m_one-log(m_pi / kpoints%latt%rcell_volume * length) + log2_minus_gamma)
      case(2)
        call messages_not_implemented("General Coulomb singularity in the 2D case")
 
      case(3)
        this%FF = m_zero
        length = m_one
        kvol_element = (m_one/(m_two*nk+m_one))**3*((m_two*m_pi)**3)/kpoints%latt%rcell_volume
        do istep = 1, nsteps
 
          do ikx = 0, nk
            qpoint(1) = ikx/(m_two*nk)*length
 
            do iky = -nk, nk
              qpoint(2) = iky/(m_two*nk)*length
 
              do ikz = -nk, nk
                qpoint(3) = ikz/(m_two*nk)*length
 
                if (abs(ikx) <= nk/3 .and. abs(iky) <= nk/3 .and. abs(ikz) <= nk/3) cycle
 
                this%FF = this%FF + aux_funct(qpoint)*kvol_element
              end do
            end do
          end do
          if (istep < nsteps) then
            length = length / m_three
            kvol_element = kvol_element / 27.0_real64
          end if
        end do
 
        !We have a factor two because we used the fact that f(q)=f(-q)
        !We multiply by 4*pi/((2*pi)^3)
        this%FF = this%FF*8.0_real64*m_pi/((m_two*m_pi)**3)
        !The remaining part is treated as a spherical BZ
        this%FF = this%FF + singul_cnst*(kpoints%latt%rcell_volume)**(m_twothird)/m_pi/kpoints%latt%rcell_volume*length
 
      end select
 
    else if (this%coulomb_singularity == singularity_gygi) then
      !See Eq. (7) of PRB 34, 4405 (1986)
      !Here we use the fact that the fcc volume is a^3/4
      this%FF = 4.423758_real64*(kpoints%latt%rcell_volume*m_four)**(m_twothird)/m_pi/kpoints%latt%rcell_volume
 
    else
      !The constant is 4*pi*(3/(4*pi))^1/3
      !We multiply by 4*pi/(2*pi^3)
      this%FF = singul_cnst*(kpoints%latt%rcell_volume)**(m_twothird)/m_pi/kpoints%latt%rcell_volume
    end if
 
 
    this%energy = m_zero
    do ik = st%d%kpt%start, st%d%kpt%end
      this%energy = this%energy + this%Fk(ik)*st%kweights(ik)
    end do
 
    if (st%d%kpt%parallel) then
      call comm_allreduce(st%d%kpt%mpi_grp, this%energy)
    end if
 
    ! In the case, the number of k-point is doubled, so the energy is twice larger
    if (st%d%ispin == spin_polarized) this%energy = m_half * this%energy
 
    this%energy = (this%energy-this%FF)*st%qtot
 
    write(message(1), '(a,f12.6,a,a,a)') 'Debug: Singularity energy ', &
      units_from_atomic(units_out%energy, this%energy), &
      ' [',trim(units_abbrev(units_out%energy)),']'
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    call profiling_out("COULOMB_SINGULARITY")
    pop_sub(singularity_correction)
 
  contains
 
    real(real64) function aux_funct(qq) result(ff)
      real(real64),            intent(in) :: qq(3)
 
      real(real64) :: half_a, qq_abs(space%dim)
 
      ! no PUSH/POP as called too often
 
      if (this%coulomb_singularity == singularity_general) then
        select case(space%dim)
        case(1)
          !The constant is -log(2) + gamma, where gamma is the Euler-Mascheroni constant
          ff =  -log(abs(sin(qq(1)*m_pi))*kpoints%latt%klattice(1,1)/(m_two*m_pi)) + 0.11593151565841244881_real64
        case(3)
          !See Eq. (16) of PRB 75, 205126 (2007)
          ff = (m_two*m_pi)**2/(m_two*(                                                             &
            (m_two*sin(qq(1)*m_pi)*sin(qq(1)*m_pi)*dot_product(kpoints%latt%klattice(:,1), kpoints%latt%klattice(:,1)) &
            + sin(qq(1)*m_two*m_pi)*sin(qq(2)*m_two*m_pi)*dot_product(kpoints%latt%klattice(:,1), kpoints%latt%klattice(:,2))) &
            + (m_two*sin(qq(2)*m_pi)*sin(qq(2)*m_pi)*dot_product(kpoints%latt%klattice(:,2), kpoints%latt%klattice(:,2)) &
            + sin(qq(2)*m_two*m_pi)*sin(qq(3)*m_two*m_pi)*dot_product(kpoints%latt%klattice(:,2), kpoints%latt%klattice(:,3))) &
            + (m_two*sin(qq(3)*m_pi)*sin(qq(3)*m_pi)*dot_product(kpoints%latt%klattice(:,3), kpoints%latt%klattice(:,3)) &
            + sin(qq(3)*m_two*m_pi)*sin(qq(1)*m_two*m_pi)*dot_product(kpoints%latt%klattice(:,3), kpoints%latt%klattice(:,1)))))
        end select
      else
        half_a = m_half*(kpoints%latt%rcell_volume*m_four)**(m_third)
        call kpoints_to_absolute(kpoints%latt, qq, qq_abs)
        !See Eq. (6) of PRB 34, 4405 (1986)
        ff = (half_a)**2/(m_three-cos(qq_abs(1)*half_a)*cos(qq_abs(2)*half_a) &
          -cos(qq_abs(1)*half_a)*cos(qq_abs(3)*half_a)         &
          -cos(qq_abs(3)*half_a)*cos(qq_abs(2)*half_a))
      end if
 
    end function aux_funct
 
  end subroutine singularity_correction
 
end module singularity_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: