How to run periodic calculations?#
For a periodic system the two-body dispersion energy is a lattice sum, which has to be either truncated in real space or transformed into reciprocal space. This page shows the three available setups on the same system, the cubic phase of ammonia with r²SCAN-D3(BJ).
Note
The background on the truncation error and how to size a cutoff for a target accuracy is covered in Periodic systems: cutoffs and summation accuracy.
To test the examples 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
Remove the truncation error with reciprocal space summation#
Enabling the low-rank expansion of the reference \(C_6\) coefficients makes the two-body energy separable, so it can be summed over the reciprocal lattice instead. The damped pair potential is bounded at the origin, which lets the reciprocal sum converge exponentially without a real space complement, and the two-body cutoff drops out of the calculation entirely.
This requires three-dimensional periodic boundary conditions and a damping function with a known reciprocal space representation, currently the rational and the zero damping function.
ewald.f90#
program pbc_ewald
use, intrinsic :: iso_fortran_env, only : r8 => real64
use dftd3, only: d3_model, d3_param, rational_damping_param, get_rational_damping, &
& new_rational_damping, new_d3_model, get_dispersion, realspace_cutoff, &
& d3_lowrank_config
use mctc_env, only: error_type
use mctc_io, only: structure_type, new
implicit none
character(len=:), allocatable :: method
type(structure_type) :: mol
type(error_type), allocatable :: error
integer, allocatable :: num(:)
real(r8), allocatable :: xyz(:, :)
real(r8) :: lattice(3, 3), energy
type(d3_model) :: disp
type(d3_param) :: inp
type(rational_damping_param) :: param
method = "r2SCAN"
num = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7]
xyz = reshape([ & ! coordinates in Bohr
& 3.5141339317522400_r8, 2.5620901186651399_r8, 0.9887044918760900_r8, &
& 2.5620901186651399_r8, 0.9887044918760900_r8, 3.5141339317522400_r8, &
& 0.9887044918760900_r8, 3.5141339317522400_r8, 2.5620901186651399_r8, &
& 8.2027323927494997_r8, 2.1265083423321101_r8, 8.3884924301184203_r8, &
& 7.2504996070913199_r8, 3.6998939691211699_r8, 5.8628740176711904_r8, &
& 5.6773029528733403_r8, 1.1742755566739300_r8, 6.8151068033293702_r8, &
& 6.8151068033293702_r8, 5.6773029528733403_r8, 1.1742755566739300_r8, &
& 8.3884924301184203_r8, 8.2027323927494997_r8, 2.1265083423321101_r8, &
& 3.6998939691211699_r8, 5.8628740176711904_r8, 7.2504996070913199_r8, &
& 1.1742755566739300_r8, 6.8151068033293702_r8, 5.6773029528733403_r8, &
& 2.1265083423321101_r8, 8.3884924301184203_r8, 8.2027323927494997_r8, &
& 5.8628740176711904_r8, 7.2504996070913199_r8, 3.6998939691211699_r8, &
& 1.9543543300807500_r8, 1.9543543300807500_r8, 1.9543543300807500_r8, &
& 6.6427638185069302_r8, 2.7342441309165002_r8, 7.4228425919137502_r8, &
& 7.4228425919137502_r8, 6.6427638185069302_r8, 2.7342441309165002_r8, &
& 2.7342441309165002_r8, 7.4228425919137502_r8, 6.6427638185069302_r8], &
& [3, size(num)])
lattice = reshape([ & ! lattice vectors in Bohr
& 9.3771283249512098_r8, 0.0000000000000000_r8, 0.0000000000000000_r8, &
& 0.0000000000000000_r8, 9.3771283249512098_r8, 0.0000000000000000_r8, &
& 0.0000000000000000_r8, 0.0000000000000000_r8, 9.3771283249512098_r8], &
& [3, 3])
call new(mol, num, xyz, lattice=lattice, 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)
! the low-rank expansion of the reference C6 coefficients enables the
! reciprocal space summation, an empty config uses the defaults
call new_d3_model(disp, mol, lowrank=d3_lowrank_config())
! the two-body cutoff is unused here, the remaining cutoffs still apply
call get_dispersion(error, mol, disp, param, realspace_cutoff(), energy)
if (allocated(error)) then
print "(2a)", "Error: ", error%message
return
end if
print "(3a, f13.10, a)", "Dispersion energy for ", method, "-D3(BJ) is ", energy, " Hartree"
end program pbc_ewald
ewald.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 = 16;
int num[16] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7};
double xyz[48] = { // coordinates in Bohr
3.5141339317522400, 2.5620901186651399, 0.9887044918760900,
2.5620901186651399, 0.9887044918760900, 3.5141339317522400,
0.9887044918760900, 3.5141339317522400, 2.5620901186651399,
8.2027323927494997, 2.1265083423321101, 8.3884924301184203,
7.2504996070913199, 3.6998939691211699, 5.8628740176711904,
5.6773029528733403, 1.1742755566739300, 6.8151068033293702,
6.8151068033293702, 5.6773029528733403, 1.1742755566739300,
8.3884924301184203, 8.2027323927494997, 2.1265083423321101,
3.6998939691211699, 5.8628740176711904, 7.2504996070913199,
1.1742755566739300, 6.8151068033293702, 5.6773029528733403,
2.1265083423321101, 8.3884924301184203, 8.2027323927494997,
5.8628740176711904, 7.2504996070913199, 3.6998939691211699,
1.9543543300807500, 1.9543543300807500, 1.9543543300807500,
6.6427638185069302, 2.7342441309165002, 7.4228425919137502,
7.4228425919137502, 6.6427638185069302, 2.7342441309165002,
2.7342441309165002, 7.4228425919137502, 6.6427638185069302};
double lattice[9] = { // lattice vectors in Bohr
9.3771283249512098, 0.0000000000000000, 0.0000000000000000,
0.0000000000000000, 9.3771283249512098, 0.0000000000000000,
0.0000000000000000, 0.0000000000000000, 9.3771283249512098};
mol = dftd3_new_structure(error, nat, num, xyz, lattice, NULL);
if (dftd3_check_error(error)) goto handle_error;
char method[7] = "r2SCAN";
param = dftd3_load_rational_damping(error, method, true);
if (dftd3_check_error(error)) goto handle_error;
d3 = dftd3_new_d3_model(error, mol);
if (dftd3_check_error(error)) goto handle_error;
// the zero arguments select the default rank and reciprocal space cutoff,
// the trailing zero sums over the reciprocal lattice directly
dftd3_set_model_ewald(error, d3, 0, 0.0, 0.0, 0);
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;
}
}
ewald.py#
import numpy as np
from dftd3.interface import RationalDampingParam, DispersionModel
num = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7])
xyz = np.array( # coordinates in Bohr
[
[3.5141339317522400, 2.5620901186651399, 0.9887044918760900],
[2.5620901186651399, 0.9887044918760900, 3.5141339317522400],
[0.9887044918760900, 3.5141339317522400, 2.5620901186651399],
[8.2027323927494997, 2.1265083423321101, 8.3884924301184203],
[7.2504996070913199, 3.6998939691211699, 5.8628740176711904],
[5.6773029528733403, 1.1742755566739300, 6.8151068033293702],
[6.8151068033293702, 5.6773029528733403, 1.1742755566739300],
[8.3884924301184203, 8.2027323927494997, 2.1265083423321101],
[3.6998939691211699, 5.8628740176711904, 7.2504996070913199],
[1.1742755566739300, 6.8151068033293702, 5.6773029528733403],
[2.1265083423321101, 8.3884924301184203, 8.2027323927494997],
[5.8628740176711904, 7.2504996070913199, 3.6998939691211699],
[1.9543543300807500, 1.9543543300807500, 1.9543543300807500],
[6.6427638185069302, 2.7342441309165002, 7.4228425919137502],
[7.4228425919137502, 6.6427638185069302, 2.7342441309165002],
[2.7342441309165002, 7.4228425919137502, 6.6427638185069302],
]
)
lattice = np.array( # lattice vectors in Bohr
[
[9.3771283249512098, 0.0000000000000000, 0.0000000000000000],
[0.0000000000000000, 9.3771283249512098, 0.0000000000000000],
[0.0000000000000000, 0.0000000000000000, 9.3771283249512098],
]
)
method = "r2SCAN"
model = DispersionModel(num, xyz, lattice)
model.set_ewald_summation()
res = model.get_dispersion(RationalDampingParam(method=method, atm=True), grad=False)
print(f"Dispersion energy for {method}-D3(BJ) is {res['energy']:13.10f} Hartree")
❯ $FC ewald.f90 $(pkg-config s-dftd3 mctc-lib --cflags --libs) && ./a.out
Dispersion energy for r2SCAN-D3(BJ) is -0.0122603293 Hartree
❯ $CC ewald.c $(pkg-config s-dftd3 --cflags --libs) && ./a.out
Dispersion energy for r2SCAN-D3(BJ) is -0.0122603293 Hartree
❯ python ewald.py
Dispersion energy for r2SCAN-D3(BJ) is -0.0122603293 Hartree
The same is available from the command line with
❯ s-dftd3 run --ewald --bj r2scan structure.poscar
Note that the pairwise decomposition and the analytical Hessian are only implemented for the real space summation and are rejected for such a model.
Choose how the reciprocal sum is evaluated#
By default the structure factors are interpolated on a mesh with cardinal B-splines and the sum is evaluated with a fast Fourier transform, which costs \(O(N \log N)\). Summing over the reciprocal lattice directly instead visits every atom at every wave vector and costs \(O(N^2)\).
The mesh is derived from the reciprocal cutoff unless it is given explicitly:
! explicit mesh, rounded up to a product of powers of two and three
call new_d3_model(disp, mol, lowrank=d3_lowrank_config(mesh=32))
! sum over the reciprocal lattice directly
call new_d3_model(disp, mol, lowrank=d3_lowrank_config(mesh=-1))
/* explicit mesh, rounded up to a product of powers of two and three */
dftd3_set_model_ewald(error, d3, 0, 0.0, 0.0, 32);
/* sum over the reciprocal lattice directly */
dftd3_set_model_ewald(error, d3, 0, 0.0, 0.0, -1);
# explicit mesh, rounded up to a product of powers of two and three
model.set_ewald_summation(mesh=32)
# sum over the reciprocal lattice directly
model.set_ewald_summation(mesh=-1)
Both reach the same energy for the cell above
Dispersion energy for r2SCAN-D3(BJ) is -0.0122603293 Hartree
From the command line the mesh is selected with --ewald-mesh and the direct summation with --ewald-direct
❯ s-dftd3 run --ewald-mesh 32 --bj r2scan structure.poscar
❯ s-dftd3 run --ewald-direct --bj r2scan structure.poscar
See Periodic systems: cutoffs and summation accuracy for how to choose the mesh.
Smooth the real space cutoff#
A hard cutoff makes the energy jump whenever an atom pair crosses it, which shows up as noise in a geometry optimization or a lattice scan. Giving the cutoff a smoothing width switches the pair contributions off over the last fraction of a Bohr instead, so the energy stays continuous. The width is the trailing argument for the two- and three-body cutoff and only changes the energy by the size of the contributions it damps.
smooth.f90#
program pbc_smooth
use, intrinsic :: iso_fortran_env, only : r8 => real64
use dftd3, only: d3_model, d3_param, rational_damping_param, get_rational_damping, &
& new_rational_damping, new_d3_model, get_dispersion, realspace_cutoff
use mctc_env, only: error_type
use mctc_io, only: structure_type, new
implicit none
character(len=:), allocatable :: method
type(structure_type) :: mol
type(error_type), allocatable :: error
integer, allocatable :: num(:)
real(r8), allocatable :: xyz(:, :)
real(r8) :: lattice(3, 3), energy
type(d3_model) :: disp
type(d3_param) :: inp
type(rational_damping_param) :: param
method = "r2SCAN"
num = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7]
xyz = reshape([ & ! coordinates in Bohr
& 3.5141339317522400_r8, 2.5620901186651399_r8, 0.9887044918760900_r8, &
& 2.5620901186651399_r8, 0.9887044918760900_r8, 3.5141339317522400_r8, &
& 0.9887044918760900_r8, 3.5141339317522400_r8, 2.5620901186651399_r8, &
& 8.2027323927494997_r8, 2.1265083423321101_r8, 8.3884924301184203_r8, &
& 7.2504996070913199_r8, 3.6998939691211699_r8, 5.8628740176711904_r8, &
& 5.6773029528733403_r8, 1.1742755566739300_r8, 6.8151068033293702_r8, &
& 6.8151068033293702_r8, 5.6773029528733403_r8, 1.1742755566739300_r8, &
& 8.3884924301184203_r8, 8.2027323927494997_r8, 2.1265083423321101_r8, &
& 3.6998939691211699_r8, 5.8628740176711904_r8, 7.2504996070913199_r8, &
& 1.1742755566739300_r8, 6.8151068033293702_r8, 5.6773029528733403_r8, &
& 2.1265083423321101_r8, 8.3884924301184203_r8, 8.2027323927494997_r8, &
& 5.8628740176711904_r8, 7.2504996070913199_r8, 3.6998939691211699_r8, &
& 1.9543543300807500_r8, 1.9543543300807500_r8, 1.9543543300807500_r8, &
& 6.6427638185069302_r8, 2.7342441309165002_r8, 7.4228425919137502_r8, &
& 7.4228425919137502_r8, 6.6427638185069302_r8, 2.7342441309165002_r8, &
& 2.7342441309165002_r8, 7.4228425919137502_r8, 6.6427638185069302_r8], &
& [3, size(num)])
lattice = reshape([ & ! lattice vectors in Bohr
& 9.3771283249512098_r8, 0.0000000000000000_r8, 0.0000000000000000_r8, &
& 0.0000000000000000_r8, 9.3771283249512098_r8, 0.0000000000000000_r8, &
& 0.0000000000000000_r8, 0.0000000000000000_r8, 9.3771283249512098_r8], &
& [3, 3])
call new(mol, num, xyz, lattice=lattice, 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)
! the widths switch the pair contributions off over the last 0.05 Bohr
! before the respective cutoff
call get_dispersion(mol, disp, param, &
& realspace_cutoff(width2=0.05_r8, width3=0.05_r8), energy)
print "(3a, f13.10, a)", "Dispersion energy for ", method, "-D3(BJ) is ", energy, " Hartree"
end program pbc_smooth
smooth.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 = 16;
int num[16] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7};
double xyz[48] = { // coordinates in Bohr
3.5141339317522400, 2.5620901186651399, 0.9887044918760900,
2.5620901186651399, 0.9887044918760900, 3.5141339317522400,
0.9887044918760900, 3.5141339317522400, 2.5620901186651399,
8.2027323927494997, 2.1265083423321101, 8.3884924301184203,
7.2504996070913199, 3.6998939691211699, 5.8628740176711904,
5.6773029528733403, 1.1742755566739300, 6.8151068033293702,
6.8151068033293702, 5.6773029528733403, 1.1742755566739300,
8.3884924301184203, 8.2027323927494997, 2.1265083423321101,
3.6998939691211699, 5.8628740176711904, 7.2504996070913199,
1.1742755566739300, 6.8151068033293702, 5.6773029528733403,
2.1265083423321101, 8.3884924301184203, 8.2027323927494997,
5.8628740176711904, 7.2504996070913199, 3.6998939691211699,
1.9543543300807500, 1.9543543300807500, 1.9543543300807500,
6.6427638185069302, 2.7342441309165002, 7.4228425919137502,
7.4228425919137502, 6.6427638185069302, 2.7342441309165002,
2.7342441309165002, 7.4228425919137502, 6.6427638185069302};
double lattice[9] = { // lattice vectors in Bohr
9.3771283249512098, 0.0000000000000000, 0.0000000000000000,
0.0000000000000000, 9.3771283249512098, 0.0000000000000000,
0.0000000000000000, 0.0000000000000000, 9.3771283249512098};
mol = dftd3_new_structure(error, nat, num, xyz, lattice, NULL);
if (dftd3_check_error(error)) goto handle_error;
char method[7] = "r2SCAN";
param = dftd3_load_rational_damping(error, method, true);
if (dftd3_check_error(error)) goto handle_error;
d3 = dftd3_new_d3_model(error, mol);
if (dftd3_check_error(error)) goto handle_error;
// the two trailing arguments are the smoothing widths, all values in Bohr
dftd3_set_model_realspace_cutoff_smooth(error, d3, 60.0, 40.0, 40.0, 0.05, 0.05);
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;
}
}
smooth.py#
import numpy as np
from dftd3.interface import RationalDampingParam, DispersionModel
num = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7])
xyz = np.array( # coordinates in Bohr
[
[3.5141339317522400, 2.5620901186651399, 0.9887044918760900],
[2.5620901186651399, 0.9887044918760900, 3.5141339317522400],
[0.9887044918760900, 3.5141339317522400, 2.5620901186651399],
[8.2027323927494997, 2.1265083423321101, 8.3884924301184203],
[7.2504996070913199, 3.6998939691211699, 5.8628740176711904],
[5.6773029528733403, 1.1742755566739300, 6.8151068033293702],
[6.8151068033293702, 5.6773029528733403, 1.1742755566739300],
[8.3884924301184203, 8.2027323927494997, 2.1265083423321101],
[3.6998939691211699, 5.8628740176711904, 7.2504996070913199],
[1.1742755566739300, 6.8151068033293702, 5.6773029528733403],
[2.1265083423321101, 8.3884924301184203, 8.2027323927494997],
[5.8628740176711904, 7.2504996070913199, 3.6998939691211699],
[1.9543543300807500, 1.9543543300807500, 1.9543543300807500],
[6.6427638185069302, 2.7342441309165002, 7.4228425919137502],
[7.4228425919137502, 6.6427638185069302, 2.7342441309165002],
[2.7342441309165002, 7.4228425919137502, 6.6427638185069302],
]
)
lattice = np.array( # lattice vectors in Bohr
[
[9.3771283249512098, 0.0000000000000000, 0.0000000000000000],
[0.0000000000000000, 9.3771283249512098, 0.0000000000000000],
[0.0000000000000000, 0.0000000000000000, 9.3771283249512098],
]
)
method = "r2SCAN"
model = DispersionModel(num, xyz, lattice)
model.set_realspace_cutoff(60.0, 40.0, 40.0, width2=0.05, width3=0.05)
res = model.get_dispersion(RationalDampingParam(method=method, atm=True), grad=False)
print(f"Dispersion energy for {method}-D3(BJ) is {res['energy']:13.10f} Hartree")
❯ $FC smooth.f90 $(pkg-config s-dftd3 mctc-lib --cflags --libs) && ./a.out
Dispersion energy for r2SCAN-D3(BJ) is -0.0122446362 Hartree
❯ $CC smooth.c $(pkg-config s-dftd3 --cflags --libs) && ./a.out
Dispersion energy for r2SCAN-D3(BJ) is -0.0122446362 Hartree
❯ python smooth.py
Dispersion energy for r2SCAN-D3(BJ) is -0.0122446362 Hartree
Smoothing removes the discontinuity, it does not remove the truncation error. The result still differs from the reciprocal space value by about 1.6 · 10⁻⁵ Hartree for this cell.
Truncate the real space sum#
Without any further setup the lattice sum is truncated at the default cutoffs, 60 Bohr for the two-body term and 40 Bohr for the three-body term and the coordination number. This is the cheapest option and the one used by every non-periodic interface.
energy.f90#
program pbc_energy
use, intrinsic :: iso_fortran_env, only : r8 => real64
use dftd3, only: d3_model, d3_param, rational_damping_param, get_rational_damping, &
& new_rational_damping, new_d3_model, get_dispersion, realspace_cutoff
use mctc_env, only: error_type
use mctc_io, only: structure_type, new
implicit none
character(len=:), allocatable :: method
type(structure_type) :: mol
type(error_type), allocatable :: error
integer, allocatable :: num(:)
real(r8), allocatable :: xyz(:, :)
real(r8) :: lattice(3, 3), energy
type(d3_model) :: disp
type(d3_param) :: inp
type(rational_damping_param) :: param
method = "r2SCAN"
num = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7]
xyz = reshape([ & ! coordinates in Bohr
& 3.5141339317522400_r8, 2.5620901186651399_r8, 0.9887044918760900_r8, &
& 2.5620901186651399_r8, 0.9887044918760900_r8, 3.5141339317522400_r8, &
& 0.9887044918760900_r8, 3.5141339317522400_r8, 2.5620901186651399_r8, &
& 8.2027323927494997_r8, 2.1265083423321101_r8, 8.3884924301184203_r8, &
& 7.2504996070913199_r8, 3.6998939691211699_r8, 5.8628740176711904_r8, &
& 5.6773029528733403_r8, 1.1742755566739300_r8, 6.8151068033293702_r8, &
& 6.8151068033293702_r8, 5.6773029528733403_r8, 1.1742755566739300_r8, &
& 8.3884924301184203_r8, 8.2027323927494997_r8, 2.1265083423321101_r8, &
& 3.6998939691211699_r8, 5.8628740176711904_r8, 7.2504996070913199_r8, &
& 1.1742755566739300_r8, 6.8151068033293702_r8, 5.6773029528733403_r8, &
& 2.1265083423321101_r8, 8.3884924301184203_r8, 8.2027323927494997_r8, &
& 5.8628740176711904_r8, 7.2504996070913199_r8, 3.6998939691211699_r8, &
& 1.9543543300807500_r8, 1.9543543300807500_r8, 1.9543543300807500_r8, &
& 6.6427638185069302_r8, 2.7342441309165002_r8, 7.4228425919137502_r8, &
& 7.4228425919137502_r8, 6.6427638185069302_r8, 2.7342441309165002_r8, &
& 2.7342441309165002_r8, 7.4228425919137502_r8, 6.6427638185069302_r8], &
& [3, size(num)])
lattice = reshape([ & ! lattice vectors in Bohr
& 9.3771283249512098_r8, 0.0000000000000000_r8, 0.0000000000000000_r8, &
& 0.0000000000000000_r8, 9.3771283249512098_r8, 0.0000000000000000_r8, &
& 0.0000000000000000_r8, 0.0000000000000000_r8, 9.3771283249512098_r8], &
& [3, 3])
call new(mol, num, xyz, lattice=lattice, 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)
! the defaults truncate the two-body sum at 60 Bohr and the coordination
! number at 40 Bohr, all cutoffs are given in Bohr
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 pbc_energy
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 = 16;
int num[16] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7};
double xyz[48] = { // coordinates in Bohr
3.5141339317522400, 2.5620901186651399, 0.9887044918760900,
2.5620901186651399, 0.9887044918760900, 3.5141339317522400,
0.9887044918760900, 3.5141339317522400, 2.5620901186651399,
8.2027323927494997, 2.1265083423321101, 8.3884924301184203,
7.2504996070913199, 3.6998939691211699, 5.8628740176711904,
5.6773029528733403, 1.1742755566739300, 6.8151068033293702,
6.8151068033293702, 5.6773029528733403, 1.1742755566739300,
8.3884924301184203, 8.2027323927494997, 2.1265083423321101,
3.6998939691211699, 5.8628740176711904, 7.2504996070913199,
1.1742755566739300, 6.8151068033293702, 5.6773029528733403,
2.1265083423321101, 8.3884924301184203, 8.2027323927494997,
5.8628740176711904, 7.2504996070913199, 3.6998939691211699,
1.9543543300807500, 1.9543543300807500, 1.9543543300807500,
6.6427638185069302, 2.7342441309165002, 7.4228425919137502,
7.4228425919137502, 6.6427638185069302, 2.7342441309165002,
2.7342441309165002, 7.4228425919137502, 6.6427638185069302};
double lattice[9] = { // lattice vectors in Bohr
9.3771283249512098, 0.0000000000000000, 0.0000000000000000,
0.0000000000000000, 9.3771283249512098, 0.0000000000000000,
0.0000000000000000, 0.0000000000000000, 9.3771283249512098};
mol = dftd3_new_structure(error, nat, num, xyz, lattice, NULL);
if (dftd3_check_error(error)) goto handle_error;
char method[7] = "r2SCAN";
param = dftd3_load_rational_damping(error, method, true);
if (dftd3_check_error(error)) goto handle_error;
d3 = dftd3_new_d3_model(error, mol);
if (dftd3_check_error(error)) goto handle_error;
// the defaults truncate the two-body sum at 60 Bohr and the coordination
// number at 40 Bohr, all cutoffs are given in Bohr
dftd3_set_model_realspace_cutoff(error, d3, 60.0, 40.0, 40.0);
if (dftd3_check_error(error)) goto handle_error;
// to get the gradient and virial, provide arrays for the output
double energy, gradient[48], virial[9];
dftd3_get_dispersion(error, mol, d3, param, &energy, gradient, virial);
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([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7])
xyz = np.array( # coordinates in Bohr
[
[3.5141339317522400, 2.5620901186651399, 0.9887044918760900],
[2.5620901186651399, 0.9887044918760900, 3.5141339317522400],
[0.9887044918760900, 3.5141339317522400, 2.5620901186651399],
[8.2027323927494997, 2.1265083423321101, 8.3884924301184203],
[7.2504996070913199, 3.6998939691211699, 5.8628740176711904],
[5.6773029528733403, 1.1742755566739300, 6.8151068033293702],
[6.8151068033293702, 5.6773029528733403, 1.1742755566739300],
[8.3884924301184203, 8.2027323927494997, 2.1265083423321101],
[3.6998939691211699, 5.8628740176711904, 7.2504996070913199],
[1.1742755566739300, 6.8151068033293702, 5.6773029528733403],
[2.1265083423321101, 8.3884924301184203, 8.2027323927494997],
[5.8628740176711904, 7.2504996070913199, 3.6998939691211699],
[1.9543543300807500, 1.9543543300807500, 1.9543543300807500],
[6.6427638185069302, 2.7342441309165002, 7.4228425919137502],
[7.4228425919137502, 6.6427638185069302, 2.7342441309165002],
[2.7342441309165002, 7.4228425919137502, 6.6427638185069302],
]
)
lattice = np.array(
[
[9.3771283249512098, 0.0000000000000000, 0.0000000000000000],
[0.0000000000000000, 9.3771283249512098, 0.0000000000000000],
[0.0000000000000000, 0.0000000000000000, 9.3771283249512098],
]
)
method = "r2SCAN"
model = DispersionModel(num, xyz, lattice)
res = model.get_dispersion(RationalDampingParam(method=method, atm=True), grad=False)
print(f"Dispersion energy for {method}-D3(BJ) is {res['energy']:13.10f} Hartree")
❯ $FC energy.f90 $(pkg-config s-dftd3 mctc-lib --cflags --libs) && ./a.out
Dispersion energy for r2SCAN-D3(BJ) is -0.0122446595 Hartree
❯ $CC energy.c $(pkg-config s-dftd3 --cflags --libs) && ./a.out
Dispersion energy for r2SCAN-D3(BJ) is -0.0122446595 Hartree
❯ python energy.py
Dispersion energy for r2SCAN-D3(BJ) is -0.0122446595 Hartree
Raising the two-body cutoff reduces the remaining error, see Periodic systems: cutoffs and summation accuracy for choosing a cutoff for a target accuracy.