Fortran API

Fortran API#

The dftd4 library seamlessly integrates with other Fortran projects via module interfaces,

Note

Generally, all quantities used in the library are stored in atomic units.

Handling of geometries and structure#

The basic infrastructure to handle molecular and periodic structures is provided by the modular computation tool chain library. The library provides a structure type which is used to represent all geometry related informations in dftd4. A structure type can be constructed from arrays or read from a file.

The constructor is provided with the generic interface new and takes an array of atomic numbers (integer) or element symbols (character(len=*)) as well as the cartesian coordinates in Bohr. Additionally, the molecular charge and the number of unpaired electrons can be provided the charge and uhf keyword, respectively. To create a periodic structure the lattice parameters can be passed as 3 by 3 matrix with the lattice keyword.

An example for using the constructor is given here

subroutine example
   use mctc_env, only : wp
   use mctc_io, only : structure_type, new
   implicit none
   type(structure_type) :: mol
   real(wp), allocatable :: xyz(:, :)
   integer, allocatable :: num(:)

   num = [6, 1, 1, 1, 1]
   xyz = reshape([ &
     &  0.00000000000000_wp, -0.00000000000000_wp,  0.00000000000000_wp, &
     & -1.19220800552211_wp,  1.19220800552211_wp,  1.19220800552211_wp, &
     &  1.19220800552211_wp, -1.19220800552211_wp,  1.19220800552211_wp, &
     & -1.19220800552211_wp, -1.19220800552211_wp, -1.19220800552211_wp, &
     &  1.19220800552211_wp,  1.19220800552211_wp, -1.19220800552211_wp],&
     & [3, size(num)])

   call new(mol, num, xyz, charge=0.0_wp, uhf=0)

   ! ...
end subroutine example

To interact with common input file formats for structures the read_structure procedure is available. The file type is inferred from the name of the file automatically or if a file type hint is provided directly from the enumerator of available file types. The read_structure routine can also use an already opened unit, but in this case the file type hint is mandatory to select the correct format to read from.

subroutine example
   use mctc_env, only : error_type
   use mctc_io, only : structure_type, read_structure, file_type
   implicit none
   type(structure_type) :: mol
   type(error_type), allocatable :: error
   character(len=:), allocatable :: input

   input = "struc.xyz"

   call read_structure(mol, input, error, file_type%xyz)
   if (allocated(error)) then
      print '(a)', error%message
      stop 1
   end if

   ! ...
end subroutine example

The structure type as well as the error type contain only allocatable members and can therefore be used without requiring explicit deconstruction.

Certain members of the structure type should be considered immutable, like the number of atoms (nat), the identifiers for unique atoms (id) and the boundary conditions (periodic). To change those specific structure parameters the structure type and all dependent objects should be reconstructed to ensure a consistent setup. Other properties, like the geometry (xyz), molecular charge (charge), number of unpaired electrons (uhf) and lattice parameters (lattice) can be changed without requiring to reconstruct dependent objects like calculators or restart data.

Error handling#

The basic error handler is an allocatable derived type, available from mctc_env as error_type, which signals an error by its allocation status.

use mctc_env, only : error_type, fatal_error
implicit none
type(error_type), allocatable :: error

call always_ok(error)
if (allocated(error)) then
   print '(a)', "Unexpected failure:", error%message
end if

call always_failed(error)
if (allocated(error)) then
   print '(a)', "Error:", error%message
end if

contains
   subroutine always_ok(error)
      type(error_type), allocatable, intent(out) :: error
   end subroutine always_ok

   subroutine always_failed(error)
      type(error_type), allocatable, intent(out) :: error

      call fatal_error(error, "Message associated with this error")
   end subroutine always_failed
end

An unhandled error might get dropped by the next procedure call.

Performing calculations#

An example for performing a calculation with DFT-D4 is shown below. The realspace_cutoff constructor also accepts optional width2 and width3 values to enable smooth cutoffs for two- and three-body dispersion contributions.

app/main.f90#
use mctc_env, only : wp, error_type, fatal_error
use mctc_io, only : new, structure_type
use dftd4, only : damping_param, get_rational_damping, &
   & get_dispersion, realspace_cutoff, &
   & d4_model, new_d4_model

 1subroutine calc_dftd4(error, mol, method, energy, gradient, sigma)
 2   !> Error handling
 3   type(error_type), allocatable, intent(out) :: error
 4   !> Molecular structure data
 5   type(structure_type), intent(in) :: mol
 6   !> Method name for which parameters should be used
 7   character(len=*), intent(in) :: method
 8   !> Dispersion energy
 9   real(wp), intent(out) :: energy
10   !> Dispersion gradient
11   real(wp), intent(out), contiguous, optional :: gradient(:, :)
12   !> Dispersion virial
13   real(wp), intent(out), contiguous, optional :: sigma(:, :)
14
15   class(damping_param), allocatable :: param
16   type(d4_model) :: disp
17
18   ! Get D4 damping parameters for given method
19   call get_rational_damping(method, param, s9=1.0_wp)
20   if (.not.allocated(param)) then
21      call fatal_error(error, "No parameters for '"//method//"' available")
22      return
23   end if
24
25   ! Initialize D4 model
26   call new_d4_model(error, disp, mol)
27   if (allocated(error)) return
28
29   call get_dispersion(mol, disp, param, realspace_cutoff(), energy, &
30      & gradient, sigma)
31
32end subroutine calc_dftd4
app/main.f90#
use mctc_env, only : wp, error_type, fatal_error
use mctc_io, only : new, structure_type
use dftd4, only : damping_param, get_rational_damping, &
   & get_dispersion, realspace_cutoff, &
   & dispersion_model, new_dispersion_model

 1subroutine calc_dftd4(error, mol, method, energy, model, gradient, sigma)
 2   !> Error handling
 3   type(error_type), allocatable, intent(out) :: error
 4   !> Molecular structure data
 5   type(structure_type), intent(in) :: mol
 6   !> Method name for which parameters should be used
 7   character(len=*), intent(in) :: method
 8   !> Dispersion energy
 9   real(wp), intent(out) :: energy
10   !> Dispersion model (D4 or D4S)
11   character(len=*), intent(in), optional :: model
12   !> Dispersion gradient
13   real(wp), intent(out), contiguous, optional :: gradient(:, :)
14   !> Dispersion virial
15   real(wp), intent(out), contiguous, optional :: sigma(:, :)
16
17   class(damping_param), allocatable :: param
18   class(dispersion_model), allocatable :: disp
19
20   ! Get D4 damping parameters for given method
21   call get_rational_damping(method, param, s9=1.0_wp)
22   if (.not.allocated(param)) then
23      call fatal_error(error, "No parameters for '"//method//"' available")
24      return
25   end if
26
27   ! Initialize D4/D4S model
28   call new_dispersion_model(error, disp, mol, model=model)
29   if (allocated(error)) return
30
31   call get_dispersion(mol, disp, param, realspace_cutoff(), energy, &
32      & gradient, sigma)
33
34end subroutine calc_dftd4

Complete Example#

A minimal program using the snippets from above could look like this:

app/main.f90#
 1program demo
 2   use, intrinsic :: iso_fortran_env, only : error_unit
 3   use mctc_env, only : wp, error_type, fatal_error
 4   use mctc_io, only : new, structure_type
 5   use dftd4, only : damping_param, get_rational_damping, &
 6      & get_dispersion, realspace_cutoff, &
 7      & d4_model, new_d4_model
 8   implicit none
 9
10   type(error_type), allocatable :: error
11   type(structure_type) :: mol
12   real(wp) :: energy
13   real(wp), allocatable :: xyz(:, :), gradient(:, :), sigma(:, :)
14   integer, allocatable :: num(:)
15   integer, parameter :: nat = 5
16
17   allocate(gradient(3, nat), sigma(3, 3))
18   allocate(num(nat), xyz(3, nat))
19
20   ! Initialize molecular structure
21   num(:) = [6, 1, 1, 1, 1]
22   xyz(:, :) = reshape([ &
23     &  0.00000000000000_wp, -0.00000000000000_wp,  0.00000000000000_wp, &
24     & -1.19220800552211_wp,  1.19220800552211_wp,  1.19220800552211_wp, &
25     &  1.19220800552211_wp, -1.19220800552211_wp,  1.19220800552211_wp, &
26     & -1.19220800552211_wp, -1.19220800552211_wp, -1.19220800552211_wp, &
27     &  1.19220800552211_wp,  1.19220800552211_wp, -1.19220800552211_wp], &
28     & [3, size(num)])
29
30   call new(mol, num, xyz, charge=0.0_wp, uhf=0)
31
32   ! Run calculation, check for errors, and print results
33   call calc_dftd4(error, mol, "pbe", energy, gradient=gradient, sigma=sigma)
34
35   if (allocated(error)) then
36      write(error_unit, '("[Error]:", 1x, a)') error%message
37      error stop
38   end if
39
40   write (*,'(a,f18.12)') "D4 dispersion energy (Hartree): ", energy
41
42contains
43
44subroutine calc_dftd4(error, mol, method, energy, gradient, sigma)
45   !> Error handling
46   type(error_type), allocatable, intent(out) :: error
47   !> Molecular structure data
48   type(structure_type), intent(in) :: mol
49   !> Method name for which parameters should be used
50   character(len=*), intent(in) :: method
51   !> Dispersion energy
52   real(wp), intent(out) :: energy
53   !> Dispersion gradient
54   real(wp), intent(out), contiguous, optional :: gradient(:, :)
55   !> Dispersion virial
56   real(wp), intent(out), contiguous, optional :: sigma(:, :)
57
58   class(damping_param), allocatable :: param
59   type(d4_model) :: disp
60
61   ! Get D4 damping parameters for given method
62   call get_rational_damping(method, param, s9=1.0_wp)
63   if (.not.allocated(param)) then
64      call fatal_error(error, "No parameters for '"//method//"' available")
65      return
66   end if
67
68   ! Initialize D4 model
69   call new_d4_model(error, disp, mol)
70   if (allocated(error)) return
71
72   call get_dispersion(mol, disp, param, realspace_cutoff(), energy, &
73      & gradient, sigma)
74
75end subroutine calc_dftd4
76
77end program demo
app/main.f90#
 1program demo
 2   use, intrinsic :: iso_fortran_env, only : error_unit
 3   use mctc_env, only : wp, error_type, fatal_error
 4   use mctc_io, only : new, structure_type
 5   use dftd4, only : damping_param, get_rational_damping, &
 6      & get_dispersion, realspace_cutoff, &
 7      & dispersion_model, new_dispersion_model
 8   implicit none
 9
10   type(error_type), allocatable :: error
11   type(structure_type) :: mol
12   real(wp) :: energy
13   real(wp), allocatable :: xyz(:, :), gradient(:, :), sigma(:, :)
14   integer, allocatable :: num(:)
15   integer, parameter :: nat = 5
16
17   allocate(gradient(3, nat), sigma(3, 3))
18   allocate(num(nat), xyz(3, nat))
19
20   ! Initialize molecular structure
21   num(:) = [6, 1, 1, 1, 1]
22   xyz(:, :) = reshape([ &
23     &  0.00000000000000_wp, -0.00000000000000_wp,  0.00000000000000_wp, &
24     & -1.19220800552211_wp,  1.19220800552211_wp,  1.19220800552211_wp, &
25     &  1.19220800552211_wp, -1.19220800552211_wp,  1.19220800552211_wp, &
26     & -1.19220800552211_wp, -1.19220800552211_wp, -1.19220800552211_wp, &
27     &  1.19220800552211_wp,  1.19220800552211_wp, -1.19220800552211_wp], &
28     & [3, size(num)])
29
30   call new(mol, num, xyz, charge=0.0_wp, uhf=0)
31
32   ! Run calculation, check for errors, and print results
33   call calc_dftd4(error, mol, "pbe", energy, gradient=gradient, &
34      & sigma=sigma, model="D4")
35
36   if (allocated(error)) then
37      write(error_unit, '("[Error]:", 1x, a)') error%message
38      error stop
39   end if
40
41   write (*,'(a,f18.12)') "D4 dispersion energy (Hartree): ", energy
42
43contains
44
45subroutine calc_dftd4(error, mol, method, energy, model, gradient, sigma)
46   !> Error handling
47   type(error_type), allocatable, intent(out) :: error
48   !> Molecular structure data
49   type(structure_type), intent(in) :: mol
50   !> Method name for which parameters should be used
51   character(len=*), intent(in) :: method
52   !> Dispersion energy
53   real(wp), intent(out) :: energy
54   !> Dispersion model (D4 or D4S)
55   character(len=*), intent(in), optional :: model
56   !> Dispersion gradient
57   real(wp), intent(out), contiguous, optional :: gradient(:, :)
58   !> Dispersion virial
59   real(wp), intent(out), contiguous, optional :: sigma(:, :)
60
61   class(damping_param), allocatable :: param
62   class(dispersion_model), allocatable :: disp
63
64   ! Get D4 damping parameters for given method
65   call get_rational_damping(method, param, s9=1.0_wp)
66   if (.not.allocated(param)) then
67      call fatal_error(error, "No parameters for '"//method//"' available")
68      return
69   end if
70
71   ! Initialize D4/D4S model
72   call new_dispersion_model(error, disp, mol, model=model)
73   if (allocated(error)) return
74
75   call get_dispersion(mol, disp, param, realspace_cutoff(), energy, &
76      & gradient, sigma)
77
78end subroutine calc_dftd4
79
80end program demo

The program can be compiled using the following minimal fpm.toml.

name = "api-minimal-3_7_0"

[dependencies]
dftd4.git = "https://github.com/dftd4/dftd4"
dftd4.tag = "v3.7.0"
multicharge.git = "https://github.com/grimme-lab/multicharge.git"
multicharge.tag = "v0.3.0"

The examples can also be found in the assets/examples directory of the repository.