Fix transverse refinement and transition-state acceptance

RELEASING.md

@@ -1,3 +1,5 @@
+Note that this file is going to be AI Generated. If you're working on this, have your agent do the job. 
+
 # Releasing pele
 
 Release tested commits from `main` using matching version tags and GitHub Releases.
@@ -60,6 +62,54 @@ rattler-build build -r conda-recipe/recipe.yaml -c conda-forge \
   --variant c_stdlib=sysroot --variant c_stdlib_version=2.17
 ```
 
+## Deployment gotchas
+
+### Floating point arithmetic differs across platforms
+
+A fixed random seed does not guarantee identical numerical results across
+platforms. CPU instructions, compilers and math libraries can round intermediate
+operations differently. Non-convex optimization can amplify differences in the
+last bits and reach different stationary points, even when the initial coordinates
+and random-generator state agree exactly.
+
+For example, the tested macOS ARM NumPy 2.5.4 wheel for CPython 3.14 uses a
+fused multiply-add to compute
+`x1*x1 + x2*x2` in its legacy Gaussian conversion. Fusing the multiply and add
+rounds that combined operation once; separate operations round the product first.
+Replaying this arithmetic reproduced all 54 recorded macOS Gaussian values,
+including six that differed from Linux despite identical generator states.
+[NumPy's conversion source](https://github.com/numpy/numpy/blob/main/numpy/random/src/legacy/legacy-distributions.c#L18-L38)
+shows the expression involved.
+
+The [seeded LJ18 investigation](https://github.com/pele-python/pele/actions/runs/38177436370)
+also recorded a difference of one unit in the last place (1 ULP) in the first
+native minimization step. The resulting macOS saddle had exactly one unstable physical
+mode and a force RMS of `8.08e-5`, satisfying the `1e-4` stopping tolerance.
+Its transverse curvature search nevertheless exhausted 100 iterations with an
+RMS of `1.96e-5`, above its separate `1e-5` tolerance. At the same fixed coordinates,
+that search converged after 103 iterations when given a larger budget. A small
+iteration cap can therefore reject a valid candidate on one platform while
+passing on another. The [trace on PR #183](https://github.com/pele-python/pele/actions/runs/38178661821)
+reproduced this ARM result; its [Intel macOS wheel trace](https://github.com/pele-python/pele/actions/runs/38178661825)
+found the same cause, with the fixed-coordinate searches needing up to 117 iterations.
+
+A later [Intel macOS Python 3.12 trajectory trace](https://github.com/pele-python/pele/actions/runs/38184087581)
+reused one wheel across two runners and repeated the original seeded accounting
+fixture. All 44 starting configurations agreed bit for bit, but four searches
+failed: three after 21 genuine positive-curvature decisions, and one after
+exhausting 100 outer iterations with two unstable modes and a large force RMS.
+None of the 1,064 recorded primary snapshots misclassified the curvature sign,
+and evaluation accounting remained exact. In one positive-curvature failure,
+independent replay confirmed a fresh force RMS of `4.85e-4`, above `1e-4`, and
+48 positive projected physical curvatures, with the lowest `+0.10566`.
+These were legitimate failed searches, unlike the valid saddle rejected by the
+100-iteration cap. The trace reproduces the abort pattern; the earlier CI
+failure's exact coordinates were not recorded. Accounting tests therefore start
+near the existing known saddle, with a physical perturbation that still requires
+refinement, and retain their success and exact-count assertions. A loosely
+quenched random configuration does not guarantee a successful global
+transition-state search within a finite iteration budget.
+
 ## Configure publishing accounts
 
 Create accounts on [PyPI](https://pypi.org/) and [TestPyPI](https://test.pypi.org/).

pele/landscape/tests/test_double_ended_connect.py

@@ -45,10 +45,24 @@ def test2(self):
         # Allow enough search iterations across compiler and runtime variants.
         connect = DoubleEndedConnect(
             m1, m2, pot, mindist, db, verbosity=10,
-            local_connect_params={"NEBparams": {"iter_density": 100}},
+            local_connect_params={"NEBparams": {"iter_density": 100},
+                                  "tsSearchParams": {"index_seed": 0}},
         )
         connect.connect()
         self.assertTrue(connect.success())
+        for transition_state in db.transition_states():
+            coords = transition_state.coords
+            xyz = coords.reshape(-1, 3)
+            xyz = xyz - xyz.mean(axis=0)
+            rigid = np.column_stack(
+                [np.tile(axis, len(xyz)) for axis in np.eye(3)]
+                + [np.cross(axis, xyz).ravel() for axis in np.eye(3)]
+            )
+            physical = np.linalg.qr(rigid, mode="complete")[0][:, 6:]
+            hessian = pot.getHessian(coords)
+            spectrum = np.linalg.eigvalsh(physical.T @ hessian @ physical)
+            self.assertLess(spectrum[0], 0.0)
+            self.assertGreater(spectrum[1], 0.0)
 
         path = connect.returnPath()
 

pele/systems/tests/test_frozen_lj.py

@@ -1,12 +1,16 @@
 import unittest
 
+import numpy as np
+
 from pele.systems.ljcluster_frozen import LJClusterFrozen
 from pele.systems import LJCluster
 from pele.thermodynamics import get_thermodynamic_information
 
 
 class TestLJCluster(unittest.TestCase):
     def setUp(self):
+        self.addCleanup(np.random.set_state, np.random.get_state())
+        np.random.seed(0)
         self.prepare_system(frozen_atoms=[0, 2, 4, 6])
 
     def prepare_system(self, frozen_atoms=None):
@@ -56,6 +60,7 @@ def test_basinhopping(self):
             bh.run(1)
 
     def test_connect(self):
+        self.system.params.double_ended_connect.local_connect_params.tsSearchParams.index_seed = 0
         db = self.make_database()
         min1, min2 = db.minima()[:2]
         connect = self.system.get_double_ended_connect(min1, min2, db)

pele/transition_states/_transition_state_refinement.py

@@ -4,6 +4,7 @@
 from pele.optimize import Result
 from pele.optimize import mylbfgs
 from pele.transition_states import FindLowestEigenVector
+from pele.transition_states._orthogopt import orthogopt
 from pele.transition_states._dimer_translator import _DimerTranslator
 from pele.transition_states._transverse_walker import _TransverseWalker
 from pele.utils.hessian import get_smallest_eig
@@ -48,6 +49,9 @@ class FindTransitionState:
                                 rotational symmetry
             orthogZeroEigs=None : the vector is unchanged
 
+    index_seed : int or None
+        Seed for the private random generator used by the additional-mode screen.
+        None uses fresh entropy without consuming the global NumPy random stream.
     lowestEigenvectorQuenchParams : dict
         these parameters are passed to the quench routine for he lowest
         eigenvector search
@@ -96,6 +100,12 @@ class FindTransitionState:
     if the tolerance for tangent space search is lower than the total tolerance,
     then it will never finish
 
+    A converged candidate is screened for additional unstable physical modes using
+    two iterative searches from independent random starts. This uses linear memory;
+    a passing screen is not an exhaustive certificate of Hessian index. The default
+    secondary-search tolerance bounds the gradient norm rather than a dimension-
+    dependent RMS. An explicitly supplied quench tolerance retains its RMS meaning.
+
     See Also
     --------
     findTransitionState : function wrapper for this class
@@ -126,8 +136,11 @@ def __init__(
         check_negative=True,
         invert_gradient=False,
         hessian_diagonalization=False,
+        index_seed=None,
     ):
         self.pot = pot
+        self.index_seed = index_seed
+        self._index_random = np.random.default_rng(index_seed)
         self.coords = np.copy(coords)
         self.nfev = 0
         self.tol = tol
@@ -224,6 +237,7 @@ def params(cls, obj=None):
         params[
             "lowestEigenvectorQuenchParams"
         ] = obj.lowestEigenvectorQuenchParams.copy()
+        params["index_seed"] = obj.index_seed
         params["tol"] = obj.tol
         params["nsteps"] = obj.nsteps
         params["nfail_max"] = obj.nfail_max
@@ -285,6 +299,7 @@ def run(self):
 
         self._compute_gradients(coords)
         iend = 0
+        index_screened = False
         for i in range(self.nsteps):
             iend = i
             # get the lowest eigenvalue and eigenvector
@@ -385,7 +400,29 @@ def run(self):
                     stepnum=i,
                 )
             if rms < self.tol:
-                break
+                self._getLowestEigenVector(coords, i, gradient=grad)
+                if (not self.leig_result.success
+                        or not np.isfinite(self.eigenval) or self.eigenval >= 0):
+                    break
+                transverse = self._get_lowest_transverse_eigenvector(coords, grad)
+                if not transverse.success:
+                    res.message.append("could not complete transverse stability screen")
+                    break
+                if transverse.eigenval > 0:
+                    index_screened = True
+                    break
+                res.message.append("additional negative curvature found")
+                if i + 1 == self.nsteps:
+                    break
+                direction = transverse.eigenvec
+                if np.dot(grad, direction) > 0:
+                    direction = -direction
+                coords = coords + min(0.01, self.maxstep_tangent) * direction
+                self._compute_gradients(coords)
+                if self._transverse_walker is not None:
+                    self.nfev += self._transverse_walker.get_result().nfev
+                    self._transverse_walker = None
+                continue
             if self.nfail >= self.nfail_max:
                 logger.warning(
                     "stopping findTransitionState.  too many failures in eigenvector search %s",
@@ -410,8 +447,12 @@ def run(self):
                     )
                     break
 
-        # done.  do one last eigenvector search because coords may have changed
-        self._getLowestEigenVector(coords, iend)
+        # A candidate passing the screen already has its final eigenpair checked.
+        if not index_screened:
+            self._getLowestEigenVector(coords, iend)
+        E = self.get_energy()
+        grad = self.get_gradient()
+        rms = np.linalg.norm(grad) * self.rmsnorm
 
         # print some data
         if self.verbosity > 0 or self.iprint > 0:
@@ -424,7 +465,7 @@ def run(self):
                 self.eigenval,
             )
 
-        success = True
+        success = index_screened
         # check if results make sense
         if self.eigenval >= 0.0:
             if self.verbosity > 2:
@@ -459,6 +500,89 @@ def run(self):
         res.nfev = self.nfev
         return res
 
+    def _get_lowest_transverse_eigenvector(self, coords, gradient):
+        """Screen for another unstable physical mode without forming a Hessian."""
+        orthogonalize = self.orthogZeroEigs
+        if orthogonalize == 0:
+            orthogonalize = orthogopt
+
+        def physical(vector):
+            vector = vector.copy()
+            if orthogonalize is not None:
+                vector = orthogonalize(vector, coords)
+            return vector
+
+        lowest = physical(self.eigenvec)
+        norm = np.linalg.norm(lowest)
+        if not np.isfinite(norm) or norm < 1e-12:
+            return Result(success=False)
+        lowest /= norm
+
+        def project(vector, position):
+            vector = physical(vector)
+            return vector - np.dot(vector, lowest) * lowest
+
+        params = self.lowestEigenvectorQuenchParams.copy()
+        # Unit eigenvectors need a dimension-independent gradient-norm threshold.
+        params.setdefault("tol", 1e-5 / np.sqrt(coords.size))
+        step = params.get("dx", 1e-6)
+
+        def hessian_product(vector):
+            _, plus = self.pot.getEnergyGradient(coords + step * vector)
+            _, minus = self.pot.getEnergyGradient(coords - step * vector)
+            self.nfev += 2
+            return (plus - minus) / (2 * step)
+
+        # Central differences also distinguish a negative mode from rigid-mode leakage.
+        primary_product = hessian_product(lowest)
+        primary_curvature = float(np.dot(lowest, primary_product))
+        primary_residual = np.linalg.norm(physical(primary_product) - primary_curvature * lowest)
+        if (not np.isfinite(primary_residual)
+                or not np.isfinite(primary_curvature)
+                or primary_curvature >= -primary_residual):
+            return Result(success=False)
+
+        # shortcut: two iterative starts can miss modes; use exhaustive spectral validation if certification is required.
+        # A local generator preserves the connection search's random sequence.
+        random = self._index_random
+        for _ in range(2):
+            guess = project(random.normal(size=coords.size), coords)
+            if np.linalg.norm(guess) < 1e-12:
+                for axis in range(coords.size):
+                    guess = np.zeros(coords.size)
+                    guess[axis] = 1.0
+                    guess = project(guess, coords)
+                    if np.linalg.norm(guess) >= 1e-12:
+                        break
+                else:
+                    # Only the checked unstable physical direction remains.
+                    return Result(success=True, eigenval=np.inf)
+            optimizer = FindLowestEigenVector(
+                coords, self.pot, eigenvec0=guess, orthogZeroEigs=project,
+                gradient=gradient, **params,
+            )
+            result = optimizer.run()
+            self.nfev += result.nfev
+            if not result.success or not np.isfinite(result.eigenval):
+                result.success = False
+                return result
+
+            direction = result.eigenvec
+            product = hessian_product(direction)
+            curvature = float(np.dot(direction, product))
+            residual = np.linalg.norm(project(product, coords) - curvature * direction)
+            coupling = 0.5 * (np.dot(lowest, product) + np.dot(direction, primary_product))
+            restricted = np.array([[primary_curvature, coupling], [coupling, curvature]])
+            result.eigenval = curvature
+            if (not np.all(np.isfinite(restricted)) or not np.isfinite(residual)
+                    or abs(curvature) <= residual):
+                result.success = False
+                return result
+            if curvature < 0:
+                result.success = np.linalg.eigvalsh(restricted)[-1] < 0
+                return result
+        return result
+
     def _get_lowest_eigenvector_RR(self, coords, gradient=None):
         """get the lowest eigenvector using Reyleigh Ritz minimization"""
         if "nsteps" in self.lowestEigenvectorQuenchParams:
@@ -614,6 +738,8 @@ def _minimizeTangentSpace(self, coords, energy=None, gradient=None):
 
     def _update_max_uphill_step(self, Fold, stepsize):
         """use a trust radius to update the maximum uphill step size"""
+        if stepsize == 0:
+            return
         Fnew = np.dot(self.eigenvec, self.get_gradient())
         # EPER=MIN(DABS(1.0D0-(FOBNEW-FOB)/(PSTEP*EVALMIN)),DABS(1.0D0-(-FOBNEW-FOB)/(PSTEP*EVALMIN)))
         a1 = 1.0 - (Fnew - Fold) / (stepsize * self.eigenval)

pele/transition_states/_transverse_walker.py

@@ -111,10 +111,17 @@ def __init__(
             gradient=transverse_gradient,
             **minimizer_kwargs
         )
+        self._initial_H0 = self.walker.H0
 
     def update_eigenvec(self, eigenvec, eigenval):
-        """update the vecotr"""
+        """update the vector"""
+        old_vector = self.tspot.vector
         self.tspot.update_vector(eigenvec)
+        if not (np.array_equal(old_vector, self.tspot.vector) or
+                np.array_equal(old_vector, -self.tspot.vector)):
+            # Secant history from another projection can move along the excluded mode.
+            self.walker.reset()
+            self.walker.H0 = self._initial_H0
 
     def update_maxstep(self, maxstep):
         """update the maximum step size of the minimizer"""

pele/transition_states/tests/lj13_index_two.xyz

@@ -0,0 +1,15 @@
+13
+Recorded index-two candidate from the connection test
+LJ 1.2122346883612989 -0.6929443307790788 0.45087668210298265
+LJ 1.3493575684234813 0.30563112715902663 -0.035345239198378668
+LJ -0.99532473343649819 -0.57767505504430305 0.1432786388422807
+LJ 1.5784817299229652 -0.96002378359696505 -0.57925284369424734
+LJ 0.6186264596851152 0.62569254249252071 -0.82016815175178082
+LJ 1.2218868343284028 -0.16768392951581501 -1.3429024291359211
+LJ -0.26761773809639516 -0.96910258148931672 -0.6026477818204794
+LJ 0.65753316464245781 -0.41046513866176954 -0.45514205783651035
+LJ 0.64174869020119463 -1.4557093976098245 -0.14256864354711837
+LJ -0.27316555555765437 0.13206557898505247 -0.32581341182977808
+LJ 0.076108603335255084 -0.19908656303161248 -1.3538815467639136
+LJ 0.08024908140762875 -0.65300286334907787 0.4309591072748527
+LJ 2.1375687167827615 -0.01370930555883952 -0.73458586264199499