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
 
  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 :: end => xc_vdw_end
  end type xc_vdw_t
 
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_d3 3
    !% The DFT-D3 scheme of S. Grimme, J. Antony, S. Ehrlich, and
    !% S. Krieg, J. Chem. Phys. 132, 154104 (2010).
    !%End
    call parse_variable(namespace, 'VDWCorrection', option__vdwcorrection__none, this%vdw_correction)
 
    if (this%vdw_correction /= option__vdwcorrection__none) then
      call messages_experimental('VDWCorrection')
 
      select case (this%vdw_correction)
      case (option__vdwcorrection__vdw_ts)
 
        !%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)
 
      case (option__vdwcorrection__vdw_d3)
        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 = 4, tz = .false.)
 
      case default
        assert(.false.)
      end select
 
    else
      this%self_consistent = .false.
    end if
 
    pop_sub(xc_vdw_init)
  end subroutine xc_vdw_init
 
  ! ---------------------------------------------------------
  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),                      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))
 
    select case (this%vdw_correction)
 
    case (option__vdwcorrection__vdw_ts)
      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
 
    case (option__vdwcorrection__vdw_d3)
 
      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)
 
    case default
      assert(.false.)
    end select
 
    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: