add basis gradient tests

lgtm and coveralls (hopefully)

more fdiff

fix

adcc/backends/psi4.py

@@ -85,7 +85,6 @@ def correlated_gradient(self, g1_ao, w_ao, g2_ao_1, g2_ao_2):
             string = "TEI" + str(atom)
             deriv2 = self.mints.ao_tei_deriv1(atom)
             for p in range(3):
-                map_key = string + cart[p]
                 deriv2_np = np.asarray(deriv2[p])
                 Gradient["TEI"][atom, p] += np.einsum(
                     'pqrs,prqs->', g2_ao_1, deriv2_np, optimize=True

adcc/gradients/orbital_response.py

@@ -123,12 +123,12 @@ def __matmul__(self, l_ov):
             - einsum("ibja,jb->ia", self.hf.ovov, l_ov)
         )
         # TODO: generalize once other solvent methods are available
-        if "pe_induction_elec" in self.hf.operators.density_dependent_operators:
+        if self.hf.environment == "pe":
             # PE contribution to the orbital Hessian
             ops = self.hf.operators
             dm = OneParticleOperator(self.hf, is_symmetric=True)
             dm.ov = l_ov
-            ret += ops.density_dependent_operators["pe_induction_elec"](dm).ov
+            ret += ops.pe_induction_elec(dm).ov
         return evaluate(ret)
 
 

adcc/gradients/test_functionality_gradients.py

@@ -0,0 +1,75 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2020 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+import unittest
+import itertools
+import adcc
+import adcc.backends
+
+from numpy.testing import assert_allclose
+
+import pytest
+
+from ..misc import expand_test_templates
+from adcc.backends.testing import cached_backend_hf
+from adcc.testdata.cache import gradient_data
+
+
+backends = [b for b in adcc.backends.available()
+            if b not in ["molsturm", "veloxchem"]]
+molecules = ["h2o"]
+basissets = ["sto3g", "ccpvdz"]
+methods = ["mp2", "adc1", "adc2"]
+combinations = list(itertools.product(molecules, basissets, methods, backends))
+
+
+@pytest.mark.skipif(len(backends) == 0, reason="No backend found.")
+@expand_test_templates(combinations)
+class TestNuclearGradients(unittest.TestCase):
+    def template_nuclear_gradient(self, molecule, basis, method, backend):
+        grad_ref = gradient_data[molecule][basis][method]
+
+        energy_ref = grad_ref["energy"]
+        grad_fdiff = grad_ref["gradient"]
+
+        scfres = cached_backend_hf(backend, molecule, basis, conv_tol=1e-13)
+        if "adc" in method:
+            # TODO: convergence needs to be very very tight...
+            # so we want to make sure all vectors are tightly converged
+            n_limit = 5
+            state = adcc.run_adc(scfres, method=method,
+                                 n_singlets=10, conv_tol=1e-11)
+            for ee in state.excitations[:n_limit]:
+                grad = adcc.nuclear_gradient(ee)
+                assert_allclose(energy_ref[ee.index], ee.total_energy, atol=1e-10)
+                assert_allclose(
+                    grad_fdiff[ee.index], grad["Total"], atol=1e-7
+                )
+        else:
+            # MP2 gradients
+            refstate = adcc.ReferenceState(scfres)
+            mp = adcc.LazyMp(refstate)
+            grad = adcc.nuclear_gradient(mp)
+            assert_allclose(energy_ref, mp.energy(2), atol=1e-8)
+            assert_allclose(
+                grad_fdiff, grad["Total"], atol=1e-8
+            )

adcc/testdata/cache.py

@@ -258,3 +258,4 @@ def read_yaml_data(fname):
 
 qchem_data = read_yaml_data("qchem_dump.yml")
 tmole_data = read_yaml_data("tmole_dump.yml")
+gradient_data = read_yaml_data("grad_dump.yml")

adcc/testdata/dump_fdiff_gradient.py

@@ -0,0 +1,139 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2021 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+import itertools
+import adcc
+import numpy as np
+import yaml
+from tqdm import tqdm
+
+from pyscf import gto
+
+from static_data import xyz
+
+
+prefactors_5p = np.array([1.0, -8.0, 8.0, -1.0]) / 12.0
+multipliers_5p = [-2, -1, 1, 2]
+coords_label = ["x", "y", "z"]
+
+
+def _molstring(elems, coords):
+    s = ""
+    for kk, (i, c) in enumerate(zip(elems, coords)):
+        s += f"{i} {c[0]} {c[1]} {c[2]}"
+        if kk != len(elems) - 1:
+            s += "\n"
+    return s
+
+
+def adc_energy(scfres, method, **kwargs):
+    state = adcc.run_adc(method=method, data_or_matrix=scfres,
+                         output=None, **kwargs)
+    return state.total_energy
+
+
+def mp_energy(scfres, method, **kwargs):
+    level = {
+        "mp2": 2,
+        "mp3": 3,
+    }
+    refstate = adcc.ReferenceState(scfres)
+    return adcc.LazyMp(refstate).energy(level[method])
+
+
+def fdiff_gradient(molstring, method, basis, step=1e-4, **kwargs):
+    m = gto.M(atom=molstring, unit='Bohr', basis=basis)
+    coords = m.atom_coords().copy()
+    elements = m.elements.copy()
+
+    n_grads = kwargs.get("n_singlets", 1)
+    conv_tol = kwargs.get("conv_tol", 1e-10) / 10
+
+    # run unperturbed system
+    scfres = adcc.backends.run_hf(
+        'pyscf', molstring, basis, conv_tol=conv_tol, conv_tol_grad=conv_tol
+    )
+    if "adc" in method:
+        en = adc_energy(scfres, method, **kwargs)
+    else:
+        en = mp_energy(scfres, method, **kwargs)
+
+    natoms = len(elements)
+    grad = np.zeros((n_grads, natoms, 3))
+    at_c = list(itertools.product(range(natoms), range(3)))
+    for i, c in tqdm(at_c):
+        for f, p in zip(multipliers_5p, prefactors_5p):
+            coords_p = coords.copy()
+            coords_p[i, c] += f * step
+            geom_p = _molstring(elements, coords_p)
+            scfres = adcc.backends.run_hf(
+                'pyscf', geom_p, basis, conv_tol=conv_tol, conv_tol_grad=conv_tol
+            )
+            if "adc" in method:
+                en_pert = adc_energy(scfres, method, **kwargs)
+            else:
+                en_pert = mp_energy(scfres, method, **kwargs)
+            grad[:, i, c] += p * en_pert / step
+    return en, grad
+
+
+def main():
+    config_excited = {
+        "n_singlets": 5,
+    }
+    basissets = [
+        "sto3g",
+        "ccpvdz",
+    ]
+    methods = [
+        "mp2",
+        "adc1",
+        "adc2",
+    ]
+    molecules = ["h2o", "hf", "formaldehyde"]
+    ret = {}
+    for molecule in molecules:
+        ret[molecule] = {}
+        for basis in basissets:
+            ret[molecule][basis] = {}
+            for method in methods:
+                kwargs = {
+                    "conv_tol": 1e-8,
+                }
+                if "adc" in method:
+                    kwargs.update(config_excited)
+                basename = f"{molecule}_{basis}_{method}"
+                print(f"Evaluating finite difference gradient for {basename}.")
+                en, grad = fdiff_gradient(xyz[molecule], method, basis, **kwargs)
+                if isinstance(en, np.ndarray):
+                    en = en.tolist()
+                cont = {
+                    "energy": en,
+                    "gradient": np.squeeze(grad).tolist(),
+                }
+                ret[molecule][basis][method] = cont
+    with open("grad_dump.yml", "w") as yamlout:
+        yaml.safe_dump(ret, yamlout)
+
+
+if __name__ == "__main__":
+    main()