Skip to content

Commit dc75007

Browse files
committed
Support PySCF
1 parent 4a95d4e commit dc75007

File tree

7 files changed

+361
-8
lines changed

7 files changed

+361
-8
lines changed

doc/api/python.rst

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -80,6 +80,12 @@ ASE support
8080
.. automodule:: dftd3.ase
8181

8282

83+
PySCF support
84+
------------
85+
86+
.. automodule:: dftd3.pyscf
87+
88+
8389
Literature
8490
----------
8591

python/README.rst

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -70,6 +70,58 @@ If the QCElemental package is installed the ``dftd3.qcschema`` module becomes im
7070
# => -0.0004204244108151285
7171
7272
73+
ASE Integration
74+
---------------
75+
76+
To integrate with `ASE <https://wiki.fysik.dtu.dk/ase/>`_ this interface implements an ASE Calculator.
77+
The ``DFTD3`` calculator becomes importable if an ASE installation is available.
78+
79+
.. code:: python
80+
81+
>>> from ase.build import molecule
82+
>>> from dftd3.ase import DFTD3
83+
>>> atoms = molecule('H2O')
84+
>>> atoms.calc = DFTD3(method="TPSS", damping="d3bj")
85+
>>> atoms.get_potential_energy()
86+
-0.0114416338147162
87+
>>> atoms.calc.set(method="PBE")
88+
{'method': 'PBE'}
89+
>>> atoms.get_potential_energy()
90+
-0.009781913226281063
91+
>>> atoms.get_forces()
92+
array([[-0.00000000e+00 -0.00000000e+00 9.56568982e-05]
93+
[-0.00000000e+00 -4.06046858e-05 -4.78284491e-05]
94+
[-0.00000000e+00 4.06046858e-05 -4.78284491e-05]])
95+
96+
To use the ``DFTD3`` calculator as dispersion correction the calculator can be combined using the `SumCalculator <https://wiki.fysik.dtu.dk/ase/ase/calculators/mixing.html>`_ from the ``ase.calculators.mixing`` module.
97+
98+
.. code:: python
99+
100+
>>> from ase.build import molecule
101+
>>> from ase.calculators.mixing import SumCalculator
102+
>>> from ase.calculators.nwchem import NWChem
103+
>>> from dftd3.ase import DFTD3
104+
>>> atoms = molecule('H2O')
105+
>>> atoms.calc = SumCalculator([DFTD3(method="PBE", damping="d3bj"), NWChem(xc="PBE")])
106+
107+
For convenience ``DFTD3`` allows to combine itself with another calculator by using the ``add_calculator`` method which returns a SumCalculator:
108+
109+
.. code:: python
110+
111+
>>> from ase.build import molecule
112+
>>> from ase.calculators.emt import EMT
113+
>>> from dftd4.ase import DFTD3
114+
>>> atoms = molecule("C60")
115+
>>> atoms.calc = DFTD3(method="pbe", damping="d3bj").add_calculator(EMT())
116+
>>> atoms.get_potential_energy()
117+
7.513593999944228
118+
>>> [calc.get_potential_energy() for calc in atoms.calc.calcs]
119+
[-4.850025823367818, 12.363619823312046]
120+
121+
The individual contributions are available by iterating over the list of calculators in ``calc.calcs``.
122+
Note that ``DFTD3`` will always place itself as first calculator in the list.
123+
124+
73125
Installing
74126
----------
75127

python/dftd3/meson.build

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -62,11 +62,13 @@ if install
6262
'library.py',
6363
'interface.py',
6464
'parameters.py',
65+
'pyscf.py',
6566
'qcschema.py',
6667
'test_ase.py',
6768
'test_library.py',
6869
'test_interface.py',
6970
'test_parameters.py',
71+
'test_pyscf.py',
7072
'test_qcschema.py',
7173
subdir: 'dftd3',
7274
)

python/dftd3/pyscf.py

Lines changed: 198 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,198 @@
1+
# This file is part of s-dftd3.
2+
# SPDX-Identifier: LGPL-3.0-or-later
3+
#
4+
# s-dftd3 is free software: you can redistribute it and/or modify it under
5+
# the terms of the Lesser GNU General Public License as published by
6+
# the Free Software Foundation, either version 3 of the License, or
7+
# (at your option) any later version.
8+
#
9+
# s-dftd3 is distributed in the hope that it will be useful,
10+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
11+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12+
# Lesser GNU General Public License for more details.
13+
#
14+
# You should have received a copy of the Lesser GNU General Public License
15+
# along with s-dftd3. If not, see <https://www.gnu.org/licenses/>.
16+
"""
17+
Compatibility layer for supporting DFT-D3 in `pyscf <https://pyscf.org/>`_.
18+
19+
Example
20+
-------
21+
>>> from pyscf import gto, scf
22+
>>> from pyscf import scf
23+
>>> import dftd3.pyscf as disp
24+
>>> mol = gto.Mole()
25+
>>> mol.atom = ''' O 0.00000000 0.00000000 -0.11081188
26+
... H -0.00000000 -0.84695236 0.59109389
27+
... H -0.00000000 0.89830571 0.52404783 '''
28+
>>> mol.basis = 'cc-pvdz'
29+
>>> _ = mol.build()
30+
>>> mf = disp.dftd3(scf.RHF(mol))
31+
>>> print(mf.kernel())
32+
-75.99396273778923
33+
"""
34+
35+
import numpy as np
36+
37+
try:
38+
from pyscf import lib, gto
39+
except ModuleNotFoundError:
40+
raise ModuleNotFoundError("This submodule requires pyscf installed")
41+
42+
from .interface import (
43+
DispersionModel,
44+
RationalDampingParam,
45+
ZeroDampingParam,
46+
ModifiedRationalDampingParam,
47+
ModifiedZeroDampingParam,
48+
OptimizedPowerDampingParam,
49+
)
50+
51+
_damping_param = {
52+
"d3bj": RationalDampingParam,
53+
"d3zero": ZeroDampingParam,
54+
"d3bjm": ModifiedRationalDampingParam,
55+
"d3mbj": ModifiedRationalDampingParam,
56+
"d3zerom": ModifiedZeroDampingParam,
57+
"d3mzero": ModifiedZeroDampingParam,
58+
"d3op": OptimizedPowerDampingParam,
59+
}
60+
61+
62+
def dftd3(mf):
63+
"""Apply DFT-D3 corrections to SCF or MCSCF methods"""
64+
65+
from pyscf.scf import hf
66+
from pyscf.mcscf import casci
67+
68+
if not isinstance(mf, (hf.SCF, casci.CASCI)):
69+
raise TypeError("mf must be an instance of SCF or CASCI")
70+
71+
with_dftd3 = DFTD3Dispersion(
72+
mf.mol,
73+
xc="hf"
74+
if isinstance(mf, casci.CASCI)
75+
else getattr(mf, "xc", "HF").upper().replace(" ", ""),
76+
)
77+
78+
if isinstance(mf, _DFTD3):
79+
mf.with_dftd3 = with_dftd3
80+
return mf
81+
82+
class DFTD3(_DFTD3, mf.__class__):
83+
def __init__(self, method, with_dftd3):
84+
self.__dict__.update(method.__dict__)
85+
self.with_dftd3 = with_dftd3
86+
self._keys.update(["with_dftd3"])
87+
88+
def dump_flags(self, verbose=None):
89+
mf.__class__.dump_flags(self, verbose)
90+
if self.with_dftd3:
91+
self.with_dftd3.dump_flags(verbose)
92+
return self
93+
94+
def energy_nuc(self):
95+
enuc = mf.__class__.energy_nuc(self)
96+
if self.with_dftd3:
97+
enuc += self.with_dftd3.kernel()[0]
98+
return enuc
99+
100+
def reset(self, mol=None):
101+
self.with_dftd3.reset(mol)
102+
return mf.__class__.reset(self, mol)
103+
104+
def nuc_grad_method(self):
105+
scf_grad = mf.__class__.nuc_grad_method(self)
106+
return grad(scf_grad)
107+
108+
Gradients = lib.alias(nuc_grad_method, alias_name="Gradients")
109+
110+
return DFTD3(mf, with_dftd3)
111+
112+
113+
def grad(scf_grad):
114+
"""
115+
Apply DFT-D3 corrections to SCF or MCSCF nuclear gradients methods
116+
"""
117+
from pyscf.grad import rhf as rhf_grad
118+
119+
if not isinstance(scf_grad, rhf_grad.Gradients):
120+
raise TypeError("scf_grad must be an instance of Gradients")
121+
122+
# Ensure that the zeroth order results include DFTD3 corrections
123+
if not getattr(scf_grad.base, "with_dftd3", None):
124+
scf_grad.base = dftd3(scf_grad.base)
125+
126+
class DFTD3Grad(_DFTD3Grad, scf_grad.__class__):
127+
def grad_nuc(self, mol=None, atmlst=None):
128+
nuc_g = scf_grad.__class__.grad_nuc(self, mol, atmlst)
129+
with_dftd3 = getattr(self.base, "with_dftd3", None)
130+
if with_dftd3:
131+
disp_g = with_dftd3.kernel()[1]
132+
if atmlst is not None:
133+
disp_g = disp_g[atmlst]
134+
nuc_g += disp_g
135+
return nuc_g
136+
137+
mfgrad = DFTD3Grad.__new__(DFTD3Grad)
138+
mfgrad.__dict__.update(scf_grad.__dict__)
139+
return mfgrad
140+
141+
142+
class DFTD3Dispersion(lib.StreamObject):
143+
def __init__(self, mol, xc="hf", version="d3bj", atm=False):
144+
self.mol = mol
145+
self.verbose = mol.verbose
146+
self.xc = xc
147+
self.atm = atm
148+
self.version = version
149+
self.edisp = None
150+
self.grads = None
151+
152+
def dump_flags(self, verbose=None):
153+
lib.logger.info(self, "** DFTD3 parameter **")
154+
lib.logger.info(self, "func %s", self.xc)
155+
lib.logger.info(
156+
self, "version %s", self.version + "-atm" if self.atm else self.version
157+
)
158+
return self
159+
160+
def kernel(self):
161+
mol = self.mol
162+
163+
disp = DispersionModel(
164+
np.array([gto.charge(mol.atom_symbol(ia)) for ia in range(mol.natm)]),
165+
mol.atom_coords(),
166+
)
167+
168+
param = _damping_param[self.version](
169+
method=self.xc,
170+
atm=self.atm,
171+
)
172+
173+
res = disp.get_dispersion(param=param, grad=True)
174+
175+
self.edisp = res.get("energy")
176+
self.grads = res.get("gradient")
177+
return self.edisp, self.grads
178+
179+
def reset(self, mol):
180+
"""Reset mol and clean up relevant attributes for scanner mode"""
181+
self.mol = mol
182+
return self
183+
184+
185+
class _DFTD3:
186+
"""
187+
Stub class used to identify instances of the `DFTD3` class
188+
"""
189+
190+
pass
191+
192+
193+
class _DFTD3Grad:
194+
"""
195+
Stub class used to identify instances of the `DFTD3Grad` class
196+
"""
197+
198+
pass

python/dftd3/test_ase.py

Lines changed: 12 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -14,13 +14,20 @@
1414
# You should have received a copy of the Lesser GNU General Public License
1515
# along with s-dftd3. If not, see <https://www.gnu.org/licenses/>.
1616

17-
from dftd3.ase import DFTD3
18-
from ase.build import molecule
19-
from ase.calculators.emt import EMT
20-
from pytest import approx, raises
2117
import numpy as np
18+
import pytest
19+
from pytest import approx, raises
20+
21+
try:
22+
import ase
23+
from dftd3.ase import DFTD3
24+
from ase.build import molecule
25+
from ase.calculators.emt import EMT
26+
except ModuleNotFoundError:
27+
ase = None
2228

2329

30+
@pytest.mark.skipif(ase is None, reason="requires ase")
2431
def test_ase_scand4():
2532
thr = 1.0e-6
2633

@@ -51,6 +58,7 @@ def test_ase_scand4():
5158
assert energies == approx([-0.03880921894019244, 3.684105315180033], abs=thr)
5259

5360

61+
@pytest.mark.skipif(ase is None, reason="requires ase")
5462
def test_ase_tpssd4():
5563
thr = 1.0e-6
5664

python/dftd3/test_pyscf.py

Lines changed: 75 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,75 @@
1+
# This file is part of s-dftd3.
2+
# SPDX-Identifier: LGPL-3.0-or-later
3+
#
4+
# s-dftd3 is free software: you can redistribute it and/or modify it under
5+
# the terms of the Lesser GNU General Public License as published by
6+
# the Free Software Foundation, either version 3 of the License, or
7+
# (at your option) any later version.
8+
#
9+
# s-dftd3 is distributed in the hope that it will be useful,
10+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
11+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12+
# Lesser GNU General Public License for more details.
13+
#
14+
# You should have received a copy of the Lesser GNU General Public License
15+
# along with s-dftd3. If not, see <https://www.gnu.org/licenses/>.
16+
17+
import numpy as np
18+
import pytest
19+
from pytest import approx, raises
20+
21+
try:
22+
import pyscf
23+
from pyscf import lib, gto, scf
24+
import dftd3.pyscf as disp
25+
except ModuleNotFoundError:
26+
pyscf = None
27+
28+
29+
@pytest.mark.skipif(pyscf is None, reason="requires pyscf")
30+
def test_dftd3():
31+
32+
_ch4 = gto.M(
33+
atom="""
34+
H 0.000000000000 0.000000000000 0.000000000000
35+
C 0.000000000000 0.000000000000 1.087900000000
36+
H 1.025681956337 0.000000000000 1.450533333333
37+
H -0.512840978169 0.888266630391 1.450533333333
38+
H -0.512840978169 -0.888266630391 1.450533333333
39+
"""
40+
)
41+
d3 = disp.DFTD3Dispersion(_ch4)
42+
d3.xc = "B3LYP"
43+
assert d3.kernel()[0] == approx(-0.0019136221730972761, abs=1.0e-9)
44+
45+
46+
@pytest.mark.skipif(pyscf is None, reason="requires pyscf")
47+
def test_dftd3_scf():
48+
_h2o = gto.M(
49+
atom="""
50+
O 0.000 0.000 0.000
51+
H 0.000 -0.757 0.587
52+
H 0.000 0.757 0.587
53+
""",
54+
symmetry=True,
55+
)
56+
57+
mf = disp.dftd3(scf.RHF(_h2o))
58+
assert mf.kernel() == approx(-74.96757204541478, abs=1.0e-8)
59+
60+
61+
@pytest.mark.skipif(pyscf is None, reason="requires pyscf")
62+
def test_dftd3_scf_grad():
63+
_h2o = gto.M(
64+
atom="O 0. 0. 0.000; H 0. -0.757 0.587; H 0. 0.757 0.587",
65+
symmetry=True,
66+
)
67+
68+
mf = disp.dftd3(scf.RHF(_h2o)).run()
69+
mfs = mf.as_scanner()
70+
e1 = mfs("O 0. 0. 0.0001; H 0. -0.757 0.587; H 0. 0.757 0.587")
71+
e2 = mfs("O 0. 0. -0.0001; H 0. -0.757 0.587; H 0. 0.757 0.587")
72+
ref = (e1 - e2) / 0.0002 * lib.param.BOHR
73+
g = mf.nuc_grad_method().kernel()
74+
75+
assert ref == approx(g[0, 2], abs=1.0e-6)

0 commit comments

Comments
 (0)