Skip to content

Commit

Permalink
mod: implemented get_all_energies for Gaussian16 and added a test
Browse files Browse the repository at this point in the history
now all OverlapCalculators return the GS energy when root is None
  • Loading branch information
eljost committed Dec 19, 2023
1 parent c376419 commit 7261cf6
Show file tree
Hide file tree
Showing 7 changed files with 95 additions and 44 deletions.
3 changes: 3 additions & 0 deletions pysisyphus/Geometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
57 changes: 42 additions & 15 deletions pysisyphus/calculators/Gaussian16.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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,
Expand All @@ -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}"
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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",
Expand Down
11 changes: 10 additions & 1 deletion pysisyphus/calculators/ORCA.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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}

Expand Down
23 changes: 10 additions & 13 deletions pysisyphus/calculators/PySCF.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@ class PySCF(OverlapCalculator):
multisteps = {
"scf": ("scf",),
"tdhf": ("scf", "tddft"),
"tdahf": ("scf", "tda"),
"dft": ("dft",),
"mp2": ("scf", "mp2"),
"tddft": ("dft", "tddft"),
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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,
Expand All @@ -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")
Expand Down
9 changes: 5 additions & 4 deletions pysisyphus/calculators/Turbomole.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"]
Expand Down
15 changes: 15 additions & 0 deletions tests/test_es_capabilities/test_es_capabilities.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand All @@ -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])
21 changes: 10 additions & 11 deletions tests/test_overlap_calculator/test_calculators.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
}


Expand All @@ -35,7 +36,7 @@
PySCF,
{
"basis": "sto3g",
"method": "tddft",
"method": "tdhf",
"nroots": 3,
},
marks=(using("pyscf")),
Expand Down Expand Up @@ -63,9 +64,7 @@
),
pytest.param(
Gaussian16,
{
"route": "hf sto-3g td=(nstates=3)"
},
{"route": "hf sto-3g td=(nstates=3)"},
marks=(using("gaussian16")),
),
),
Expand Down

0 comments on commit 7261cf6

Please sign in to comment.