Merge pull request eljost#208 from eljost/dev

Dev

pysisyphus/Geometry.py

@@ -312,7 +312,7 @@ def __eq__(self, other):
 
     def __sub__(self, other):
         self.assert_compatibility(other)
-        if self.coord_type == "cart":
+        if self.coord_type in ("cart", "cartesian"):
             diff = self.coords - other.coords
         elif self.coord_type in ("redund", "dlc"):
             # Take periodicity of dihedrals into account by calling

pysisyphus/calculators/EnergyMin.py

@@ -1,43 +1,139 @@
 # [1] https://doi.org/10.1063/5.0021923
 #     Extending NEB method to reaction pathways involving multiple spin states
 #     Zhao et al, 2020
+# [2] https://doi.org/10.1021/jp0761618
+#     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 AU2KCALPERMOL, AU2KJPERMOL
+from pysisyphus.constants import AU2KJPERMOL
+from pysisyphus.helpers import norm_max_rms
 
 
 class EnergyMin(Calculator):
-    def __init__(self, calculator1, calculator2, **kwargs):
+    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 = False
+        self.mix = mix
+        self.recalc_in = self.check_after
+        self.fixed_calc = None
 
     def do_calculations(self, name, atoms, coords, **prepare_kwargs):
-        results1 = getattr(self.calc1, name)(atoms, coords, **prepare_kwargs)
-        results2 = getattr(self.calc2, 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:
-            alpha = 0.02 / AU2KCALPERMOL
-            sigma = 3.5 / AU2KCALPERMOL
             # 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*AU2KJPERMOL:.2f} kJ mol⁻¹)"
+                f"Mix mode, ΔE={energy_diff:.6f} au ({energy_diff_kJ:.2f} kJ mol⁻¹)"
             )
-            energy_diff_sq = energy_diff ** 2
-            denom = energy_diff + alpha
-            energy = (en_i + en_j) / 2 + sigma * energy_diff_sq / denom
+            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,
@@ -49,23 +145,53 @@ def do_calculations(self, name, atoms, coords, **prepare_kwargs):
                 forces_i, forces_j = (
                     (forces2, forces1) if (energy1 < energy2) else (forces1, forces2)
                 )
-                forces = (forces_i + forces_j) / 2 + sigma * (
-                    energy_diff_sq + 2 * alpha * energy_diff
-                ) / denom ** 2 * (forces_j - forces_i)
+                # 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
-                self.log(f"norm(mixed forces)={np.linalg.norm(forces):.6f} au a0⁻¹")
+                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
-        en_diff_kJ = abs(energy1 - energy2) * AU2KJPERMOL
+        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}, {en_diff_kJ:.2f} kJ mol⁻¹ lower."
+            f"results for {en1_or_en2}, {energy_diff_kJ: >10.2f} kJ mol⁻¹ lower."
         )
+        self.calc_counter += 1
         return results
 
     def get_energy(self, atoms, coords, **prepare_kwargs):
@@ -76,3 +202,23 @@ def get_forces(self, atoms, coords, **prepare_kwargs):
 
     def get_hessian(self, atoms, coords, **prepare_kwargs):
         return self.do_calculations("get_hessian", atoms, coords, **prepare_kwargs)
+
+    def get_chkfiles(self) -> dict:
+        chkfiles = {}
+        for key in ("calc1", "calc2"):
+            try:
+                chkfiles[key] = getattr(self, key).get_chkfiles()
+            except AttributeError:
+                pass
+        return chkfiles
+
+    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)

pysisyphus/calculators/MultiCalc.py

@@ -18,14 +18,18 @@
     pass
 
 
-def calcs_from_dict(calc_dict, base_name, charge, mult, pal, mem):
+def calcs_from_dict(calc_dict, base_name, calc_number, charge, mult, pal, mem):
     keys_calcs = {}
     calc_kwargs_ = {
+        "calc_number": calc_number,
         "charge": charge,
         "mult": mult,
         "pal": pal,
         "mem": mem,
     }
+    # Try to distinguish the calculators according to their base_name. If no
+    # base_name is given we rely on the calc_number.
+    base_name = base_name if base_name != "" else f"{calc_number:03d}"
     for i, (key, kwargs) in enumerate(calc_dict.items()):
         calc_kwargs = calc_kwargs_.copy()
         calc_kwargs["base_name"] = f"{base_name}_{key}"
@@ -46,7 +50,13 @@ def __init__(self, calcs, **kwargs):
             key: calc_dict.pop("do_hess", False) for key, calc_dict in calcs.items()
         }
         self.keys_calcs = calcs_from_dict(
-            calcs, self.base_name, self.charge, self.mult, self.pal, self.mem
+            calcs,
+            self.base_name,
+            self.calc_number,
+            self.charge,
+            self.mult,
+            self.pal,
+            self.mem,
         )
         max_len = max([len(key) for key in self.keys_calcs.keys()])
         self.max_fmt = f" >{max_len+2}"

pysisyphus/calculators/ORCA.py

@@ -17,7 +17,7 @@ def make_sym_mat(table_block):
     mat_size = int(table_block[1])
     # Orca prints blocks of 5 columns
     arr = np.array(table_block[2:], dtype=float)
-    assert arr.size == mat_size**2
+    assert arr.size == mat_size ** 2
     block_size = 5 * mat_size
     cbs = [
         arr[i * block_size : (i + 1) * block_size].reshape(mat_size, -1)
@@ -556,13 +556,13 @@ def parse_gbw(self, gbw_fn):
             # print('Number of Operators: {}'.format(operators))
             # print('Basis Dimension: {}'.format(dimension))
 
-            coeffs_fmt = "<" + dimension**2 * "d"
+            coeffs_fmt = "<" + dimension ** 2 * "d"
 
             assert operators == 1, "Unrestricted case is not implemented!"
 
             for i in range(operators):
                 # print('\nOperator: {}'.format(i))
-                coeffs = struct.unpack(coeffs_fmt, handle.read(8 * dimension**2))
+                coeffs = struct.unpack(coeffs_fmt, handle.read(8 * dimension ** 2))
                 occupations = struct.iter_unpack("<d", handle.read(8 * dimension))
                 energies = struct.iter_unpack("<d", handle.read(8 * dimension))
                 irreps = struct.iter_unpack("<i", handle.read(4 * dimension))
@@ -605,7 +605,7 @@ def set_mo_coeffs_in_gbw(in_gbw_fn, out_gbw_fn, mo_coeffs):
 
         tot_offset = offset + 4 + 4
         start = gbw_bytes[:tot_offset]
-        end = gbw_bytes[tot_offset + 8 * dimension**2 :]
+        end = gbw_bytes[tot_offset + 8 * dimension ** 2 :]
         # Construct new gbw content by replacing the MO coefficients in the middle
         mod_gbw_bytes = start + mo_coeffs.T.tobytes() + end
 

pysisyphus/cos/ChainOfStates.py

@@ -16,7 +16,7 @@
 
 class ChainOfStates:
     logger = logging.getLogger("cos")
-    valid_coord_types = "cart dlc".split()
+    valid_coord_types = ("cart", "cartesian", "dlc")
 
     def __init__(
         self,
@@ -87,6 +87,14 @@ def __init__(
         except AttributeError:
             self.typed_prims = None
 
+    @property
+    def calculator(self):
+        try:
+            calc = self.images[0].calculator
+        except IndexError:
+            calc = None
+        return calc
+
     def log(self, message):
         self.logger.debug(f"Counter {self.counter+1:03d}, {message}")
 
@@ -186,7 +194,7 @@ def set_coords_at(self, i, coords):
         Then tries to set cartesian coordinate as self.images[i].coords
         which will raise an error when coord_type != "cart".
         """
-        assert self.images[i].coord_type == "cart", (
+        assert self.images[i].coord_type in ("cart", "cartesian"), (
             "ChainOfStates.set_coords_at() has to be reworked to support "
             "internal coordiantes. Try to set 'align: False' in the 'opt' "
             "section of the .yaml input file."
@@ -259,7 +267,6 @@ def calculate_forces(self):
             all_energies = np.array([image.all_energies for image in self.images])
             energy_diffs = np.diff(all_energies, axis=1).flatten()
             calc_inds = all_energies.argmin(axis=1)
-            print("calc_inds", calc_inds, ";", len(calc_inds), "images")
             mix_at = []
             for i, calc_ind in enumerate(calc_inds[:-1]):
                 next_ind = calc_inds[i + 1]
@@ -274,7 +281,9 @@ def calculate_forces(self):
             for ind in mix_at:
                 self.images[ind].calculator.mix = True
                 # Recalculate correct energy and forces
-                print(f"Switch after calc_ind={calc_ind} at index {ind}. Recalculating.")
+                print(
+                    f"Switch after calc_ind={calc_ind} at index {ind}. Recalculating."
+                )
                 self.images[ind].calc_energy_and_forces()
                 self.org_forces_indices.append(ind)
                 calc_ind = calc_inds[ind]
@@ -355,9 +364,22 @@ def set_zero_forces_for_fixed_images(self):
             self.images[-1].cart_forces = zero_forces
             self.log("Zeroed forces on fixed last image.")
 
-    def get_tangent(self, i, kind="upwinding", lanczos_guess=None):
+    def get_tangent(
+        self, i, kind="upwinding", lanczos_guess=None, disable_lanczos=False
+    ):
         """[1] Equations (8) - (11)"""
 
+        # Converge to lowest curvature mode at the climbing image.
+        # In the current implementation the given kind may be overwritten when
+        # Lanczos iterations are enabled and there are climbing images. By
+        # setting 'disable_lanczos=True' the provided kind is never overwritten.
+        if (
+            not disable_lanczos
+            and self.started_climbing_lanczos
+            and (i in self.get_climbing_indices())
+        ):
+            kind = "lanczos"
+
         tangent_kinds = ("upwinding", "simple", "bisect", "lanczos")
         assert kind in tangent_kinds, "Invalid kind! Valid kinds are: {tangent_kinds}"
         prev_index = max(i - 1, 0)