mod: implemented get_all_energies for Gaussian16 and added a test

now all OverlapCalculators return the GS energy when root is None

Author: Johannes Steinmetzer

pysisyphus/Geometry.py

@@ -873,6 +873,9 @@ def all_energies(self):
             self.set_results(results)
         return self._all_energies
 
+    def get_root_energy(self, root):
+        return self.all_energies[root]
+
     def has_all_energies(self):
         return self._all_energies is not None
 

pysisyphus/calculators/Gaussian16.py

@@ -55,6 +55,7 @@ def __init__(
         }
         exc_keyword = [key for key in "td tda cis".split() if key in keywords]
         self.nstates = self.nroots
+        self.exc_args = None
         if exc_keyword:
             self.exc_key = exc_keyword[0]
             exc_dict = keywords[self.exc_key]
@@ -101,6 +102,7 @@ def __init__(
         self.gaussian_input = textwrap.dedent(self.gaussian_input)
 
         self.parser_funcs = {
+            "energy": self.parse_energy,
             "force": self.parse_force,
             "hessian": self.parse_hessian,
             "stable": self.parse_stable,
@@ -118,8 +120,7 @@ def __init__(
         )
 
     def make_exc_str(self):
-        # Ground state calculation
-        if not self.track and (self.root is None):
+        if self.exc_args is None:
             return ""
         root = self.root if self.root is not None else 1
         root_str = f"root={root}"
@@ -341,7 +342,25 @@ def store_and_track(self, results, func, atoms, coords, **prepare_kwargs):
         return results
 
     def get_energy(self, atoms, coords, **prepare_kwargs):
-        return self.get_forces(atoms, coords, **prepare_kwargs)
+        did_stable = False
+        if self.stable:
+            is_stable = self.run_stable(atoms, coords)
+            self.log(f"Wavefunction is now stable: {is_stable}")
+            did_stable = True
+        inp = self.prepare_input(
+            atoms, coords, "sp", did_stable=did_stable, **prepare_kwargs
+        )
+        kwargs = {
+            "calc": "energy",
+        }
+        results = self.run(inp, **kwargs)
+        results = self.store_and_track(
+            results, self.get_energy, atoms, coords, **prepare_kwargs
+        )
+        return results
+
+    def get_all_energies(self, atoms, coords, **prepare_kwargs):
+        return self.get_energy(atoms, coords, **prepare_kwargs)
 
     def get_forces(self, atoms, coords, **prepare_kwargs):
         did_stable = False
@@ -574,33 +593,41 @@ def prepare_overlap_data(self, path):
         all_energies[1:] += exc_energies
         return C, X, Y, all_energies
 
-    def parse_force(self, path):
+    def parse_energy(self, path):
         results = {}
-        keys = ("SCF Energy", "Total Energy", "Cartesian Gradient")
+        keys = ("SCF Energy", "Total Energy")
         fchk_path = Path(path) / f"{self.fn_base}.fchk"
         fchk_dict = self.parse_fchk(fchk_path, keys)
         results["energy"] = fchk_dict["SCF Energy"]
-        results["forces"] = -fchk_dict["Cartesian Gradient"]
 
         if self.nstates:
             # This sets the proper excited state energy in the
             # results dict and also stores all energies.
             exc_energies = self.parse_tddft(path)
             # G16 root input is 1 based, so we substract 1 to get
             # the right index here.
-            root = self.root if self.root is not None else 1
-            root_exc_en = exc_energies[root - 1]
             gs_energy = fchk_dict["SCF Energy"]
-            # Add excitation energy to ground state energy.
-            results["energy"] += root_exc_en
-            # Create a new array including the ground state energy
-            # to save the energies of all states.
-            all_ens = np.full(len(exc_energies) + 1, gs_energy)
-            all_ens[1:] += exc_energies
-            results["all_energies"] = all_ens
+            # Modify "energy" when a root is selected
+            if self.root is not None:
+                root_exc_en = exc_energies[self.root - 1]
+                # Add excitation energy to ground state energy.
+                results["energy"] += root_exc_en
+                # Create a new array including the ground state energy
+                # to save the energies of all states.
+                all_ens = np.full(len(exc_energies) + 1, gs_energy)
+                all_ens[1:] += exc_energies
+                results["all_energies"] = all_ens
 
         return results
 
+    def parse_force(self, path):
+        results = self.parse_energy(path)
+        keys = ("Cartesian Gradient",)
+        fchk_path = Path(path) / f"{self.fn_base}.fchk"
+        fchk_dict = self.parse_fchk(fchk_path, keys)
+        results["forces"] = -fchk_dict["Cartesian Gradient"]
+        return results
+
     def parse_hessian(self, path):
         keys = (
             "Total Energy",

pysisyphus/calculators/ORCA.py

@@ -616,6 +616,7 @@ def __init__(
 
         self.do_tddft = bool(es_block_header_set & td_blocks)
         self.do_ice = bool(es_block_header_set & ice_blocks)
+        self.do_es = any((self.do_tddft, self.do_ice))
         # There can be at most on ES block at a time
         assert not (self.do_tddft and self.do_ice)
         if self.es_block_header:
@@ -948,7 +949,15 @@ def parse_energy(self, path):
         log_fn = log_fn[0]
         with open(log_fn) as handle:
             text = handle.read()
-        mobj = re.search(r"FINAL SINGLE POINT ENERGY\s+([\d\-\.]+)", text)
+        # By default reports the total energy of the first ES, when ES were calculated.
+        # But we are interested in the GS energy, when self.root is None ...
+        if not self.do_es or self.root is not None:
+            en_re = r"FINAL SINGLE POINT ENERGY\s+([\d\-\.]+)"
+        # ... so we look at the energy that was reported after the SCF.
+        else:
+            en_re = "Total Energy\s+:\s+([\d\-\.]+) Eh"
+        en_re = re.compile(en_re)
+        mobj = en_re.search(text)
         energy = float(mobj[1])
         return {"energy": energy}
 

pysisyphus/calculators/PySCF.py

@@ -23,6 +23,7 @@ class PySCF(OverlapCalculator):
     multisteps = {
         "scf": ("scf",),
         "tdhf": ("scf", "tddft"),
+        "tdahf": ("scf", "tda"),
         "dft": ("dft",),
         "mp2": ("scf", "mp2"),
         "tddft": ("dft", "tddft"),
@@ -53,12 +54,13 @@ def __init__(
         self.basis = basis
         self.xc = xc
         self.method = method.lower()
-        if self.method in ("tda", "tddft") and self.xc is None:
-            self.multisteps[self.method] = ("scf", self.method)
-        if self.xc and self.method != "tddft":
+        self.do_es = self.method in ("tda", "tddft", "tdhf", "tdahf")
+        # Set method to DFT when an XC-functional was selected, but DFT wasn't explicitly
+        # enabled. When an ES method was chosen it is kept.
+        if self.xc and self.method not in ("tddft", "tda"):
             self.method = "dft"
-        if self.method == "tddft":
-            assert self.nroots, "nroots must be set with method='tddft'!"
+        if self.do_es:
+            assert self.nroots, f"nroots must be set with method='{self.method}'!"
         self.auxbasis = auxbasis
         self.keep_chk = keep_chk
         self.verbose = int(verbose)
@@ -164,12 +166,11 @@ def get_energy(self, atoms, coords, **prepare_kwargs):
 
         mol = self.prepare_input(atoms, coords)
         mf = self.run(mol, point_charges=point_charges)
-        energy = mf.e_tot
-        root = 0 if self.root is None else self.root
-        try:
-            energy = energy[root]
-        except (IndexError, TypeError):
-            pass
+        all_energies = self.parse_all_energies()
+        if self.root is not None:
+            energy = all_energies[self.root]
+        else:
+            energy = all_energies[0]
 
         results = {
             "energy": energy,
@@ -191,7 +192,7 @@ def get_forces(self, atoms, coords, **prepare_kwargs):
         mol = self.prepare_input(atoms, coords)
         mf = self.run(mol, point_charges=point_charges)
         grad_driver = mf.Gradients()
-        if self.root:
+        if self.root is not None:
             grad_driver.state = self.root
         gradient = grad_driver.kernel()
         self.log("Completed gradient step")

pysisyphus/calculators/Turbomole.py

@@ -698,10 +698,11 @@ def parse_energy(self, path):
         en_regex = re.compile(r"Total energy\s*:?\s*=?\s*([\d\-\.]+)", re.IGNORECASE)
         tot_ens = en_regex.findall(text)
 
-        if self.td:
-            # Drop ground state energy that is repeated
-            root = self.root if self.root is not None else 1
-            tot_en = tot_ens[1:][root]
+        # Only modify energy when self.root is set; otherwise stick with the GS energy.
+        if self.td and self.root is not None:
+            # Drop ground state energy that is repeated. That is why we don't subtract
+            # 1 from self.root.
+            tot_en = tot_ens[1:][self.root]
         elif self.ricc2 and self.ricc2_opt:
             results = parse_turbo_gradient(path)
             tot_en = results["energy"]

tests/test_es_capabilities/test_es_capabilities.py

@@ -57,6 +57,13 @@
             },
             marks=using("PySCF"),
         ),
+        pytest.param(
+            Gaussian16,
+            {
+                "route": "hf def2svp td(nstates=3)",
+            },
+            marks=using("gaussian16"),
+        ),
     ),
 )
 def test_h2o_all_energies(mult, ref_energies, calc_cls, calc_kwargs, this_dir):
@@ -72,3 +79,11 @@ def test_h2o_all_energies(mult, ref_energies, calc_cls, calc_kwargs, this_dir):
     # PySCF and Turbomole agree extermely well, at least in the restricted calculation.
     # ORCA deviates up to 5e-5 Eh.
     np.testing.assert_allclose(all_energies, ref_energies, atol=5e-5)
+
+    # As we did not set any root the geometries energy should correspond to the GS energy
+    energy = geom.energy
+    assert energy == pytest.approx(ref_energies[0])
+
+    for root in range(4):
+        root_en = geom.get_root_energy(root)
+        assert root_en == pytest.approx(ref_energies[root])

tests/test_overlap_calculator/test_calculators.py

@@ -6,17 +6,18 @@
 
 
 _ROOT_REF_ENERGIES = {
+    # -74.9607 au is the GS
     # -74.4893 au is the S1, not the GS
-    (ORCA, None): -74.4893,
+    (ORCA, None): -74.9607,
     (ORCA, 2): -74.3984,
-    (PySCF, None): -74.4893,
-    (PySCF, 2): -74.3714,
-    (DFTBp, None): -4.077751,
-    (DFTBp, 2): -3.33313,
-    (Turbomole, None): -74.48927,
+    (PySCF, None): -74.9607,
+    (PySCF, 2): -74.3984,  # ???
+    (Turbomole, None): -74.9607,
     (Turbomole, 2): -74.39837,
-    (Gaussian16, None): -74.48927,
+    (Gaussian16, None): -74.9607,
     (Gaussian16, 2): -74.39837,
+    (DFTBp, None): -4.077751,
+    (DFTBp, 2): -3.33313,
 }
 
 
@@ -35,7 +36,7 @@
             PySCF,
             {
                 "basis": "sto3g",
-                "method": "tddft",
+                "method": "tdhf",
                 "nroots": 3,
             },
             marks=(using("pyscf")),
@@ -63,9 +64,7 @@
         ),
         pytest.param(
             Gaussian16,
-            {
-                "route": "hf sto-3g td=(nstates=3)"
-            },
+            {"route": "hf sto-3g td=(nstates=3)"},
             marks=(using("gaussian16")),
         ),
     ),