Installation

pip install tad-dftd3

conda install tad-dftd3

git clone https://github.com/dftd3/tad-dftd3
cd tad-dftd3

mamba env create -n torch -f environment.yml
mamba activate torch

pip install .

Development

For development, additionally install the following tools in your environment.

mamba install black covdefaults mypy pre-commit pylint pytest pytest-cov pytest-xdist tox
pip install pytest-random-order

With pip, add the option -e for installing in development mode, and add [dev] for the development dependencies

pip install -e .[dev]

The pre-commit hooks are initialized by running the following command in the root of the repository.

pre-commit install

For testing all Python environments, simply run tox.

tox

Note that this randomizes the order of tests but skips “large” tests. To modify this behavior, tox has to skip the optional posargs.

tox -- test

Testing against the s-dftd3 reference

Some tests compare tad-dftd3’s output against reference values computed at test time by calling the s-dftd3 Fortran implementation through its Python bindings, rather than against numbers stored in this repository. The dftd3 package is therefore a required test dependency, not an optional one – it is pulled in by the [dev] extra, and there is no offline fallback if it is missing.

Examples

All examples can also be found in the examples directory.

• single.py

• batch.py

• forces.py

• hessian.py

The following example shows how to calculate the DFT-D3 dispersion energy for a single structure.

import torch
import tad_dftd3 as d3
import tad_mctc as mctc

numbers = mctc.convert.symbol_to_number(symbols="C C C C N C S H H H H H".split())
positions = torch.tensor(
    [
        [-2.56745685564671, -0.02509985979910, 0.00000000000000],
        [-1.39177582455797, +2.27696188880014, 0.00000000000000],
        [+1.27784995624894, +2.45107479759386, 0.00000000000000],
        [+2.62801937615793, +0.25927727028120, 0.00000000000000],
        [+1.41097033661123, -1.99890996077412, 0.00000000000000],
        [-1.17186102298849, -2.34220576284180, 0.00000000000000],
        [-2.39505990368378, -5.22635838332362, 0.00000000000000],
        [+2.41961980455457, -3.62158019253045, 0.00000000000000],
        [-2.51744374846065, +3.98181713686746, 0.00000000000000],
        [+2.24269048384775, +4.24389473203647, 0.00000000000000],
        [+4.66488984573956, +0.17907568006409, 0.00000000000000],
        [-4.60044244782237, -0.17794734637413, 0.00000000000000],
    ]
)
param = {
    "a1": torch.tensor(0.49484001),
    "s8": torch.tensor(0.78981345),
    "a2": torch.tensor(5.73083694),
}

energy = d3.dftd3(numbers, positions, param)

torch.set_printoptions(precision=10)
print(energy)
# tensor([-0.0004075971, -0.0003940886, -0.0003817684, -0.0003949536,
#         -0.0003577212, -0.0004110279, -0.0005385976, -0.0001808242,
#         -0.0001563670, -0.0001503394, -0.0001577045, -0.0001764488])

The next example shows the calculation of dispersion energies for a batch of structures, while retaining access to all intermediates used for calculating the dispersion energy.

import torch
import tad_dftd3 as d3
import tad_mctc as mctc

sample1 = dict(
    numbers=mctc.convert.symbol_to_number("Pb H H H H Bi H H H".split()),
    positions=torch.tensor(
        [
            [-0.00000020988889, -4.98043478877778, +0.00000000000000],
            [+3.06964045311111, -6.06324400177778, +0.00000000000000],
            [-1.53482054188889, -6.06324400177778, -2.65838526500000],
            [-1.53482054188889, -6.06324400177778, +2.65838526500000],
            [-0.00000020988889, -1.72196703577778, +0.00000000000000],
            [-0.00000020988889, +4.77334244722222, +0.00000000000000],
            [+1.35700257511111, +6.70626379422222, -2.35039772300000],
            [-2.71400388988889, +6.70626379422222, +0.00000000000000],
            [+1.35700257511111, +6.70626379422222, +2.35039772300000],
        ]
    ),
)
sample2 = dict(
    numbers=mctc.convert.symbol_to_number(
        "C C C C C C I H H H H H S H C H H H".split(" ")
    ),
    positions=torch.tensor(
        [
            [-1.42754169820131, -1.50508961850828, -1.93430551124333],
            [+1.19860572924150, -1.66299114873979, -2.03189643761298],
            [+2.65876001301880, +0.37736955363609, -1.23426391650599],
            [+1.50963368042358, +2.57230374419743, -0.34128058818180],
            [-1.12092277855371, +2.71045691257517, -0.25246348639234],
            [-2.60071517756218, +0.67879949508239, -1.04550707592673],
            [-2.86169588073340, +5.99660765711210, +1.08394899986031],
            [+2.09930989272956, -3.36144811062374, -2.72237695164263],
            [+2.64405246349916, +4.15317840474646, +0.27856972788526],
            [+4.69864865613751, +0.26922271535391, -1.30274048619151],
            [-4.63786461351839, +0.79856258572808, -0.96906659938432],
            [-2.57447518692275, -3.08132039046931, -2.54875517521577],
            [-5.88211879210329, 11.88491819358157, +2.31866455902233],
            [-8.18022701418703, 10.95619984550779, +1.83940856333092],
            [-5.08172874482867, 12.66714386256482, -0.92419491629867],
            [-3.18311711399702, 13.44626574330220, -0.86977613647871],
            [-5.07177399637298, 10.99164969235585, -2.10739192258756],
            [-6.35955320518616, 14.08073002965080, -1.68204314084441],
        ]
    ),
)
numbers = mctc.batch.pack(
    (
        sample1["numbers"],
        sample2["numbers"],
    )
)
positions = mctc.batch.pack(
    (
        sample1["positions"],
        sample2["positions"],
    )
)
ref = d3.reference.Reference()
rcov = mctc.data.COV_D3()[numbers]
rvdw = mctc.data.VDW_PAIRWISE()[numbers.unsqueeze(-1), numbers.unsqueeze(-2)]
r4r2 = d3.data.R4R2()[numbers]
param = {
    "a1": torch.tensor(0.49484001),
    "s8": torch.tensor(0.78981345),
    "a2": torch.tensor(5.73083694),
}

cn = mctc.ncoord.cn_d3(
    numbers, positions, counting_function=mctc.ncoord.exp_count, rcov=rcov
)
weights = d3.model.weight_references(numbers, cn, ref, d3.model.gaussian_weight)
c6 = d3.model.atomic_c6(numbers, weights, ref)
energy = d3.disp.dispersion(
    numbers,
    positions,
    param,
    c6,
    rvdw,
    r4r2,
    d3.disp.rational_damping,
)

torch.set_printoptions(precision=10)
print(torch.sum(energy, dim=-1))
# tensor([-0.0014092578, -0.0057840119])

Since the dispersion energy is differentiable with respect to the atomic positions, the D3 contribution to the gradient (and hence the forces) is obtained from a simple backward pass.

import torch
import tad_dftd3 as d3
import tad_mctc as mctc

numbers = mctc.convert.symbol_to_number(symbols="C C C C N C S H H H H H".split())
positions = torch.tensor(
    [
        [-2.56745685564671, -0.02509985979910, 0.00000000000000],
        [-1.39177582455797, +2.27696188880014, 0.00000000000000],
        [+1.27784995624894, +2.45107479759386, 0.00000000000000],
        [+2.62801937615793, +0.25927727028120, 0.00000000000000],
        [+1.41097033661123, -1.99890996077412, 0.00000000000000],
        [-1.17186102298849, -2.34220576284180, 0.00000000000000],
        [-2.39505990368378, -5.22635838332362, 0.00000000000000],
        [+2.41961980455457, -3.62158019253045, 0.00000000000000],
        [-2.51744374846065, +3.98181713686746, 0.00000000000000],
        [+2.24269048384775, +4.24389473203647, 0.00000000000000],
        [+4.66488984573956, +0.17907568006409, 0.00000000000000],
        [-4.60044244782237, -0.17794734637413, 0.00000000000000],
    ],
    requires_grad=True,
)
param = {
    "a1": torch.tensor(0.49484001),
    "s8": torch.tensor(0.78981345),
    "a2": torch.tensor(5.73083694),
}

energy = d3.dftd3(numbers, positions, param)

(grad,) = torch.autograd.grad(energy.sum(), positions)
forces = -grad

The full example, including a comparison against numerical gradients, is available in forces.py. Second derivatives, i.e., the D3 contribution to the Hessian, are obtained by applying reverse-mode automatic differentiation twice (see hessian.py).

Limitations

The current implementation only works for molecular structures. Periodic boundary conditions are not implemented, i.e., no stress tensor or lattice gradient is available.

The code is fully vectorized for maximum efficiency. Therefore, all quantities are stored as full tensors, which makes calculations rather memory intensive. Especially, the ATM term can become limiting as it requires a 3D tensor of dimension (n_atoms, n_atoms, n_atoms).

Contributing

This is a volunteer open source projects and contributions are always welcome. Please, take a moment to read the contributing guidelines.

License

Licensed under the Apache License, Version 2.0 (the “License”); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an “as is” basis, without warranties or conditions of any kind, either express or implied. See the License for the specific language governing permissions and limitations under the License.

Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in this project by you, as defined in the Apache-2.0 license, shall be licensed as above, without any additional terms or conditions.