How to use this library?

How to use this library?#

This section contains a few self-contained examples on how to use D3.

Compute energy with rational damping#

This example shows how to compute the dispersion energy with the rational damping function.

energy.f90#
program test_simple_d3
   use, intrinsic :: iso_fortran_env, only : r8 => real64
   use mctc_env, only: error_type
   use mctc_io, only: structure_type, new
   use dftd3, only: d3_model, d3_param, rational_damping_param, get_rational_damping, &
      & new_rational_damping, new_d3_model, get_dispersion, realspace_cutoff
   implicit none

   character(len=:), allocatable :: method
   type(structure_type) :: mol
   type(error_type), allocatable :: error
   integer, allocatable :: num(:)
   real(r8), allocatable :: xyz(:, :)
   real(r8) :: energy

   type(d3_model) :: disp
   type(d3_param) :: inp
   type(rational_damping_param) :: param

   method = 'PBE0'
   num = [6, 1, 1, 1, 1]
   xyz = reshape([ &  ! coordinates in Bohr
     &  0.0000000_r8, -0.0000000_r8,  0.0000000_r8, &
     & -1.1922080_r8,  1.1922080_r8,  1.1922080_r8, &
     &  1.1922080_r8, -1.1922080_r8,  1.1922080_r8, &
     & -1.1922080_r8, -1.1922080_r8, -1.1922080_r8, &
     &  1.1922080_r8,  1.1922080_r8, -1.1922080_r8],&
     & [3, size(num)])
   call new(mol, num, xyz, charge=0.0_r8, uhf=0)

   call get_rational_damping(inp, method, error, s9=1.0_r8)
   if (allocated(error)) then
      print '(2a)', "Error: ", error%message
      return
   end if
   call new_rational_damping(param, inp)
   call new_d3_model(disp, mol)

   call get_dispersion(mol, disp, param, realspace_cutoff(), energy)
   print '(3a, f13.10, a)', 'Dispersion energy for ', method, '-D3(BJ) is ', energy, ' Hartree'

end program test_simple_d3
energy.c#
#include <stdio.h>
#include <stdbool.h>

#include "dftd3.h"

int main(void)
{
  dftd3_error error = dftd3_new_error();
  dftd3_structure mol = NULL;
  dftd3_model d3 = NULL;
  dftd3_param param = NULL;

  int nat = 5;
  int num[5] = {6, 1, 1, 1, 1};
  double xyz[15] = {  // coordinates in Bohr
     0.0000000, -0.0000000,  0.0000000,
    -1.1922080,  1.1922080,  1.1922080,
     1.1922080, -1.1922080,  1.1922080,
    -1.1922080, -1.1922080, -1.1922080,
     1.1922080,  1.1922080, -1.1922080};
  mol = dftd3_new_structure(error, nat, num, xyz, NULL, NULL);
  if (dftd3_check_error(error)) goto handle_error;

  char method[5] = "PBE0";
  param = dftd3_load_rational_damping(error, method, false);
  if (dftd3_check_error(error)) goto handle_error;

  d3 = dftd3_new_d3_model(error, mol);
  if (dftd3_check_error(error)) goto handle_error;

  double energy;
  dftd3_get_dispersion(error, mol, d3, param, &energy, NULL, NULL);
  if (dftd3_check_error(error)) goto handle_error;

  printf("Dispersion energy for %s-D3(BJ) is %13.10lf Hartree\n", method, energy);

  dftd3_delete(error);
  dftd3_delete(mol);
  dftd3_delete(d3);
  dftd3_delete(param);
  return 0;

handle_error:
  char msg[512];
  dftd3_get_error(error, msg, NULL);
  printf("Error: %s\n", msg);

  dftd3_delete(error);
  dftd3_delete(mol);
  dftd3_delete(d3);
  dftd3_delete(param);
  return 1;
}
energy.py#
import numpy as np
from dftd3.interface import RationalDampingParam, DispersionModel

num = np.array([6, 1, 1, 1, 1])
xyz = np.array(  # coordinates in Bohr
    [
        [ 0.0000000, -0.0000000,  0.0000000],
        [-1.1922080,  1.1922080,  1.1922080],
        [ 1.1922080, -1.1922080,  1.1922080],
        [-1.1922080, -1.1922080, -1.1922080],
        [ 1.1922080,  1.1922080, -1.1922080],
    ]
)
method = "PBE0"

model = DispersionModel(num, xyz)
res = model.get_dispersion(RationalDampingParam(method=method), grad=False)
print(f"Dispersion energy for {method}-D3(BJ) is {res['energy']:13.10f} Hartree")

To test this example you can install the dependencies with

mamba create d3 simple-dftd3 fortran-compiler pkg-config
mamba activate d3

mamba create d3 simple-dftd3 c-compiler pkg-config
mamba activate d3

mamba create d3 dftd3-python
mamba activate d3

You can run the example code with

❯ $FC energy.f90 $(pkg-config s-dftd3 mctc-lib --cflags --libs) && ./a.out
Dispersion energy for PBE0-D3(BJ) is -0.0009218696 Hartree

❯ $CC energy.c $(pkg-config s-dftd3 mctc-lib --cflags --libs) && ./a.out
Dispersion energy for PBE0-D3(BJ) is -0.0009218696 Hartree

❯ python energy.py
Dispersion energy for PBE0-D3(BJ) is -0.0009219059 Hartree

Disable selected atoms with ghost indices#

The same idea works for all public interfaces. Use the atom indices of the fragment that should be excluded from the dispersion contribution, then pass them to the model or CLI entry point:

The Fortran constructor uses 1-based atom indices, while the C and Python interfaces use 0-based indices. For the CLI, the same selection is available as

❯ s-dftd3 run --ghost 4,5,6 structure.xyz