xc_vdw.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2023 . 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_vdw_oct_m
  use atom_oct_m
  use debug_oct_m
  use global_oct_m
  use grid_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use lattice_vectors_oct_m
  use libvdwxc_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use species_oct_m
  use species_factory_oct_m
  use states_elec_oct_m
  use vdw_ts_oct_m
  use xc_oct_m
  use xc_interaction_oct_m
  use xc_functional_oct_m
 
  ! from the dftd3 library
  use dftd3_api
 
  implicit none
 
  private
  public ::             &
    xc_vdw_t, &
    d3_lib_options
 
  type xc_vdw_t
    private
    integer,                  public :: vdw_correction = -1
    logical                          :: self_consistent = .false.
    type(vdw_ts_t),           public :: vdw_ts
    type(dftd3_calc)                 :: vdw_d3
    real(real64), public,       allocatable :: forces(:, :)
 
    real(real64), public                    :: stress(3, 3)
 
  contains
    procedure :: init => xc_vdw_init
    procedure :: calc => xc_vdw_calc
    procedure :: d3_version => xc_vdw_octopus_input_to_version
    procedure :: end => xc_vdw_end
  end type xc_vdw_t
 
  ! No reason for this to be int64, but options are appended with this
  ! precision when parsed.
  integer(int64), parameter :: d3_lib_options(5) = [&
    option__vdwcorrection__vdw_d2, &
    option__vdwcorrection__vdw_d3_zero_damping, &
    option__vdwcorrection__vdw_d3_bj, &
    option__vdwcorrection__vdw_d3m_zero_damping, &
    option__vdwcorrection__vdw_d3m_bj &
    ]
 
contains
 
  subroutine xc_vdw_init(this, namespace, space, gr, xc, ions, x_id, c_id)
    class(xc_vdw_t),         intent(inout) :: this
    type(namespace_t),       intent(in)    :: namespace
    class(space_t),          intent(in)    :: space
    type(grid_t),            intent(in)    :: gr
    type(xc_t),              intent(inout) :: xc
    type(ions_t),            intent(in)    :: ions
    integer,                 intent(in)    :: x_id, c_id
 
    integer :: iatom
    type(dftd3_input) :: d3_input
    character(len=20) :: d3func_def, d3func
 
    push_sub(xc_vdw_init)
 
    if (in_family(xc%family, [xc_family_libvdwxc])) then
      call libvdwxc_set_geometry(xc%functional(func_c,1)%libvdwxc, namespace, space, gr)
    end if
 
    !%Variable VDWCorrection
    !%Type integer
    !%Default no
    !%Section Hamiltonian::XC
    !%Description
    !% (Experimental) This variable selects which van der Waals
    !% correction to apply to the correlation functional.
    !%Option none 0
    !% No correction is applied.
    !%Option vdw_ts 1
    !% The scheme of Tkatchenko and Scheffler, Phys. Rev. Lett. 102
    !% 073005 (2009).
    !%Option vdw_d2 2
    !% DFT-D2 corrections, equivalent to specifying version = 2 in the DFTD3 library.
    !%Option vdw_d3_zero_damping 3
    !% DFT-D3 corrections, damping the dispersion contribution to zero for short ranges.
    !% This is equivalent to specifying version = 3 in the DFTD3 library.
    !% For more details, see [J. Chem. Phys. 132, 154104 (2010)](https://doi.org/10.1063/1.3382344).
    !%Option vdw_d3_bj 4
    !% DFT-D3 with Becke-Johnson finite-damping, variant 2. In Octopus 15.0 and older, this was only option with the DFT-D3 library, and was
    !% specified with "vdw_d3".
    !% This is equivalent to specifying version = 4 in the DFTD3 library.
    !% For more details, see [Effect of the damping function in dispersion corrected density functional theory](https://doi.org/10.1002/jcc.21759).
    !%Option vdw_d3m_zero_damping 5
    !% DFT-D3 corrections, damping the dispersion contribution to zero for short ranges. This uses the modified damping function and parameters
    !% of Sherrill and coworkers. For more details see [J. Phys. Chem. Lett. 2016, 7, 12, 2197–2203](https://doi.org/10.1021/acs.jpclett.6b00780)
    !% This is equivalent to specifying version = 5 in the DFTD3 library.
    !%Option vdw_d3m_bj 6
    !% DFT-D3 with Becke-Johnson finite-damping, and the modified parameters of Sherrill and coworkers.
    !% For more details see [J. Phys. Chem. Lett. 2016, 7, 12, 2197–2203](https://doi.org/10.1021/acs.jpclett.6b00780)
    !% This is equivalent to specifying version = 6 in the DFTD3 library.
    !%End
    call parse_variable(namespace, 'VDWCorrection', option__vdwcorrection__none, this%vdw_correction)
 
    if (this%vdw_correction == option__vdwcorrection__none) then
      this%self_consistent = .false.
      pop_sub(xc_vdw_init)
      return
    endif
 
    call messages_experimental('VDWCorrection')
 
    if( this%vdw_correction == option__vdwcorrection__vdw_ts) then
      !%Variable VDWSelfConsistent
      !%Type logical
      !%Default yes
      !%Section Hamiltonian::XC
      !%Description
      !% This variable controls whether the VDW correction is applied
      !% self-consistently, the default, or just as a correction to
      !% the total energy. This option only works with vdw_ts.
      !%End
      call parse_variable(namespace, 'VDWSelfConsistent', .true., this%self_consistent)
 
      call vdw_ts_init(this%vdw_ts, namespace, ions)
 
      ! DFTD3 library
    else
      this%self_consistent = .false.
 
      if (space%dim /= 3) then
        call messages_write('xc_vdw_d3 can only be used in 3-dimensional systems')
        call messages_fatal(namespace=namespace)
      end if
 
      do iatom = 1, ions%natoms
        if (.not. ions%atom(iatom)%species%represents_real_atom()) then
          call messages_write('xc_vdw_d3 is not implemented when non-atomic species are present')
          call messages_fatal(namespace=namespace)
        end if
      end do
 
      d3func_def = ''
 
      ! The list of valid values can be found in 'external_libs/dftd3/core.f90'.
      ! For the moment I include the most common ones.
      if (x_id == option__xcfunctional__gga_x_b88 .and. c_id*libxc_c_index == option__xcfunctional__gga_c_lyp) then
        d3func_def = 'b-lyp'
      end if
      if (x_id == option__xcfunctional__gga_x_pbe .and. c_id*libxc_c_index == option__xcfunctional__gga_c_pbe) then
        d3func_def = 'pbe'
      end if
      if (x_id == option__xcfunctional__gga_x_pbe_sol .and. c_id*libxc_c_index == option__xcfunctional__gga_c_pbe_sol) then
        d3func_def = 'pbesol'
      end if
      if (c_id*libxc_c_index == option__xcfunctional__hyb_gga_xc_b3lyp) then
        d3func_def = 'b3-lyp'
      end if
      if (c_id*libxc_c_index == option__xcfunctional__hyb_gga_xc_pbeh) then
        d3func_def = 'pbe0'
      end if
 
      !%Variable VDWD3Functional
      !%Type string
      !%Section Hamiltonian::XC
      !%Description
      !% (Experimental) You can use this variable to override the
      !% parametrization used by the DFT-D3 van deer Waals
      !% correction. Normally you need not set this variable, as the
      !% proper value will be selected by Octopus (if available).
      !%
      !% This variable takes a string value, the valid values can
      !% be found in the source file 'external_libs/dftd3/core.f90'.
      !% For example you can use:
      !%
      !%  VDWD3Functional = 'pbe'
      !%
      !%End
      if (parse_is_defined(namespace, 'VDWD3Functional')) call messages_experimental('VDWD3Functional')
      call parse_variable(namespace, 'VDWD3Functional', d3func_def, d3func)
 
      if (d3func == '') then
        call messages_write('Cannot find  a matching parametrization  of DFT-D3 for the current')
        call messages_new_line()
        call messages_write('XCFunctional.  Please select a different XCFunctional, or select a')
        call messages_new_line()
        call messages_write('functional for DFT-D3 using the <tt>VDWD3Functional</tt> variable.')
        call messages_fatal(namespace=namespace)
      end if
 
      if (space%periodic_dim /= 0 .and. space%periodic_dim /= 3) then
        call messages_write('For partially periodic systems,  the xc_vdw_d3 interaction is assumed')
        call messages_new_line()
        call messages_write('to be periodic in three dimensions.')
        call messages_warning(namespace=namespace)
      end if
 
      call dftd3_init(this%vdw_d3, d3_input, trim(conf%share)//'/dftd3/pars.dat')
      call dftd3_set_functional(this%vdw_d3, func = d3func, version = this%d3_version(), tz = .false.)
    endif
 
    pop_sub(xc_vdw_init)
  end subroutine xc_vdw_init
 
 
  function xc_vdw_octopus_input_to_version(this) result(d3_version)
    class(xc_vdw_t), intent(in) :: this
 
    integer                     :: d3_version
 
    push_sub(xc_vdw_octopus_input_to_version)
 
    select case (this%vdw_correction)
    case (option__vdwcorrection__vdw_d2)
      d3_version = 2
    case (option__vdwcorrection__vdw_d3_zero_damping)
      d3_version = 3
    case (option__vdwcorrection__vdw_d3_bj)
      d3_version = 4
    case (option__vdwcorrection__vdw_d3m_zero_damping)
      d3_version = 5
    case (option__vdwcorrection__vdw_d3m_bj)
      d3_version = 6
    case default
      assert(.false.)
    end select
 
    pop_sub(xc_vdw_octopus_input_to_version)
 
  end function xc_vdw_octopus_input_to_version
 
 
  ! ---------------------------------------------------------
  subroutine xc_vdw_end(this)
    class(xc_vdw_t),   intent(inout) :: this
 
    push_sub(xc_vdw_end)
 
    select case (this%vdw_correction)
    case (option__vdwcorrection__vdw_ts)
      call vdw_ts_end(this%vdw_ts)
    end select
 
    pop_sub(xc_vdw_end)
  end subroutine xc_vdw_end
 
  ! ---------------------------------------------------------
  subroutine xc_vdw_calc(this, namespace, space, latt, atom, natoms, pos, gr, st, energy, vxc)
    class(xc_vdw_t),                   intent(inout) :: this
    type(namespace_t),                 intent(in)    :: namespace
    type(space_t),                     intent(in)    :: space
    type(lattice_vectors_t),           intent(in)    :: latt
    type(atom_t),                      intent(in)    :: atom(:)
    integer,                           intent(in)    :: natoms
    real(real64),                      intent(in)    :: pos(1:space%dim,1:natoms)
    type(grid_t),                      intent(in)    :: gr
    type(states_elec_t),               intent(inout) :: st
    real(real64),                      intent(out)   :: energy
    real(real64), contiguous,          intent(inout) :: vxc(:,:)
 
    integer :: ispin, iatom, idir
    real(real64), allocatable :: vxc_vdw(:)
    integer, allocatable :: atnum(:)
    real(real64) :: rlattice(3, 3)
    ! For non-periodic directions, we make the length of the box very large
    real(real64), parameter :: aperiodic_scaling = 1000.0_real64
 
    energy = m_zero
    if (this%vdw_correction == option__vdwcorrection__none) return
 
    push_sub(xc_vdw_calc)
 
    assert(space%dim == 3)
 
    safe_allocate(vxc_vdw(1:gr%np))
    safe_allocate(this%forces(1:space%dim, 1:natoms))
 
    if(this%vdw_correction == option__vdwcorrection__vdw_ts) then
      vxc_vdw = m_zero
      call vdw_ts_calculate(this%vdw_ts, namespace, space, latt, atom, natoms, pos, &
        gr, st%d%nspin, st%rho, energy, vxc_vdw, this%forces)
      this%stress = m_zero
 
    elseif(any(this%vdw_correction == d3_lib_options)) then
 
      safe_allocate(atnum(1:natoms))
 
      atnum = [(nint(atom(iatom)%species%get_z()), iatom=1, natoms)]
 
      if (space%is_periodic()) then
        ! DFTD3 treats interactions as 3D periodic. In the case of
        ! mixed-periodicity, we set the length of the periodic lattice along
        ! the non-periodic dimensions to be much larger than the system
        ! size, such that for all practical purposes the system is treated
        ! as isolated.
        rlattice = latt%rlattice(1:3, 1:3)
        if (space%periodic_dim < 3) then
          message(1) = "Info: Using DFT-D3 van der Walls corrections for a system with mixed-periodicity,"
          message(2) = "      but DFTD3 treats system as fully periodic. Octopus will set the DFT-D3"
          message(3) = "      lattice vectors along non-periodic dimensions to a suitably large value."
          call messages_info(3, namespace=namespace)
 
          do idir = space%periodic_dim + 1, 3
            rlattice(idir,idir) = (maxval(abs(pos(idir,:))) + m_one)*aperiodic_scaling
          end do
        end if
        call dftd3_pbc_dispersion(this%vdw_d3, pos, atnum, rlattice, energy, this%forces, &
          this%stress)
      else
        call dftd3_dispersion(this%vdw_d3, pos, atnum, energy, this%forces)
        this%stress = m_zero
      end if
 
      safe_deallocate_a(atnum)
 
    else
      assert(.false.)
    endif
 
    if (this%self_consistent) then
      do ispin = 1, st%d%nspin
        call lalg_axpy(gr%np, m_one, vxc_vdw, vxc(:, ispin))
      end do
    end if
 
    safe_deallocate_a(vxc_vdw)
 
    pop_sub(xc_vdw_calc)
  end subroutine xc_vdw_calc
 
end module xc_vdw_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: