Pith. sign in

REVIEW 4 major objections 6 minor 7 references

Gaussian process force uncertainty makes instanton path refinement cost independent of the number of beads, while tunneling rates stay within 20% of exact.

Reviewed by Pith at T0; open to challenge. T0 means a machine referee read the full paper against a public rubric. the ladder, T0–T4 →

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

T0 review · deepseek-v4-flash

2026-08-02 22:22 UTC pith:FN6IILUZ

load-bearing objection Useful acceleration for instanton rates with a sensible selective-Hessian idea, but the bead-independence claim rests on an unvalidated uncertainty estimate and the manuscript shows clear signs of being unfinished. the 4 major comments →

arxiv 2602.16962 v2 pith:FN6IILUZ submitted 2026-02-18 physics.chem-ph

Towards Efficient Instanton Rate Calculations using Machine Learning Surrogates

Chenghao Zhang , Amke Nimmrich , Axel Gomez , Munira Khalil , Niranjan Govind This is my paper
classification physics.chem-ph
keywords instanton theoryring polymerGaussian process regressiontunneling ratestunneling splittingline integral string methodselective Hessian trainingproton transfer
verification ladder T0 review T1 audit T2 compute T3 formal T4 reserved

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

This paper is trying to establish that Gaussian process regression can remove the traditional trade-off between the number of beads in a ring polymer instanton calculation and the cost of converging the tunneling path. Using the GPR posterior force variance as a stopping criterion, the authors show that increasing the bead count from 10 to 80 requires only about 8 additional force evaluations, rather than hundreds. They further show that a selective Hessian training strategy, which treats flexible proton-coupled modes differently from rigid modes, cuts force evaluations by 40–62.5% while keeping instanton rates within 20% of the rigorous DFT result. If this works broadly, it makes high-accuracy tunneling calculations much more affordable for molecular systems.

Core claim

On its own terms, the paper claims that the GPR-enhanced line integral string method converges instanton paths with a number of force evaluations that is effectively independent of the number of beads: malonaldehyde needs 64 force evaluations at N=10 beads and only 8 additional at N=80, while Z-3-aminopropenal needs 70 plus 10 additional. The convergence criterion is the GPR-predicted force variance rather than a fixed bead count. For rate calculations, selective Hessian training—using only 3 beads for rigid modes and 10–20 beads for flexible modes—reduces force evaluations by 40–62.5% while keeping predicted rates within 20% of exact instanton rates, and often within 6%. GPU-accelerated Bla

What carries the argument

The central object is a Gaussian process regression surrogate of the potential energy surface that includes potential, gradient, and Hessian observations, and whose posterior force variance (the trace of the force covariance in Cartesian coordinates) serves as a convergence check for path optimization. The line integral string method provides the path optimizer, while BBMM reduces GPR hyperparameter training from cubic to quadratic complexity, with additional GPU acceleration. Selective Hessian training splits internal modes into an active subspace of flexible proton-coupled modes (modeled with GPR) and a null subspace of rigid modes (modeled with linear regression), enabling large reduction

Load-bearing premise

The central claim collapses if the Gaussian process uncertainty estimate is not a trustworthy measure of how wrong its force predictions really are, and the paper never compares predicted uncertainties against actual force errors.

What would settle it

During a typical malonaldehyde instanton run, compute the GPR-predicted force variance at held-out beads together with the actual error |F_ML - F_DFT|. If the ratio of predicted to actual error is far from 1, or if the required number of force evaluations changes systematically with the initial bead count, the stopping criterion is not calibrated and the bead-independence claim fails.

Watch this falsifier. Get emailed when new claim-graph text bears on it.

If this is right

  • Bead count can be chosen purely for discretization accuracy, since refining from N=10 to N=80 costs only about 8–10 extra electronic-structure force evaluations.
  • Selective Hessian training achieves tunneling rates within 20% of the rigorous instanton rate while cutting force evaluations by 40–62.5%, with errors often near 6%.
  • The accuracy of surrogate-predicted rates improves when the underlying DFT calculations use tighter convergence thresholds and finer grids, indicating that electronic-structure noise is a key error source.
  • Cubic spline interpolation is a simple and accurate alternative for Hessian fitting when noise is low, but its error grows at higher temperatures when the instanton path shortens and interpolation points cluster.
  • The same GPR-LI-String machinery locates tunneling-splitting instanton paths in formic acid dimer and malonaldehyde with roughly 100 potential/force evaluations, yielding splittings in reasonable agreement with experiment and high-level theory.
  • The paper flags that the Z-matrix descriptor is not permutation invariant, so tunneling-splitting calculations require a manual symmetry operation; a permutation-invariant descriptor would remove this step.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • If the posterior force variance is truly calibrated, the same bead-independence should extend to other chain-of-states methods, to larger molecules, and to potentials with stronger anharmonicity; a natural test is to hold the bead count fixed and verify that the stopping criterion does not depend on the initial number of beads.
  • The active/null mode split is selected manually with no reproducible cutoff; automating the split by ranking each mode's coupling to the path tangent would make the 40–62.5% cost savings reproducible without user judgment.
  • The observed dominance of electronic-structure noise suggests that combining these surrogate instanton paths with transfer-learned or higher-level corrections to the potential could push rate errors well below the current 20% envelope.
  • The failure of cubic splines under high noise suggests a hybrid scheme that uses GPR where local uncertainty is high and spline interpolation where data are clean could combine robustness with low cost.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit.

Referee Report

4 major / 6 minor

Summary. The paper develops a Gaussian process regression (GPR) enhanced Line Integral String (LI-String) method for ring-polymer instanton calculations of tunneling rates and tunneling splittings. The central methodological claims are: (i) using the GPR posterior force uncertainty as a stopping criterion makes the number of ab initio force evaluations effectively independent of the number of beads used to discretize the instanton path (Section 3.1 and Fig. 1); (ii) GPU-accelerated Blackbox Matrix-Matrix Multiplication (BBMM) reduces GPR hyperparameter training time by roughly an order of magnitude (Section 3.2 and Fig. 2); (iii) a selective Hessian training strategy, which separates flexible and rigid internal modes, reduces Hessian evaluation cost by 40--62.5% while keeping instanton rates within 20% of direct DFT reference rates (Section 3.3 and Tables 1--2). The method is applied to malonaldehyde and Z-3-aminopropenal for rates, and to formic acid dimer and malonaldehyde for tunneling splittings (Tables 3--4). The paper also compares cubic-spline Hessian interpolation with GPR Hessian surrogates, finding spline interpolation less robust at higher temperatures and for noisier data.

Significance. If the claims are substantiated, the paper would provide a practically useful acceleration of on-the-fly instanton rate and tunneling splitting calculations, which are currently expensive because each bead requires forces and Hessians. The rate benchmarks in Tables 1 and 2 report explicit percent errors against independently computed DFT reference rates, which is a strength: the surrogate rates are not evaluated against quantities used to fit the surrogate. The reported GPU speedups and selective-Hessian cost reductions are also concrete and potentially valuable. However, the central bead-independence claim rests on a single run per system and on treating the GPR posterior force variance as a calibrated stopping signal without any comparison between predicted and actual force errors. Several numerical details needed for reproducibility (GPR convergence thresholds, flexible/rigid mode assignment criteria) are not stated. The 'within 20%' accuracy claim is also not uniformly supported by the tables. The significance is therefore conditional on calibration and reproducibility checks.

major comments (4)
  1. [Section 3.1, Fig. 1, Eq. (12)] The headline claim that the number of force evaluations becomes effectively independent of bead count is based on one run per system with no error bars or repeated-initialization statistics. The GPR-enhanced path search depends on the initial training set, GP hyperparameters, and the stopping threshold; a single trajectory cannot establish 'effectively independent.' The GPR force-uncertainty stopping threshold is not stated anywhere, and Eq. (12) only gives the general LI-String convergence criteria. Please report the actual threshold used for the posterior force variance, and provide statistics over multiple independent runs with different initial training configurations.
  2. [Appendix B, Eqs. (15)--(17), Section 3.1] The load-bearing assumption is that the GPR posterior force variance, Eq. (17), is a faithful estimate of the true force error. The paper never compares predicted uncertainties with actual errors (e.g., |F_GPR - F_DFT| at converged beads). If the GP hyperparameters or the active/null mode decomposition are misspecified, the variance can be overconfident (premature stopping on an unconverged path) or underconfident (destroying the claimed cost reduction). The flexible/rigid mode assignment is also selected manually, with no reproducible cutoff or algorithm; the numbers of flexible modes vary (6, 8, 12, 13) across systems and temperatures. A calibration plot and a reproducible mode-selection criterion are needed before the cost claims can be considered supported.
  3. [Tables 1 and 2, Section 4.1] The abstract and conclusions state that tunneling rates are predicted 'within 20%' of exact values, but several entries in Tables 1 and 2 exceed 20%: Table 1, T=275 K, standard, case (a), N=40 gives 20.7%; Table 1, T=275 K, tight, case (b), N=40 gives 21.6%; Table 2, T=250 K, standard, case (b), N=320 gives 20.0%; Table 2, T=160 K, standard, case (b), N=40 gives 21.6%. Please either soften the claim to reflect the actual spread, define the statistical meaning of 'within 20%,' or explain why these specific entries are consistent with the stated accuracy claim.
  4. [Section 4.1, Tables 1--2] The reported cost reductions ('reduces the number of force evaluations by 44%, 52%, 62.5%') are not defined in terms of a quantitative cost model. For example, in Table 1 at T=275 K, case (b) uses 10 Hessian data points along 8 flexible modes and 3 Hessian data points along 13 rigid modes, while case (a) uses 10 full Hessian data points. It is unclear why adding partial Hessian evaluations (13 total vs. 10 full) reduces the number of force evaluations by 44%. If the reduction is due to the lower cost of partial Hessians for rigid modes, the cost model and the actual number of electronic-structure force evaluations should be stated explicitly.
minor comments (6)
  1. [Fig. 2 caption] The caption contains the placeholder text 'NEED to update this plot for the arxiv submission.' This must be removed before publication.
  2. [General] The 'within 20%' claim in the abstract should be qualified by the fact that some Table 1 and Table 2 entries exceed 20%, and the text should state whether this refers to a typical, maximum, or average error.
  3. [Tables 3 and 4] The text says the instanton path can be represented with 20 beads for formic acid dimer, but the tunneling-splitting tables list bead numbers 320, 640, and 1280. Please clarify the relationship between the bead count used for path optimization and the bead count used for the splitting evaluation.
  4. [Section 4.2.2] For malonaldehyde, the computed splitting of 61.43 cm-1 (B3LYP/def2-TZVPP) is a factor of ~2.8 larger than the experimental value of 21.583 cm-1. Calling this 'reasonable agreement' is generous; the discussion should more explicitly quantify the deviation and its sensitivity to the DFT barrier.
  5. [General] There are several typographical errors, e.g., 'implmentation' at the end of Section 2, 'togehter' in Section 4.1.2, and 'Z3-aminopropenal' in the Table 2 header. A careful proofread is needed.
  6. [Reproducibility] No data or code availability statement is provided. Given the method's reliance on GP hyperparameters, convergence thresholds, and mode-assignment choices, a reproducibility statement or a documented implementation would substantially strengthen the paper.

Circularity Check

0 steps flagged

No significant circularity: surrogate rates are benchmarked against independent DFT instanton rates; bead-independence is a stated design choice, not a hidden fit.

full rationale

The central rate predictions in Tables 1 and 2 are not circular: surrogate rates are compared against rigorous DFT instanton rates computed on the same potential, not against quantities used to fit the GPR model, and tunneling splittings are compared with experiment and independent fitted/neural-network PES results. The self-citations (Ref. 40 for the adaptive regression strategy and constrained-dynamics temperature, Refs. 41/42 for LI-String) are ordinary methodological citations and are not used to forbid alternatives or to prove a uniqueness claim; the key algorithm components (GPR, BBMM, selective Hessian training) are attributed to external references and are concretely benchmarked here. The bead-independence claim is a consequence of using the GPR posterior force variance (Eq. 17) as the stopping criterion, and that criterion is a known design choice rather than an output fitted to the reported force-evaluation counts; the empirical Fig. 1 observation could in principle have failed if GPR uncertainty at newly inserted beads had remained high. The main weakness is therefore a validity/reproducibility gap, not circularity: the paper does not calibrate Eq. 17 against actual GPR-vs-DFT force errors and does not state the stopping threshold, so the 8-10 additional evaluation figures are not independently auditable. The 'NEED to update this plot for the arxiv submission' note in the Fig. 2 caption is an incompleteness artifact and does not affect the circularity assessment. Overall, the derivation chain is self-contained and externally benchmarked; no step reduces by construction to its own input.

Axiom & Free-Parameter Ledger

5 free parameters · 5 axioms · 0 invented entities

The paper introduces no new physical entities. Its central claims rest on computational surrogates (GPR, linear regression, spline interpolation) and on the domain assumptions of instanton theory and DFT accuracy. The main unquantified burden is the manual selection of flexible/rigid modes and the unvalidated use of GPR uncertainty as a convergence proxy.

free parameters (5)
  • GPR hyperparameters (kernel lengthscales, likelihood noise) = Not reported per system; optimized by minimizing negative log marginal likelihood (Eq. 2)
    Central to surrogate accuracy; fitted automatically but the manuscript does not provide the values or kernel details needed for reproduction.
  • Number of rigid-mode Hessian beads (3) and flexible-mode Hessian beads (10/20) = 3 rigid; 10-20 flexible depending on system and temperature
    Manual choices in the selective Hessian strategy; directly control the reported 40-62.5% cost reduction.
  • Flexible/rigid mode assignment (number of flexible modes: 6, 8, 12, 13) = Varies by system and temperature
    The classification of internal modes into active/null subspace is load-bearing for the adaptive regression and selective Hessian methods, but no reproducible cutoff is given.
  • Force-minimum fitting cutoffs fc and Vc = fc = 1e-3 a.u., Vc = 0.01 eV
    Hand-chosen constants used to replace GPR force predictions near reactant/product minima.
  • GPR convergence thresholds in Eq. 12 = Not reported
    The reported force-evaluation counts depend on these thresholds, but the manuscript does not state their numerical values.
axioms (5)
  • domain assumption Ring-polymer instanton theory with harmonic fluctuations (Eq. 1) provides a valid approximation for the tunneling rates and splittings studied.
    The entire paper builds on this semiclassical approximation; it is cited from the literature rather than re-derived.
  • domain assumption DFT potential energy surfaces (B3LYP, wB97X, PBE0) are accurate enough that surrogate errors can be benchmarked against direct DFT calculations and conclusions about tunneling are meaningful.
    The paper compares GPR rates to DFT rates, but the physics conclusions (e.g., comparison to experiment) inherit DFT limitations, which the authors acknowledge for malonaldehyde splittings.
  • domain assumption The GPR posterior force variance is a calibrated measure of true force error and can be used as a stopping criterion.
    The bead-independence claim depends on this assumption, but the paper does not validate predicted uncertainties against actual force errors.
  • domain assumption Rigid modes can be accurately modeled with linear regression in the null subspace, while only flexible modes require GPR.
    The selective Hessian strategy assumes this separation is valid and that the manual mode assignment is correct for each temperature.
  • domain assumption The Z-matrix descriptor, after symmetry correction, is adequate for constructing the GPR model.
    The paper notes Z-matrix descriptors are not permutation-invariant and applies a symmetry operation to fix the resulting asymmetry; this is a workaround rather than an invariant descriptor.

pith-pipeline@v1.3.0-alltime-deepseek · 16939 in / 10952 out tokens · 107652 ms · 2026-08-02T22:22:03.561310+00:00 · methodology

0 comments
read the original abstract

We develop a Gaussian process regression enhanced line integral string method to accelerate ring polymer instanton calculations of tunneling rates in molecular proton transfer reactions. By exploiting uncertainty estimates from the surrogate modeling, we show that the number of force evaluations required to converge an instanton path becomes effectively independent of the number of beads used to discretize the pathway. To reduce the computational overhead associated with training, particularly when Hessian information is included, we implement an efficient training strategy by combining physical GPR prior, Hessian-free GPR training and graphics processing unit accelerated black box matrix matrix multiplication, achieving an order of magnitude speedups relative to standard implementations. For rate calculations, we introduce a selective Hessian training strategy that distinguishes flexible modes strongly coupled to the transferring proton from more rigid modes weakly coupled to the reaction coordinate. This enables the construction of accurate surrogate potential energy surfaces with reduced Hessian evaluations. We apply both cubic spline interpolation method and Gaussian Process Regression to approximate the instanton rate for the prototypical systems, malonaldehyde, Z-3-aminopropenal and 7,9-dinitro-10-hydroxybenzo[h]quinoline. In our numerical test, the spline interpolation emerges as a simple and computationally efficient approach for the instanton rate calculations.

Figures

Figures reproduced from arXiv: 2602.16962 by Amke Nimmrich, Axel Gomez, Chenghao Zhang, Munira Khalil, Niranjan Govind.

Figure 1
Figure 1. Figure 1: Comparison of the number of force evaluations required to achieve convergence as a [PITH_FULL_IMAGE:figures/full_fig_p010_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: GPR hyperparameter training time for three different molecules using the pseudo [PITH_FULL_IMAGE:figures/full_fig_p012_2.png] view at source ↗

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

7 extracted references · 2 linked inside Pith

  1. [1]

    The multi-configurational time-dependent Hartree approach

    (1) Meyer, H.-D.; Manthe, U.; Cederbaum, L. The multi-configurational time-dependent Hartree approach. Chemical Physics Letters1990, 165, 73–78. (2) Beck, M.; J¨ ackle, A.; Worth, G.; Meyer, H.-D. The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets. Physics Reports2000, 324, 1–105. (3) Wan...

    Pith/arXiv arXiv 1906

  2. [10]

    (37) ´Asgeirsson, V.; Arnaldsson, A.; J´ onsson, H

    2012; pp 45–55. (37) ´Asgeirsson, V.; Arnaldsson, A.; J´ onsson, H. Efficient evaluation of atom tunneling combined with electronic structure calculations.The Journal of Chemical Physics2018, 148, 102334. (38) Henkelman, G.; Uberuaga, B. P.; J´ onsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. The J...

    2012

  3. [27]

    Minimum Mode Saddle Point Searches Using Gaussian Process Regression with Inverse-Distance Covariance Function

    (54) Koistinen, O.-P.; ´Asgeirsson, V.; Vehtari, A.; J´ onsson, H. Minimum Mode Saddle Point Searches Using Gaussian Process Regression with Inverse-Distance Covariance Function. Journal of Chemical Theory and Computation2020, 16, 499–509, PMID: 31801018. (55) Denzel, A.; K¨ astner, J. Hessian Matrix Update Scheme for Transition State Search Based on Gaus...

    2020

  4. [31]

    N.; Richardson, J

    32 (49) Beyer, A. N.; Richardson, J. O.; Knowles, P. J.; Rommel, J.; Althorpe, S. C. Quantum Tunneling Rates of Gas-Phase Reactions from On-the-Fly Instanton Calculations. The Journal of Physical Chemistry Letters2016, 7, 4374–4379, PMID: 27775889. (50) Bitzek, E.; Koskinen, P.; G¨ ahler, F.; Moseler, M.; Gumbsch, P. Structural Relaxation Made Simple. Phy...

    2006

  5. [1269]

    P.; Kondor, R.; Cs´ anyi, G

    (73) Bart´ ok, A. P.; Kondor, R.; Cs´ anyi, G. On representing chemical environments.Phys. Rev. B2013, 87, 184115. (74) Li, F.; Yang, X.; Liu, X.; Cao, J.; Bian, W. An Ab Initio Neural Network Potential En- 35 ergy Surface for the Dimer of Formic Acid and Further Quantum Tunneling Dynamics. ACS Omega2023, 8, 17296–17303. (75) Zhang, Y.; Li, W.; Luo, W.; Z...

    1991

  6. [1999]

    A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines

    (60) Hutchinson, M. A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Communications in Statistics - Simulation and Computation1990, 19, 433–450. (61) Meyer, R. A.; Musco, C.; Musco, C.; Woodruff, D. P. Hutch++: Optimal Stochastic Trace Estimation. 2021;https://arxiv.org/abs/2010.09649. (62) Litman, Y.; Richardso...

    Pith/arXiv arXiv 2021

  7. [6897]

    O.; Althorpe, S

    (64) Richardson, J. O.; Althorpe, S. C. Ring-polymer instanton method for calculating tun- neling splittings. The Journal of Chemical Physics2011, 134, 054109. (65) Kawatsu, T.; Miura, S. Efficient algorithms for semiclassical instanton calculations based on discretized path integrals.The Journal of Chemical Physics2014, 141, 024101. 34 (66) Richardson, J...

    2017