# Force-free identification of minimum-energy pathways and transition states
for stochastic electronic structure theories

Gopal R. Iyer  Department of Chemistry, Brown University, Providence, Rhode
Island 02912, USA  Noah Whelpley  Department of Chemistry, Brown University,
Providence, Rhode Island 02912, USA  Juha Tiihonen  Department of Physics,
Nanoscience Center, University of Jyväskylä, Jyväskylä, Finland  Paul R. C.
Kent  Computational Sciences and Engineering Division, Oak Ridge National
Laboratory, Oak Ridge, Tennessee 37831, USA  Jaron T. Krogel  Materials
Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge,
Tennessee 37831, USA<EMAIL_ADDRESS>Brenda M. Rubenstein  Department of
Chemistry, Brown University, Providence, Rhode Island 02912, USA
<EMAIL_ADDRESS>
###### Abstract

The accurate mapping of potential energy surfaces (PESs) is crucial to our
understanding of the numerous physical and chemical processes mediated by
atomic rearrangements - such as conformational changes and chemical reactions
- and the thermodynamic and kinetic feasibility of these processes. Stochastic
electronic structure theories, e.g., Quantum Monte Carlo (QMC) methods, enable
highly accurate total energy calculations which in principle can be used to
construct the PES. However, their stochastic nature poses a challenge to the
computation and use of forces and Hessians, which are typically required in
algorithms for minimum-energy pathway (MEP) and transition state (TS)
identification, such as the nudged-elastic band (NEB) algorithm and its
climbing image formulation. Here, we present strategies that utilize the
surrogate Hessian line-search method - previously developed for QMC structural
optimization \- to efficiently identify MEP and TS structures without
requiring force calculations at the level of the stochastic electronic
structure theory. By modifying the surrogate Hessian algorithm to operate in
path-orthogonal subspaces and on saddle points, we show that it is possible to
identify MEPs and TSs using a force-free QMC approach. We demonstrate these
strategies via two examples, the inversion of the ammonia (NH3) molecule and
the nucleophilic substitution (S${}_{\textrm{N}}$2) reaction F${}^{-}+$ CH3F
$\rightarrow$ FCH3 \+ F-. We validate our results using Density Functional
Theory (DFT)- and coupled cluster-based NEB calculations. We then introduce a
hybrid DFT-QMC approach to compute thermodynamic and kinetic quantities - free
energy differences, rate constants, and equilibrium constants - that
incorporates stochastically-optimized structures and their energies, and show
that this scheme improves upon DFT accuracy. Our methods generalize
straightforwardly to other systems and other high-accuracy theories that
similarly face challenges computing energy gradients, paving the way for
highly accurate PES mapping, transition state determination, and thermodynamic
and kinetic calculations, at significantly reduced computational expense.

##  1 Introduction

One of the crowning achievements of modern quantum chemistry is the ability to
compute potential energy surfaces (PESs) entirely from first principles. 1, 2
With knowledge of these surfaces, researchers can predict the reactivity,
dynamics, and kinetics of a wide range of technologically important and
scientifically interesting atomistic systems. 3, 4 For instance, the stable
structures of the reactants and products of a catalytic reaction can be
determined by identifying the minima on its PES using a variety of geometry
optimization algorithms. 5, 6, 7 Reaction mechanisms can, moreover, be
predicted by determining the minimum-energy pathway (MEP) that connects the
reactant and product states of a reactive system using, for instance, chain-
of-states methods8, 9 or the nudged elastic band (NEB) method10 and its
variations11, 12. These methods can often be straightforwardly extended -
e.g., using the climbing-image nudged elastic band (CI-NEB) method13 \- to
identify the (usually) first-order saddle point of the PES along the MEP,
which corresponds to the transition state. Reaction mechanisms predicted in
this fashion provide invaluable atomistic information about reaction pathways,
granting researchers key insights into how reactions can be redesigned and
potentially controlled to minimize energy, material, and other resource
requirements. Indeed, the ability to computationally predict reaction
mechanisms has become an indispensable tool in the fields of combustion
analysis,14 catalysis,15, 16 and chemical reactor design. 17

Nonetheless, while virtually all quantum chemical methods can predict the
energy of a chemical system given the geometries of the species it contains,
it can be challenging for some methods to do the inverse, i.e., to very
accurately predict equilibrium (reactant/product) and transition states, and
reaction pathways given the energies of structures sampled in some local
region of the PES. This is because being able to efficiently determine
equilibrium structures, transition states, and MEPs requires knowledge of
energy derivatives, including forces (first derivatives) or Hessian matrices
(second derivatives). These quantities then guide the dynamics of initial
structures to the true equilibrium or saddle point structures that are
consistent with the given level of ab initio theory. Computing such
derivatives in methods like Density Functional Theory (DFT) can be done in a
relatively straightforward fashion using the Hellman-Feynman Theorem,18 making
DFT one of the most popular methods for determining molecular structures.
However, computing these derivatives in more expensive many-body theories,
including Coupled Cluster Theory (particularly with the inclusion of single,
double, and perturbative triple excitations, CCSD(T)),19, 20 the Random Phase
Approximation (RPA),21, 22 or Quantum Monte Carlo (QMC)23, 24, 25, 26, 27, 28,
29, 30, 31, 32 methods can be either approximate or highly costly. In the case
of stochastic many-body theories such as QMC, this also requires formal
treatment of the energy derivatives in the presence of inherent statistical
errors in the energies. 23 As a result, instead of determining these
geometries or paths using highly accurate many-body methods, many resort to
using potentially less accurate theories like DFT to determine structures and
then higher-accuracy theories to determine their corresponding energies. 33
Although DFT geometries can be sufficiently accurate in the many cases in
which electron correlation does not play a critical role in dictating
geometries, computing properties using many-body theories based on single-
reference geometries can lead to a fundamental mismatch in theories that can
result in inaccurate or inconsistent predictions. For instance, as the authors
have previously shown, lower-accuracy DFT geometries of 2D materials that can
be readily stretched/compressed along particular axes can yield misleading
predictions of magnetic and other electronic properties when studied using
more accurate quantum chemical methods. 34, 35 It would thus be ideal to use
the same high-accuracy method to both predict structures and their energies
and other properties.

Recently, a variety of hybrid approaches for circumventing and reducing the
cost of directly computing energy gradients using relatively expensive high-
accuracy methods have been proposed. Several approaches, particularly machine
learning (ML)-based methods for mapping PESs and computing QMC forces, have
been demonstrated successfully. 36, 37, 38, 39 While such ML methods open up
exciting new avenues for the scalable modeling of large molecules and complex
chemical reactions, 39 they are typically highly reliant on the generation of
large training sets. This can limit their practicality for systems for which
it is difficult to assemble sufficient data, curtailing their generalizability
beyond the systems on which they have been trained. Furthermore, the accuracy
of ML methods can be significantly improved if they are trained with both
energies _and_ forces,37, 40 which may be challenging to accomplish if forces
are not readily provided by the methods used to furnish training data.

An alternative set of approaches that do not necessitate large training sets
or forces are so-called surrogate approaches, which use information easier to
obtain from less accurate, but also less expensive, theories to guide more
accurate, more expensive theories. Some of the first ideas along these lines
were proposed in Ref. 41. More recently, Ref. 42 introduced the surrogate
Hessian line-search method, which locates energy minima to within some tunable
structural error tolerance via a series of one-dimensional line searches that
account for and leverage statistical errors in QMC energy calculations. This
method was then successfully applied to the structural optimization of both
molecular42 and material systems. 34, 35

Building upon the success of this approach for identifying equilibrium
geometries, in this work, we develop new formalisms and algorithms for
generalizing this approach to the determination of transition state geometries
and other structures along the minimum energy pathways of chemical processes.
Key to generalizing these techniques is the realization that reaction pathway
identification tasks can be viewed as iterative energy optimization steps.
Transition states (saddle points) can be identified by minimizing the energy
along all directions except one that has negative curvature, along which the
energy is maximized instead. MEPs can be identified by minimizing the energy
along all directions in the subspace orthogonal to the tangent to the path at
each structure along some guess discretized path. We note here that the idea
of reformulating MEP identification as a series of energy optimization
problems instead of force-based dynamical problems has been suggested before,
notably in Ref. 43. Additionally, Saccani et al. (Ref. 44) demonstrated the
application of the NEB algorithm to QMC-based MEP identification using
correlated sampling. 26 However, to the best of our knowledge, ours is the
first approach that generalizes MEP identification to stochastic electronic
structure methods by utilizing Hessian information from lower-level surrogate
theories, eliminating the need for computing derivatives at the QMC level. We
demonstrate the accuracy and utility of this formalism by determining the
transition states and minimum energy pathways of two paradigmatic chemical
processes: (i) ammonia (NH3) inversion, which involves two structural degrees
of freedom, and (ii) and an S${}_{\textrm{N}}$2 nucleophilic substitution
reaction, F- \+ CH3F $\rightarrow$ FCH3 \+ F-,45 that involves four degrees of
freedom, using Diffusion Quantum Monte Carlo (DMC) informed by DFT Hessian
computations. We compare the minimum energy pathways we compute using DMC with
those obtained using DFT- and Coupled Cluster-based NEB calculations,
demonstrating that our surrogate techniques can provide at least post-DFT
quality pathways at only a small factor of $M\times(N-1)$ greater cost than
traditional single-point DMC energy calculations, where $M$ is the number of
points along which a discrete one-dimensional line-search grid is generated
and $N$ is the number of structural parameters used to parameterize the
system. We subsequently show that we can use a combined DFT-QMC approach to
more accurately estimate thermodynamic quantities, such as free energies of
reactions, and kinetic quantities, such as rate constants, which can often be
more directly compared to experiments. Overall, this work provides clear
avenues for determining reaction pathways and transition states with chemical
(or greater) accuracy at substantially lower computational cost than direct
approaches, but without the need for the large training datasets that
statistical learning-based methods require. We believe that these techniques
will provide new tools for the study of reactions and processes that involve
significant amounts of electron correlation and the design of novel
homogeneous and heterogeneous catalysts that can accelerate these reactions.

This manuscript is organized as follows. In Section 2, we introduce the
theoretical foundations of the surrogate Hessian line-search method and then
show how this formalism can be generalized to transition states (Section
2.1.1) and minimum energy pathways (Section 2.1.2). We also describe how
thermodynamic corrections to obtain free energies and equilibrium constants
can be performed using the hybrid DFT-QMC approach (Section 2.1.3). We then
present our results and benchmarks on our two illustrative chemical processes
in Section 3. Lastly, in Section 4, we conclude with the impacts and potential
further generalizations of this method.

##  2 Methods

###  2.1 Theoretical Background

####  2.1.1 Surrogate Hessian Line-Search Method for Equilibrium and
Transition States

The original surrogate Hessian line-search method for identifying equilibrium
structures (energy minima) has previously been introduced in Ref. 42. For
brevity, we only summarize the method here while noting the specific
modifications necessary to extend the method to transition state
identification.

In the surrogate Hessian line-search method, a surrogate method (such as DFT)
that can more rapidly compute energy gradients is used to calculate Hessians
that can then be used to determine the directions along which line searches
for minima should be performed using energies from a higher-accuracy and
costlier theory. This method relies on the assumption that the local curvature
of the target, higher-accuracy PES near a critical (0-derivative) point -
energy minimum or saddle point - is approximated reasonably well by the PES
calculated using a surrogate theory. This allows one to compute (a) the
optimal directions along which to perform the line-search, and (b) the maximum
tolerable statistical noise (if a stochastic method is employed) at each point
on the line-search grid in order to be able to determine an optimized geometry
within some structural tolerance. This information is then used to perform a
series of one-dimensional line-searches to locate the critical point within
some confidence interval.

While the surrogate and higher accuracy approaches that can be employed are
highly generalizable, throughout the rest of this manuscript, we focus on the
specific choices of Density Functional Theory (DFT) as our surrogate method
and DMC46, 47 as our higher accuracy method. This choice of methods is natural
since DFT calculations are often used as inputs into DMC calculations. In this
context, the key steps of the surrogate Hessian line-search method can be
summarized as follows:

  1. 1.

Initialization and parameterization: To provide an initial seed structure to
the equilibrium or saddle point line-search, we use an optimized equilibrium
or transition state geometry obtained using the surrogate theory (DFT). This
may be the result of a typical structural relaxation or CI-NEB calculation,
respectively. In systems where some degrees of freedom are correlated due to
structural symmetries, it is often desirable to define a few structural
parameters that capture these symmetries and only describe the relevant
regions of configuration space instead of operating in the full space of $3N$
Cartesian coordinates (where $N$ is the number of atoms). We achieve this by
defining a mapping function
$\mathrm{F^{I}}:\mathbb{R}^{3N}\rightarrow\mathbb{R}^{d}$ (given $d$
structural parameters). It is defined to ensure that all structures that share
these symmetries map to the same representation in $\mathbb{R}^{d}$, while its
inverse ($(\mathrm{F^{I}})^{-1}$) maps a structure in $\mathbb{R}^{d}$ to a
single structure in $\mathbb{R}^{3N}$. When working with solids, it is also
reasonable to include the dimensions of the unit cell as additional structural
parameters to allow for relaxations of the overall lattice. The structural
parameterizations of the NH3 inversion and S${}_{\textrm{N}}$2 reaction
systems studied in this manuscript are described in Section 2.2.

  2. 2.

Hessian evaluation: Given a reduced structural representation of the initial
state, we then compute the Hessian matrix of the PES in the vicinity of that
structure using the surrogate theory. Here, since we only define a small
number of degrees of freedom, we evaluate the Hessian using a finite-
difference scheme. For more complex systems, it may be beneficial to evaluate
it in generalized coordinates and then map it to the reduced space. Notably,
when dealing with equilibrium geometries, we expect all the Hessian
eigenvalues to be positive, while we expect $n$ eigenvalues to be negative
when dealing with $n^{\textrm{th}}$-order saddle points. This distinction will
become important when defining the objective function of the line-search.
Furthermore, the Hessian eigenvectors provide near-optimal directions along
which to perform the line-search using the stochastic theory.

  3. 3.

Tolerance optimization and line-search: The Hessian eigenvectors now provide
us with directions along which we can perform a series of one-dimensional
line-search operations to determine the desired equilibrium or saddle point.
This requires evaluating energies using the stochastic theory at a few points
along each line-search direction. To determine these points, it is necessary
to identify the optimal (minimal) degree of stochastic energy sampling
required to be able to specify the structural parameters with a certain pre-
specified level of precision ($\delta P$). For each line-search direction,
this is equivalent to evaluating the maximum point-wise energy noise tolerated
along each line-search direction, $\sigma^{\textrm{max}}$, and the
corresponding structural window along which the energy grid is sampled,
$L^{\textrm{max}}$, both of which are functions of $\delta P$. Finally, along
each line-search direction $x$, the objective of the line-search for
determining the equilibrium geometry can be expressed simply as an energy
minimization criterion

$x_{\textrm{min}}^{*}=\operatorname*{arg\,min}E(x).$ (1)

The modification to this objective for locating a saddle point is then
relatively straightforward. If a certain line-search direction corresponds to
a negative eigenvalue, instead of minimizing the energy along that direction,
we maximize it. This satisfies the definition of a saddle point as lying at a
local maximum on the PES along some (small number of) directions while lying
at a local minimum along all others. The line-search objective is then
modified as

$x_{\textrm{saddle}}^{*}=\begin{cases}\operatorname*{arg\,min}E(x)&\text{if
}\lambda_{x}>0\\\ \operatorname*{arg\,max}E(x)&\text{if
}\lambda_{x}<0\end{cases},$ (2)

where $\lambda_{x}$ is the eigenvalue corresponding to the Hessian eigenvector
$x$.

We direct the reader to Ref. 42 for further details regarding the surrogate
Hessian line-search method.

####  2.1.2 Path-Orthogonal Subspace Optimization

We now describe how the surrogate Hessian line-search method can be adapted to
identify MEPs. We begin by defining an MEP according to Fukui’s formulation of
the intrinsic reaction coordinate (see Ref. 48) as a one-dimensional sub-
manifold in configuration space that satisfies the following conditions: (1)
connects two equilibrium states; (2) passes through a saddle point (assumed to
be first-order for simplicity); and (3) consists exclusively of points that
lie in local energy minima in all directions except along the tangent to the
path at those points.

MEPs are typically identified using chain-of-states methods such as NEB (Fig.
1(a)). 10 The NEB algorithm begins with the construction of an initial guess
pathway, usually a linear interpolation of $M$ points between two equilibrium
states (satisfying condition 1). For each of these pathway structures
(images), tangents are calculated along the path. Then, two types of forces
acting on (each atom of) each image are computed: (1) the “true” force, i.e.,
the negative gradient of the potential energy calculated according to the
prescribed level of theory, and (2) an artificial “spring” force calculated as
a function of the excess distance between adjacent images compared to some
equilibrium distance, usually set to be $\mathcal{D}/(M+1)$, where
$\mathcal{D}$ is the Euclidean distance between the initial and final states.
(Note that the result of an NEB calculation ought to be agnostic to the spring
constant, provided it is set to a reasonable value, typically $\sim 0.1$
eV/Å2.) These two forces are then resolved into components parallel and
orthogonal to the path-tangent at each point, after which the component of the
true force along the tangent and the component of the spring force orthogonal
to the tangent are set to zero. Finally, the remaining forces are used to
dynamically propagate the system according to some dynamical rule, such as the
BFGS or velocity Verlet algorithms. This process is performed iteratively
until a self-consistent path is obtained, i.e., the force is converged to
within some tolerance (satisfying condition 3). Since NEB identifies a
discrete path, condition 2 is not immediately satisfied. However, CI-NEB
addresses this by maximizing the energy of the highest-energy state along the
converged NEB path. This is done by setting the path-orthogonal force
components to zero and only moving the image along the path.

((a))

((b))

Figure 1: Schematic illustration of the (a) NEB and (b) surrogate Hessian
subspace optimization algorithms on a contour map of the true PES. Yellow and
red contours indicate local energy minima and maxima, respectively.

The optimization of an NEB path based on force convergence criteria can be
understood as a constrained geometry optimization of each path image
orthogonal to the path. Therefore, to identify the MEP given a discrete
initial path between two equilibrium structures, we only require a geometry
optimization algorithm capable of operating in the subspace orthogonal to the
tangent at each point of the path.

In the case of the surrogate Hessian line-search algorithm, this translates to
calculating the surrogate Hessian matrix in the path-orthogonal subspace at
each point along an initial guess pathway and performing the line-search
within this subspace. Borrowing a DFT-NEB pathway as an initial guess thus
permits us to use DMC to obtain a more accurate, corrected MEP. Hereafter, we
refer to this as the subspace optimization method. The key steps of subspace
optimization for identifying MEPs are the following:

  1. 1.

Initialization and parameterization: We begin with a seed MEP that is obtained
using an NEB or CI-NEB calculation within the surrogate theory, i.e., DFT in
our case. Suppose this path is given by $\\{\mathbf{R}_{0}\equiv
e_{A},\mathbf{R}_{1},\mathbf{R}_{2},\ldots,\mathbf{R}_{M},\mathbf{R}_{M+1}\equiv
e_{B}\\}$, where $e_{A}$ and $e_{B}$ are the initial and final equilibria
(e.g., the reactants and products or intermediates for a chemical reaction
pathway). Each intermediate image between the initial and final equilibria is
then parameterized using the same scheme
$\mathrm{F^{I}}:\mathbb{R}^{3N}\rightarrow\mathbb{R}^{d}$ that is used to
parameterize the equilibrium and transition states to give
$\\{P_{0},P_{1},P_{2},\ldots,P_{M},P_{M+1}\\}$, where
$P_{i}\equiv\mathrm{F^{I}}(\mathbf{R}_{i})$.

  2. 2.

Tangent evaluation: For each intermediate structure along the pathway, we
evaluate the tangent to the path at that structure. There are several ways to
accomplish this. Here, we take the vector difference between the current
structure and the adjacent structure that is closer to the saddle point. If
the current structure is itself the closest structure to the saddle point, we
take the mean of the normalized vector differences between the current image
and each of its two adjacent images. The key difference with the standard
implementation is that we compute the vector differences in structural
parameter space rather than in Cartesian coordinates:

$\tau_{i}=\begin{cases}P_{i+1}-P_{i}&\text{if }i<i_{\textrm{saddle}}\\\
\frac{P_{i+1}-P_{i}}{|P_{i+1}-P_{i}|}+\frac{P_{i}-P_{i-1}}{|P_{i}-P_{i-1}|}&\text{if
}i=i_{\textrm{saddle}}\\\ P_{i}-P_{i-1}&\text{if
}i>i_{\textrm{saddle}}\end{cases},$ (3)

where $i_{\textrm{saddle}}$ is the index of the structure closest to the
saddle point.111We make the simplifying assumption here that the surrogate CI-
NEB MEP (correctly) identifies a path whose energy increases monotonically
until the saddle point and then decreases monotonically. This allows us to
directly index the structure closest to the saddle point using the surrogate
MEP energies. If we don’t make this assumption, the tangent equation would
have to be modified slightly to account for the energy of each image relative
to its neighbors’ energies in any given iteration of the algorithm.

  3. 3.

Subspace vectors and Hessian evaluation: Having computed the normalized path-
tangents, we then use the Gram-Schmidt process to construct unit vectors that
are orthogonal to each tangent. For a given tangent $\tau_{i}$, these unit
vectors $\mathbf{v}_{i}$ then define the path-orthogonal subspace of the
overall structural parameter space. At this point, we define secondary
mappings $\mathrm{F^{II}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d-1}$ that map
a point from structural parameter space to the path-orthogonal subspace at
each intermediate structure, which is set to be the origin of the subspace.
This mapping, $\mathcal{P}_{i}\equiv\mathrm{F^{II}}(P_{i}|\mathbf{v}_{i})$, is
then used as the basis for computing a finite-difference Hessian matrix in the
subspace with the surrogate theory. We also construct an inverse mapping
$(\mathrm{F^{II}})^{-1}:\mathbb{R}^{d-1}\rightarrow\mathbb{R}^{d}$ to keep
track of line-search optimized structures in the original space of structural
parameters.

  4. 4.

Tolerance optimization and line-search: Finally, for each intermediate
structure, we optimize the line-search parameters
$(L_{i}^{\textrm{max}},\sigma_{i}^{\textrm{max}})$ in the subspace and perform
the usual line-search, as described in Section 2.1.1, to locate the local
energy minimum within the subspace. We then map each optimized structure,
$\mathcal{P}_{i}^{\textrm{opt}}$, back to the original space of structural
parameters as $(\mathrm{F^{II}})^{-1}(\mathcal{P}_{i}^{\textrm{opt}})$, update
each intermediate structure to be the new line-search optimum, and repeat from
Step 2 onwards.

We summarize the subspace optimization method in Algorithm 1 and show a
schematic diagram in Fig. 1(b). We note that, since the tangents are pre-
computed for all the intermediate structures in a given iteration of subspace
optimization, the aforementioned operations - corresponding to the inner for
loop in Algorithm 1 \- can be performed on each structure in parallel.

Algorithm 1 Surrogate Hessian subspace optimization algorithm

Consider an $N$-atom system with two equilibrium states, $e_{A}$ and $e_{B}$,
for which an $M$-image MEP is to be determined with structural accuracy
$\delta P=\\{\delta p_{1},\delta p_{2},\ldots,\delta p_{M}\\}$.

start

$\\{\mathbf{R}_{0}=e_{A},\mathbf{R}_{1},\mathbf{R}_{2},\ldots,\mathbf{R}_{M},\mathbf{R}_{M+1}=e_{B}\\}\leftarrow$
CI-NEB($e_{A},e_{B}|M$)

$P\leftarrow$ $\\{$$\mathrm{F^{I}}$($\mathbf{R}_{0}$),
$\mathrm{F^{I}}$($\mathbf{R}_{2}$),
$\ldots,$$\mathrm{F^{I}}$($\mathbf{R}_{M+1}$)$\\}$

$\triangleright$ $\mathrm{F^{I}}:\mathbb{R}^{3N}\rightarrow\mathbb{R}^{d}$
maps a structure from real space to structural parameter space

for $j\leftarrow 1,2,\ldots,N_{\textrm{iter}}$ do

for $i\leftarrow 1,2,\ldots,M$ do

$\tau_{i}\leftarrow$ Tangent($P_{i}|P_{i-1},P_{i+1}$)

$\triangleright$ Computes the path-tangent at image $i$ given its neighboring
images

$\mathbf{v}_{i}\leftarrow$ GramSchmidt($\tau_{i}$)

$\triangleright$ Uses the Gram-Schmidt process to identify $(d-1)$ vectors
orthogonal to $\tau_{i}$

$\mathcal{P}_{i}\leftarrow$ $\mathrm{F^{II}}$($P_{i}|\mathbf{v}_{i}$)

$\triangleright$ $\mathrm{F^{II}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d-1}$
maps a structure from structural parameter space to a $(d-1)$-dimensional
path-orthogonal subspace

$\triangleright$ Now, we perform the usual surrogate Hessian line-search as
described in Ref. 42

$H_{i}\leftarrow$ Hessian($\mathcal{P}_{i}$) =
$\nabla^{2}E_{\textrm{surrogate}}(\mathcal{P}_{i})$

$(L_{i}^{\textrm{max}},\sigma_{i}^{\textrm{max}})\leftarrow$
Optimize($H_{i}|\delta P$)

$\mathcal{P}_{i}^{\textrm{opt}}(\Delta\mathcal{P}_{i}^{\textrm{opt}})\leftarrow$
LineSearch($L_{i}^{\textrm{max}},\sigma_{i}^{\textrm{max}}$)

end for

$\\{P_{1},P_{2},\ldots,P_{M}\\}\leftarrow$
$\\{$$(\mathrm{F^{II}})^{-1}$($\mathcal{P}_{1}^{\textrm{opt}}$),
$(\mathrm{F^{II}})^{-1}$($\mathcal{P}_{2}^{\textrm{opt}}$), $\ldots$,
$(\mathrm{F^{II}})^{-1}$($\mathcal{P}_{M}^{\textrm{opt}}$)$\\}$

end for

end

Note that while Algorithm 1 depicts running the outer loop for a predefined
number of iterations ($N_{\textrm{iter}}$), this would ideally be replaced
with a path convergence-based loop. In practice, however, we find that
beginning with a CI-NEB initial guess allows us to obtain accurate paths after
just a single iteration, due to which we retain the for loop format shown in
the algorithm.

The pathway obtained from this algorithm also presents an alternative approach
for locating the saddle point. Instead of computing the surrogate Hessian and
initializing the line-search at the highest-energy structure of the DFT-NEB
pathway as described in Section 2.1.1, we only compute the Hessian at the
highest-energy DFT-NEB structure and initialize the line-search at the
highest-energy structure on the DMC subspace optimization pathway. The utility
of this approach is that it permits us to initialize the saddle point line-
search at a point on a more accurate (DMC) pathway while only borrowing the
search directions (Hessian eigenvectors) from the surrogate theory, thereby
providing a more self-consistent protocol for transition state identification.
We hereafter refer to this as the self-consistent line-search approach for
transition states and discuss the application of this approach in Section 3.1.

####  2.1.3 Hybrid DFT-QMC Approach for Thermodynamics and Kinetics

Having obtained equilibrium and transition states, and path images, it is
desirable to calculate thermodynamic and kinetic quantities using our DMC-
optimized geometries and energies. In particular, let us suppose we are
interested in computing the free energy changes and equilibrium and rate
constants associated with a molecular transformation along an MEP using data
from a 0 K ground state ab initio theory. This can be done in thermochemistry
packages using molecular statistical mechanical principles by evaluating the
translational, rotational, and vibrational energies of the structures of
interest that compose the enthalpic and entropic contributions to the free
energy, given the temperature and the pressure. The equilibrium constant can
then be evaluated as

$K_{\textrm{eq}}=\exp\left(-\frac{\Delta G^{\circ}}{k_{B}T}\right),$ (4)

where $\Delta G^{\circ}$ is the free energy difference between the final and
initial equilibrium states. 49

Given the transition state, the rate constant for the forward reaction can be
obtained using transition state theory as

$k_{\rightarrow}=\frac{k_{B}T}{h}\exp\left(-\frac{\Delta
G^{\oplus}}{k_{B}T}\right),$ (5)

where we use $\Delta G^{\oplus}$ to refer to $(\Delta H^{\oplus}-T\Delta
S^{\oplus})$. 49 $\Delta S^{\oplus}$ is the entropy change in transforming
from the initial equilibrium state to the transition state, and $\Delta
H^{\oplus}$ is the activation energy of the forward reaction.
$K_{\textrm{eq}}$ and $k_{\rightarrow}$ then automatically determine the
backward reaction rate constant as
$k_{\leftarrow}=k_{\rightarrow}/K_{\textrm{eq}}$.

Performing these calculations using DMC, however, is not a straightforward
task. This is because, in practice, obtaining the vibrational energies
requires the evaluation of a finite-difference Hessian matrix in the full
$\mathbb{R}^{3N}$ Cartesian representation of the $N$-atom system. In addition
to requiring a number of expensive DMC calculations, the intrinsic statistical
nature of QMC methods introduces difficulty in stably and reliably
diagonalizing the stochastically sampled Hessian matrix, which is indeed the
original motivation for using the surrogate Hessian line-search method.

Surrogate Hessian line-search and subspace optimization are practicable
methods due to the - in most cases - correct assumption that the local
curvature of the surrogate and stochastic PESs are similar near the structures
of interest to a good approximation. It is then reasonable to suggest that we
may rely on DFT to provide only the vibrational energies while obtaining the
non-vibrational contributions to the enthalpy and entropy from DMC.

In the ideal gas limit, the contributions to the enthalpy may be decomposed as

$\begin{split}H(T)&=E_{\textrm{pot}}+E_{\textrm{ZPE}}+\int_{0}^{T}C_{P}dT\\\
&=E_{\textrm{pot}}+E_{\textrm{ZPE}}+\int_{0}^{T}(k_{B}+C_{V}^{\textrm{trans}}+C_{V}^{\textrm{rot}}+C_{V}^{\textrm{vib}}+C_{V}^{\textrm{elec}})dT,\end{split}$
(6)

where $E_{\textrm{pot}}$ is the potential energy, $E_{\textrm{ZPE}}$ is the
zero-point energy, and $C_{P}$ is the heat capacity at constant pressure,
which can be further broken down into the sum of the translational
($C_{V}^{\textrm{trans}}$), rotational ($C_{V}^{\textrm{rot}}$), vibrational
($C_{V}^{\textrm{vib}}$), and electronic ($C_{V}^{\textrm{elec}}$) heat
capacities at constant volume. 1 We suppose that the zero-point energy and the
vibrational component of the heat capacity can be obtained from DFT. Ignoring
the electronic heat capacity, it is now clear that the potential energy,
$E_{\textrm{pot}}$, is the only other component of the enthalpy that is
sensitive to the choice of ab initio method used to optimize the structure and
calculate its energy. We may then make an approximation to the enthalpy that
leverages the improved accuracy of DMC over DFT:

$\begin{split}H_{\textrm{hybrid}}(T)=\underbrace{E_{\textrm{pot}}^{\textrm{DMC}}}_{\text{calculated
at
}\langle\textbf{R}_{\textrm{DMC}}^{*}\rangle}+\underbrace{E_{\textrm{ZPE}}^{\textrm{DFT}}+\int_{0}^{T}C_{V,\textrm{vib}}dT}_{\text{calculated
at }\textbf{R}_{\textrm{DFT}}^{*}}\\\
+\int_{0}^{T}(k_{B}+C_{V}^{\textrm{trans}}+C_{V}^{\textrm{rot}})dT,\end{split}$
(7)

with the statistical error in the DMC potential energy now affecting the
enthalpy directly as

$\sigma(H_{\textrm{hybrid}}(T))=\sigma(E_{\textrm{pot}}^{\textrm{DMC}}),$ (8)

which is just the standard error of the DMC energy.

It is important to note that $E_{\textrm{pot}}^{\textrm{DMC}}$ and
$(E_{\textrm{ZPE}}^{\textrm{DFT}},C_{V}^{\textrm{vib}}|_{\textrm{DFT}})$ are
calculated using two slightly different geometries in Eq. 7.
$E_{\textrm{pot}}^{\textrm{DMC}}$ is taken to be the potential energy of the
mean surrogate Hessian-optimized DMC geometry
($\langle\textbf{R}_{\textrm{DMC}}^{*}\rangle$), while the vibrational
energies and all ensuing contributions are calculated using the DFT-optimized
geometry ($\textbf{R}_{\textrm{DFT}}^{*}$). This is crucial for the following
reason. Assume we are interested in an equilibrium structure. Then, the local
features of the PES near this structure will only be reflected in the DFT
Hessian matrix if the matrix is calculated around the DFT-equilibrium
structure. Since this may be slightly displaced from the DMC-equilibrium
structure, constructing a DFT Hessian matrix centered at the latter may
capture the curvature of an irrelevant region of the DFT PES. Instead, we
combine the curvature (Hessian matrix) near the equilibrium of the surrogate,
i.e., DFT, PES with the potential energy near the equilibrium of the more
accurate, stochastic PES (DMC). The same principle applies to saddle point
structures. This strategy is valid because the notion of a critical point
within some sufficiently local region of a PES is independent of the exact
structural parameters to which it corresponds and the specific level-of-theory
that was used to obtain it. This allows us to combine two slightly different
structures optimized using different methods because they theoretically refer
to the same region of the PES.

Now, let us turn our attention to the entropy, which can be decomposed as

$S(T,P)=S_{\textrm{trans}}(T,P^{\circ})+S_{\textrm{rot}}(T)+S_{\textrm{vib}}(T)+S_{\textrm{elec}}-k_{B}\ln\frac{P}{P^{\circ}},$
(9)

where $P^{\circ}$ is the pressure at the standard state and
$S_{\textrm{trans}}$, $S_{\textrm{rot}}$, $S_{\textrm{vib}}$, and
$S_{\textrm{elec}}$ refer to the translational, rotational, vibrational, and
electronic entropies respectively. 1 Again, $S_{\textrm{trans}}$ and
$S_{\textrm{elec}}$ are unaffected by slight differences between structures
optimized by different ab initio methods. As earlier, we can rely on DFT to
compute $S_{\textrm{vib}}$. Computing $S_{\textrm{rot}}$ analytically,
however, is a non-trivial task since it involves the calculation of molecular
moments-of-inertia - given the uncertainty in the structural parameters of a
DMC-optimized molecule, it is challenging to predict the propagation of these
errors to the moments-of-inertia and the rotational entropy. Not only are the
errors in structural parameters correlated, but the propagation of this error
also becomes hard to generalize for differently parameterized systems.

To overcome this, we employ a Monte Carlo scheme to estimate the rotational
entropy stochastically. First, we sample a large number of structures using
the mean and variance of the DMC-optimized structural parameters. After
converting these structures back to their Cartesian representation, we
calculate the rotational entropy for each of these sampled structures
independently, and finally report the average entropy and variance:

$\begin{split}S_{\textrm{hybrid}}(T,P)=\underbrace{S_{\textrm{rot}}(T)}_{\text{Monte
Carlo sampling around
}\langle\textbf{R}_{\textrm{DMC}}^{*}\rangle}+\underbrace{S_{\textrm{vib}}(T)}_{\text{calculated
at }\textbf{R}_{\textrm{DFT}}^{*}}\\\
+S_{\textrm{trans}}(T,P^{\circ})+S_{\textrm{elec}}-k_{B}\ln\frac{P}{P^{\circ}},\end{split}$
(10)

so that

$\sigma(S_{\textrm{hybrid}}(T))=\sqrt{\frac{\sum_{i=1}^{N_{\textrm{MC}}}\left(S_{\textrm{rot}}^{i}(T)-\langle
S_{\textrm{rot}}(T)\rangle\right)^{2}}{N_{\textrm{MC}}}},$ (11)

where $S_{\textrm{rot}}^{i}(T)$ is the rotational entropy of a structural
sample $i$ around $\langle\mathbf{R}_{\textrm{DMC}}^{*}\rangle$ and
$N_{\textrm{MC}}$ is the number of Monte Carlo samples.

From here, it is easy to calculate the free energy changes and equilibrium and
rate constants from Eqs. 4 and 5.

###  2.2 Illustrative Test Cases

In this work, we focus on two illustrative chemical processes as test cases:
(1) ammonia (NH3) inversion and (2) the nucleophilic substitution
(S${}_{\textrm{N}}$2) reaction F- \+ CH3F $\rightarrow$ FCH3 \+ F-. Ammonia
inversion is a process that occurs $\sim$30 billion times per second50 at room
temperature in which the three hydrogen atoms that form a trigonal pyramid
around nitrogen flip from one side of the N to the other. In the
S${}_{\textrm{N}}$2 reaction, one fluorine substitutes for another in CH3F in
a concerted (hence the 2) fashion, resulting in an inversion of the methyl
group to accommodate the bond rearrangement. This and other
S${}_{\textrm{N}}$2-type systems have been studied extensively using ab initio
methods; see, for example, Ref. 45. We now describe examples of minimal
structural parameterizations of the NH3 inversion and S${}_{\textrm{N}}$2
reaction systems that describe the transformations (inversion and bond
rearrangement, respectively) that these systems undergo.
We define the geometry of the NH3 system using two structural parameters (Fig.
2(a)): the mean N-H bond length in Å ($p_{0}$) and the mean angle in radians
($p_{1}$) formed by the N-H bonds with the negative $x$-axis (labeled
$x^{\prime}$). Without loss of generality, we assume that the centroid of the
three H atoms lies on the $x$-axis in all configurations during the inversion
of the molecule. We set the N atom to always coincide with the origin and
constrain one of the H atoms (labeled H1) to lie in the $x-y$ plane. We
similarly parameterize the S${}_{\textrm{N}}$2 reaction system (Fig. 2(b))
using the mean C-H bond length ($p_{0}$) and the mean
$\measuredangle$HC$x^{\prime}$ bond angle ($p_{1}$). Additionally, we include
the two C-F bond lengths, C-F1 ($p_{2}$) and C-F2 ($p_{3}$). We set the C atom
to coincide with the origin, the H1 atom to lie in the $x-y$ plane, and the
C-F1 bond to lie along the $x^{\prime}$ axis. It follows that
$\measuredangle$HC$x^{\prime}$ = $\measuredangle$HCF1. Explicit expressions
for both parameterizations are provided in Section LABEL:supporting:mapping of
the Supporting Information.

In both systems, the initial and final equilibrium geometries are identical
and share $C_{3v}$ symmetry. However, grounding the parameterization in a
Cartesian reference frame enables us to distinguish between the initial and
final states in our analysis.

((a))

((b))

Figure 2: Structural parameterization of (a) the NH3 inversion system, and (b)
the S${}_{\textrm{N}}$2 reaction system. All bond lengths are measured in Å
and bond angles in radians.

###  2.3 Computational Details

For the NH3 system, we set error tolerances of 0.01 Å for the bond length
($\delta p_{0}$) and 0.05 rad for the bond angle ($\delta p_{1}$). For the
S${}_{\textrm{N}}$2 reaction system, we set the following tolerances: $(\delta
p_{0}=0.01$ Å$,\delta p_{1}=0.05$ rad$,\delta p_{2}=0.05$ Å, and $\delta
p_{3}=0.05$ Å$)$. We choose higher tolerances for the two C-F bond lengths
because they vary to a greater extent over the course of the reaction, which
enables their variation to be captured accurately with lower precision and
hence lower computational cost.

We perform all DFT and Coupled Cluster with Single and Double excitations
(CCSD) calculations in PySCF. 51, 52 We use correlation-consistent effective
core potentials (ccECPs) and basis sets (cc-pVTZ-ccECP). 53, 54 We use the
PBE55 functional for our DFT calculations and use identical DFT settings for
computing energies, trial wave functions, and DFT-NEB paths.

For each QMC calculation, following DFT trial wave function generation, we
optimize one- and two-body Jastrow factors via 8 optimization cycles of
100,000 samples. We then further optimize the trial wave function with a VMC
calculation by taking 1000 samples. Finally, we perform a DMC calculation with
a timestep of 0.005 a.u. and run it for the estimated number of steps required
to compute the energies within the stipulated tolerances. We note that
translating error tolerances in structural parameters ($\delta p$) into the
number of DMC steps required to saturate these tolerances involves a series of
approximations (see Ref. 42) and can result in higher or lower structural
errors than intended. To mitigate this, we extend the number of DMC steps in
some of our calculations by a factor of 16 or 64 (structure-wise values can be
found in Section LABEL:supporting:extension). We perform all QMC calculations
using QMCPACK56, 57 and assemble the DFT-VMC-DMC pipeline using the Nexus
workflow management system. 58

We perform all NEB calculations in the Atomic Simulation Environment (ASE)59
using the BFGS59 algorithm and a force tolerance of 0.01 eV/Å. For the NH3
system, we include 3 intermediate images between the initial and final
equilibrium states, whereas for the S${}_{\textrm{N}}$2 system, we include 5
intermediate images. To locate the saddle point, we perform a climbing image
calculation on the structure with the highest energy along the MEP. We also
use the BFGS algorithm for DFT and CCSD equilibrium geometry optimization. In
this case, we set the force tolerance to 0.001 eV/Å.

##  3 Results and Discussion

###  3.1 Equilibrium and Transition State Geometries

First, we report the results of the surrogate Hessian line-search method
applied to the optimization of the equilibrium and transition state geometries
of our two test systems. Tables 1 (NH3) and 2 (S${}_{\textrm{N}}$2) compare
the DFT and CCSD structural relaxations and CI-NEB results for the transition
states to those obtained using the surrogate Hessian line-search approach
using DMC as described in Section 2.1.1 for the NH3 inversion and
S${}_{\textrm{N}}$2 reaction systems. Additionally, to illustrate that our
approaches are self-consistent, we include the results for saddle point
identification using the self-consistent approach described in Section 2.1.2.
For this approach, we compute the line-search directions as the eigenvectors
of the surrogate Hessian at the DFT-NEB saddle point, while the line-search
itself is initialized at a high-energy image along the DMC MEP obtained using
subspace optimization. Given the structural symmetries of the NH3 and
S${}_{\textrm{N}}$2 systems and the odd number of intermediate images we have
considered here, the highest-energy image already coincides with the saddle
point in both cases. Thus, to make for more challenging test cases, we instead
seed the self-consistent saddle point line-search at a point adjacent to the
highest-energy image along each path. This corresponds to DMC path image 3 in
Fig. 4 for the NH3 case, and DMC path image 2 in Fig. 6 for the
S${}_{\textrm{N}}2$ case.

| Equilibrium A | Transition State | Equilibrium B
---|---|---|---
$p_{0}$ (Å) |  $p_{1}$ (rad) |  $p_{0}$ (Å) |  $p_{1}$ (rad) |  $p_{0}$ (Å) |  $p_{1}$ (rad)
DFT | 1.019 | 1.958 | 0.998 | 1.571 | 1.019 | 1.184
CCSD | 1.010 | 1.959 | 0.990 | 1.571 | 1.010 | 1.183
DMC (LS) | 1.006(1) | 1.957(5) | 0.991(2) | 1.569(14) | 1.007(2) | 1.202(10)
DMC (LS-SC) | - | - | 0.992(2) | 1.574(34) | - | -
Table 1: NH3 structural parameters obtained using: [Row 1] DFT structural relaxation for the equilibrium states and CI-NEB for the transition state; [Row 2] CCSD structural relaxation for the equilibrium states and CI-NEB for the transition state; [Row 3] DMC surrogate Hessian line-search; and [Row 4] DMC self-consistent saddle point line-search initialized at path image 3 in Fig. 4. | Equilibrium A
---|---
$p_{0}$ (Å) |  $p_{1}$ (rad) |  $p_{2}$ (Å) |  $p_{3}$ (Å)
DFT | 1.094 | 1.218 | 2.575 | 1.434
CCSD | 1.079 | 1.243 | 2.575 | 1.434
DMC (LS) | 1.079(3) | 1.244(16) | 2.580(17) | 1.412(11)
| Transition State
$p_{0}$ (Å) |  $p_{1}$ (rad) |  $p_{2}$ (Å) |  $p_{3}$ (Å)
DFT | 1.076 | 1.571 | 1.863 | 1.863
CCSD | 1.063 | 1.571 | 1.849 | 1.849
DMC (LS) | 1.066(3) | 1.568(14) | 1.824(13) | 1.816(14)
DMC (LS-SC) | 1.068(6) | 1.565(20) | 1.851(30) | 1.794(28)
| Equilibrium B
$p_{0}$ (Å) |  $p_{1}$ (rad) |  $p_{2}$ (Å) |  $p_{3}$ (Å)
DFT | 1.094 | 1.923 | 1.434 | 2.575
CCSD | 1.079 | 1.899 | 1.434 | 2.575
DMC (LS) | 1.078(2) | 1.894(8) | 1.412(10) | 2.580(14)
Table 2: S${}_{\textrm{N}}$2 reaction structural parameters obtained using:
[Row 1] DFT structural relaxation for equilibrium states and CI-NEB for the
transition state; [Row 2] CCSD structural relaxation for equilibrium states
and CI-NEB for the transition state; [Row 3] DMC surrogate Hessian line-
search; [Row 4 in Transition State block] DMC self-consistent saddle point
line-search initialized at path image 2 in Fig. 6.

From Tables 1 and 2, it is clear that the DMC line-search yields geometries
that are in excellent agreement with CCSD geometries obtained using BFGS
structural relaxation, with differences in bond lengths of less than
$\sim$0.02 Å and differences in bond angles of less than $\sim$0.02 rad (or
1.15∘). Additionally, agreement between the direct DMC line-search (DMC (LS))
and self-consistent DMC line-search (DMC (LS-SC)), within statistical
variation, indicates the capability of the surrogate Hessian approach to
determine transition state structures both in the presence and absence of
initial guess transition state geometries obtained from surrogate theories
such as DFT. We note that, while there are slightly larger ($\sim$0.05 Å)
differences between DMC and CCSD in the two C-F bond lengths ($p_{2}$ and
$p_{3}$) in the S${}_{\textrm{N}}$2 example, particularly in the transition
state structure, these are a consequence of the very low sensitivity of the
potential energy to these structural parameters. Since the surrogate Hessian
line-search method is based on non-linearly combining structural tolerances
into an energetic tolerance, locating a critical point to within some
energetic tolerance can incur - based on how strongly (or weakly) they are
coupled - commensurate variability in the structural parameters. This will
become clearer when we compare the energy barriers along the MEP obtained by
each method. For a more detailed analysis of the coupling of energetic and
structural tolerances, we refer the reader to Ref. 42.

###  3.2 Minimum-Energy Pathways

We now discuss the use of the subspace optimization algorithm to determine the
MEPs of the NH3 inversion and the S${}_{\textrm{N}}$2 reaction systems. We
begin by defining two equilibrium states, i.e., the DFT-optimized states
labeled Equilibrium A and Equilibrium B in Tables 1 and 2. For example, in the
NH3 system, Equilibrium B refers to the molecule shown in Fig. 2(a) with the
angle $\measuredangle$HN$x^{\prime}$ being acute, while Equilibrium A refers
to the inverted structure of this molecule, where the angle
$\measuredangle$HN$x^{\prime}$ becomes obtuse. In the S${}_{\textrm{N}}$2
case, both the angle $\measuredangle$HCF1 and the C-F1 and C-F2 bond lengths
change in going from Equilibrium A to Equilibrium B. For consistency of MEP
comparison between the various methods, we use these DFT-optimized
configurations as the fixed endpoints for all MEP calculations.

Starting with a linear interpolation between these two equilibria, we
independently perform NEB calculations using DFT and CCSD. Then, we use the
DFT-NEB states to compute the subspace surrogate Hessian and its eigenvectors
for each intermediate structure along the pathway. To compare a DFT-based
surrogate Hessian subspace optimization with the DFT-NEB pathway, we displace
the DFT-NEB images by a small fixed amount of 0.05 units along each of the
subspace eigenvectors. (It is difficult to associate this shift with a
particular unit, such as Å or radian, because the subspace is constructed by
taking combinations of vectors along the individual structural parameters, due
to which a certain vector in the path-orthogonal subspace has a component
along more than one structural parameter with no unique unit of measurement.)
We then use these displaced states to initialize the DFT-based subspace
optimization. We then proceed to identify the DMC MEP by performing a DMC-
based subspace optimization using the DFT-NEB pathway as an initial guess with
no displacement along the subspace eigenvectors. To illustrate the differences
between these pathways, we plot the NH3 inversion pathway in structural
parameter space overlaid on an energy heat map calculated using a fine grid of
CCSD calculations in Fig. 3.

Figure 3: MEPs for NH3 inversion obtained using subspace optimization and NEB
for different methods overlaid on a CCSD energy heat map. Energies are
measured in eV. The linear interpolation path is used as an initial guess for
the DFT-NEB and CCSD-NEB calculations. The displaced DFT guess path (displaced
from the fully converged DFT-NEB path) is used to seed the DFT subspace
optimization calculation. The x\- and y-axes indicate the NH3 structural
parameters, $p_{0}$ and $p_{1}$, respectively.

It is clear from Fig. 3 that the subspace optimization algorithm correctly
identifies the MEP due to the close agreement between the DFT-NEB and DFT-
subspace optimization pathways. Furthermore, the DMC subspace optimization
pathway shows better structural agreement with the CCSD-NEB pathway than the
DFT-based pathways, indicating the potential for more accurate (compared to
DFT) and inexpensive (compared to CCSD) MEP identification using surrogate
Hessian subspace optimization.

At this point, we note that, despite good agreement with NEB, subspace
optimization and NEB are not strictly equivalent algorithms and therefore
cannot always be expected to yield fully identical MEPs. The following are a
few notable distinctions between typical implementations of NEB and subspace
optimization:

  1. 1.

Tangent calculation: During an NEB calculation, tangents to the MEP at each
intermediate structure are calculated using vector differences between
adjacent structures in Cartesian coordinates. Subspace optimization instead
computes these tangents in the space of structural parameters. Since the
operations that map Cartesian coordinates to reduced structural parameters may
be non-distributive over vector differences between structures, the tangents -
and thereby, the path-orthogonal directions - computed using NEB and subspace
optimization may not be identical. Neither approach is necessarily “more
correct” since the choice of coordinates is a matter of representational
convenience. An alternative representation-invariant approach (not addressed
here) would be to construct an internal coordinate representation that
accounts for the implicit metric change that occurs when computing tangents at
points located in different regions of the PES60.

  2. 2.

Hessian evaluation: In NEB calculations that use force-based dynamical
algorithms such as BFGS to propagate the chain of states toward the MEP, the
Hessian matrix is evaluated and updated at each step of the dynamics. This
results in the possibility of multiple Hessian evaluations even within a
single NEB iteration. By contrast, the subspace optimization algorithm is
currently implemented to only perform a single evaluation of the subspace
Hessian in a given iteration, which is then used to perform multiple line-
search iterations. It is for this reason that providing good starting
structures for Hessian evaluation is necessary for this algorithm to operate
efficiently. While multiple Hessian evaluations can be implemented, we do not
include this in our approach.

  3. 3.

Use of springs: Since NEB is based on the resolution of forces along and
orthogonal to path-tangents, it is reasonable to introduce artificial spring
forces to maintain an even spacing of images along the path. However, since
such resolution of components is not possible in a purely energy-based
optimization scheme, subspace optimization does not include spring forces and
instead relies on evenly spaced starting images.

  4. 4.

Termination/convergence criterion: Implementations of NEB are typically based
on force-based convergence, where the method is set to iterate until the force
on each atom falls below some threshold. Subspace optimization, on the other
hand, may be implemented using a path-based convergence criterion, where the
structural difference between consecutive paths must fall below a certain
threshold for convergence.

Having obtained MEPs from each of the methods, we now consider the energy
barriers along the MEPs obtained with each method. We plot this for the NH3
system in Fig. 4, with the inclusion of a Gaussian process regression fit to
the DMC data points to depict the smooth variation of the mean energy and
standard energy error along the MEP. It must be noted that, while we plot each
energy barrier with a common $x$-axis labeled ‘Path image,’ these do not
correspond to the same reaction coordinate, since the MEP obtained using each
method is slightly different from the rest.

Figure 4: Energy vs. NH3 inversion path image for different methods. The DFT-
NEB and DFT subspace optimization points overlap almost completely.

It is clear from Fig. 4 that DMC-based subspace optimization is capable of
resolving energy barriers to within a few meV and is in very good agreement
with the CCSD-NEB energy barrier. We also note that in Fig. 4, the DFT-NEB and
DFT subspace optimization barriers are so close that the DFT-NEB data points
are not visible.

We now turn to the S${}_{\textrm{N}}$2 reaction system, which possesses more
degrees of freedom and is expected to exhibit greater electron correlation,
leading to more significant differences between the MEPs. Since it is
difficult to show the variation of all structural parameters in the same plot,
we show the variation of each structural parameter with respect to the path
image in Fig. 5. Again, it is evident that subspace optimization accurately
captures the trends in the variation of structural parameters across the MEP
using both DFT and DMC. We compare the energy barriers obtained by the
different methods in Fig. 6. Here, we see that DMC subspace optimization is in
excellent agreement with CCSD along the MEP, while the DFT-NEB and DFT
subspace optimization differ by less than 2.5 meV, or $\sim$ 0.1 mHa, from
each other. Note that the DFT and DMC MEP barrier heights differ by nearly 20
meV from one another, illustrating that correlation assumes a more significant
role in this test case.

Figure 5: Structural parameters vs. path image of the S${}_{\textrm{N}}$2
reaction system for different methods (see top right panel for legend). The
linear interpolation path is used as an initial guess for the DFT-NEB and
CCSD-NEB calculations. The displaced DFT guess path (displaced from the fully
converged DFT-NEB path) is used to seed the DFT subspace optimization
calculation. Figure 6: Energy vs. S${}_{\textrm{N}}$2 reaction path image for
different methods.

We note that, due to the distinctions between NEB and subspace optimization
mentioned above, the two methods do not yield perfectly identical MEP
structures at the DFT level. In our method, the structures along the MEP are
determined as a consequence of path-orthogonal energy optimization. Given that
several structures can be energetically very similar, the energy optimization
algorithm may rightly identify any one of these structures as lying on the
MEP, although it may not necessarily coincide with the structure identified by
a force-minimization algorithm such as NEB. Due to these practically
insignificant energy differences, paths obtained by both methods can be
considered to consist of representative points along the MEP despite small
structural dissimilarities.

###  3.3 Free Energies and Kinetic Parameters

We now discuss the evaluation of thermodynamic and kinetic quantities obtained
using the hybrid DFT-QMC approach described in Section 2.1.3. In particular,
we report the free energy difference between the saddle point and equilibrium
structures ($\Delta G^{\oplus}$ in Eq. 4) and between the two equilibrium
states ($\Delta G^{\circ}$ in Eq. 5). We also report the associated kinetic
parameters, $k_{\rightarrow}$ and $K_{\textrm{eq}}$. We benchmark these
results against DFT and CCSD calculations using the standard approach, in
which all energetic and entropic contributions are computed deterministically
using the same geometry. 59, 1

During the Monte Carlo structure sampling to estimate the rotational entropy
($S_{\textrm{rot}}(T)$ in Eq. 10), we sample 2500 structures each for the
saddle point and equilibrium states of NH3 using the mean and standard error
of the structural parameters identified using DMC in Table 1. Similarly, in
the case of the S${}_{\textrm{N}}$2 reaction, we sample 10,000 structures
around the mean and standard error values shown in Table 2.

|  $\Delta G^{\oplus}$ (eV) |  $k_{\rightarrow}(\times 10^{16}$ s${}^{-1})$ |  $\Delta G^{\circ}$ (eV) | $K_{\textrm{eq}}$
---|---|---|---|---
DFT | -0.199 | 1.436 | 0.000 | 1.000
Hybrid DFT-QMC | -0.200(0) | 1.483(24) | -0.000(1) | 1.018(33)
CCSD | -0.222 | 3.514 | 0.000 | 1.000
Table 3: Free energy differences and equilibrium and rate constants for the NH3 inversion system. |  $\Delta G^{\oplus}$ (eV) |  $k_{\rightarrow}(\times 10^{12}$ s${}^{-1})$ |  $\Delta G^{\circ}$ (eV) | $K_{\textrm{eq}}$
---|---|---|---|---
DFT | 0.004 | 5.317 | 0.000 | 1.000
Hybrid DFT-QMC | 0.019(1) | 2.918(122) | 0.000(0) | 0.989(19)
CCSD | 0.020 | 2.852 | 0.000 | 1.000
Table 4: Free energy differences and equilibrium and rate constants for the
S${}_{\textrm{N}}$2 reaction system.

The free energy differences and equilibrium and rate constants are shown in
Table 3 for the NH3 inversion system and in Table 4 for the
S${}_{\textrm{N}}$2 reaction system. Tables
LABEL:tab:NH3_thermodynamic_breakdown and
LABEL:tab:SN2_thermodynamic_breakdown show the free energy contributions of
the various ab initio method-dependent terms in Eqs. 6, 7, 9, and 10. As can
be seen in Table LABEL:tab:NH3_thermodynamic_breakdown, the dominant
contribution to $\Delta G^{\oplus}$ in the NH3 system is the difference in
vibrational entropies of the transition state and the equilibrium state, which
DFT underestimates by $\sim$0.022 eV compared to CCSD. Since the hybrid DFT-
QMC approach directly borrows the vibrational entropy from the DFT
calculation, this explains the high similarity between the DFT and hybrid DFT-
QMC values of $\Delta G^{\oplus}$ and $k_{\rightarrow}$ for NH3 inversion. In
the S${}_{\textrm{N}}$2 system, however, as shown in Table
LABEL:tab:SN2_thermodynamic_breakdown, the dominant contribution now becomes
the difference in potential energies of the transition state and equilibrium
state, which DFT predicts to be $\sim$0.008 eV, whereas DMC and CCSD predict
to be $\sim$0.022 eV. Since DMC computes the potential energies more
accurately than DFT, the hybrid DFT-QMC approach is able to obtain $\Delta
G^{\oplus}$ and $k_{\rightarrow}$ values in much better agreement with the
CCSD result.

##  4 Conclusions

In this work, we have demonstrated the extension of the surrogate Hessian
line-search method previously employed to determine equilibrium geometries to
the identification of transition states and minimum-energy pathways. Modifying
the objective function of the line-search to maximize the energy along the
negative-eigenvalue directions, as shown in Eq. 2, permits us to identify
saddle points of arbitrary order. Furthermore, performing the line-search
within restricted, path-orthogonal subspaces yields points along the MEP that
can be specified to arbitrary structural accuracy, leveraging the controllable
computational cost of DMC and other stochastic ab initio theories. Finally, we
have shown that it is possible to combine stochastic potential energies and
Monte Carlo-sampled rotational entropies obtained using DMC with vibrational
contributions from DFT to compute corrected free energy changes and kinetic
constants that can improve upon DFT results and, potentially, approach coupled
cluster accuracy using only energetic information from DMC calculations
without the need for gradient evaluation at the DMC level.

The surrogate Hessian line-search and subspace optimization methods are both
primarily based on energy optimization. Imposing structural tolerances on the
geometries that are identified using these methods, therefore, directly
influences the precision to which the energy is specified. In practice, the
structural tolerance must be set large enough so as not to incur
disproportionate computational expense while still meeting the level of
precision required for the study. If an accurately designed surrogate Hessian-
based calculation results in structures that show slight deviations from more
accurate methods such as CCSD or CCSD(T), this is likely an indication of the
“flatness” of the PES near those structures, which, in turn, reveals
interesting information about the physics of the system.

Subspace optimization, as noted in Section 3.2, features a few methodological
departures from NEB and other chain-of-states methods. However, both classes
of methods accomplish the same task and identify slightly different - yet
theoretically consistent - points along the MEP. Steps can be taken to
eliminate the differences between the two approaches, such as standardizing
the tangent calculation method, repeated reevaluation of the surrogate Hessian
at every step taken by the line-search, etc. We leave these modifications for
future work.

High-accuracy PES calculations play a significant role in several areas of
computational chemistry and materials science research. At large system sizes,
$\mathcal{O}(N^{6-7})$\- or higher-scaling coupled cluster- and full
configuration interaction-based methods remain largely inaccessible.
Stochastic electronic structure theories such as QMC - accelerated using
surrogate theories such as DFT - present opportunities to unravel the physics
and chemistry of large, correlated systems at high accuracy and for
potentially lower computational cost, given their $\mathcal{O}(N^{3-4})$
scaling. The work presented here also provides a path toward obtaining
critical information about reaction pathways without the explicit use of
energy gradients, which can be challenging to obtain in a variety of otherwise
promising electronic structure methods. 22 In the context of DMC, our approach
will have particular value for large systems or systems containing heavy
atoms, for which force calculations can be especially cumbersome.

Overall, our work represents a step further in the expanding suite of
capabilities of QMC methods, enabling the accurate description of systems that
undergo dynamical transitions along minimum-energy pathways, such as diffusion
processes,61 chemical reactions,44 including catalytic reactions,62, 33
conformational changes, and phase transitions. 63

## Supporting Information

The Supporting Information document contains the structural parameterization
equations for the NH3 and S${}_{\textrm{N}}$2 systems, DMC calculation
extension factors, and the breakdown of the free energy contributions. All
code associated with this work can be found at
https://github.com/gopal-iyer/surrogate_hessian_reaction_pathway_search.

G.R.I. thanks Shubham Sharma, Benjamin Foulon, and Bjarne Kreitz for useful
discussions. G.R.I. (derivation, implementation, and analysis of the
algorithms presented, drafting of the manuscript) was funded by AFOSR Award
Number FA9550-19-1-9999 and the Brown University Chemistry Department Vince
Wernig Fellowship. J.T. (mentorship) was self-supported. J.T.K. (concept,
mentorship), P.R.C.K. (concept), and B.R. (concept, mentorship, manuscript
writing) were supported by the U.S. Department of Energy, Office of Science,
Basic Energy Sciences, Materials Sciences and Engineering Division, as part of
the Computational Materials Sciences Program and Center for Predictive
Simulation of Functional Materials. This research was conducted using
computational resources and services at the Center for Computation and
Visualization, Brown University.

This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-
AC05-00OR22725 with the U.S. Department of Energy. The United States
Government retains and the publisher, by accepting the article for
publication, acknowledges that the United States Government retains a non-
exclusive, paid-up, irrevocable, world-wide license to publish or reproduce
the published form of this manuscript, or allow others to do so, for United
States Government purposes.

## References

  * Cramer 2013 Cramer, C. _Essentials of Computational Chemistry: Theories and Models_ ; Wiley, 2013
  * Czakó et al. 2014 Czakó, G.; Szabó, I.; Telekes, H. On the Choice of the Ab Initio Level of Theory for Potential Energy Surface Developments. _The Journal of Physical Chemistry A_ 2014, _118_ , 646–654
  * Bowman et al. 2011 Bowman, J. M.; Czakó, G.; Fu, B. High-dimensional ab initio potential energy surfaces for reaction dynamics calculations. _Phys. Chem. Chem. Phys._ 2011, _13_ , 8094–8111
  * Bowman et al. 2011 Bowman, J. M.; Czakó, G.; Fu, B. High-dimensional ab initio potential energy surfaces for reaction dynamics calculations. _Phys. Chem. Chem. Phys._ 2011, _13_ , 8094–8111
  * Schlegel 2011 Schlegel, H. B. Geometry optimization. _WIREs Computational Molecular Science_ 2011, _1_ , 790–809
  * Garijo del Río et al. 2019 Garijo del Río, E.; Mortensen, J. J.; Jacobsen, K. W. Local Bayesian optimizer for atomic structures. _Phys. Rev. B_ 2019, _100_ , 104103
  * Bitzek et al. 2006 Bitzek, E.; Koskinen, P.; Gähler, F.; Moseler, M.; Gumbsch, P. Structural Relaxation Made Simple. _Phys. Rev. Lett._ 2006, _97_ , 170201
  * Pratt 1986 Pratt, L. R. A statistical method for identifying transition states in high dimensional problems. _The Journal of Chemical Physics_ 1986, _85_ , 5045–5048
  * Dellago et al. 1998 Dellago, C.; Bolhuis, P. G.; Csajka, F. S.; Chandler, D. Transition path sampling and the calculation of rate constants. _The Journal of Chemical Physics_ 1998, _108_ , 1964–1977
  * 10 JÓNSSON, H.; MILLS, G.; JACOBSEN, K. W. _Classical and Quantum Dynamics in Condensed Phase Simulations_ ; pp 385–404
  * Sheppard et al. 2012 Sheppard, D.; Xiao, P.; Chemelewski, W.; Johnson, D. D.; Henkelman, G. A generalized solid-state nudged elastic band method. _The Journal of Chemical Physics_ 2012, _136_ , 074103
  * Kolsbjerg et al. 2016 Kolsbjerg, E. L.; Groves, M. N.; Hammer, B. An automated nudged elastic band method. _The Journal of Chemical Physics_ 2016, _145_ , 094107
  * Henkelman et al. 2000 Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. _The Journal of Chemical Physics_ 2000, _113_ , 9901–9904
  * Alaghemandi et al. 2017 Alaghemandi, M.; Newcomb, L. B.; Green, J. R. Ignition in an Atomistic Model of Hydrogen Oxidation. _The Journal of Physical Chemistry A_ 2017, _121_ , 1686–1692, PMID: 28169533
  * Thomas 2009 Thomas, J. Handbook Of Heterogeneous Catalysis. 2., completely revised and enlarged Edition. Vol. 1–8. Edited by G. Ertl, H. Knözinger, F. Schüth, and J. Weitkamp. _Angewandte Chemie International Edition_ 2009, _48_
  * Nørskov et al. 2014 Nørskov, J. K.; Studt, F.; Abild-Pedersen, F.; Bligaard, T. _Fundamental Concepts in Heterogeneous Catalysis_ ; John Wiley & Sons, Inc., 2014; First published: 29 August 2014
  * Wang et al. 2014 Wang, L.-P.; Titov, A.; McGibbon, R.; Liu, F.; Pande, V. S.; Martínez, T. J. Discovering chemistry with an ab initio nanoreactor. _Nature Chemistry_ 2014, _6_ , 1044–1048
  * Politzer and Murray 2018 Politzer, P.; Murray, J. S. The Hellmann-Feynman theorem: a perspective. _Journal of Molecular Modeling_ 2018, _24_ , 266
  * Salter et al. 1989 Salter, E. A.; Trucks, G. W.; Bartlett, R. J. Analytic energy derivatives in many‐body methods. I. First derivatives. _The Journal of Chemical Physics_ 1989, _90_ , 1752–1766
  * Stanton and Gauss 2000 Stanton, J. F.; Gauss, J. Analytic second derivatives in high-order many-body perturbation and coupled-cluster theories: Computational considerations and applications. _International Reviews in Physical Chemistry_ 2000, _19_ , 61–95
  * Burow et al. 2014 Burow, A. M.; Bates, J. E.; Furche, F.; Eshuis, H. Analytical First-Order Molecular Properties and Forces within the Adiabatic Connection Random Phase Approximation. _Journal of Chemical Theory and Computation_ 2014, _10_ , 180–194, _eprint: https://doi.org/10.1021/ct4008553
  * Ramberger et al. 2017 Ramberger, B.; Schäfer, T.; Kresse, G. Analytic Interatomic Forces in the Random Phase Approximation. _Phys. Rev. Lett._ 2017, _118_ , 106403
  * Assaraf and Caffarel 2000 Assaraf, R.; Caffarel, M. Computing forces with quantum Monte Carlo. _The Journal of Chemical Physics_ 2000, _113_ , 4028–4034
  * Badinski et al. 2010 Badinski, A.; Haynes, P. D.; Trail, J. R.; Needs, R. J. Methods for calculating forces within quantum Monte Carlo simulations. _Journal of Physics: Condensed Matter_ 2010, _22_ , 074202
  * Chiesa et al. 2005 Chiesa, S.; Ceperley, D. M.; Zhang, S. Accurate, Efficient, and Simple Forces Computed with Quantum Monte Carlo Methods. _Phys. Rev. Lett._ 2005, _94_ , 036404
  * Filippi and Umrigar 2000 Filippi, C.; Umrigar, C. J. Correlated sampling in quantum Monte Carlo: A route to forces. _Phys. Rev. B_ 2000, _61_ , R16291–R16294
  * Assaraf and Caffarel 2003 Assaraf, R.; Caffarel, M. Zero-variance zero-bias principle for observables in quantum Monte Carlo: Application to forces. _The Journal of Chemical Physics_ 2003, _119_ , 10536–10552
  * Ríos and Conduit 2019 Ríos, P. L.; Conduit, G. J. Tail-regression estimator for heavy-tailed distributions of known tail indices and its application to continuum quantum Monte Carlo data. _Phys. Rev. E_ 2019, _99_ , 063312
  * Nakano et al. 2022 Nakano, K.; Raghav, A.; Sorella, S. Space-warp coordinate transformation for efficient ionic force calculations in quantum Monte Carlo. _The Journal of Chemical Physics_ 2022, _156_ , 034101
  * Moroni et al. 2014 Moroni, S.; Saccani, S.; Filippi, C. Practical Schemes for Accurate Forces in Quantum Monte Carlo. _Journal of Chemical Theory and Computation_ 2014, _10_ , 4823–4829, _eprint: https://doi.org/10.1021/ct500780r
  * Jiang et al. 2022 Jiang, T.; Fang, W.; Alavi, A.; Chen, J. General Analytical Nuclear Forces and Molecular Potential Energy Surface from Full Configuration Interaction Quantum Monte Carlo. _Journal of Chemical Theory and Computation_ 2022, _18_ , 7233–7242, PMID: 36326847
  * Chen and Zhang 2022 Chen, S.; Zhang, S. A structural optimization algorithm with stochastic forces and stresses. _Nature Computational Science_ 2022, _2_ , 736–744
  * Iyer and Rubenstein 2022 Iyer, G. R.; Rubenstein, B. M. Finite-Size Error Cancellation in Diffusion Monte Carlo Calculations of Surface Chemistry. _The Journal of Physical Chemistry A_ 2022, _126_ , 4636–4646
  * Staros et al. 2022 Staros, D.; Hu, G.; Tiihonen, J.; Nanguneri, R.; Krogel, J.; Bennett, M. C.; Heinonen, O.; Ganesh, P.; Rubenstein, B. A combined first principles study of the structural, magnetic, and phonon properties of monolayer CrI3. _The Journal of Chemical Physics_ 2022, _156_ , 014707
  * Shin et al. 2021 Shin, H.; Krogel, J. T.; Gasperich, K.; Kent, P. R. C.; Benali, A.; Heinonen, O. Optimized structure and electronic band gap of monolayer GeSe from quantum Monte Carlo methods. _Phys. Rev. Mater._ 2021, _5_ , 024002
  * Ryczko et al. 2022 Ryczko, K.; Krogel, J. T.; Tamblyn, I. Machine Learning Diffusion Monte Carlo Energies. _Journal of Chemical Theory and Computation_ 2022, _18_ , 7695–7701, PMID: 36317712
  * Huang and Rubenstein 2023 Huang, C.; Rubenstein, B. M. Machine Learning Diffusion Monte Carlo Forces. _The Journal of Physical Chemistry A_ 2023, _127_ , 339–355, PMID: 36576803
  * Chen et al. 2023 Chen, M. S.; Lee, J.; Ye, H.-Z.; Berkelbach, T. C.; Reichman, D. R.; Markland, T. E. Data-Efficient Machine Learning Potentials from Transfer Learning of Periodic Correlated Electronic Structure Methods: Liquid Water at AFQMC, CCSD, and CCSD(T) Accuracy. _Journal of Chemical Theory and Computation_ 2023, _19_ , 4510–4519, PMID: 36730728
  * Archibald et al. 2018 Archibald, R.; Krogel, J. T.; Kent, P. R. C. Gaussian process based optimization of molecular geometries using statistically sampled energy surfaces from quantum Monte Carlo. _The Journal of Chemical Physics_ 2018, _149_ , 164116
  * Ceperley et al. 2023 Ceperley, D. M.; Jensen, S.; Yang, Y.; Niu, H.; Pierleoni, C.; Holzmann, M. Training models using forces computed by stochastic electronic structure methods. 2023
  * Grossman and Mitas 2005 Grossman, J. C.; Mitas, L. Efficient Quantum Monte Carlo Energies for Molecular Dynamics Simulations. _Phys. Rev. Lett._ 2005, _94_ , 056403
  * Tiihonen et al. 2022 Tiihonen, J.; Kent, P. R. C.; Krogel, J. T. Surrogate Hessian accelerated structural optimization for stochastic electronic structure theories. _The Journal of Chemical Physics_ 2022, _156_ , 054104
  * Samanta and E 2013 Samanta, A.; E, W. Optimization-Based String Method for Finding Minimum Energy Path. _Communications in Computational Physics_ 2013, _14_ , 265–275
  * Saccani et al. 2013 Saccani, S.; Filippi, C.; Moroni, S. Minimum energy pathways via quantum Monte Carlo. _The Journal of Chemical Physics_ 2013, _138_ , 084109
  * Gonzales et al. 2001 Gonzales, J. M.; Cox, R. S.; Brown, S. T.; Allen, W. D.; Schaefer, H. F. Assessment of Density Functional Theory for Model SN2 Reactions: CH3X + F- (X = F, Cl, CN, OH, SH, NH2, PH2). _The Journal of Physical Chemistry A_ 2001, _105_ , 11327–11346
  * Foulkes et al. 2001 Foulkes, W. M. C.; Mitas, L.; Needs, R. J.; Rajagopal, G. Quantum Monte Carlo simulations of solids. _Rev. Mod. Phys._ 2001, _73_ , 33–83
  * Ceperley 2010 Ceperley, D. M. In _Theoretical and Computational Methods in Mineral Physics_ ; Wentzcovitch, R. M., Stixrude, L., Eds.; De Gruyter: Berlin, Boston, 2010; pp 129–136
  * Fukui 1970 Fukui, K. Formulation of the reaction coordinate. _The Journal of Physical Chemistry_ 1970, _74_ , 4161–4163
  * Motagamwala and Dumesic 2021 Motagamwala, A. H.; Dumesic, J. A. Microkinetic Modeling: A Tool for Rational Catalyst Design. _Chemical Reviews_ 2021, _121_ , 1049–1076
  * Anderson 1972 Anderson, P. W. More Is Different. _Science_ 1972, _177_ , 393–396
  * Sun et al. 2018 Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.; McClain, J. D.; Sayfutyarova, E. R.; Sharma, S.; Wouters, S.; Chan, G. K.-L. PySCF: the Python-based simulations of chemistry framework. _WIREs Computational Molecular Science_ 2018, _8_ , e1340
  * Sun et al. 2020 Sun, Q. et al. Recent developments in the PySCF program package. _The Journal of Chemical Physics_ 2020, _153_ , 024109
  * Annaberdiyev et al. 2018 Annaberdiyev, A.; Wang, G.; Melton, C. A.; Bennett, M. C.; Shulenburger, L.; Mitas, L. A new generation of effective core potentials from correlated calculations: 3d transition metal series. _The Journal of Chemical Physics_ 2018, _149_ , 134108
  * Wang et al. 2019 Wang, G.; Annaberdiyev, A.; Melton, C. A.; Bennett, M. C.; Shulenburger, L.; Mitas, L. A new generation of effective core potentials from correlated calculations: 4s and 4p main group elements and first row additions. _The Journal of Chemical Physics_ 2019, _151_ , 144110
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 1996, _77_ , 3865–3868
  * Kim et al. 2018 Kim, J. et al. QMCPACK: an open source ab initio quantum Monte Carlo package for the electronic structure of atoms, molecules and solids. _Journal of Physics: Condensed Matter_ 2018, _30_ , 195901
  * Kent et al. 2020 Kent, P. R. C. et al. QMCPACK: Advances in the development, efficiency, and application of auxiliary field and real-space variational and diffusion quantum Monte Carlo. _The Journal of Chemical Physics_ 2020, _152_ , 174105
  * Krogel 2016 Krogel, J. T. Nexus: A modular workflow management system for quantum simulation codes. _Computer Physics Communications_ 2016, _198_ , 154–168
  * Larsen et al. 2017 Larsen, A. H. et al. The atomic simulation environment—a Python library for working with atoms. _Journal of Physics: Condensed Matter_ 2017, _29_ , 273002
  * Zhu et al. 2019 Zhu, X.; Thompson, K. C.; Martínez, T. J. Geodesic interpolation for reaction pathways. _The Journal of Chemical Physics_ 2019, _150_ , 164103
  * Ma et al. 2011 Ma, J.; Michaelides, A.; Alfè, D.; Schimka, L.; Kresse, G.; Wang, E. Adsorption and diffusion of water on graphene from first principles. _Phys. Rev. B_ 2011, _84_ , 033402
  * Shi et al. 2023 Shi, B. X.; Zen, A.; Kapil, V.; Nagy, P. R.; Grüneis, A.; Michaelides, A. Many-Body Methods for Surface Chemistry Come of Age: Achieving Consensus with Experiments. _Journal of the American Chemical Society_ 2023, _145_ , 25372–25381, PMID: 37948071
  * Wines et al. 2023 Wines, D.; Tiihonen, J.; Saritas, K.; Krogel, J. T.; Ataca, C. A Quantum Monte Carlo Study of the Structural, Energetic, and Magnetic Properties of Two-Dimensional H and T Phase VSe2. _The Journal of Physical Chemistry Letters_ 2023, _14_ , 3553–3560, PMID: 37017431