From 7261cf634c38573372bc2a9173b653719b534fde Mon Sep 17 00:00:00 2001 From: Johannes Steinmetzer Date: Tue, 19 Dec 2023 14:26:09 +0100 Subject: [PATCH] mod: implemented get_all_energies for Gaussian16 and added a test now all OverlapCalculators return the GS energy when root is None --- pysisyphus/Geometry.py | 3 + pysisyphus/calculators/Gaussian16.py | 57 ++++++++++++++----- pysisyphus/calculators/ORCA.py | 11 +++- pysisyphus/calculators/PySCF.py | 23 ++++---- pysisyphus/calculators/Turbomole.py | 9 +-- .../test_es_capabilities.py | 15 +++++ .../test_calculators.py | 21 ++++--- 7 files changed, 95 insertions(+), 44 deletions(-) diff --git a/pysisyphus/Geometry.py b/pysisyphus/Geometry.py index 83c492c2b..7f0807f73 100644 --- a/pysisyphus/Geometry.py +++ b/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 diff --git a/pysisyphus/calculators/Gaussian16.py b/pysisyphus/calculators/Gaussian16.py index 8bb739d1c..472be0fee 100644 --- a/pysisyphus/calculators/Gaussian16.py +++ b/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,13 +593,12 @@ 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 @@ -588,19 +606,28 @@ def parse_force(self, path): 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", diff --git a/pysisyphus/calculators/ORCA.py b/pysisyphus/calculators/ORCA.py index 2a5320ddb..a9b71c563 100644 --- a/pysisyphus/calculators/ORCA.py +++ b/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} diff --git a/pysisyphus/calculators/PySCF.py b/pysisyphus/calculators/PySCF.py index 00e7c2b8b..2ddad117c 100644 --- a/pysisyphus/calculators/PySCF.py +++ b/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,9 @@ 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.method = "dft" - if self.method == "tddft": - assert self.nroots, "nroots must be set with method='tddft'!" + self.do_es = self.method in ("tda", "tddft", "tdhf", "tdahf") + 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 +162,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 +188,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") diff --git a/pysisyphus/calculators/Turbomole.py b/pysisyphus/calculators/Turbomole.py index 86abfb4ec..de0fe0487 100644 --- a/pysisyphus/calculators/Turbomole.py +++ b/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"] diff --git a/tests/test_es_capabilities/test_es_capabilities.py b/tests/test_es_capabilities/test_es_capabilities.py index 84519cde8..ebc946021 100644 --- a/tests/test_es_capabilities/test_es_capabilities.py +++ b/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]) diff --git a/tests/test_overlap_calculator/test_calculators.py b/tests/test_overlap_calculator/test_calculators.py index b8bac1f05..085b46f5d 100644 --- a/tests/test_overlap_calculator/test_calculators.py +++ b/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")), ), ),