Source code for pysisyphus.calculators.EnergyMin

# [1]
#     Extending NEB method to reaction pathways involving multiple spin states
#     Zhao et al, 2020
# [2]
#     Optimizing Conical Intersections without Derivative Coupling Vectors:
#     Application to Multistate Multireference Second-Order Perturbation Theory
#     Levine, Coe, Martinez 2008

import numpy as np

from pysisyphus.calculators.Calculator import Calculator
from pysisyphus.constants import AU2KJPERMOL
from pysisyphus.helpers import norm_max_rms

[docs] class EnergyMin(Calculator):
[docs] def __init__( self, calculator1: Calculator, calculator2: Calculator, mix: bool = False, alpha: float = 0.02, # Hartree sigma: float = 3.5, # Unitless; default for ethene case in [1] min_energy_diff: float = 0.0, check_after: int = 0, **kwargs, ): """ Use energy and derivatives of the calculator with lower energy. This calculators carries out two calculations with different settings and returns the results of the lower energy one. This can be used to consider flips between a singlet and a triplet PES etc. Parameters ---------- calculator1 Wrapped QC calculator that provides energies and its derivatives. calculator2 Wrapped QC calculator that provides energies and its derivatives. mix Enable mixing of both forces, according to the approach outlined in [2]. Can be used to optimize guesses for MECPs. Pass alpha Smoothing parameter in Hartree. See [2] for a discussion. sigma Unitless gap size parameter. The final gap becomes smaller for bigga sigmas. Has to be adapted for each case. See [2] for a discussion (p. 407 right column and p. 408 left column.) min_energy_diff Energy difference in Hartree. When set to a value != 0 and the energy difference between both calculators drops below this value, execution of both calculations is diabled for 'check_after' cycles. In these cycles the calculator choice remains fixed. After 'check_after' cycles, both energies will be calculated and it is checked, if the previous calculator choice remains valid. In conjunction with 'check_after' both arguments can be used to save computational ressources. check_after Amount of cycles in which the calculator choice remains fixed. Other Parameters ---------------- **kwargs Keyword arguments passed to the Calculator baseclass. """ super().__init__(**kwargs) self.calc1 = calculator1 self.calc2 = calculator2 self.alpha = alpha self.sigma = sigma self.min_energy_diff = float(min_energy_diff) self.check_after = int(check_after) assert self.check_after >= 0, "'check_after' must not be negative!" self.mix = mix self.recalc_in = self.check_after self.fixed_calc = None
[docs] def do_calculations(self, name, atoms, coords, **prepare_kwargs): def run_calculation(calc, name=name): results = getattr(calc, name)(atoms, coords, **prepare_kwargs) return results if (self.fixed_calc is not None) and self.recalc_in > 0: self.log( f"Used fixed calculator '{self.fixed_calc}'. Re-checking both " f"calculators in {self.recalc_in} cycles." ) results = run_calculation(self.fixed_calc) self.recalc_in -= 1 self.calc_counter += 1 return results elif (self.fixed_calc is not None) and self.recalc_in == 0: self.log(f"Unset fixed calculator {self.calc1}.") self.fixed_calc = None # Avoid unnecessary costly Hessian calculations; decide which Hessian # to calculate from previous energy calculation. if name == "get_hessian": tmp_energy1 = run_calculation(self.calc1, "get_energy")["energy"] tmp_energy2 = run_calculation(self.calc2, "get_energy")["energy"] if tmp_energy1 <= tmp_energy2: results1 = run_calculation(self.calc1) results2 = {"energy": tmp_energy2} self.log("Calculated Hessian for calc1, skipped it for calc2.") else: results1 = {"energy": tmp_energy1} results2 = run_calculation(self.calc2) self.log("Calculated Hessian for calc2, skipped it for calc1.") # Do both full calculation otherwise else: results1 = run_calculation(self.calc1, name) results2 = run_calculation(self.calc2, name) energy1 = results1["energy"] energy2 = results2["energy"] all_energies = np.array((energy1, energy2)) # Mixed forces to optimize crossing points if self.mix: # Must be positive, so substract lower energy from higher energy. # I is the higher state, j the lower one. en_i, en_j = ( (energy2, energy1) if (energy1 < energy2) else (energy1, energy2) ) energy_diff = en_i - en_j energy_diff_kJ = energy_diff * AU2KJPERMOL assert energy_diff > 0.0 self.log( f"Mix mode, ΔE={energy_diff:.6f} au ({energy_diff_kJ:.2f} kJ mol⁻¹)" ) energy_diff_sq = energy_diff**2 denom = energy_diff + self.alpha energy = (en_i + en_j) / 2 + self.sigma * energy_diff_sq / denom results = { "energy": energy, "all_energies": all_energies, } self.log(f"Mixed energy: {energy:.6f} au") if name == "get_forces": forces1 = results1["forces"] forces2 = results2["forces"] forces_i, forces_j = ( (forces2, forces1) if (energy1 < energy2) else (forces1, forces2) ) # There seems to be a typo in Eq. (8) in [1]; the final term # should be (dV_J - dV_I) instead of the (dV_I - dV_J). # The formula is correct though, in the original publication # (Eq. (7) in [2]). forces = (forces_i + forces_j) / 2 + self.sigma * ( energy_diff_sq + 2 * self.alpha * energy_diff ) / denom**2 * (forces_i - forces_j) results["forces"] = forces for key, forces in ( ("mixed forces", forces), ("forces1", forces1), ("forces2", forces2), ): norm, max_, rms_ = norm_max_rms(forces) self.log( f"{key: >14s}: (norm={norm:.4f}, max(|{key: >14s}|)={max_:.4f}, " f"rms={rms_:.4f}) au a0⁻¹" ) self.calc_counter += 1 return results # Mixed forces end min_ind = [1, 0][int(energy1 < energy2)] en1_or_en2 = ("calc1", "calc2")[min_ind] energy_diff = energy1 - energy2 # Try to fix calculator, if requested if (self.min_energy_diff and self.check_after) and ( # When the actual difference is above to minimum differences # or # no calculator is fixed yet (energy_diff > self.min_energy_diff) or (self.fixed_calc is None) ): self.fixed_calc = (self.calc1, self.calc2)[min_ind] self.recalc_in = self.check_after self.log( f"Fixed calculator choice ({en1_or_en2}) for {self.recalc_in} cycles." ) results = (results1, results2)[min_ind] results["all_energies"] = all_energies energy_diff_kJ = abs(energy_diff) * AU2KJPERMOL self.log( f"energy_calc1={energy1:.6f} au, energy_calc2={energy2:.6f} au, returning " f"results for {en1_or_en2}, {energy_diff_kJ: >10.2f} kJ mol⁻¹ lower." ) self.calc_counter += 1 return results
[docs] def get_energy(self, atoms, coords, **prepare_kwargs): return self.do_calculations("get_energy", atoms, coords, **prepare_kwargs)
[docs] def get_forces(self, atoms, coords, **prepare_kwargs): return self.do_calculations("get_forces", atoms, coords, **prepare_kwargs)
[docs] def get_hessian(self, atoms, coords, **prepare_kwargs): return self.do_calculations("get_hessian", atoms, coords, **prepare_kwargs)
[docs] def get_chkfiles(self) -> dict: chkfiles = {} for key in ("calc1", "calc2"): try: chkfiles[key] = getattr(self, key).get_chkfiles() except AttributeError: pass return chkfiles
[docs] def set_chkfiles(self, chkfiles: dict): for key in ("calc1", "calc2"): try: getattr(self, key).set_chkfiles(chkfiles[key]) msg = f"Set chkfile on '{key}'" except KeyError: msg = f"Found no information for '{key}' in chkfiles!" except AttributeError: msg = f"Setting chkfiles is not supported by '{key}'" self.log(msg)