BUG?: Anomalous convergence for benchmark system

[HaoZeke]

I'd like to report a reproduction error. The first system of the SI described in [1], namely singlet_000 does not seem to converge to the nearest first order saddle point. The attached slides and reproduction document details everything needed to reproduce the error. Basically it seems like Sella "misses" the nearest saddle and instead converges to a trivial rotation.

For clarity all the details are in this issue, but they can also be validated offline with the following files:

sella_anomaly_slides.pdf

sella_repro.pdf

Anomaly and background

Baseline setup

To ensure maximal reproduciblity, we assume a pixi installation ¹, which may be obtained via:

# Linux & macOS
curl -fsSL https://pixi.sh/install.sh | bash

Now we can setup the directory structure, and get the supplementary information from here (with a web-browser, also attached):

• https://pubs.acs.org/doi/suppl/10.1021/acs.jctc.2c00395/suppl_file/ct2c00395_si_001.zip

    mkdir -p sella_reproducer
    cd sella_reproducer
    # Extract this, can be done manually as well
    unzip ct2c00395_si_001

Now we can setup the environment.

    cd sella_reproducer
    pixi init
    pixi add ase nwchem
    pixi add --pypi sella

Reproduction

The extracted script from the supplementary information has all the details needed, so it can be run almost as-is. The primary change is to ensure that the NWCHEM location is taken from the environment. So apply this patch.

    --- calc_orig.py  2024-11-09 16:37:11.341264575 +0000
    +++ calc.py  2024-11-09 16:37:13.914580395 +0000
    @@ -16,10 +16,10 @@
     run_path = 'runs'
    
     # Path to local scratch directory for NWChem scratch files
    -scratch_path = '/scratch/ehermes'
    +scratch_path = '/tmp'
    
     # Full path to your nwchem executable
    -nwchem_path = '/home/ehermes/build/nwchem/bin/LINUX64/nwchem'
    +nwchem_path = os.environ["NWCHEM_COMMAND"]
    
     # Memory for NWChem to allocate. You probably don't need to change this.
     memory = '2 gb'

With the standard approach, that is:

    cd sella_reproducer/sella_si
    # save the patch into patch_os.diff (also attached)
    patch < patch_os.diff

Now we’re ready to run this.

    cd sella_reproducer/sella_si
    pixi shell
    export NWCHEM_COMMAND=$(which nwchem)
    python calc.py

This should take around four minutes, 44 steps, and 108 frames in the
trajectory. By the end of this, the figures in the slides can be extracted by
visualizing the generated trajectory with ASE.

    cd sella_reproducer/sella_si
    pixi shell
    ase gui runs/singlets/000/singlets_000.traj
    # energy values are also present

The Hessian confirmation checks use the standard ase.vibrations module, and
are available on request.

Footnotes

¹ If you already have a working sella setup you may skip this and use that instead.

1. Hermes, E. D., Sargsyan, K., Najm, H. N. & Zádor, J. Sella, an Open-Source Automation-Friendly Molecular Saddle Point Optimizer. J. Chem. Theory Comput. 18, 6974–6988 (2022).

[juditzador]

Thank you for the detailed submission. We will take a look at this.

[levineds]

Update after rerunning the supplied singlets/000.xyz case directly:

I reproduced the reported behavior on current master (475f02c) using Python 3.12, ASE 3.29.0, and conda-forge NWChem 7.3.1, with the original HF/3-21G, internal=True, fmax=0.01, and 200-step settings. Sella converged in 48 optimizer steps and 104 force evaluations to fmax = 0.00699 eV/Å. The energy change was -1.32663 eV, matching the -1.32 eV shown in the attached slides. The hydrogen transfer also occurs along the trajectory.

An analytic NWChem frequency calculation at the final geometry gives one non-rigid imaginary frequency, 230.79i cm^-1, plus the expected six translational/rotational zero modes. The unstable eigenvector has a 0.970 normalized overlap with rigid rotation of the newly formed methyl group about its adjacent C-C bond (0.998 when restricted to the three methyl hydrogens). This confirms that the final stationary point is the same kind of rotation-dominated first-order saddle described in the report.

Therefore, the reported outcome remains reproducible on current master; the subsequent rigid-body-mode and geodesic-integration fixes did not eliminate it. I would currently describe this as a mode-selection or robustness limitation rather than a demonstrated correctness bug: without a supplied v0, Sella does not guarantee the geometrically nearest saddle or the same saddle selected by another optimizer. The most useful next test is to seed v0 along the intended H-exchange mode. If Sella then loses that mode and returns to the rotational saddle, that would point much more directly to a mode-following defect.

[HaoZeke]

I'm not sure I agree @levineds , the issue is more that the trajectory taken by sella includes going near a saddle.

Basically if it just reached the rotational saddle without going near the proximal one it's ok to say it's because it's not guaranteed but it's weirder to deterministically fail to converge to a saddle it passes along the optimization trajectory.

IMO it suggests there's something off about the effective surface sella is using .

[levineds]

In our current-master rerun, however, the H-transfer midpoint still has fmax = 3.14 eV/Å, and the frame closest to the reported eOn saddle energy has fmax = 1.89 eV/Å. Thus, the available trajectory shows a similar bonding motif, but does not yet establish that Sella passes close to the stationary point.

Could you share the final eOn saddle geometry and, ideally, its unstable mode or trajectory? We can then quantify the closest approach and restart Sella nearby using Cartesian and internal coordinates, with and without a mode-aligned v0. If only the internal-coordinate calculation systematically leaves that saddle, that would provide strong evidence of a genuine mode-following or coordinate-model defect.