xc_fbe.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2023 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 xc_fbe_oct_m
  use comm_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use electron_space_oct_m
  use exchange_operator_oct_m
  use global_oct_m
  use grid_oct_m
  use kpoints_oct_m
  use lalg_basic_oct_m
  use messages_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use nl_operator_oct_m
  use parser_oct_m
  use profiling_oct_m
  use poisson_oct_m
  use ring_pattern_oct_m
  use solvers_oct_m
  use space_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use varinfo_oct_m
  use xc_functional_oct_m
 
  implicit none
 
  private
  public ::               &
    x_fbe_calc,           &
    lda_c_fbe,            &
    fbe_c_lda_sl
 
  type(grid_t), pointer :: gr_aux  => null()
  real(real64), pointer :: rho_aux(:) => null()
  real(real64), allocatable  :: diag_lapl(:)
 
 
contains
 
  ! -------------------------------------------------------------------------------------
  subroutine x_fbe_calc (id, namespace, psolver, gr, st, space, ex, vxc)
    integer,                     intent(in)    :: id
    type(namespace_t),           intent(in)    :: namespace
    type(poisson_t),             intent(in)    :: psolver
    type(grid_t),                intent(in)    :: gr
    type(states_elec_t),         intent(inout) :: st
    type(space_t),               intent(in)    :: space
    real(real64),                intent(inout) :: ex
    real(real64), optional,      intent(inout) :: vxc(:,:)
 
    real(real64), allocatable :: fxc(:,:,:), internal_vxc(:,:)
 
    push_sub(x_fbe_calc)
 
    select case(id)
    case(xc_oep_x_fbe)
      if (states_are_real(st)) then
        call dx_fbe_calc(namespace, psolver, gr, gr%der, st, ex, vxc=vxc)
      else
        call zx_fbe_calc(namespace, psolver, gr, gr%der, st, ex, vxc=vxc)
      end if
    case(xc_oep_x_fbe_sl)
      safe_allocate(fxc(1:gr%np_part, 1:gr%box%dim, 1:st%d%spin_channels))
      safe_allocate(internal_vxc(1:gr%np, 1:st%d%spin_channels))
      internal_vxc = m_zero
      ! We first compute the force density
      if (states_are_real(st)) then
        call dx_fbe_calc(namespace, psolver, gr, gr%der, st, ex, vxc=internal_vxc, fxc=fxc)
      else
        call zx_fbe_calc(namespace, psolver, gr, gr%der, st, ex, vxc=internal_vxc, fxc=fxc)
      end if
 
      ! We solve the Sturm-Liouville equation
      if (present(vxc)) then
        call solve_sturm_liouville(namespace, gr, st, space, fxc, internal_vxc)
      end if
 
      ! Get the energy from the virial relation
      ex = get_virial_energy(gr, st%d%spin_channels, fxc)
 
      ! Adds the calculated potential
      call lalg_axpy(gr%np, st%d%spin_channels, m_one, internal_vxc, vxc)
 
      safe_deallocate_a(fxc)
      safe_deallocate_a(internal_vxc)
    case default
      assert(.false.)
    end select
 
    pop_sub(x_fbe_calc)
  end subroutine x_fbe_calc
 
  ! -------------------------------------------------------------------------------------
  subroutine solve_sturm_liouville(namespace, gr, st, space, fxc, vxc)
    type(namespace_t),            intent(in)     :: namespace
    type(grid_t), target,         intent(in)     :: gr
    type(states_elec_t), target,  intent(in)     :: st
    type(space_t),                intent(in)     :: space
    real(real64),  contiguous,    intent(inout)  :: fxc(:,:,:)
    real(real64),                 intent(inout)  :: vxc(:,:)
 
    real(real64), allocatable :: rhs(:)
    integer :: iter, ispin
    real(real64) :: res
    real(real64), parameter :: threshold = 1e-7_real64
    character(len=32) :: name
 
    type(nl_operator_t) :: op(1)
 
    push_sub(solve_sturm_liouville)
 
    assert(ubound(fxc, dim=1) >= gr%np_part)
 
    gr_aux => gr
    call mesh_init_mesh_aux(gr)
 
    ! the smoothing is performed uing the same stencil as the Laplacian
    name = 'FBE preconditioner'
    call derivatives_get_lapl(gr%der, namespace, op, space, name, 1)
    safe_allocate(diag_lapl(1:op(1)%np))
    call dnl_operator_operate_diag(op(1), diag_lapl)
    call nl_operator_end(op(1))
 
    safe_allocate(rhs(1:gr%np))
 
    do ispin = 1, st%d%spin_channels
      call dderivatives_div(gr%der, fxc(:, :, ispin), rhs)
      rho_aux => st%rho(:, ispin)
 
      iter = 500
      call dqmr_sym_gen_dotu(gr%np, vxc(:, ispin), rhs, &
        sl_operator, dmf_dotu_aux, dmf_nrm2_aux, preconditioner, &
        iter, residue = res, threshold = threshold, showprogress = .false.)
 
      write(message(1), '(a, i6, a)') "Info: Sturm-Liouville solver converged in  ", iter, " iterations."
      write(message(2), '(a, es14.6)') "Info: The residue is ", res
      call messages_info(2, namespace=namespace)
    end do
 
    safe_deallocate_a(rhs)
 
    safe_deallocate_a(diag_lapl)
 
    nullify(rho_aux)
    nullify(gr_aux)
 
    pop_sub(solve_sturm_liouville)
  contains
    !----------------------------------------------------------------
    subroutine sl_operator(x, hx)
      real(real64), contiguous, intent(in)  :: x(:)
      real(real64), contiguous, intent(out) :: hx(:)
 
      integer :: ip, idir
      real(real64), allocatable :: vxc(:)
      real(real64), allocatable :: grad_vxc(:,:)
 
      safe_allocate(vxc(1:gr_aux%np_part))
      safe_allocate(grad_vxc(1:gr_aux%np_part, 1:gr_aux%box%dim))
      call lalg_copy(gr_aux%np, x, vxc)
 
      call dderivatives_grad(gr_aux%der, vxc, grad_vxc)
 
      !$omp parallel
      do idir = 1, gr_aux%box%dim
        !$omp do
        do ip = 1, gr_aux%np
          grad_vxc(ip, idir) = grad_vxc(ip, idir)*rho_aux(ip)
        end do
        !$omp end do nowait
      end do
      !$omp end parallel
 
      call dderivatives_div(gr_aux%der, grad_vxc, hx)
 
      safe_deallocate_a(vxc)
      safe_deallocate_a(grad_vxc)
    end subroutine sl_operator
 
    !----------------------------------------------------------------
    subroutine preconditioner(x, hx)
      real(real64), contiguous, intent(in)  :: x(:)
      real(real64), contiguous, intent(out) :: hx(:)
 
      integer :: ip
 
      !$omp parallel do
      do ip = 1, gr_aux%np
        hx(ip) = x(ip) / (max(rho_aux(ip), 1d-12) * diag_lapl(ip))
      end do
 
    end subroutine preconditioner
 
  end subroutine solve_sturm_liouville
 
  ! -------------------------------------------------------------------------------------
  real(real64) function get_virial_energy(gr, nspin, fxc) result(exc)
    type(grid_t),  intent(in) :: gr
    integer,       intent(in) :: nspin
    real(real64),  intent(in) :: fxc(:,:,:)
 
    integer :: isp, idir, ip
    real(real64), allocatable :: rfxc(:)
    real(real64) :: xx(gr%box%dim), rr
 
    push_sub(get_virial_energy)
 
    exc = m_zero
    do isp = 1, nspin
      safe_allocate(rfxc(1:gr%np))
      do ip = 1, gr%np
        rfxc(ip) = m_zero
        call mesh_r(gr, ip, rr, coords=xx)
        do idir = 1, gr%box%dim
          rfxc(ip) = rfxc(ip) + fxc(ip, idir, isp) * xx(idir)
        end do
      end do
      exc = exc + dmf_integrate(gr, rfxc)
      safe_deallocate_a(rfxc)
    end do
 
    pop_sub(get_virial_energy)
  end function get_virial_energy
 
 
  ! -------------------------------------------------------------------------------------
  subroutine lda_c_fbe (st, n_blocks, l_dens, l_dedd, l_zk)
    type(states_elec_t),         intent(in)    :: st
    integer,                     intent(in)    :: n_blocks
    real(real64),                intent(in)    :: l_dens(:,:)
    real(real64),                intent(inout) :: l_dedd(:,:)
    real(real64), optional,      intent(inout) :: l_zk(:)
 
    integer :: ip, ispin
    real(real64) :: rho, beta, beta2, e_c
    real(real64) :: q
 
    push_sub(lda_c_fbe)
 
    ! Set q such that we get the leading order of the r_s->0 limit for the HEG
    q = ((5.0_real64*sqrt(m_pi)**5)/(m_three*(m_one-log(m_two))))**(m_third)
    if (present(l_zk)) l_zk = m_zero
 
    do ip = 1, n_blocks
      rho = sum(l_dens(1:st%d%spin_channels, ip))
      if (rho < 1e-20_real64) then
        l_dedd(1:st%d%spin_channels, ip) = m_zero
        cycle
      end if
      rho = max(rho, 1e-12_real64)
      beta = q*rho**m_third
      beta2 = beta**2
 
      ! Potential
      ! First part of the potential
      l_dedd(1:st%d%spin_channels, ip) = (m_pi/(q**3))*((sqrt(m_pi)*beta/(m_one+sqrt(m_pi)*beta))**2 -m_one) * beta
      ! Second part of the potential
      l_dedd(1:st%d%spin_channels, ip) = l_dedd(1:st%d%spin_channels, ip) &
        - (5.0_real64*sqrt(m_pi))/(m_three*q**3)*(log(m_one+sqrt(m_pi)*beta) &
        -m_half/(m_one+sqrt(m_pi)*beta)**2 + m_two/(m_one+sqrt(m_pi)*beta)) + (5.0_real64*sqrt(m_pi))/(m_two*q**3)
 
      if (st%d%nspin == 1 .and. present(l_zk)) then
        ! Energy density
        ! First part of the energy density
        e_c = (9.0_real64*q**3)/m_two/beta &
          - m_two*q**3*sqrt(m_pi) &
          - 12.0_real64/beta2*(q**3/sqrt(m_pi)) &
          + m_three/(m_pi*rho)*(m_one/(m_one+sqrt(m_pi)*beta) - m_one &
          + 5.0_real64*log(m_one+sqrt(m_pi)*beta))
 
        ! Second part of the energy density
        e_c = e_c - 5.0_real64/6.0_real64*( &
          7.0_real64*q**3/beta &
          + m_three/(m_pi*rho*(m_one+sqrt(m_pi)*beta)) &
          - 17.0_real64*q**3/sqrt(m_pi)/beta2 &
          - 11.0_real64*q**3*sqrt(m_pi)/(m_three) &
          + (20.0_real64/(m_pi*rho) + m_two*sqrt(m_pi)*q**3)*log(m_one+sqrt(m_pi)*beta) &
          - m_three/(m_pi*rho))
        e_c = e_c/(q**6)
        l_zk(ip) = e_c
      else if(st%d%nspin == 2) then
        ! Here we have no energy density, so leave the potential unchanged
        ! This is the approximate potential that we implement here
        do ispin = 1, st%d%spin_channels
          l_dedd(ispin, ip) = l_dedd(ispin, ip) * m_two * l_dens(-ispin+3, ip) / rho
        end do
      end if
    end do
 
    pop_sub(lda_c_fbe)
  end subroutine lda_c_fbe
 
  ! -------------------------------------------------------------------------------------
  subroutine fbe_c_lda_sl (namespace, psolver, gr, st, space, ec, vxc)
    type(namespace_t),           intent(in)    :: namespace
    type(poisson_t),             intent(in)    :: psolver
    type(grid_t),                intent(in)    :: gr
    type(states_elec_t),         intent(inout) :: st
    type(space_t),               intent(in)    :: space
    real(real64),                intent(inout) :: ec
    real(real64), optional,             intent(inout) :: vxc(:,:)
 
    integer :: idir, ip, ispin
    real(real64), allocatable :: fxc(:,:,:), internal_vxc(:,:), grad_rho(:,:,:), tmp1(:,:), tmp2(:,:)
    real(real64) :: q, beta, rho, l_gdens
 
    push_sub(fbe_c_lda_sl)
 
    safe_allocate(internal_vxc(1:gr%np, 1:st%d%spin_channels))
 
    ! Needed to get the initial guess for the iterative solution of the Sturm-Liouville equation
    safe_allocate(tmp1(1:st%d%spin_channels, 1:gr%np))
    safe_allocate(tmp2(1:st%d%spin_channels, 1:gr%np))
    tmp1 = transpose(st%rho(1:gr%np, 1:st%d%spin_channels))
    call lda_c_fbe(st, gr%np, tmp1, tmp2)
    internal_vxc = transpose(tmp2)
    safe_deallocate_a(tmp1)
    safe_deallocate_a(tmp2)
 
    ! Set q such that we get the leading order of the r_s->0 limit for the HEG
    q = ((5.0_real64*sqrt(m_pi)**5)/(m_three*(m_one-log(m_two))))**(m_third)
 
    safe_allocate(fxc(1:gr%np_part, 1:gr%box%dim, 1:st%d%spin_channels))
    safe_allocate(grad_rho(1:gr%np, 1:gr%box%dim, 1:st%d%spin_channels))
    do ispin = 1, st%d%spin_channels
      call dderivatives_grad(gr%der, st%rho(:, ispin), grad_rho(:, :, ispin))
    end do
 
    do ispin = 1, st%d%spin_channels
      do idir = 1, gr%box%dim
        do ip = 1, gr%np
          rho = sum(st%rho(ip, 1:st%d%spin_channels))
          if (st%rho(ip, ispin) < 1e-20_real64) then
            fxc(ip, idir, ispin) = m_zero
            cycle
          end if
          rho = max(rho, 1e-12_real64)
          beta = rho**m_third * q
 
          l_gdens = sum(grad_rho(ip, idir, 1:st%d%spin_channels))
 
          if (st%d%spin_channels == 1) then
            fxc(ip, idir, ispin) = l_gdens * &
              ( m_pi * beta**2/((m_one + sqrt(m_pi)*beta)**2) - m_one &
              + m_third * m_pi * beta**2 / ((m_one + sqrt(m_pi)*beta)**3) )
          else
            fxc(ip, idir, ispin) = m_two * (grad_rho(ip, idir, 3-ispin) * &
              (m_pi * beta**2/((m_one + sqrt(m_pi)*beta)**2) - m_one ) &
              + l_gdens * (m_third * m_pi * beta**2 / ((m_one + sqrt(m_pi)*beta)**3) ) &
              * st%rho(ip, 3-ispin) / rho)
          end if
 
          fxc(ip, idir, ispin) = fxc(ip, idir, ispin) * m_pi/(m_three*beta**2)  * st%rho(ip, ispin)
        end do
      end do
    end do
 
    ! We solve the Sturm-Liouville equation
    if (present(vxc)) then
      call solve_sturm_liouville(namespace, gr, st, space, fxc, internal_vxc)
    end if
 
    ! Get the energy from the virial relation
    ec = get_virial_energy(gr, st%d%spin_channels, fxc)
 
    ! Adds the calculated potential
    call lalg_axpy(gr%np, st%d%spin_channels, m_one, internal_vxc, vxc)
 
    safe_deallocate_a(fxc)
 
    pop_sub(fbe_c_lda_sl)
  end subroutine fbe_c_lda_sl
 
 
#include "undef.F90"
#include "real.F90"
#include "xc_fbe_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "xc_fbe_inc.F90"
 
end module xc_fbe_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: