From 6e004b15e05922c242d6bf1f58b5cb2e11f198ee Mon Sep 17 00:00:00 2001 From: Praharsh Suryadevara Date: Sun, 11 Oct 2026 16:24:13 -0500 Subject: [PATCH 01/11] Fix transverse refinement and saddle acceptance --- .../tests/test_double_ended_connect.py | 13 ++ .../_transition_state_refinement.py | 120 +++++++++++++++++- pele/transition_states/_transverse_walker.py | 2 + .../tests/lj13_index_two.xyz | 15 +++ .../tests/test_find_transition_state_hef.py | 22 +++- .../tests/test_transition_state_index.py | 110 ++++++++++++++++ .../tests/test_transition_state_search.py | 49 ++++++- .../tests/test_transverse_walker.py | 34 ++++- 8 files changed, 358 insertions(+), 7 deletions(-) create mode 100644 pele/transition_states/tests/lj13_index_two.xyz create mode 100644 pele/transition_states/tests/test_transition_state_index.py diff --git a/pele/landscape/tests/test_double_ended_connect.py b/pele/landscape/tests/test_double_ended_connect.py index 3c1a83c49..793f2687b 100644 --- a/pele/landscape/tests/test_double_ended_connect.py +++ b/pele/landscape/tests/test_double_ended_connect.py @@ -49,6 +49,19 @@ def test2(self): ) 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() diff --git a/pele/transition_states/_transition_state_refinement.py b/pele/transition_states/_transition_state_refinement.py index 3d19f1783..544fa548e 100644 --- a/pele/transition_states/_transition_state_refinement.py +++ b/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 @@ -285,6 +286,7 @@ def run(self): self._compute_gradients(coords) iend = 0 + index_validated = False for i in range(self.nsteps): iend = i # get the lowest eigenvalue and eigenvector @@ -385,7 +387,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 verify transverse stability") + break + if transverse.eigenval > 0: + index_validated = 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 +434,12 @@ def run(self): ) break - # done. do one last eigenvector search because coords may have changed - self._getLowestEigenVector(coords, iend) + # A successful candidate already has its final eigenpair checked. + if not index_validated: + 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 +452,7 @@ def run(self): self.eigenval, ) - success = True + success = index_validated # check if results make sense if self.eigenval >= 0.0: if self.verbosity > 2: @@ -459,6 +487,88 @@ def run(self): res.nfev = self.nfev return res + def _get_lowest_transverse_eigenvector(self, coords, gradient): + """Check 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() + params.setdefault("nsteps", 100) + 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 = np.random.default_rng(0) + 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 +724,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) diff --git a/pele/transition_states/_transverse_walker.py b/pele/transition_states/_transverse_walker.py index c5503ea37..ae5b15986 100644 --- a/pele/transition_states/_transverse_walker.py +++ b/pele/transition_states/_transverse_walker.py @@ -115,6 +115,8 @@ def __init__( def update_eigenvec(self, eigenvec, eigenval): """update the vecotr""" self.tspot.update_vector(eigenvec) + # Secant history from another projection can move along the excluded mode. + self.walker.reset() def update_maxstep(self, maxstep): """update the maximum step size of the minimizer""" diff --git a/pele/transition_states/tests/lj13_index_two.xyz b/pele/transition_states/tests/lj13_index_two.xyz new file mode 100644 index 000000000..bf4aba4a1 --- /dev/null +++ b/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 diff --git a/pele/transition_states/tests/test_find_transition_state_hef.py b/pele/transition_states/tests/test_find_transition_state_hef.py index 7966ec1dd..3240b49c1 100644 --- a/pele/transition_states/tests/test_find_transition_state_hef.py +++ b/pele/transition_states/tests/test_find_transition_state_hef.py @@ -40,13 +40,32 @@ def test_params(self): params = FindTransitionState.params() def test1(self): - x = self.system.get_random_configuration() + path = os.path.join(os.path.dirname(__file__), "lj18_ts.xyz") + with open(path) as source: + x = read_xyz(source).coords.ravel() + perturbation = np.sin(np.arange(x.size) + 1.0) + perturbation = self.system.get_orthogonalize_to_zero_eigenvectors()(perturbation, x) + x = x + 0.01 * perturbation / np.linalg.norm(perturbation) dimer = self.make_dimer(x) res = dimer.run() self.assertTrue(res.success) self.assertEqual(res.nfev, self.pot.nfev) self.assertGreater(res.nfev, 1) + # Check physical curvature independently of the iterative index validation. + _, _, hessian = self.pot.pot.getEnergyGradientHessian(res.coords) + centered = res.coords.reshape(-1, 3) + centered = centered - centered.mean(axis=0) + rigid = np.column_stack([ + np.tile(axis, len(centered)) for axis in np.eye(3) + ] + [np.cross(np.tile(axis, (len(centered), 1)), centered).ravel() + for axis in np.eye(3)]) + self.assertEqual(np.linalg.matrix_rank(rigid), 6) + basis = np.linalg.qr(rigid, mode="complete")[0][:, 6:] + eigenvalues = np.linalg.eigvalsh(basis.T @ hessian @ basis) + self.assertLess(eigenvalues[0], -1e-5) + self.assertGreater(eigenvalues[1], 1e-5) + def test2(self): # get the path of the file directory path = os.path.dirname(os.path.abspath(__file__)) @@ -77,6 +96,7 @@ def make_dimer(self, x, **kwargs): self.pot, orthogZeroEigs=self.system.get_orthogonalize_to_zero_eigenvectors(), invert_gradient=True, + **kwargs, ) diff --git a/pele/transition_states/tests/test_transition_state_index.py b/pele/transition_states/tests/test_transition_state_index.py new file mode 100644 index 000000000..e78e9ef4b --- /dev/null +++ b/pele/transition_states/tests/test_transition_state_index.py @@ -0,0 +1,110 @@ +import unittest +from unittest import mock + +import numpy as np + +from pele.potentials import BasePotential +from pele.transition_states import FindTransitionState + + +class _QuadraticPotential(BasePotential): + def __init__(self, hessian): + self.hessian = np.asarray(hessian, dtype=float) + self.nfev = 0 + + def getEnergyGradient(self, coords): + self.nfev += 1 + gradient = self.hessian @ coords + return 0.5 * np.dot(coords, gradient), gradient + + def getHessian(self, coords): + raise AssertionError("Default index validation must not form a dense Hessian") + + +class TestTransitionStateIndexEdges(unittest.TestCase): + def test_rejects_rotated_index_two_with_stable_initial_rr_direction(self): + # Eigenvalues are -2, -1, +1; the first transverse probe points along +1. + pot = _QuadraticPotential([ + [-2.0, 0.0, 0.0], + [0.0, -0.918372375740774, -0.3957172974185572], + [0.0, -0.3957172974185572, 0.918372375740774], + ]) + ret = FindTransitionState( + np.zeros(3), pot, eigenvec0=np.array([1.0, 0.0, 0.0]), + orthogZeroEigs=None, nsteps=1, verbosity=0, + ).run() + self.assertFalse(ret.success) + + def test_accepts_index_one_when_zero_modes_leave_no_transverse_direction(self): + def remove_zero_mode(vector, coords): + return np.array([vector[0], 0.0]) + + ret = FindTransitionState( + np.zeros(2), _QuadraticPotential([[-2.0, 0.0], [0.0, 0.0]]), + eigenvec0=np.array([1.0, 0.0]), orthogZeroEigs=remove_zero_mode, + nsteps=1, verbosity=0, + ).run() + self.assertTrue(ret.success) + self.assertLess(ret.rms, 1e-4) + + def test_accepts_one_dimensional_index_one_stationary_point(self): + ret = FindTransitionState( + np.zeros(1), _QuadraticPotential([[-2.0]]), + eigenvec0=np.ones(1), orthogZeroEigs=None, nsteps=1, verbosity=0, + ).run() + self.assertTrue(ret.success) + self.assertAlmostEqual(ret.eigenval, -2.0) + self.assertLess(ret.rms, 1e-4) + + def test_success_requires_a_converged_returned_primary_eigenpair(self): + search = FindTransitionState( + np.zeros(1), _QuadraticPotential([[-2.0]]), + eigenvec0=np.ones(1), orthogZeroEigs=None, nsteps=1, verbosity=0, + ) + find_primary = search._get_lowest_eigenvector_RR + calls = 0 + + def fail_redundant_refresh(*args, **kwargs): + nonlocal calls + result = find_primary(*args, **kwargs) + calls += 1 + if calls >= 3: + result.success = False + return result + + with mock.patch.object(search, "_get_lowest_eigenvector_RR", + side_effect=fail_redundant_refresh): + ret = search.run() + # Skipping the redundant refresh is valid; accepting its failed result is not. + self.assertTrue(not ret.success or search.leig_result.success) + + def test_recovery_failure_returns_energy_and_gradient_at_returned_coords(self): + pot = _QuadraticPotential(np.diag([-2.0, -1.0, 1.0])) + search = FindTransitionState( + np.zeros(3), pot, eigenvec0=np.array([1.0, 0.0, 0.0]), + orthogZeroEigs=None, demand_initial_negative_vec=True, + nsteps=2, verbosity=0, + ) + search.npositive_max = 0 + find_primary = search._getLowestEigenVector + + def abort_after_recovery(coords, iteration, gradient=None): + overlap = find_primary(coords, iteration, gradient=gradient) + if iteration > 0: + search.eigenval = 1.0 + return overlap + + with mock.patch.object(search, "_getLowestEigenVector", + side_effect=abort_after_recovery): + ret = search.run() + self.assertFalse(ret.success) + self.assertGreater(np.linalg.norm(ret.coords), 1e-5) + self.assertEqual(ret.nfev, pot.nfev) + energy, gradient = pot.getEnergyGradient(ret.coords) + self.assertAlmostEqual(ret.energy, energy, places=12) + np.testing.assert_allclose(ret.grad, gradient, rtol=1e-12, atol=1e-12) + self.assertAlmostEqual(ret.rms, np.linalg.norm(gradient) / np.sqrt(3.0), places=12) + + +if __name__ == "__main__": + unittest.main() diff --git a/pele/transition_states/tests/test_transition_state_search.py b/pele/transition_states/tests/test_transition_state_search.py index f6aae047e..e67af21fa 100644 --- a/pele/transition_states/tests/test_transition_state_search.py +++ b/pele/transition_states/tests/test_transition_state_search.py @@ -1,11 +1,13 @@ import unittest +from pathlib import Path import numpy as np import matplotlib.pyplot as plt from pele.transition_states import FindTransitionState, findTransitionState -from pele.potentials import BasePotential +from pele.potentials import BasePotential, LJ +from pele.utils.xyz import read_xyz from pele.potentials.tests import _base_test from pele.potentials.tests._base_test import assert_arrays_almost_equal @@ -176,5 +178,50 @@ def test1(self): print(ret) +class IndexTwoPotential(BasePotential): + def getEnergyGradient(self, coords): + curvature = np.array([-2.0, -1.0, 1.0]) + return 0.5 * np.dot(curvature * coords, coords), curvature * coords + + def getHessian(self, coords): + raise AssertionError("Default index validation must not form a dense Hessian") + + +class RecoverableIndexTwoPotential(BasePotential): + def getEnergyGradient(self, coords): + x, y, z = coords + return (-x*x + 0.25*(y*y - 1.0)**2 + 0.5*z*z, + np.array([-2.0*x, y*(y*y - 1.0), z])) + + +class TestTransitionStateIndex(unittest.TestCase): + def test_rejects_index_two_stationary_point(self): + ret = FindTransitionState( + np.zeros(3), IndexTwoPotential(), eigenvec0=np.array([1., 0., 0.]), + orthogZeroEigs=None, nsteps=1, verbosity=0, + ).run() + self.assertFalse(ret.success) + + def test_rejects_recorded_lj13_index_two_saddle(self): + with Path(__file__).with_name("lj13_index_two.xyz").open() as source: + coords = read_xyz(source).coords.ravel() + pot = LJ() + values, vectors = np.linalg.eigh(pot.getHessian(coords)) + self.assertLess(values[1], -1.0) + ret = FindTransitionState( + coords, pot, eigenvec0=vectors[:, 0], nsteps=1, verbosity=0, + ).run() + self.assertFalse(ret.success) + + def test_escapes_extra_unstable_mode(self): + ret = FindTransitionState( + np.zeros(3), RecoverableIndexTwoPotential(), + eigenvec0=np.array([1., 0., 0.]), orthogZeroEigs=None, verbosity=0, + ).run() + self.assertTrue(ret.success) + self.assertAlmostEqual(abs(ret.coords[1]), 1.0, delta=1e-3) + self.assertLess(ret.rms, 1e-4) + + if __name__ == "__main__": unittest.main() diff --git a/pele/transition_states/tests/test_transverse_walker.py b/pele/transition_states/tests/test_transverse_walker.py index 1aa8ea4ea..6a5287360 100644 --- a/pele/transition_states/tests/test_transverse_walker.py +++ b/pele/transition_states/tests/test_transverse_walker.py @@ -23,7 +23,39 @@ def test(self): xnew = ret.coords dx = xnew - x0 - self.assertLess(np.dot(dx, self.evec), 1e-6) + self.assertLess(abs(np.dot(dx, self.evec)), 1e-6) + + +class TestTransverseWalkerHistory(unittest.TestCase): + def test_motion_stays_in_updated_plane(self): + class QuadraticPotential: + hessian = np.array([[3.0, 1.0, 0.5], + [1.0, 2.0, 0.25], + [0.5, 0.25, 1.0]]) + + def getEnergyGradient(self, coords): + gradient = self.hessian @ coords + return 0.5 * np.dot(coords, gradient), gradient + + pot = QuadraticPotential() + opt = _TransverseWalker( + np.array([1.0, -0.7, 0.4]), pot, np.array([0.0, 0.0, 1.0]), + tol=1e-12, + ) + # Build valid secant history in the original transverse plane. + opt.run(3) + axis = np.array([1.0, 0.0, 0.0]) + coords = np.array([-0.5, 1.0, 1.0]) + energy, gradient = pot.getEnergyGradient(coords) + opt.update_eigenvec(axis, -1.0) + opt.update_coords(coords, energy, gradient) + + start = coords.copy() + for _ in range(5): + previous = coords.copy() + coords = opt.run(1).coords.copy() + self.assertLess(abs(np.dot(coords - previous, axis)), 1e-14) + self.assertGreater(np.linalg.norm(coords - start), 1e-5) class TestTransverseWalker_NFEV(unittest.TestCase): From 3cb19393964aec7dc2343b2b77ebc933f1e87c1a Mon Sep 17 00:00:00 2001 From: Praharsh Suryadevara Date: Sun, 11 Oct 2026 17:00:18 -0500 Subject: [PATCH 02/11] Trace seeded transition-state failures in PR CI --- .github/workflows/test.yml | 1 + .../tests/test_find_transition_state_hef.py | 100 +++++++++++++++++- 2 files changed, 100 insertions(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index a89f61601..09c30559f 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -9,6 +9,7 @@ on: env: BUILD_TYPE: Release + PYTEST_ADDOPTS: --capture=tee-sys # a new push to a branch/PR cancels its older, still running workflow concurrency: diff --git a/pele/transition_states/tests/test_find_transition_state_hef.py b/pele/transition_states/tests/test_find_transition_state_hef.py index 3240b49c1..d6c1210d3 100644 --- a/pele/transition_states/tests/test_find_transition_state_hef.py +++ b/pele/transition_states/tests/test_find_transition_state_hef.py @@ -1,5 +1,10 @@ import unittest +from unittest import mock import os +import hashlib +import json +import platform +import sys import numpy as np from pele.systems import LJCluster @@ -110,9 +115,102 @@ def setUp(self): self.x = self.system.get_random_minimized_configuration(tol=10.0).coords def do_check(self, **kwargs): + # Temporary PR diagnostic: observe the original search without changing it. + from pele.transition_states import FindLowestEigenVector + from pele.utils import rotations + + def trace(event, **data): + print("HEF_TRACE " + json.dumps(dict(test=self.id(), event=event, **data)), flush=True) + + def digest(array): + return hashlib.sha256(array.tobytes()).hexdigest() + + def summary(result): + data = {} + for name in ("success", "nsteps", "nfev", "rms", "eigenval", "message"): + value = getattr(result, name, None) + if isinstance(value, np.generic): + value = value.item() + if isinstance(value, float) and not np.isfinite(value): + value = None + data[name] = value + return data + + try: + from threadpoolctl import threadpool_info + pools = threadpool_info() + except ImportError: + pools = None + state = np.random.get_state() + rng_hash = hashlib.sha256( + state[1].tobytes() + repr((state[0], state[2], state[3], state[4])).encode() + ).hexdigest() + trace("setup", python=sys.version, numpy=np.__version__, platform=platform.platform(), + machine=platform.machine(), threadpools=pools, + environment={name: os.environ.get(name) for name in + ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "VECLIB_MAXIMUM_THREADS")}, + coords_hash=digest(self.x), coords=self.x.tolist(), rng_hash=rng_hash, + kwargs=kwargs) + + phase = "search" + original_run = FindLowestEigenVector.run + original_random = rotations.vec_random_ndim + + def traced_random(shape): + vector = original_random(shape) + trace("random_axis", phase=phase, vector_hash=digest(vector), values=vector.tolist()) + return vector + + def traced_run(optimizer, *args, **run_kwargs): + orthogonalizer = optimizer.eigpot.orthogZeroEigs + label = "secondary" if getattr(orthogonalizer, "__name__", None) == "project" else "primary" + budget = args[0] if args else run_kwargs.get("niter", optimizer.minimizer.nsteps) + try: + result = original_run(optimizer, *args, **run_kwargs) + except Exception as error: + trace("rr_error", phase=phase, label=label, budget=budget, + error=type(error).__name__, detail=str(error)) + raise + trace("rr", phase=phase, label=label, budget=budget, + first_order=optimizer.eigpot.first_order, dx=optimizer.eigpot.diff, + **summary(result)) + return result + self.pot.nfev = 0 opt = FindTransitionState(self.x, self.pot, **kwargs) - ret = opt.run() + with mock.patch.object(rotations, "vec_random_ndim", side_effect=traced_random), \ + mock.patch.object(FindLowestEigenVector, "run", new=traced_run): + ret = opt.run() + trace("final", energy=float(ret.energy), coords=ret.coords.tolist(), + eigenvec=ret.eigenvec.tolist(), counted_nfev=self.pot.nfev, **summary(ret)) + try: + raw_potential = self.pot.pot + xyz = ret.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)] + ) + rank = int(np.linalg.matrix_rank(rigid)) + basis = np.linalg.qr(rigid, mode="complete")[0][:, rank:] + hessian = raw_potential.getHessian(ret.coords) + values = np.linalg.eigvalsh(basis.T @ hessian @ basis) + trace("physical_spectrum", rigid_rank=rank, values=values.tolist(), + negative_modes=int(np.count_nonzero(values < -1e-5))) + if not ret.success: + for budget in (100, 500): + phase = "fixed_replay_%d" % budget + params = opt.lowestEigenvectorQuenchParams.copy() + params["nsteps"] = budget + replay = FindTransitionState( + ret.coords.copy(), raw_potential, eigenvec0=ret.eigenvec.copy(), + orthogZeroEigs=opt.orthogZeroEigs, + lowestEigenvectorQuenchParams=params, verbosity=0, + ) + result = replay._get_lowest_transverse_eigenvector(ret.coords.copy(), ret.grad.copy()) + trace("fixed_replay", budget=budget, **summary(result)) + except Exception as error: + trace("posthoc_error", error=type(error).__name__, detail=str(error)) self.assertEqual(ret.nfev, self.pot.nfev) self.assertTrue(ret.success) self.assertGreater(ret.nfev, 0) From 73950dbe5135cc3bdbcf32e6a19a3ec6e6f2eef4 Mon Sep 17 00:00:00 2001 From: Praharsh Suryadevara Date: Sun, 11 Oct 2026 17:14:43 -0500 Subject: [PATCH 03/11] Document arithmetic differences during deployment validation --- RELEASING.md | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/RELEASING.md b/RELEASING.md index 278692a79..67977c49a 100644 --- a/RELEASING.md +++ b/RELEASING.md @@ -60,6 +60,49 @@ 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. + +When validating a deployment: + +- Check reported solver success, the force stopping criterion and the expected + number of unstable physical modes after excluding rigid zero modes. Compare + these properties rather than requiring identical final coordinates. +- Keep force, eigenvector-search and test-comparison tolerances distinct. A + larger iteration budget can resolve non-convergence without loosening any + tolerance; it does not guarantee convergence for every candidate. +- Record the Python, NumPy and SciPy versions, architecture, compiler, math + backend, thread settings, random state, coordinates and optimizer iteration + counts. Save a failing candidate and replay it at fixed coordinates to + distinguish arithmetic differences in the path from a convergence or + correctness defect. + ## Configure publishing accounts Create accounts on [PyPI](https://pypi.org/) and [TestPyPI](https://test.pypi.org/). From f00e9638d9bd99fc32fe93a1cba11243945b5e6f Mon Sep 17 00:00:00 2001 From: Praharsh Suryadevara Date: Sun, 11 Oct 2026 17:23:26 -0500 Subject: [PATCH 04/11] Use the default budget for transverse stability checks --- .github/workflows/test.yml | 1 - RELEASING.md | 4 +- .../_transition_state_refinement.py | 1 - .../tests/test_find_transition_state_hef.py | 100 +----------------- .../tests/test_transition_state_index.py | 13 +++ 5 files changed, 17 insertions(+), 102 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 09c30559f..a89f61601 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -9,7 +9,6 @@ on: env: BUILD_TYPE: Release - PYTEST_ADDOPTS: --capture=tee-sys # a new push to a branch/PR cancels its older, still running workflow concurrency: diff --git a/RELEASING.md b/RELEASING.md index 67977c49a..6814dfd15 100644 --- a/RELEASING.md +++ b/RELEASING.md @@ -87,7 +87,9 @@ 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. +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. When validating a deployment: diff --git a/pele/transition_states/_transition_state_refinement.py b/pele/transition_states/_transition_state_refinement.py index 544fa548e..6bcc8af52 100644 --- a/pele/transition_states/_transition_state_refinement.py +++ b/pele/transition_states/_transition_state_refinement.py @@ -510,7 +510,6 @@ def project(vector, position): return vector - np.dot(vector, lowest) * lowest params = self.lowestEigenvectorQuenchParams.copy() - params.setdefault("nsteps", 100) step = params.get("dx", 1e-6) def hessian_product(vector): diff --git a/pele/transition_states/tests/test_find_transition_state_hef.py b/pele/transition_states/tests/test_find_transition_state_hef.py index d6c1210d3..3240b49c1 100644 --- a/pele/transition_states/tests/test_find_transition_state_hef.py +++ b/pele/transition_states/tests/test_find_transition_state_hef.py @@ -1,10 +1,5 @@ import unittest -from unittest import mock import os -import hashlib -import json -import platform -import sys import numpy as np from pele.systems import LJCluster @@ -115,102 +110,9 @@ def setUp(self): self.x = self.system.get_random_minimized_configuration(tol=10.0).coords def do_check(self, **kwargs): - # Temporary PR diagnostic: observe the original search without changing it. - from pele.transition_states import FindLowestEigenVector - from pele.utils import rotations - - def trace(event, **data): - print("HEF_TRACE " + json.dumps(dict(test=self.id(), event=event, **data)), flush=True) - - def digest(array): - return hashlib.sha256(array.tobytes()).hexdigest() - - def summary(result): - data = {} - for name in ("success", "nsteps", "nfev", "rms", "eigenval", "message"): - value = getattr(result, name, None) - if isinstance(value, np.generic): - value = value.item() - if isinstance(value, float) and not np.isfinite(value): - value = None - data[name] = value - return data - - try: - from threadpoolctl import threadpool_info - pools = threadpool_info() - except ImportError: - pools = None - state = np.random.get_state() - rng_hash = hashlib.sha256( - state[1].tobytes() + repr((state[0], state[2], state[3], state[4])).encode() - ).hexdigest() - trace("setup", python=sys.version, numpy=np.__version__, platform=platform.platform(), - machine=platform.machine(), threadpools=pools, - environment={name: os.environ.get(name) for name in - ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "VECLIB_MAXIMUM_THREADS")}, - coords_hash=digest(self.x), coords=self.x.tolist(), rng_hash=rng_hash, - kwargs=kwargs) - - phase = "search" - original_run = FindLowestEigenVector.run - original_random = rotations.vec_random_ndim - - def traced_random(shape): - vector = original_random(shape) - trace("random_axis", phase=phase, vector_hash=digest(vector), values=vector.tolist()) - return vector - - def traced_run(optimizer, *args, **run_kwargs): - orthogonalizer = optimizer.eigpot.orthogZeroEigs - label = "secondary" if getattr(orthogonalizer, "__name__", None) == "project" else "primary" - budget = args[0] if args else run_kwargs.get("niter", optimizer.minimizer.nsteps) - try: - result = original_run(optimizer, *args, **run_kwargs) - except Exception as error: - trace("rr_error", phase=phase, label=label, budget=budget, - error=type(error).__name__, detail=str(error)) - raise - trace("rr", phase=phase, label=label, budget=budget, - first_order=optimizer.eigpot.first_order, dx=optimizer.eigpot.diff, - **summary(result)) - return result - self.pot.nfev = 0 opt = FindTransitionState(self.x, self.pot, **kwargs) - with mock.patch.object(rotations, "vec_random_ndim", side_effect=traced_random), \ - mock.patch.object(FindLowestEigenVector, "run", new=traced_run): - ret = opt.run() - trace("final", energy=float(ret.energy), coords=ret.coords.tolist(), - eigenvec=ret.eigenvec.tolist(), counted_nfev=self.pot.nfev, **summary(ret)) - try: - raw_potential = self.pot.pot - xyz = ret.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)] - ) - rank = int(np.linalg.matrix_rank(rigid)) - basis = np.linalg.qr(rigid, mode="complete")[0][:, rank:] - hessian = raw_potential.getHessian(ret.coords) - values = np.linalg.eigvalsh(basis.T @ hessian @ basis) - trace("physical_spectrum", rigid_rank=rank, values=values.tolist(), - negative_modes=int(np.count_nonzero(values < -1e-5))) - if not ret.success: - for budget in (100, 500): - phase = "fixed_replay_%d" % budget - params = opt.lowestEigenvectorQuenchParams.copy() - params["nsteps"] = budget - replay = FindTransitionState( - ret.coords.copy(), raw_potential, eigenvec0=ret.eigenvec.copy(), - orthogZeroEigs=opt.orthogZeroEigs, - lowestEigenvectorQuenchParams=params, verbosity=0, - ) - result = replay._get_lowest_transverse_eigenvector(ret.coords.copy(), ret.grad.copy()) - trace("fixed_replay", budget=budget, **summary(result)) - except Exception as error: - trace("posthoc_error", error=type(error).__name__, detail=str(error)) + ret = opt.run() self.assertEqual(ret.nfev, self.pot.nfev) self.assertTrue(ret.success) self.assertGreater(ret.nfev, 0) diff --git a/pele/transition_states/tests/test_transition_state_index.py b/pele/transition_states/tests/test_transition_state_index.py index e78e9ef4b..9e51e037b 100644 --- a/pele/transition_states/tests/test_transition_state_index.py +++ b/pele/transition_states/tests/test_transition_state_index.py @@ -56,6 +56,19 @@ def test_accepts_one_dimensional_index_one_stationary_point(self): self.assertAlmostEqual(ret.eigenval, -2.0) self.assertLess(ret.rms, 1e-4) + def test_accepts_stiff_index_one_stationary_point(self): + # One negative mode and eleven positive modes with disparate stiffnesses. + curvatures = np.r_[-2.0, np.geomspace(1.0, 1000.0, 11)] + pot = _QuadraticPotential(np.diag(curvatures)) + primary = np.zeros(curvatures.size) + primary[0] = 1.0 + ret = FindTransitionState( + np.zeros(curvatures.size), pot, eigenvec0=primary, + orthogZeroEigs=None, nsteps=1, verbosity=0, + ).run() + self.assertEqual(ret.nfev, pot.nfev) + self.assertTrue(ret.success, msg=repr(ret)) + def test_success_requires_a_converged_returned_primary_eigenpair(self): search = FindTransitionState( np.zeros(1), _QuadraticPotential([[-2.0]]), From 93fb8556cea1029de194dc299335bca5d3e2f851 Mon Sep 17 00:00:00 2001 From: Praharsh Suryadevara Date: Sun, 11 Oct 2026 17:30:43 -0500 Subject: [PATCH 05/11] Add AI generation note to RELEASING.md Added note about AI-generated content for clarity. --- RELEASING.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/RELEASING.md b/RELEASING.md index 6814dfd15..337c9188a 100644 --- a/RELEASING.md +++ b/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. From 7bb721a2686196f5289a9cb4675f79ada34f67c3 Mon Sep 17 00:00:00 2001 From: Praharsh Suryadevara Date: Sun, 11 Oct 2026 18:00:39 -0500 Subject: [PATCH 06/11] Trace remaining Intel Python 3.12 convergence failure --- .github/workflows/release.yml | 7 ++-- .../tests/test_find_transition_state_hef.py | 32 ++++++++++++++++++- 2 files changed, 33 insertions(+), 6 deletions(-) diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index ac2003496..fbf651467 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -97,14 +97,11 @@ jobs: strategy: fail-fast: false matrix: - python: [cp311, cp312, cp313, cp314] + # Temporary: isolate the remaining Intel Python 3.12 failure. + python: [cp312] platform: - - os: ubuntu-24.04 - arch: x86_64 - os: macos-15-intel arch: x86_64 - - os: macos-15 - arch: arm64 steps: - uses: actions/setup-python@v6 with: diff --git a/pele/transition_states/tests/test_find_transition_state_hef.py b/pele/transition_states/tests/test_find_transition_state_hef.py index 3240b49c1..caf3edb40 100644 --- a/pele/transition_states/tests/test_find_transition_state_hef.py +++ b/pele/transition_states/tests/test_find_transition_state_hef.py @@ -114,7 +114,37 @@ def do_check(self, **kwargs): opt = FindTransitionState(self.x, self.pot, **kwargs) ret = opt.run() self.assertEqual(ret.nfev, self.pot.nfev) - self.assertTrue(ret.success) + if not ret.success: + import json + import platform + import sys + + _, raw_gradient, hessian = self.pot.pot.getEnergyGradientHessian(ret.coords) + centered = ret.coords.reshape(-1, 3) - ret.coords.reshape(-1, 3).mean(axis=0) + rigid = np.column_stack([ + np.tile(axis, len(centered)) for axis in np.eye(3) + ] + [np.cross(np.tile(axis, (len(centered), 1)), centered).ravel() + for axis in np.eye(3)]) + rank = np.linalg.matrix_rank(rigid) + basis = np.linalg.qr(rigid, mode="complete")[0][:, rank:] + spectrum = np.linalg.eigvalsh(basis.T @ hessian @ basis) + primary_curvature = float(ret.eigenvec @ hessian @ ret.eigenvec) + print("HEF_TRACE " + json.dumps({ + "test": self.id(), "python": sys.version, "numpy": np.__version__, + "platform": platform.platform(), "coords": ret.coords.tolist(), + "energy": float(ret.energy), "rms": float(ret.rms), + "raw_rms": float(np.linalg.norm(raw_gradient) / np.sqrt(ret.coords.size)), + "nfev": int(ret.nfev), "counted_nfev": int(self.pot.nfev), + "nsteps": int(ret.nsteps), "message": ret.message, + "eigenvec": ret.eigenvec.tolist(), "eigenval": float(ret.eigenval), + "primary": {key: opt.leig_result.get(key) + for key in ("success", "nsteps", "rms", "eigenval")}, + "primary_curvature": primary_curvature, + "primary_residual": float(np.linalg.norm( + basis.T @ (hessian @ ret.eigenvec - primary_curvature * ret.eigenvec))), + "physical_spectrum": spectrum.tolist(), + }, default=lambda value: value.item() if isinstance(value, np.generic) else str(value))) + self.assertTrue(ret.success, msg=repr(ret)) self.assertGreater(ret.nfev, 0) def test(self): From 6ad7421f57ff65e36222f1daae040a083f9e66e1 Mon Sep 17 00:00:00 2001 From: Praharsh Suryadevara Date: Sun, 11 Oct 2026 18:23:10 -0500 Subject: [PATCH 07/11] Trace Intel HEF trajectories using the existing wheel --- .github/scripts/trace-hef.py | 85 +++++++++++++++++++++++++++++++++ .github/workflows/trace-hef.yml | 50 +++++++++++++++++++ 2 files changed, 135 insertions(+) create mode 100644 .github/scripts/trace-hef.py create mode 100644 .github/workflows/trace-hef.yml diff --git a/.github/scripts/trace-hef.py b/.github/scripts/trace-hef.py new file mode 100644 index 000000000..e40516c59 --- /dev/null +++ b/.github/scripts/trace-hef.py @@ -0,0 +1,85 @@ +"""Temporary observer for the installed Intel wheel; no solver settings change.""" +import json +import platform +import sys +from pathlib import Path +from unittest import mock + +import numpy as np +import pytest +import scipy +from threadpoolctl import threadpool_info + +import pele +from pele.transition_states.tests import test_find_transition_state_hef as tests + +output = open(sys.argv[1], 'w') +original_search = tests.FindTransitionState +active_case = None + + +def record(kind, **values): + output.write(json.dumps(dict(kind=kind, case=active_case, **values), + default=lambda value: value.item() if isinstance(value, np.generic) else str(value)) + '\n') + output.flush() + + +class TraceSearch(original_search): + def _getLowestEigenVector(self, coords, iteration, gradient=None): + overlap = super()._getLowestEigenVector(coords, iteration, gradient=gradient) + energy, raw_gradient, hessian = self.pot.pot.getEnergyGradientHessian(coords) + centered = coords.reshape(-1, 3) - coords.reshape(-1, 3).mean(axis=0) + rigid = np.column_stack([ + np.tile(axis, len(centered)) for axis in np.eye(3) + ] + [np.cross(np.tile(axis, (len(centered), 1)), centered).ravel() + for axis in np.eye(3)]) + basis = np.linalg.qr(rigid, mode='complete')[0][:, np.linalg.matrix_rank(rigid):] + spectrum = np.linalg.eigvalsh(basis.T @ hessian @ basis) + curvature = float(self.eigenvec @ hessian @ self.eigenvec) + record('primary', iteration=iteration, npositive=self.npositive, + coords=coords.tolist(), eigenvec=self.eigenvec.tolist(), + energy=energy, raw_rms=float(np.linalg.norm(raw_gradient) / np.sqrt(coords.size)), + primary={key: self.leig_result.get(key) for key in ('success', 'nsteps', 'rms', 'eigenval')}, + primary_curvature=curvature, + primary_residual=float(np.linalg.norm(basis.T @ (hessian @ self.eigenvec - curvature * self.eigenvec))), + physical_spectrum=spectrum.tolist()) + return overlap + + def run(self): + record('start', coords=self.coords.tolist()) + result = super().run() + record('result', success=result.success, energy=result.energy, rms=result.rms, + nsteps=result.nsteps, nfev=result.nfev, counted_nfev=self.pot.nfev, + coords=result.coords.tolist(), eigenval=result.eigenval, message=result.message) + assert result.nfev == self.pot.nfev + return result + + +class TraceAccountingTests: + def pytest_runtest_setup(self, item): + global active_case + if item.cls is tests.TestFindTransitionState_NFEV: + active_case = item.nodeid + tests.FindTransitionState = TraceSearch + + def pytest_runtest_teardown(self, item): + tests.FindTransitionState = original_search + + +record('runtime', python=sys.version, numpy=np.__version__, scipy=scipy.__version__, + platform=platform.platform(), threads=threadpool_info()) +package = Path(pele.__file__).parent +suite_status = pytest.main([str(package), '--ignore=' + str(package / 'gui'), '-q'], + plugins=[TraceAccountingTests()]) +record('suite', status=int(suite_status)) +for trial in range(10): + for variant in ('test', 'test2'): + active_case = f'repeat-{trial}-{variant}' + case = tests.TestFindTransitionState_NFEV(variant) + case.setUp() + with mock.patch.object(tests, 'FindTransitionState', TraceSearch): + try: + getattr(case, variant)() + except AssertionError as error: + record('assertion', message=str(error)) +output.close() diff --git a/.github/workflows/trace-hef.yml b/.github/workflows/trace-hef.yml new file mode 100644 index 000000000..8a284e385 --- /dev/null +++ b/.github/workflows/trace-hef.yml @@ -0,0 +1,50 @@ +name: Trace Intel macOS HEF trajectory + +on: + pull_request: + branches: [main] + paths: [.github/workflows/trace-hef.yml, .github/scripts/trace-hef.py] + workflow_dispatch: + +permissions: + contents: read + actions: read + +jobs: + trace: + runs-on: macos-15-intel + timeout-minutes: 15 + strategy: + fail-fast: false + matrix: + sample: [1, 2] + steps: + - uses: actions/checkout@v6 + with: + persist-credentials: false + - uses: actions/setup-python@v6 + with: + python-version: '3.12.10' + - uses: actions/download-artifact@v6 + with: + name: pele-wheels-macos-15-intel-x86_64-cp312 + run-id: 38182653643 + github-token: ${{ github.token }} + path: ${{ runner.temp }}/wheels + - name: Install the existing portable wheel + run: python -m pip install "$RUNNER_TEMP"/wheels/*.whl numpy==2.5.4 scipy==1.18.1 pytest threadpoolctl + - name: Trace full suite and repeated accounting searches + working-directory: ${{ runner.temp }} + env: + OMP_NUM_THREADS: '1' + OPENBLAS_NUM_THREADS: '1' + PYTHONPATH: '' + run: | + sysctl -n machdep.cpu.brand_string + python "$GITHUB_WORKSPACE/.github/scripts/trace-hef.py" "$RUNNER_TEMP/hef-trace.jsonl" + - uses: actions/upload-artifact@v5 + if: always() + with: + name: hef-trajectory-${{ matrix.sample }} + path: ${{ runner.temp }}/hef-trace.jsonl + if-no-files-found: error From a276bfba2ab5c78b212ce4d97d6ef0d937dd54b6 Mon Sep 17 00:00:00 2001 From: Praharsh Suryadevara Date: Sun, 11 Oct 2026 18:40:48 -0500 Subject: [PATCH 08/11] Stabilize transition-state test inputs after completed macOS traces --- .github/scripts/trace-hef.py | 85 ------------------- .github/workflows/release.yml | 7 +- .github/workflows/trace-hef.yml | 50 ----------- RELEASING.md | 19 +++++ .../tests/test_find_transition_state_hef.py | 40 ++------- .../tests/test_generalized_dimer.py | 2 + 6 files changed, 34 insertions(+), 169 deletions(-) delete mode 100644 .github/scripts/trace-hef.py delete mode 100644 .github/workflows/trace-hef.yml diff --git a/.github/scripts/trace-hef.py b/.github/scripts/trace-hef.py deleted file mode 100644 index e40516c59..000000000 --- a/.github/scripts/trace-hef.py +++ /dev/null @@ -1,85 +0,0 @@ -"""Temporary observer for the installed Intel wheel; no solver settings change.""" -import json -import platform -import sys -from pathlib import Path -from unittest import mock - -import numpy as np -import pytest -import scipy -from threadpoolctl import threadpool_info - -import pele -from pele.transition_states.tests import test_find_transition_state_hef as tests - -output = open(sys.argv[1], 'w') -original_search = tests.FindTransitionState -active_case = None - - -def record(kind, **values): - output.write(json.dumps(dict(kind=kind, case=active_case, **values), - default=lambda value: value.item() if isinstance(value, np.generic) else str(value)) + '\n') - output.flush() - - -class TraceSearch(original_search): - def _getLowestEigenVector(self, coords, iteration, gradient=None): - overlap = super()._getLowestEigenVector(coords, iteration, gradient=gradient) - energy, raw_gradient, hessian = self.pot.pot.getEnergyGradientHessian(coords) - centered = coords.reshape(-1, 3) - coords.reshape(-1, 3).mean(axis=0) - rigid = np.column_stack([ - np.tile(axis, len(centered)) for axis in np.eye(3) - ] + [np.cross(np.tile(axis, (len(centered), 1)), centered).ravel() - for axis in np.eye(3)]) - basis = np.linalg.qr(rigid, mode='complete')[0][:, np.linalg.matrix_rank(rigid):] - spectrum = np.linalg.eigvalsh(basis.T @ hessian @ basis) - curvature = float(self.eigenvec @ hessian @ self.eigenvec) - record('primary', iteration=iteration, npositive=self.npositive, - coords=coords.tolist(), eigenvec=self.eigenvec.tolist(), - energy=energy, raw_rms=float(np.linalg.norm(raw_gradient) / np.sqrt(coords.size)), - primary={key: self.leig_result.get(key) for key in ('success', 'nsteps', 'rms', 'eigenval')}, - primary_curvature=curvature, - primary_residual=float(np.linalg.norm(basis.T @ (hessian @ self.eigenvec - curvature * self.eigenvec))), - physical_spectrum=spectrum.tolist()) - return overlap - - def run(self): - record('start', coords=self.coords.tolist()) - result = super().run() - record('result', success=result.success, energy=result.energy, rms=result.rms, - nsteps=result.nsteps, nfev=result.nfev, counted_nfev=self.pot.nfev, - coords=result.coords.tolist(), eigenval=result.eigenval, message=result.message) - assert result.nfev == self.pot.nfev - return result - - -class TraceAccountingTests: - def pytest_runtest_setup(self, item): - global active_case - if item.cls is tests.TestFindTransitionState_NFEV: - active_case = item.nodeid - tests.FindTransitionState = TraceSearch - - def pytest_runtest_teardown(self, item): - tests.FindTransitionState = original_search - - -record('runtime', python=sys.version, numpy=np.__version__, scipy=scipy.__version__, - platform=platform.platform(), threads=threadpool_info()) -package = Path(pele.__file__).parent -suite_status = pytest.main([str(package), '--ignore=' + str(package / 'gui'), '-q'], - plugins=[TraceAccountingTests()]) -record('suite', status=int(suite_status)) -for trial in range(10): - for variant in ('test', 'test2'): - active_case = f'repeat-{trial}-{variant}' - case = tests.TestFindTransitionState_NFEV(variant) - case.setUp() - with mock.patch.object(tests, 'FindTransitionState', TraceSearch): - try: - getattr(case, variant)() - except AssertionError as error: - record('assertion', message=str(error)) -output.close() diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index fbf651467..ac2003496 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -97,11 +97,14 @@ jobs: strategy: fail-fast: false matrix: - # Temporary: isolate the remaining Intel Python 3.12 failure. - python: [cp312] + python: [cp311, cp312, cp313, cp314] platform: + - os: ubuntu-24.04 + arch: x86_64 - os: macos-15-intel arch: x86_64 + - os: macos-15 + arch: arm64 steps: - uses: actions/setup-python@v6 with: diff --git a/.github/workflows/trace-hef.yml b/.github/workflows/trace-hef.yml deleted file mode 100644 index 8a284e385..000000000 --- a/.github/workflows/trace-hef.yml +++ /dev/null @@ -1,50 +0,0 @@ -name: Trace Intel macOS HEF trajectory - -on: - pull_request: - branches: [main] - paths: [.github/workflows/trace-hef.yml, .github/scripts/trace-hef.py] - workflow_dispatch: - -permissions: - contents: read - actions: read - -jobs: - trace: - runs-on: macos-15-intel - timeout-minutes: 15 - strategy: - fail-fast: false - matrix: - sample: [1, 2] - steps: - - uses: actions/checkout@v6 - with: - persist-credentials: false - - uses: actions/setup-python@v6 - with: - python-version: '3.12.10' - - uses: actions/download-artifact@v6 - with: - name: pele-wheels-macos-15-intel-x86_64-cp312 - run-id: 38182653643 - github-token: ${{ github.token }} - path: ${{ runner.temp }}/wheels - - name: Install the existing portable wheel - run: python -m pip install "$RUNNER_TEMP"/wheels/*.whl numpy==2.5.4 scipy==1.18.1 pytest threadpoolctl - - name: Trace full suite and repeated accounting searches - working-directory: ${{ runner.temp }} - env: - OMP_NUM_THREADS: '1' - OPENBLAS_NUM_THREADS: '1' - PYTHONPATH: '' - run: | - sysctl -n machdep.cpu.brand_string - python "$GITHUB_WORKSPACE/.github/scripts/trace-hef.py" "$RUNNER_TEMP/hef-trace.jsonl" - - uses: actions/upload-artifact@v5 - if: always() - with: - name: hef-trajectory-${{ matrix.sample }} - path: ${{ runner.temp }}/hef-trace.jsonl - if-no-files-found: error diff --git a/RELEASING.md b/RELEASING.md index 337c9188a..e564f8c02 100644 --- a/RELEASING.md +++ b/RELEASING.md @@ -93,11 +93,30 @@ passing on another. The [trace on PR #183](https://github.com/pele-python/pele/a 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. + When validating a deployment: - Check reported solver success, the force stopping criterion and the expected number of unstable physical modes after excluding rigid zero modes. Compare these properties rather than requiring identical final coordinates. +- Give stochastic tests their own seeds, so changing the random draws in one + fixture cannot change the initial direction or iteration count in another test. - Keep force, eigenvector-search and test-comparison tolerances distinct. A larger iteration budget can resolve non-convergence without loosening any tolerance; it does not guarantee convergence for every candidate. diff --git a/pele/transition_states/tests/test_find_transition_state_hef.py b/pele/transition_states/tests/test_find_transition_state_hef.py index caf3edb40..9442e7aa3 100644 --- a/pele/transition_states/tests/test_find_transition_state_hef.py +++ b/pele/transition_states/tests/test_find_transition_state_hef.py @@ -107,44 +107,20 @@ def setUp(self): np.random.seed(0) self.system = LJCluster(18) self.pot = _PotWrapper(self.system.get_potential()) - self.x = self.system.get_random_minimized_configuration(tol=10.0).coords + # Test evaluation accounting with actual refinement in a controlled saddle basin. + path = os.path.join(os.path.dirname(__file__), "lj18_ts.xyz") + with open(path) as source: + coords = read_xyz(source).coords.ravel() + perturbation = np.sin(np.arange(coords.size) + 1.0) + perturbation = self.system.get_orthogonalize_to_zero_eigenvectors()(perturbation, coords) + self.x = coords + 0.01 * perturbation / np.linalg.norm(perturbation) def do_check(self, **kwargs): self.pot.nfev = 0 opt = FindTransitionState(self.x, self.pot, **kwargs) ret = opt.run() self.assertEqual(ret.nfev, self.pot.nfev) - if not ret.success: - import json - import platform - import sys - - _, raw_gradient, hessian = self.pot.pot.getEnergyGradientHessian(ret.coords) - centered = ret.coords.reshape(-1, 3) - ret.coords.reshape(-1, 3).mean(axis=0) - rigid = np.column_stack([ - np.tile(axis, len(centered)) for axis in np.eye(3) - ] + [np.cross(np.tile(axis, (len(centered), 1)), centered).ravel() - for axis in np.eye(3)]) - rank = np.linalg.matrix_rank(rigid) - basis = np.linalg.qr(rigid, mode="complete")[0][:, rank:] - spectrum = np.linalg.eigvalsh(basis.T @ hessian @ basis) - primary_curvature = float(ret.eigenvec @ hessian @ ret.eigenvec) - print("HEF_TRACE " + json.dumps({ - "test": self.id(), "python": sys.version, "numpy": np.__version__, - "platform": platform.platform(), "coords": ret.coords.tolist(), - "energy": float(ret.energy), "rms": float(ret.rms), - "raw_rms": float(np.linalg.norm(raw_gradient) / np.sqrt(ret.coords.size)), - "nfev": int(ret.nfev), "counted_nfev": int(self.pot.nfev), - "nsteps": int(ret.nsteps), "message": ret.message, - "eigenvec": ret.eigenvec.tolist(), "eigenval": float(ret.eigenval), - "primary": {key: opt.leig_result.get(key) - for key in ("success", "nsteps", "rms", "eigenval")}, - "primary_curvature": primary_curvature, - "primary_residual": float(np.linalg.norm( - basis.T @ (hessian @ ret.eigenvec - primary_curvature * ret.eigenvec))), - "physical_spectrum": spectrum.tolist(), - }, default=lambda value: value.item() if isinstance(value, np.generic) else str(value))) - self.assertTrue(ret.success, msg=repr(ret)) + self.assertTrue(ret.success) self.assertGreater(ret.nfev, 0) def test(self): diff --git a/pele/transition_states/tests/test_generalized_dimer.py b/pele/transition_states/tests/test_generalized_dimer.py index 0f67e820d..3195a01df 100644 --- a/pele/transition_states/tests/test_generalized_dimer.py +++ b/pele/transition_states/tests/test_generalized_dimer.py @@ -1,5 +1,6 @@ import unittest import os +import numpy as np from pele.systems import LJCluster from pele.transition_states._generalized_dimer import GeneralizedDimer @@ -8,6 +9,7 @@ class TestGeneralizedDimer(unittest.TestCase): def setUp(self): + np.random.seed(0) self.system = LJCluster(18) self.pot = self.system.get_potential() From 800e5ec5eaca4ea0202039bbaddad911efa0c74b Mon Sep 17 00:00:00 2001 From: Praharsh Suryadevara Date: Sun, 11 Oct 2026 18:41:23 -0500 Subject: [PATCH 09/11] Clarify what independent test seeds guarantee --- RELEASING.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/RELEASING.md b/RELEASING.md index e564f8c02..9d89d19ec 100644 --- a/RELEASING.md +++ b/RELEASING.md @@ -115,8 +115,8 @@ When validating a deployment: - Check reported solver success, the force stopping criterion and the expected number of unstable physical modes after excluding rigid zero modes. Compare these properties rather than requiring identical final coordinates. -- Give stochastic tests their own seeds, so changing the random draws in one - fixture cannot change the initial direction or iteration count in another test. +- Give stochastic tests their own seeds so their starting random state is + independent of draws in earlier fixtures. - Keep force, eigenvector-search and test-comparison tolerances distinct. A larger iteration budget can resolve non-convergence without loosening any tolerance; it does not guarantee convergence for every candidate. From b2ea75df61ec2e8f91dc73e314ca6d1da673daf0 Mon Sep 17 00:00:00 2001 From: Praharsh Suryadevara Date: Sun, 11 Oct 2026 20:06:54 -0500 Subject: [PATCH 10/11] Keep saddle screening linear and preserve valid LBFGS history --- RELEASING.md | 18 +++++++ .../tests/test_double_ended_connect.py | 3 +- pele/systems/tests/test_frozen_lj.py | 5 ++ .../_transition_state_refinement.py | 31 ++++++++--- pele/transition_states/_transverse_walker.py | 11 ++-- .../tests/test_find_transition_state_hef.py | 10 ++-- .../tests/test_transition_state_index.py | 54 ++++++++++++++++--- .../tests/test_transition_state_search.py | 20 +++---- .../tests/test_transverse_walker.py | 43 +++++++++++++++ 9 files changed, 164 insertions(+), 31 deletions(-) diff --git a/RELEASING.md b/RELEASING.md index 9d89d19ec..aa92c959f 100644 --- a/RELEASING.md +++ b/RELEASING.md @@ -110,6 +110,24 @@ 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. +The default transition-state acceptance check is a probabilistic screen, not an +exhaustive index-one certificate. It retains two iterative searches and linear +workspace; no dense Hessian or growing spectral basis is allocated. Each search +uses a private generator that advances between screens. Set `index_seed` to +reproduce its starts, and seed the separate global NumPy generator when the +primary eigenvector search also needs reproducible initialization. Passing the +screen can still miss a negative mode; no per-million failure rate is established. + +For the secondary searches, the default RMS tolerance is `1e-5 / sqrt(ndof)`, +which bounds their gradient norm by `1e-5` independently of dimension. An explicit +`lowestEigenvectorQuenchParams["tol"]` retains its existing RMS meaning. This does +not change the force stopping tolerance or the 500-step default search budget. +Seeded gradient-only quadratics at 262,144 dimensions exposed two missed negative +modes among 16 cases with the old RMS rule; the scaled rule detected all 16. At +4,096 dimensions with weaker curvature, it detected all eight cases where the +old rule missed four. These are model-specific checks, not a statistical +reliability estimate. + When validating a deployment: - Check reported solver success, the force stopping criterion and the expected diff --git a/pele/landscape/tests/test_double_ended_connect.py b/pele/landscape/tests/test_double_ended_connect.py index 793f2687b..5698139dd 100644 --- a/pele/landscape/tests/test_double_ended_connect.py +++ b/pele/landscape/tests/test_double_ended_connect.py @@ -45,7 +45,8 @@ 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()) diff --git a/pele/systems/tests/test_frozen_lj.py b/pele/systems/tests/test_frozen_lj.py index 5a9dcd78e..48715260a 100644 --- a/pele/systems/tests/test_frozen_lj.py +++ b/pele/systems/tests/test_frozen_lj.py @@ -1,5 +1,7 @@ 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 @@ -7,6 +9,8 @@ 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) diff --git a/pele/transition_states/_transition_state_refinement.py b/pele/transition_states/_transition_state_refinement.py index 6bcc8af52..ad25ea70a 100644 --- a/pele/transition_states/_transition_state_refinement.py +++ b/pele/transition_states/_transition_state_refinement.py @@ -49,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 @@ -97,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 @@ -127,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 @@ -225,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 @@ -286,7 +299,7 @@ def run(self): self._compute_gradients(coords) iend = 0 - index_validated = False + index_screened = False for i in range(self.nsteps): iend = i # get the lowest eigenvalue and eigenvector @@ -393,10 +406,10 @@ def run(self): break transverse = self._get_lowest_transverse_eigenvector(coords, grad) if not transverse.success: - res.message.append("could not verify transverse stability") + res.message.append("could not complete transverse stability screen") break if transverse.eigenval > 0: - index_validated = True + index_screened = True break res.message.append("additional negative curvature found") if i + 1 == self.nsteps: @@ -434,8 +447,8 @@ def run(self): ) break - # A successful candidate already has its final eigenpair checked. - if not index_validated: + # 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() @@ -452,7 +465,7 @@ def run(self): self.eigenval, ) - success = index_validated + success = index_screened # check if results make sense if self.eigenval >= 0.0: if self.verbosity > 2: @@ -488,7 +501,7 @@ def run(self): return res def _get_lowest_transverse_eigenvector(self, coords, gradient): - """Check for another unstable physical mode without forming a Hessian.""" + """Screen for another unstable physical mode without forming a Hessian.""" orthogonalize = self.orthogZeroEigs if orthogonalize == 0: orthogonalize = orthogopt @@ -510,6 +523,8 @@ def project(vector, position): 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): @@ -529,7 +544,7 @@ def hessian_product(vector): # 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 = np.random.default_rng(0) + random = self._index_random for _ in range(2): guess = project(random.normal(size=coords.size), coords) if np.linalg.norm(guess) < 1e-12: diff --git a/pele/transition_states/_transverse_walker.py b/pele/transition_states/_transverse_walker.py index ae5b15986..c406b88bf 100644 --- a/pele/transition_states/_transverse_walker.py +++ b/pele/transition_states/_transverse_walker.py @@ -111,12 +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) - # Secant history from another projection can move along the excluded mode. - self.walker.reset() + 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""" diff --git a/pele/transition_states/tests/test_find_transition_state_hef.py b/pele/transition_states/tests/test_find_transition_state_hef.py index 9442e7aa3..14b9cf021 100644 --- a/pele/transition_states/tests/test_find_transition_state_hef.py +++ b/pele/transition_states/tests/test_find_transition_state_hef.py @@ -24,6 +24,8 @@ def getEnergyGradient(self, coords): class TestFindTransitionState(unittest.TestCase): def setUp(self): + self.addCleanup(np.random.set_state, np.random.get_state()) + np.random.seed(0) self.system = LJCluster(18) self.pot = PotWrapper(self.system.get_potential()) @@ -31,7 +33,7 @@ def make_dimer(self, x, **kwargs): return FindTransitionState( x, self.pot, - orthogZeroEigs=self.system.get_orthogonalize_to_zero_eigenvectors(), + index_seed=0, orthogZeroEigs=self.system.get_orthogonalize_to_zero_eigenvectors(), **kwargs # verbosity=10 ) @@ -52,7 +54,7 @@ def test1(self): self.assertEqual(res.nfev, self.pot.nfev) self.assertGreater(res.nfev, 1) - # Check physical curvature independently of the iterative index validation. + # Check physical curvature independently of the iterative index screen. _, _, hessian = self.pot.pot.getEnergyGradientHessian(res.coords) centered = res.coords.reshape(-1, 3) centered = centered - centered.mean(axis=0) @@ -94,7 +96,7 @@ def make_dimer(self, x, **kwargs): return FindTransitionState( x, self.pot, - orthogZeroEigs=self.system.get_orthogonalize_to_zero_eigenvectors(), + index_seed=0, orthogZeroEigs=self.system.get_orthogonalize_to_zero_eigenvectors(), invert_gradient=True, **kwargs, ) @@ -117,7 +119,7 @@ def setUp(self): def do_check(self, **kwargs): self.pot.nfev = 0 - opt = FindTransitionState(self.x, self.pot, **kwargs) + opt = FindTransitionState(self.x, self.pot, index_seed=0, **kwargs) ret = opt.run() self.assertEqual(ret.nfev, self.pot.nfev) self.assertTrue(ret.success) diff --git a/pele/transition_states/tests/test_transition_state_index.py b/pele/transition_states/tests/test_transition_state_index.py index 9e51e037b..b78d2734f 100644 --- a/pele/transition_states/tests/test_transition_state_index.py +++ b/pele/transition_states/tests/test_transition_state_index.py @@ -31,17 +31,59 @@ def test_rejects_rotated_index_two_with_stable_initial_rr_direction(self): ]) ret = FindTransitionState( np.zeros(3), pot, eigenvec0=np.array([1.0, 0.0, 0.0]), - orthogZeroEigs=None, nsteps=1, verbosity=0, + index_seed=0, orthogZeroEigs=None, nsteps=1, verbosity=0, ).run() self.assertFalse(ret.success) + def test_random_start_finds_index_two_mode_missed_by_fixed_guesses(self): + # Spectrum -2, -1, +1, +1; both old fixed probes miss the -1 mode. + pot = _QuadraticPotential([ + [-2.0, 0.0, 0.0, 0.0], + [0.0, -0.07552173356870506, -0.3733085508688484, 0.9246281380138601], + [0.0, -0.3733085508688484, 0.870426352344002, 0.3209340913541529], + [0.0, 0.9246281380138601, 0.3209340913541529, 0.20509538122470305], + ]) + np.testing.assert_allclose(np.linalg.eigvalsh(pot.hessian), + [-2.0, -1.0, 1.0, 1.0], atol=1e-14) + ret = FindTransitionState( + np.zeros(4), pot, eigenvec0=np.array([1.0, 0.0, 0.0, 0.0]), + index_seed=1, orthogZeroEigs=None, nsteps=1, verbosity=0, + ).run() + self.assertEqual(ret.nfev, pot.nfev) + self.assertFalse(ret.success) + + def test_weak_high_dimensional_mode_is_not_lost_to_rms_stopping(self): + dimension = 4096 + hidden = np.random.default_rng(314159).normal(size=dimension) + hidden[0] = 0.0 + hidden /= np.linalg.norm(hidden) + + class RankTwoPotential(BasePotential): + def getEnergyGradient(self, coords): + # Spectrum -2, -0.01, +0.01... without allocating a Hessian. + gradient = 0.01 * coords - 0.02 * np.dot(hidden, coords) * hidden + gradient[0] -= 2.01 * coords[0] + return 0.5 * np.dot(coords, gradient), gradient + + primary = np.zeros(dimension) + primary[0] = 1.0 + search = FindTransitionState( + np.zeros(dimension), RankTwoPotential(), eigenvec0=primary, + index_seed=0, orthogZeroEigs=None, verbosity=0, + ) + result = search._get_lowest_transverse_eigenvector( + np.zeros(dimension), np.zeros(dimension), + ) + self.assertTrue(result.success) + self.assertLess(result.eigenval, 0.0) + def test_accepts_index_one_when_zero_modes_leave_no_transverse_direction(self): def remove_zero_mode(vector, coords): return np.array([vector[0], 0.0]) ret = FindTransitionState( np.zeros(2), _QuadraticPotential([[-2.0, 0.0], [0.0, 0.0]]), - eigenvec0=np.array([1.0, 0.0]), orthogZeroEigs=remove_zero_mode, + eigenvec0=np.array([1.0, 0.0]), index_seed=0, orthogZeroEigs=remove_zero_mode, nsteps=1, verbosity=0, ).run() self.assertTrue(ret.success) @@ -50,7 +92,7 @@ def remove_zero_mode(vector, coords): def test_accepts_one_dimensional_index_one_stationary_point(self): ret = FindTransitionState( np.zeros(1), _QuadraticPotential([[-2.0]]), - eigenvec0=np.ones(1), orthogZeroEigs=None, nsteps=1, verbosity=0, + eigenvec0=np.ones(1), index_seed=0, orthogZeroEigs=None, nsteps=1, verbosity=0, ).run() self.assertTrue(ret.success) self.assertAlmostEqual(ret.eigenval, -2.0) @@ -64,7 +106,7 @@ def test_accepts_stiff_index_one_stationary_point(self): primary[0] = 1.0 ret = FindTransitionState( np.zeros(curvatures.size), pot, eigenvec0=primary, - orthogZeroEigs=None, nsteps=1, verbosity=0, + index_seed=0, orthogZeroEigs=None, nsteps=1, verbosity=0, ).run() self.assertEqual(ret.nfev, pot.nfev) self.assertTrue(ret.success, msg=repr(ret)) @@ -72,7 +114,7 @@ def test_accepts_stiff_index_one_stationary_point(self): def test_success_requires_a_converged_returned_primary_eigenpair(self): search = FindTransitionState( np.zeros(1), _QuadraticPotential([[-2.0]]), - eigenvec0=np.ones(1), orthogZeroEigs=None, nsteps=1, verbosity=0, + eigenvec0=np.ones(1), index_seed=0, orthogZeroEigs=None, nsteps=1, verbosity=0, ) find_primary = search._get_lowest_eigenvector_RR calls = 0 @@ -95,7 +137,7 @@ def test_recovery_failure_returns_energy_and_gradient_at_returned_coords(self): pot = _QuadraticPotential(np.diag([-2.0, -1.0, 1.0])) search = FindTransitionState( np.zeros(3), pot, eigenvec0=np.array([1.0, 0.0, 0.0]), - orthogZeroEigs=None, demand_initial_negative_vec=True, + index_seed=0, orthogZeroEigs=None, demand_initial_negative_vec=True, nsteps=2, verbosity=0, ) search.npositive_max = 0 diff --git a/pele/transition_states/tests/test_transition_state_search.py b/pele/transition_states/tests/test_transition_state_search.py index e67af21fa..bc78f4d20 100644 --- a/pele/transition_states/tests/test_transition_state_search.py +++ b/pele/transition_states/tests/test_transition_state_search.py @@ -77,6 +77,8 @@ def print_event(coords=None, **kwargs): class TestFindTransitionStateSimplePot(unittest.TestCase): def setUp(self): + self.addCleanup(np.random.set_state, np.random.get_state()) + np.random.seed(0) self.pot = SimpleTSPot() self.x0 = np.array([0.1, 0.1]) self.xts = np.zeros(self.x0.size) @@ -87,7 +89,7 @@ def test1(self): opt = FindTransitionState( self.x0, self.pot, - orthogZeroEigs=None, + index_seed=0, orthogZeroEigs=None, # iprint=1, # verbosity=10, event=print_event, # tol=1e-3, @@ -101,7 +103,7 @@ def test1(self): self.assertEqual(ret.nfev + 1, self.pot.nfev) def test_wrapper(self): - ret = findTransitionState(self.x0, self.pot, orthogZeroEigs=None) + ret = findTransitionState(self.x0, self.pot, index_seed=0, orthogZeroEigs=None) self.assertTrue(ret.success) assert_arrays_almost_equal(self, ret.coords, self.xts, places=3) @@ -114,7 +116,7 @@ def event(**kwargs): opt = FindTransitionState( self.x0, self.pot, - orthogZeroEigs=None, + index_seed=0, orthogZeroEigs=None, tangentSpaceQuenchParams=dict(maxstep=1.0), event=event, ) @@ -128,7 +130,7 @@ def test_from_near_minimum(self): opt = FindTransitionState( x0, self.pot, - orthogZeroEigs=None, + index_seed=0, orthogZeroEigs=None, iprint=1, verbosity=10, # event=print_event, # tol=1e-3, @@ -147,7 +149,7 @@ def test_from_near_minimum_demand_negative_eigenvalue(self): opt = FindTransitionState( x0, self.pot, - orthogZeroEigs=None, + index_seed=0, orthogZeroEigs=None, demand_initial_negative_vec=True, iprint=1, verbosity=10, # event=print_event, @@ -168,7 +170,7 @@ def test1(self): opt = FindTransitionState( x0, pot, - orthogZeroEigs=None, + index_seed=0, orthogZeroEigs=None, iprint=1, verbosity=10, # event=print_event, hessian_diagonalization=True, @@ -198,7 +200,7 @@ class TestTransitionStateIndex(unittest.TestCase): def test_rejects_index_two_stationary_point(self): ret = FindTransitionState( np.zeros(3), IndexTwoPotential(), eigenvec0=np.array([1., 0., 0.]), - orthogZeroEigs=None, nsteps=1, verbosity=0, + index_seed=0, orthogZeroEigs=None, nsteps=1, verbosity=0, ).run() self.assertFalse(ret.success) @@ -209,14 +211,14 @@ def test_rejects_recorded_lj13_index_two_saddle(self): values, vectors = np.linalg.eigh(pot.getHessian(coords)) self.assertLess(values[1], -1.0) ret = FindTransitionState( - coords, pot, eigenvec0=vectors[:, 0], nsteps=1, verbosity=0, + coords, pot, eigenvec0=vectors[:, 0], index_seed=0, nsteps=1, verbosity=0, ).run() self.assertFalse(ret.success) def test_escapes_extra_unstable_mode(self): ret = FindTransitionState( np.zeros(3), RecoverableIndexTwoPotential(), - eigenvec0=np.array([1., 0., 0.]), orthogZeroEigs=None, verbosity=0, + eigenvec0=np.array([1., 0., 0.]), index_seed=0, orthogZeroEigs=None, verbosity=0, ).run() self.assertTrue(ret.success) self.assertAlmostEqual(abs(ret.coords[1]), 1.0, delta=1e-3) diff --git a/pele/transition_states/tests/test_transverse_walker.py b/pele/transition_states/tests/test_transverse_walker.py index 6a5287360..ba6c72360 100644 --- a/pele/transition_states/tests/test_transverse_walker.py +++ b/pele/transition_states/tests/test_transverse_walker.py @@ -27,6 +27,49 @@ def test(self): class TestTransverseWalkerHistory(unittest.TestCase): + def test_unchanged_plane_keeps_the_same_trajectory(self): + class QuadraticPotential: + def getEnergyGradient(self, coords): + gradient = np.array([1.0, 10.0, 100.0]) * coords + return 0.5 * np.dot(coords, gradient), gradient + + pot = QuadraticPotential() + coords = np.array([1.0, -0.7, 0.4]) + axis = np.array([0.0, 0.0, 1.0]) + baseline = _TransverseWalker(coords, pot, axis, tol=1e-12) + expected = baseline.run(4).coords.copy() + for sign in (1.0, -1.0): + with self.subTest(sign=sign): + opt = _TransverseWalker(coords, pot, axis, tol=1e-12) + current = opt.run(3).coords.copy() + energy, gradient = pot.getEnergyGradient(current) + opt.update_eigenvec(sign * axis, -1.0) + opt.update_coords(current, energy, gradient) + np.testing.assert_allclose(opt.run(1).coords, expected, + rtol=0.0, atol=1e-14) + + def test_rotation_restores_configured_initial_scale(self): + class QuadraticPotential: + def getEnergyGradient(self, coords): + gradient = np.array([1.0, 2.0, 3.0]) * coords + return 0.5 * np.dot(coords, gradient), gradient + + pot = QuadraticPotential() + opt = _TransverseWalker( + np.array([1.0, -0.7, 0.4]), pot, np.array([0.0, 0.0, 1.0]), + tol=1e-12, H0=0.7, maxstep=10.0, + ) + opt.run(3) + axis = np.array([1.0, 0.0, 0.0]) + coords = np.array([-0.5, 1.0, 1.0]) + energy, gradient = pot.getEnergyGradient(coords) + opt.update_eigenvec(axis, -1.0) + opt.update_coords(coords, energy, gradient) + fresh = _TransverseWalker(coords, pot, axis, tol=1e-12, + H0=0.7, maxstep=10.0) + np.testing.assert_allclose(opt.run(1).coords, fresh.run(1).coords, + rtol=0.0, atol=1e-14) + def test_motion_stays_in_updated_plane(self): class QuadraticPotential: hessian = np.array([[3.0, 1.0, 0.5], From f97d66675a8253cc439406f8befa97741b2af9ca Mon Sep 17 00:00:00 2001 From: Praharsh Suryadevara Date: Sun, 11 Oct 2026 20:22:08 -0500 Subject: [PATCH 11/11] Clean up RELEASING.md by removing redundant details Removed detailed explanations about transition-state acceptance checks, secondary searches, and validation processes. --- RELEASING.md | 34 ---------------------------------- 1 file changed, 34 deletions(-) diff --git a/RELEASING.md b/RELEASING.md index aa92c959f..99c5113af 100644 --- a/RELEASING.md +++ b/RELEASING.md @@ -110,40 +110,6 @@ 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. -The default transition-state acceptance check is a probabilistic screen, not an -exhaustive index-one certificate. It retains two iterative searches and linear -workspace; no dense Hessian or growing spectral basis is allocated. Each search -uses a private generator that advances between screens. Set `index_seed` to -reproduce its starts, and seed the separate global NumPy generator when the -primary eigenvector search also needs reproducible initialization. Passing the -screen can still miss a negative mode; no per-million failure rate is established. - -For the secondary searches, the default RMS tolerance is `1e-5 / sqrt(ndof)`, -which bounds their gradient norm by `1e-5` independently of dimension. An explicit -`lowestEigenvectorQuenchParams["tol"]` retains its existing RMS meaning. This does -not change the force stopping tolerance or the 500-step default search budget. -Seeded gradient-only quadratics at 262,144 dimensions exposed two missed negative -modes among 16 cases with the old RMS rule; the scaled rule detected all 16. At -4,096 dimensions with weaker curvature, it detected all eight cases where the -old rule missed four. These are model-specific checks, not a statistical -reliability estimate. - -When validating a deployment: - -- Check reported solver success, the force stopping criterion and the expected - number of unstable physical modes after excluding rigid zero modes. Compare - these properties rather than requiring identical final coordinates. -- Give stochastic tests their own seeds so their starting random state is - independent of draws in earlier fixtures. -- Keep force, eigenvector-search and test-comparison tolerances distinct. A - larger iteration budget can resolve non-convergence without loosening any - tolerance; it does not guarantee convergence for every candidate. -- Record the Python, NumPy and SciPy versions, architecture, compiler, math - backend, thread settings, random state, coordinates and optimizer iteration - counts. Save a failing candidate and replay it at fixed coordinates to - distinguish arithmetic differences in the path from a convergence or - correctness defect. - ## Configure publishing accounts Create accounts on [PyPI](https://pypi.org/) and [TestPyPI](https://test.pypi.org/).