Testing D3 damping parameters

Testing D3 damping parameters#

In this tutorial we are implementing a command line tool to compute the dispersion energy for a reaction and test multiple damping functions for D3. For this we will use the Fortran API via the dftd3 module.

Using the Fortran API#

To make the dftd3 available in our project we will use the Fortran package manager (fpm) to manage the necessary dependencies. We create a minimal package manifest with the following content:

fpm.toml#

name = "param-scanner"
version = "1.0.0"

[dependencies]
s-dftd3.git = "https://github.com/dftd3/simple-dftd3"

Our library will be called param-scanner because it will be scanning through the available damping parameters for D3. We also need to declare a version, here we go with 1.0.0 but feel free to choose any version which feels appropriate. Finally, we declare our dependencies via the URL of the git repository.

Computing dispersion for reactions#

For our parameter scanner to give meaningful results, we would like to target relative energies or reaction energies. We start by defining a type to hold the reaction definition, combining an array of structure_type values and an array of stochiometric coefficients. The structure_type value is not defined by the dftd3 module but the mctc_io module, which we already include as a transient dependency.

Note

For more details on mctc_io checkout its module documentation.

We setup our reaction_type as part of our main module d3_param_scan:

src/scan.f90#

module d3_param_scan
   use mctc_env, only : wp
   use mctc_io, only : structure_type
   implicit none
   private

   public :: reaction_type

   type :: reaction_type
      type(structure_type), allocatable :: mol(:)
      real(wp), allocatable :: stochiometry(:)
   end type reaction_type

end module d3_param_scan

Note

Here we do not define the precision ourselves but use the one defined by mctc_env.

Now can define a function to evaluate the dispersion energy for a reaction. For this we need to use the get_dispersion subroutine provided in the dftd3 module. We start with checking the signature of this function:

interface
   !> Calculate scalar dispersion energy.
   subroutine get_dispersion(mol, disp, param, cutoff, energy, gradient, sigma)
      use mctc_env, only : wp
      use mctc_io, only : structure_type
      use dftd3, only : d3_model, damping_param, realspace_cutoff
      !> Molecular structure data
      class(structure_type), intent(in) :: mol
      !> Dispersion model
      class(d3_model), intent(in) :: disp
      !> Damping parameters
      class(damping_param), intent(in) :: param
      !> Realspace cutoffs
      type(realspace_cutoff), intent(in) :: cutoff
      !> Dispersion energy
      real(wp), intent(out) :: energy
      !> Dispersion gradient
      real(wp), intent(out), contiguous, optional :: gradient(:, :)
      !> Dispersion virial
      real(wp), intent(out), contiguous, optional :: sigma(:, :)
   end subroutine get_dispersion
end interface

The main inputs we need to construct are the d3_model, the damping_param and the realspace_cutoff objects. We check the documentation and find that the d3_model can be constructed from a structure_type. Similarly, for the realspace_cutoff we find that the derived type is rather simple and can be easily constructed. For the damping_param we find that there are several implementations available. Our goal is to scan through the damping functions, therefore we will have a closer look to damping_param soon. First, we use get_dispersion to define our own dispersion energy evaluator for our reaction values.

src/scan.f90#

subroutine get_dispersion_for_reaction(reaction, param, energy)
   use dftd3, only : get_dispersion, damping_param, realspace_cutoff, &
      & d3_model, new_d3_model
   type(reaction_type), intent(in) :: reaction
   class(damping_param), intent(in) :: param
   real(wp), intent(inout) :: energy

   integer :: component
   real(wp) :: energy_component
   type(d3_model) :: disp

   do component = 1, size(reaction%mol)
      call new_d3_model(disp, reaction%mol(component))
      call get_dispersion(reaction%mol(component), disp, param, realspace_cutoff(), &
         & energy_component)
      energy = energy + energy_component * reaction%stochiometry(component)
   end do

end subroutine get_dispersion_for_reaction

We keep the interface for our get_dispersion_for_reaction function simple, we focus on the two inputs, the reaction data and the damping parameters, and obtain a single output, the dispersion energy. For the energy dummy argument we choose intent(inout) in contrast to the get_dispersion interface where it is intent(out).

Creating a parameter scanner#

Next we will look into actually iterating through all the available damping schemes available for D3. Before we start we define the interface for our procedure:

interface
   !> Scan all available damping parameters for a given method by computing
   !> and comparing the dispersion energy of a reaction.
   subroutine scan_param_for_reaction(error, reaction, method, dft_energy)
      import :: reaction_type
      use mctc_env, only : wp, error_type
      !> Error handling
      type(error_type), allocatable :: error
      !> Reaction data for relative energy computation
      type(reaction_type), intent(in) :: reaction
      !> Method name to query for damping parameters
      character(*), intent(in) :: method
      !> Optional DFT energy of the reaction
      real(wp), intent(in), optional :: dft_energy
   end subroutine scan_param_for_reaction
end interface

For our interface we will be using the error_type provided by the mctc_env module. In a larger code base we would probably have our own error handling strategy, for this tutorial however we will reuse the same as used in dftd3 to keep the code simple. Also, we will not return any output back to the calling routine here but rather directly print them out.

With this design decision made we can now have a closer look into the damping schemes and parameters again. We will start with the zero damping scheme D3 was originally published with, overall there are five different damping schemes currently supported in dftd3. For this we need to create an instance of the zero_damping_param, which is an implementation of damping_param. Based on the just defined interface we implement our first block to obtain zero damping parameters for the input method:

src/scan.f90#

subroutine scan_param_for_reaction(error, reaction, method, dft_energy)
   use mctc_env, only : error_type, fatal_error
   use dftd3, only : d3_param
   type(error_type), allocatable :: error
   type(reaction_type), intent(in) :: reaction
   character(*), intent(in) :: method
   real(wp), intent(in), optional :: dft_energy

   logical :: has_param(5)
   real(wp) :: disp_energies(2, 5)
   type(d3_param) :: inp
   integer, parameter :: zero_damping_idx = 1, rational_damping_idx = 2, &
      & mzero_damping_idx = 3, mrational_damping_idx = 4, op_damping_idx = 5

   has_param(:) = .false.
   disp_energies(:, :) = 0.0_wp
   if (present(dft_energy)) disp_energies(:, :) = dft_energy

   zero_d: block
      use dftd3, only : zero_damping_param, get_zero_damping, new_zero_damping
      type(zero_damping_param) :: param
      call get_zero_damping(inp, method, error, s9=0.0_wp)
      if (allocated(error)) exit zero_d

      has_param(zero_damping_idx) = .true.
      call new_zero_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(1, zero_damping_idx))

      inp%s9 = 1.0_wp
      call new_zero_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(2, zero_damping_idx))
   end block zero_d
   if (allocated(error)) deallocate(error)
end subroutine scan_param_for_reaction

For handling the availability of the parameters we use that the get_zero_damping function fails with populating the error handler. In case we find the error handler is allocated we leave the block, but we clear the error by deallocating it, the success of loading the parameters will be tracked by the has_param array. The dispersion energy computed by our get_dispersion_for_reaction function will be stored in disp_energies. For disp_energies we initialize with the DFT energy if present or zero otherwise and use the fact that the energy dummy argument will be updated rather than overridden.

The damping parameters are retrieved in a simple derived type d3_param, which we can use to initialize the actually damping_param implementation. For zero_damping_param this is the new_zero_damping constructor. The d3_param type is transparent, we can change parameters if we want, for example to evaluate the dispersion energy onces with two-body contributions only and once with non-additive terms added we can override the s9 value and reconstruct the object.

The command line driver#

Now that we have the block for the zero damping finished, we can look into creating our main driver and test our program for the first time. For our command line driver we will go with positional arguments only for a start. For example to evaluate the PBE0-D3 energies for the forth system of S66 from the previous tutorial we would use:

param-scanner PBE0 1 dimer.xyz -1 monomer-a.xyz -1 monomer-b.xyz -0.012069608770

Again we make use of some of the features conveniently provided in mctc_env and mctc_io this time we use the possibility to retrieve command line arguments and the reader for geometry inputs. Also, we will be using the error_type handler to report errors.

With this we can create our full main program as

app/main.f90#

program param_scanner
   use mctc_env, only : wp, error_type, get_argument, fatal_error
   use mctc_io, only : structure_type, read_structure
   use d3_param_scan, only : reaction_type, scan_param_for_reaction
   implicit none

   type(reaction_type) :: reaction
   type(error_type), allocatable :: error
   character(:), allocatable :: method
   real(wp), allocatable :: dft_energy
   integer :: n_args, n_mol, iarg, imol

   n_args = command_argument_count()

   if (n_args < 3) then
      print '(a)', "Usage: param-scanner <method> <coeff1> <mol1> ... [dft energy]"
      stop 1
   end if

   n_mol = (n_args - 1) / 2

   call get_argument(1, method)

   allocate(reaction%mol(n_mol))
   allocate(reaction%stochiometry(n_mol))
   if (allocated(error)) then
      print '(a)', error%message
      stop 1
   end if

end program param_scanner

This should already be possible to run, while it does not do anything, we can already verify our error handling. Let’s test our program with no arguments:

❯ fpm run
Usage: param-scanner <method> <coeff1> <mol1> ... [dft energy]
STOP 1

As a next step we define some helper functions, we want to read a geometry from a command line argument defined by its index. We define a small subroutine which we will include as a contained procedure in our main program:

app/main.f90#

subroutine read_mol(idx, mol, error)
   integer, intent(in) :: idx
   type(structure_type), intent(out) :: mol
   type(error_type), allocatable, intent(out) :: error

   character(len=:), allocatable :: tmp

   call get_argument(idx, tmp)
   call read_structure(mol, tmp, error=error)
end subroutine read_mol

Similarly, we define a procedure to read floating point values from a command line argument index.

app/main.f90#

subroutine read_real(idx, val, error)
   integer, intent(in) :: idx
   real(wp), intent(out) :: val
   type(error_type), allocatable, intent(out) :: error

   character(len=:), allocatable :: tmp
   integer :: stat

   call get_argument(idx, tmp)
   read(tmp, *, iostat=stat) val
   if (stat /= 0) then
      call fatal_error(error, "Could not read floating point value from '"//tmp//"'")
   end if
end subroutine read_real

With this we can add a loop over all our arguments to populate the reaction:

app/main.f90#

   imol = 0
   do iarg = 1, 2 * n_mol, 2
      imol = imol + 1
      call read_real(iarg + 1, reaction%stochiometry(imol), error)
      if (allocated(error)) exit
      call read_mol(iarg + 2, reaction%mol(imol), error)
      if (allocated(error)) exit
   end do

And even already call our scanner function

app/main.f90#

   call scan_param_for_reaction(error, reaction, method, dft_energy)
   if (allocated(error)) then
      print '(a)', error%message
      stop 1
   end if

Adding output for the results#

Running the example command line should work now, but not produce any output yet. Next we are going to look into reporting our results. We add a block to printout our dispersion energies after the zero damping block:

src/scan.f90#

   block
      use mctc_io_convert, only : autokj
      integer :: ipar

      print '(a)', "Energies in kJ/mol"
      print '(66("-"))'
      print '(1x, a, t20, 2a15, a16)', "method", "E(2)", "E(2+3)", "%E(3)"
      print '(66("-"))'
      do ipar = 1, 5
         if (.not.has_param(ipar)) cycle
         print '(1x, a, t20, 3f15.3, "%")', &
            & method // "-" // trim(label(ipar)), &
            & disp_energies(:, ipar) * autokj, &
            & (disp_energies(1, ipar) - disp_energies(2, ipar)) / disp_energies(2, ipar) * 100
      end do
      print '(66("-"))'
   end block

Note

We choose kJ/mol as unit here, but any other unit you prefer can be used here as well.

Also, we need to define the label parameter for making the printout more pretty:

src/scan.f90#

   character(*), parameter :: label(5) = [character(len=20) :: &
      & "D3(0)", "D3(BJ)", "D3M(0)", "D3M(BJ)", "D3(op)"]

Full library code

src/scan.f90#

module d3_param_scan
   use mctc_env, only : wp
   use mctc_io, only : structure_type
   implicit none
   private

   public :: reaction_type
   public :: scan_param_for_reaction

   type :: reaction_type
      type(structure_type), allocatable :: mol(:)
      real(wp), allocatable :: stochiometry(:)
   end type reaction_type

contains

subroutine get_dispersion_for_reaction(reaction, param, energy)
   use dftd3, only : get_dispersion, damping_param, realspace_cutoff, &
      & d3_model, new_d3_model
   type(reaction_type), intent(in) :: reaction
   class(damping_param), intent(in) :: param
   real(wp), intent(inout) :: energy

   integer :: component
   real(wp) :: energy_component
   type(d3_model) :: disp

   do component = 1, size(reaction%mol)
      call new_d3_model(disp, reaction%mol(component))
      call get_dispersion(reaction%mol(component), disp, param, realspace_cutoff(), &
         & energy_component)
      energy = energy + energy_component * reaction%stochiometry(component)
   end do

end subroutine get_dispersion_for_reaction

subroutine scan_param_for_reaction(error, reaction, method, dft_energy)
   use mctc_env, only : error_type, fatal_error
   use dftd3, only : d3_param
   type(error_type), allocatable :: error
   type(reaction_type), intent(in) :: reaction
   character(*), intent(in) :: method
   real(wp), intent(in), optional :: dft_energy

   logical :: has_param(5)
   real(wp) :: disp_energies(2, 5)
   type(d3_param) :: inp
   integer, parameter :: zero_damping_idx = 1, rational_damping_idx = 2, &
      & mzero_damping_idx = 3, mrational_damping_idx = 4, op_damping_idx = 5
   character(*), parameter :: label(5) = [character(len=20) :: &
      & "D3(0)", "D3(BJ)", "D3M(0)", "D3M(BJ)", "D3(op)"]

   has_param(:) = .false.
   disp_energies(:, :) = 0.0_wp
   if (present(dft_energy)) disp_energies(:, :) = dft_energy

   zero_d: block
      use dftd3, only : zero_damping_param, get_zero_damping, new_zero_damping
      type(zero_damping_param) :: param
      call get_zero_damping(inp, method, error, s9=0.0_wp)
      if (allocated(error)) exit zero_d

      has_param(zero_damping_idx) = .true.
      call new_zero_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(1, zero_damping_idx))

      inp%s9 = 1.0_wp
      call new_zero_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(2, zero_damping_idx))
   end block zero_d
   if (allocated(error)) deallocate(error)

   block
      use mctc_io_convert, only : autokj
      integer :: ipar

      print '(a)', "Energies in kJ/mol"
      print '(66("-"))'
      print '(1x, a, t20, 2a15, a16)', "method", "E(2)", "E(2+3)", "%E(3)"
      print '(66("-"))'
      do ipar = 1, 5
         if (.not.has_param(ipar)) cycle
         print '(1x, a, t20, 3f15.3, "%")', &
            & method // "-" // trim(label(ipar)), &
            & disp_energies(:, ipar) * autokj, &
            & (disp_energies(1, ipar) - disp_energies(2, ipar)) / disp_energies(2, ipar) * 100
      end do
      print '(66("-"))'
   end block

end subroutine scan_param_for_reaction

Running our example command should give the following output:

❯ fpm run -- PBE0 1 dimer.xyz -1 monomer-a.xyz -1 monomer-b.xyz
Energies in kJ/mol
------------------------------------------------------------------
 method                       E(2)         E(2+3)           %E(3)
------------------------------------------------------------------
 PBE0-D3(0)                 -4.591         -4.607         -0.350%
------------------------------------------------------------------

Currently, we don’t process the DFT energy for the reaction. Let’s add this feature to our main driver by including

app/main.f90#

   if (.not.allocated(error) .and. 2 * n_mol < n_args - 1) then
      allocate(dft_energy)
      call read_real(n_args, dft_energy, error)
   end if

Now we can try our main program again with the DFT relative energy:

❯ fpm run -- PBE0 1 dimer.xyz -1 monomer-a.xyz -1 monomer-b.xyz -0.012069608770
Energies in kJ/mol
------------------------------------------------------------------
 method                       E(2)         E(2+3)           %E(3)
------------------------------------------------------------------
 PBE0-D3(0)                -36.280        -36.296         -0.044%
------------------------------------------------------------------

Supporting all damping functions#

Finally, we fill in the remaining blocks for the other damping parameters. The procedures needed for this are very similar to the ones used for the zero-damping and therefore we leave this as an exercise. The output with all damping functions should look like

❯ fpm run -- PBE0 1 dimer.xyz -1 monomer-a.xyz -1 monomer-b.xyz -0.012069608770
Energies in kJ/mol
------------------------------------------------------------------
 method                       E(2)         E(2+3)           %E(3)
------------------------------------------------------------------
 PBE0-D3(0)                -36.280        -36.296         -0.044%
 PBE0-D3(BJ)               -35.463        -35.479         -0.045%
 PBE0-D3M(BJ)              -35.677        -35.693         -0.045%
 PBE0-D3(op)               -35.616        -35.632         -0.045%
------------------------------------------------------------------

Note

We have two main categories of damping functions supported with D3. First, rational-type damping which makes the dispersion energy go to a constant value at short distances, in this case constant contributions cancel in reactions if the same pairs are present. Second, zero-type damping where the dispersion energy is fully removed at short distances. Finally, there are versions like the optimized power damping scheme which can be either of rational-type or zero-type damping depending on the power parameters.

We can compare our numbers with the PBE0-D3(BJ) calculations from the last tutorial to confirm our implementation is indeed correct. Overall, we can see a bit of spread in the results for PBE0 with D3 based on this small example case, however the dispersion energies are rougly all in the same range. For deciding on which damping function to choose for a density functional, it is always best to check the performance on the specific systems of interest.

Full library code

src/scan.f90#

module d3_param_scan
   use mctc_env, only : wp
   use mctc_io, only : structure_type
   implicit none
   private

   public :: reaction_type
   public :: scan_param_for_reaction

   type :: reaction_type
      type(structure_type), allocatable :: mol(:)
      real(wp), allocatable :: stochiometry(:)
   end type reaction_type

contains

subroutine get_dispersion_for_reaction(reaction, param, energy)
   use dftd3, only : get_dispersion, damping_param, realspace_cutoff, &
      & d3_model, new_d3_model
   type(reaction_type), intent(in) :: reaction
   class(damping_param), intent(in) :: param
   real(wp), intent(inout) :: energy

   integer :: component
   real(wp) :: energy_component
   type(d3_model) :: disp

   do component = 1, size(reaction%mol)
      call new_d3_model(disp, reaction%mol(component))
      call get_dispersion(reaction%mol(component), disp, param, realspace_cutoff(), &
         & energy_component)
      energy = energy + energy_component * reaction%stochiometry(component)
   end do

end subroutine get_dispersion_for_reaction

subroutine scan_param_for_reaction(error, reaction, method, dft_energy)
   use mctc_env, only : error_type, fatal_error
   use dftd3, only : d3_param
   type(error_type), allocatable :: error
   type(reaction_type), intent(in) :: reaction
   character(*), intent(in) :: method
   real(wp), intent(in), optional :: dft_energy

   logical :: has_param(5)
   real(wp) :: disp_energies(2, 5)
   type(d3_param) :: inp
   integer, parameter :: zero_damping_idx = 1, rational_damping_idx = 2, &
      & mzero_damping_idx = 3, mrational_damping_idx = 4, op_damping_idx = 5
   character(*), parameter :: label(5) = [character(len=20) :: &
      & "D3(0)", "D3(BJ)", "D3M(0)", "D3M(BJ)", "D3(op)"]

   has_param(:) = .false.
   disp_energies(:, :) = 0.0_wp
   if (present(dft_energy)) disp_energies(:, :) = dft_energy

   zero_d: block
      use dftd3, only : zero_damping_param, get_zero_damping, new_zero_damping
      type(zero_damping_param) :: param
      call get_zero_damping(inp, method, error, s9=0.0_wp)
      if (allocated(error)) exit zero_d

      has_param(zero_damping_idx) = .true.
      call new_zero_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(1, zero_damping_idx))

      inp%s9 = 1.0_wp
      call new_zero_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(2, zero_damping_idx))
   end block zero_d
   if (allocated(error)) deallocate(error)

   rational_d: block
      use dftd3, only : rational_damping_param, get_rational_damping, new_rational_damping
      type(rational_damping_param) :: param
      call get_rational_damping(inp, method, error, s9=0.0_wp)
      if (allocated(error)) exit rational_d

      has_param(rational_damping_idx) = .true.
      call new_rational_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(1, rational_damping_idx))

      inp%s9 = 1.0_wp
      call new_rational_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(2, rational_damping_idx))
   end block rational_d
   if (allocated(error)) deallocate(error)

   mzero_d: block
      use dftd3, only : mzero_damping_param, get_mzero_damping, new_mzero_damping
      type(mzero_damping_param) :: param
      call get_mzero_damping(inp, method, error, s9=0.0_wp)
      if (allocated(error)) exit mzero_d

      has_param(zero_damping_idx) = .true.
      call new_mzero_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(1, mzero_damping_idx))

      inp%s9 = 1.0_wp
      call new_mzero_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(2, mzero_damping_idx))
   end block mzero_d
   if (allocated(error)) deallocate(error)

   mrational_d: block
      use dftd3, only : rational_damping_param, get_mrational_damping, new_rational_damping
      type(rational_damping_param) :: param
      call get_mrational_damping(inp, method, error, s9=0.0_wp)
      if (allocated(error)) exit mrational_d

      has_param(mrational_damping_idx) = .true.
      call new_rational_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(1, mrational_damping_idx))

      inp%s9 = 1.0_wp
      call new_rational_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(2, mrational_damping_idx))
   end block mrational_d
   if (allocated(error)) deallocate(error)

   optimizedpower_d: block
      use dftd3, only : optimizedpower_damping_param, get_optimizedpower_damping, &
         & new_optimizedpower_damping
      type(optimizedpower_damping_param) :: param
      call get_optimizedpower_damping(inp, method, error, s9=0.0_wp)
      if (allocated(error)) exit optimizedpower_d

      has_param(op_damping_idx) = .true.
      call new_optimizedpower_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(1, op_damping_idx))

      inp%s9 = 1.0_wp
      call new_optimizedpower_damping(param, inp)
      call get_dispersion_for_reaction(reaction, param, disp_energies(2, op_damping_idx))
   end block optimizedpower_d
   if (allocated(error)) deallocate(error)

   if (.not.any(has_param)) then
      call fatal_error(error, "No parameters found for method '"//method//"'")
      return
   end if

   block
      use mctc_io_convert, only : autokj
      integer :: ipar

      print '(a)', "Energies in kJ/mol"
      print '(66("-"))'
      print '(1x, a, t20, 2a15, a16)', "method", "E(2)", "E(2+3)", "%E(3)"
      print '(66("-"))'
      do ipar = 1, 5
         if (.not.has_param(ipar)) cycle
         print '(1x, a, t20, 3f15.3, "%")', &
            & method // "-" // trim(label(ipar)), &
            & disp_energies(:, ipar) * autokj, &
            & (disp_energies(1, ipar) - disp_energies(2, ipar)) / disp_energies(2, ipar) * 100
      end do
      print '(66("-"))'
   end block

end subroutine scan_param_for_reaction

Summary#

This concludes our tutorial on writing a simple command line tool to scan through damping parameters for D3.

Important

In this tutorial we learned

how to use the dftd3 module and Fortran API to compute dispersion energies
create damping parameters for different damping schemes
how to handle geometry and error types used in dftd3

Testing D3 damping parameters

Contents

Testing D3 damping parameters#

Using the Fortran API#

Computing dispersion for reactions#

Creating a parameter scanner#

The command line driver#

Adding output for the results#

Supporting all damping functions#

Summary#