Examples

Examples#

Below, you find the most common usage examples for tad-dftd4.

Single Molecule#

Calculate the DFT-D4 energy of a single molecule.

 1# SPDX-Identifier: CC0-1.0
 2import tad_mctc as mctc
 3import torch
 4
 5import tad_dftd4 as d4
 6
 7numbers = mctc.convert.symbol_to_number(
 8    symbols="C C C C N C S H H H H H".split()
 9)
10
11# coordinates in Bohr
12positions = torch.tensor(
13    [
14        [-2.56745685564671, -0.02509985979910, 0.00000000000000],
15        [-1.39177582455797, +2.27696188880014, 0.00000000000000],
16        [+1.27784995624894, +2.45107479759386, 0.00000000000000],
17        [+2.62801937615793, +0.25927727028120, 0.00000000000000],
18        [+1.41097033661123, -1.99890996077412, 0.00000000000000],
19        [-1.17186102298849, -2.34220576284180, 0.00000000000000],
20        [-2.39505990368378, -5.22635838332362, 0.00000000000000],
21        [+2.41961980455457, -3.62158019253045, 0.00000000000000],
22        [-2.51744374846065, +3.98181713686746, 0.00000000000000],
23        [+2.24269048384775, +4.24389473203647, 0.00000000000000],
24        [+4.66488984573956, +0.17907568006409, 0.00000000000000],
25        [-4.60044244782237, -0.17794734637413, 0.00000000000000],
26    ]
27)
28
29# total charge of the system
30charge = torch.tensor(0.0)
31
32# TPSSh-D4-ATM parameters
33param = d4.damping.Param(
34    s6=positions.new_tensor(1.0),
35    s8=positions.new_tensor(1.85897750),
36    s9=positions.new_tensor(1.0),
37    a1=positions.new_tensor(0.44286966),
38    a2=positions.new_tensor(4.60230534),
39)
40
41# parameters can also be obtained using the functional name:
42# param = d4.get_params(method="d4", functional="tpssh")
43
44energy1 = d4.dftd4(numbers, positions, charge, param)
45
46# class-based interface
47disp = d4.dispersion.DispD4()
48energy2 = disp.calculate(numbers, positions, charge, param)
49
50torch.set_printoptions(precision=10)
51
52ref = torch.tensor(
53    [
54        -0.0020841344,
55        -0.0018971195,
56        -0.0018107513,
57        -0.0018305695,
58        -0.0021737693,
59        -0.0019484236,
60        -0.0022788253,
61        -0.0004080658,
62        -0.0004261866,
63        -0.0004199839,
64        -0.0004280768,
65        -0.0005108935,
66    ]
67)
68assert torch.allclose(energy1, ref, atol=1e-8), "Energy does not match"
69assert torch.allclose(energy2, ref, atol=1e-8), "Energy does not match"
70
71
72print(energy1)
73# tensor([-0.0020841344, -0.0018971195, -0.0018107513, -0.0018305695,
74#         -0.0021737693, -0.0019484236, -0.0022788253, -0.0004080658,
75#         -0.0004261866, -0.0004199839, -0.0004280768, -0.0005108935])

Batched Calculations#

 1# SPDX-Identifier: CC0-1.0
 2import tad_mctc as mctc
 3import torch
 4
 5import tad_dftd4 as d4
 6
 7# S22 system 4: formamide dimer
 8numbers = mctc.batch.pack(
 9    (
10        mctc.convert.symbol_to_number("C C N N H H H H H H O O".split()),
11        mctc.convert.symbol_to_number("C O N H H H".split()),
12    )
13)
14
15# coordinates in Bohr
16positions = mctc.batch.pack(
17    (
18        torch.tensor(
19            [
20                [-3.81469488143921, +0.09993441402912, 0.00000000000000],
21                [+3.81469488143921, -0.09993441402912, 0.00000000000000],
22                [-2.66030049324036, -2.15898251533508, 0.00000000000000],
23                [+2.66030049324036, +2.15898251533508, 0.00000000000000],
24                [-0.73178529739380, -2.28237795829773, 0.00000000000000],
25                [-5.89039325714111, -0.02589114569128, 0.00000000000000],
26                [-3.71254944801331, -3.73605775833130, 0.00000000000000],
27                [+3.71254944801331, +3.73605775833130, 0.00000000000000],
28                [+0.73178529739380, +2.28237795829773, 0.00000000000000],
29                [+5.89039325714111, +0.02589114569128, 0.00000000000000],
30                [-2.74426102638245, +2.16115570068359, 0.00000000000000],
31                [+2.74426102638245, -2.16115570068359, 0.00000000000000],
32            ]
33        ),
34        torch.tensor(
35            [
36                [-0.55569743203406, +1.09030425468557, 0.00000000000000],
37                [+0.51473634678469, +3.15152550263611, 0.00000000000000],
38                [+0.59869690244446, -1.16861263789477, 0.00000000000000],
39                [-0.45355203669134, -2.74568780438064, 0.00000000000000],
40                [+2.52721209544999, -1.29200800956867, 0.00000000000000],
41                [-2.63139587595376, +0.96447869452240, 0.00000000000000],
42            ]
43        ),
44    )
45)
46
47# total charge of both system
48charge = torch.tensor([0.0, 0.0])
49
50# TPSSh-D4-ATM parameters
51param = d4.Param(
52    s6=positions.new_tensor(1.0),
53    s8=positions.new_tensor(1.85897750),
54    s9=positions.new_tensor(1.0),
55    a1=positions.new_tensor(0.44286966),
56    a2=positions.new_tensor(4.60230534),
57)
58
59# calculate dispersion energy in Hartree
60energy = torch.sum(d4.dftd4(numbers, positions, charge, param), -1)
61torch.set_printoptions(precision=10)
62print(energy)
63# tensor([-0.0088341432, -0.0027013607])
64print(energy[0] - 2 * energy[1])
65# tensor(-0.0034314217)

Gradient / Forces#

 1# SPDX-Identifier: CC0-1.0
 2import tad_mctc as mctc
 3import torch
 4
 5import tad_dftd4 as d4
 6
 7dd: mctc.typing.DD = {"device": torch.device("cpu"), "dtype": torch.float64}
 8
 9#####################
10# System Parameters #
11#####################
12
13numbers = mctc.convert.symbol_to_number(
14    symbols="C C C C N C S H H H H H".split()
15)
16
17# coordinates in Bohr
18positions = torch.tensor(
19    [
20        [-2.56745685564671, -0.02509985979910, 0.00000000000000],
21        [-1.39177582455797, +2.27696188880014, 0.00000000000000],
22        [+1.27784995624894, +2.45107479759386, 0.00000000000000],
23        [+2.62801937615793, +0.25927727028120, 0.00000000000000],
24        [+1.41097033661123, -1.99890996077412, 0.00000000000000],
25        [-1.17186102298849, -2.34220576284180, 0.00000000000000],
26        [-2.39505990368378, -5.22635838332362, 0.00000000000000],
27        [+2.41961980455457, -3.62158019253045, 0.00000000000000],
28        [-2.51744374846065, +3.98181713686746, 0.00000000000000],
29        [+2.24269048384775, +4.24389473203647, 0.00000000000000],
30        [+4.66488984573956, +0.17907568006409, 0.00000000000000],
31        [-4.60044244782237, -0.17794734637413, 0.00000000000000],
32    ],
33    **dd
34)
35
36# total charge of the system
37charge = torch.tensor(0.0)
38
39# TPSSh-D4-ATM parameters
40param = d4.get_params(method="d4", functional="tpssh")
41
42
43#######################
44# Analytical Gradient #
45#######################
46
47pos = positions.clone().requires_grad_(True)
48energy = d4.dftd4(numbers, pos, charge, param)
49
50(grad,) = torch.autograd.grad(energy.sum(), pos)
51
52
53######################
54# Numerical Gradient #
55######################
56
57num_grad = torch.zeros(positions.shape, **dd)
58step = 1e-5
59
60for i in range(numbers.shape[-1]):
61    for j in range(3):
62        positions[i, j] += step
63        e1 = d4.dftd4(numbers, positions, charge, param).sum()
64
65        positions[i, j] -= 2 * step
66        e2 = d4.dftd4(numbers, positions, charge, param).sum()
67
68        positions[i, j] += step
69        num_grad[i, j] = (e1 - e2) / (2 * step)
70
71
72# Check if analytical and numerical gradients match
73assert torch.allclose(grad, num_grad, atol=1e-8), "Gradient check failed!"

D4S Model#

This example shows how to use the D4SModel.

 1# SPDX-Identifier: CC0-1.0
 2import tad_mctc as mctc
 3import torch
 4
 5import tad_dftd4 as d4
 6
 7numbers = mctc.convert.symbol_to_number(
 8    symbols="C C C C N C S H H H H H".split()
 9)
10
11# coordinates in Bohr
12positions = torch.tensor(
13    [
14        [-2.56745685564671, -0.02509985979910, 0.00000000000000],
15        [-1.39177582455797, +2.27696188880014, 0.00000000000000],
16        [+1.27784995624894, +2.45107479759386, 0.00000000000000],
17        [+2.62801937615793, +0.25927727028120, 0.00000000000000],
18        [+1.41097033661123, -1.99890996077412, 0.00000000000000],
19        [-1.17186102298849, -2.34220576284180, 0.00000000000000],
20        [-2.39505990368378, -5.22635838332362, 0.00000000000000],
21        [+2.41961980455457, -3.62158019253045, 0.00000000000000],
22        [-2.51744374846065, +3.98181713686746, 0.00000000000000],
23        [+2.24269048384775, +4.24389473203647, 0.00000000000000],
24        [+4.66488984573956, +0.17907568006409, 0.00000000000000],
25        [-4.60044244782237, -0.17794734637413, 0.00000000000000],
26    ]
27)
28
29# total charge of the system
30charge = torch.tensor(0.0)
31
32# Create the D4S model
33model = d4.model.D4SModel(numbers)
34
35param = d4.get_params(method="d4", functional="tpssh")
36energy = d4.dftd4(numbers, positions, charge, param, model=model)
37torch.set_printoptions(precision=10)
38print(energy)
39# tensor([-0.0020843975, -0.0019013016, -0.0018165035, -0.0018363572,
40#         -0.0021877293, -0.0019495023, -0.0022923108, -0.0004326892,
41#         -0.0004439871, -0.0004362087, -0.0004454589, -0.0005344027])
Contents