[1,2]Yannic Rath 1]National Physical Laboratory, Teddington, TW11 0LW, United
Kingdom

[2]George H. Booth

2]Department of Physics and Thomas Young Centre, King’s College London,
Strand, London, WC2R 2LS, United Kingdom

# Interpolating many-body wave functions for accelerated molecular dynamics on
near-exact electronic surfaces

<EMAIL_ADDRESS>[<EMAIL_ADDRESS>[

###### Abstract

While there have been many developments in computational probes of both
strongly-correlated molecular systems and machine-learning accelerated
molecular dynamics, there remains a significant gap in capabilities between
them, where it is necessary to describe the accurate electronic structure over
timescales in which atoms move. We describe a practical approach to bridge
these fields by interpolating the correlated many-electron state through
chemical space, whilst avoiding the exponential complexity of these underlying
states. With a small number of accurate correlated wave function calculations
as a training set, we demonstrate provable convergence to near-exact potential
energy surfaces for subsequent dynamics with propagation of a valid many-body
wave function and inference of its variational energy at all points, whilst
retaining a mean-field computational scaling. This represents a profoundly
different paradigm to the direct interpolation of properties through chemical
space in established machine-learning approaches. It benefits from access to
all electronic properties of interest from the same model without relying on
local descriptors, and demonstrates improved performance compared to the
direct training on energies themselves. We combine this with modern
systematically-improvable electronic structure methods to resolve the
molecular dynamics for a number of correlated electron problems, including the
proton dynamics of a Zundel cation trajectory, where we highlight the
qualitative improvement from traditional machine learning or ab initio
dynamics on mean-field surfaces.

## Introduction

The quantum fluctuations of interacting electrons represent the critical
interaction between atoms which underpin all atomic bonding, dynamics and
reactivity. Computational approaches for systems with strongly interacting
electrons have undergone a number of major developments in recent decades, as
emerging methods enable a description of correlated electronic structure for
ever larger and more realistic systems with unprecedented accuracy
[PhysRevX.7.031059, PhysRevX.10.011041, doi:10.1021/acs.jpclett.0c02621].
These modern approaches across both chemical and materials science include
those based on tensor networks [10.1063/1.4905329, 10.1063/5.0050902,
Larsson2022TheCD, doi:10.1126/science.abm2295, Li2019], stochastic methods
[Booth2012TowardsAE, 10.1063/5.0150706, Ren2023,
doi:10.1021/acs.jctc.1c01162], selected configuration interaction
[KOCH1993193, 10.1063/1.5063376, 10.1063/1.1679199] and machine-learning
inspired wave function ansatze [Choo_2020,
https://doi.org/10.48550/arxiv.2109.12606, zhaoScalableNeuralQuantum2022,
https://doi.org/10.48550/arxiv.2301.03755, Yang2020, PhysRevB.107.205119,
bortone2023impact, PhysRevResearch.2.033429, Hermann2019DeepneuralnetworkSO,
gerard2022goldstandard]. These approaches have allowed for the near-exact
solution to the quantum many-body problem in these systems, providing high-
accuracy insights for a few fixed atomic configurations, but have in general
had little or no impact on our understanding of the physics and chemistry on
the timescales of atomic and molecular motion.

The reasons for this are obvious; while a small number of single point
calculations with fixed atoms are possible, the different timescales of atomic
dynamics and electronic quantum fluctuations mean that on the order of at
least thousands of sequential electronic structure calculations are required.
This is essential to propagate the atoms in molecular systems to relevant
timescales, entailing generally prohibitive computation expense for these
high-accuracy methods. This is particularly challenging for these emerging
methods which can lack a ‘black-box’ use, requiring care to ensure reliable
convergence at each point, while often also lacking analytic atomic forces to
propagate the nuclear coordinates in time [doi:10.1021/acs.chemrev.9b00496].
For this reason ‘ab-initio Born-Oppenheimer molecular dynamics’ (AI-BOMD),
where the atoms are classically propagated according to the potential energy
surface of the electrons, is almost synonymous with a more empirical density
functional description of the electronic structure (with exceptions
[doi:10.1021/jz401931f, 10.3389/fmats.2015.00029, 10.1063/1.4941091]), which
lacks systematic improvability and has many well documented deficiencies
[marx2000ab]. These include an over-stabilisation of delocalized electronic
states, as well as often inaccurate descriptions of dispersion forces,
transition states or bond breaking among others [Cohen2012ChallengesFD]. These
are critical parts of the phase space in real chemical dynamics, and the acute
need for more reliable potential energy surfaces which build on the
developments in accurate electronic structure is clear.

The most widespread and successful resolution to this need has come from the
‘machine learning’ of force fields [Kabylda2023EfficientID, Fu2022ForcesAN,
Unke2020MachineLF, Gkeka2020MachineLF]. These interpolate across chemical
space between accurate single-point estimates of the electronic energy, based
on local descriptors of the environment of each atom
[doi:10.1021/acs.chemrev.1c00021]. While this approach to straddling the
electronic and atomic timescales has been arguably one of the most successful
contributions of machine learning to quantum-level simulations to date, it is
not without its own drawbacks [Kulik_2022]. In particular, the local nature of
the descriptors can lead to difficulty describing long-range interactions
[doi:10.1021/acs.jctc.3c00704], as well as ‘holes’ where non-variational
inferred energy estimates can lead to a collapse in the statistical sampling
of phase space to these unphysical minima [doi:10.1021/acs.jctc.1c00647]. On a
more fundamental level, since these approaches integrate out the electronic
structure, there is in general no fundamental electronic variable at each
sampled point (such as the wave function), meaning that the electronic
properties which can be extracted are limited to the ones which correspond to
the model definition. If the evolution of e.g. the dipole moment or charge
distribution across a trajectory was desired this would not be accessible from
a force field, and extensions to non-adiabatic effects are also far from
straightforward in this framework, noting however significant recent research
in these directions [Schtt2019UnifyingML, Fedik2022ExtendingML,
doi:10.1021/acs.jpclett.0c00527].

We take a different perspective and show that rather than interpolating
observables such as the potential energy, we can instead interpolate the many-
body electronic wave function itself through chemical space. Importantly,
despite the many-body wave function of each training point being in general
exponentially complex, inference of properties from the model can be achieved
in a scaling which is the same as (hybrid) density functional theory,
rendering this a practical scheme. This decouples the unfavourable scaling of
high accuracy single point electronic structure calculations from the
evaluation of the interpolated potential energy surface, and thus allows for
the use of these electronic structure methods for atomic dynamics on realistic
timescales. We show that the resulting potential energy surfaces for molecular
dynamics over relevant timescales are systematically improvable to near-
exactness via interpolation between highly accurate training configurations.
Since a valid correlated many-electron state is propagated through chemical
space, this paradigm enables all electronic properties of interest to be
simultaneously accessible within the same model, without relying on local or
low-rank descriptors. Furthermore, since the energy is computed as a rigorous
quantum expectation value over this inferred state it provides a fully
variational potential energy estimate (precluding ‘holes’) for all atomic
configurations, allows for clear evidence of systematic improvability to
exactness as the training set is enlarged, an inductive bias of the model away
from poorly described regions of phase space, and simple access to analytic
atomic forces of the model for efficient propagation of dynamics.

We combine this approach for interpolating wave functions with modern density
matrix renormalization group (DMRG) methods, allowing convergence of the
strongly correlated potential energy surfaces to near exactness
[10.1063/1.4905329, 10.1063/5.0050902]. We demonstrate that this can provide a
fully correlated electronic description of reactive molecular dynamics beyond
traditional parameterized or machine-learned force fields, over time scales
which would be inaccessible without this acceleration. We show this can result
in qualitative differences in behavior for a number of proto-typical molecular
dynamics simulations compared to both density functional and traditional
machine-learned force field approaches [10.1063/5.0160898]. Finally, we
consider a high-energy reaction trajectory on a near-exact potential surface
for the Zundel cation, a key intermediate for the Grotthuss mechanism for
hydrogen diffusion through aqueous solutions [Agmon1995TheGM,
Cukierman2006EtTG, 10.1063/1.4941091]. We find a hydrogen atom ejected from
the system over the order of a hundred femtoseconds, contradicting widely
employed density function results, and demonstrating a marked difference in
the charge distribution over the ions.

## Interpolating wave functions

We first consider how to interpolate a single many-electron wave function
between two different atomic configurations. We assume that we have an exact
(FCI) correlated many-electron state defined within an atom-centered basis set
of $L$ functions, for a specific set of atomic coordinates ${\bf R}$
[KNOWLES198975]. This wave function is a linear superposition over
exponentially many electron configurations (Slater determinants) spanning the
Hilbert space, as

$|\Psi({\bf
R})\rangle=\sum_{n_{1},n_{2},\dots,n_{L}}C_{n_{1},n_{2},\dots,n_{L}}|n_{1},n_{2},\dots,n_{L}\rangle,$
(1)

where $C_{n_{1},n_{2},\dots,n_{L}}$ is the rank-$L$ tensor of probability
amplitudes over the electronic configurations, and $n_{i}$ indexes the four
local Fock states of the $i^{\textrm{th}}$ orbital; either unoccupied, spin-
up, spin-down, or doubly occupied with electrons for each orbital. In general,
both the probability amplitudes, and the single-particle orbitals defining
each electronic configuration $|n_{1},n_{2},\dots,n_{L}\rangle$ will change
with atomic configuration ${\bf R}$. However, we aim to represent an
approximation to the correlated electronic state at a different chemical
configuration (and therefore electronic Hilbert space) with the same tensor of
probability amplitudes over these electronic states. We exploit the fact that
the properties of the exact state are invariant to orthogonal rotations of the
single-particle orbitals, but that the probability amplitudes themselves will
vary with this choice. Therefore, to enable transferrability between chemical
environments, we seek a choice of orbital representation in which the
probability amplitudes of the exact many-electron state change least between
atomic configurations of interest.

A plausible choice is a basis of local atomic-like functions, appealing to the
fact that a large portion of the electronic fluctuations among atomic-local
orbitals will remain qualitatively similar as atoms are moved by small
amounts. Similarly, regions of similar chemical bonding will also have common
features in their probability amplitudes defining e.g. covalent fluctuations
between neighboring atoms [mejutozaera2023quantum]. However, for reasons which
will become clear, we also require that the orbitals represent an orthonormal
set for all atomic configurations. To ensure this, while (in a least-squares
sense) optimally preserving this atomic-like character of the orbitals, we
symmetrically (Löwdin) orthonormalize the atomic-orbital basis (see Methods),
defining orthonormal ‘SAO’ orbital sets for each atomic configuration
[10.1063/1.1747632, https://doi.org/10.1002/qua.981]. We can then choose to
interpolate the state (and all resulting properties) between atomic
configurations without re-optimizing the many-electron state by simply
transferring the probability amplitudes, while ensuring the consistent SAO
basis definition.

This simple approach is limited, since the many-body amplitudes will in
general change as the atoms move. However, we can generalize the state while
retaining a valid wave function by linearly combining probability amplitude
tensors in this transferable SAO representation from a larger ‘training’ set
optimized at other atomic configurations. We then variationally optimize the
relative contributions of each of the $N$ training states for any test atomic
configuration. This is achieved in closed form as the diagonalization of a
generalized eigenvalue problem in the basis of the training states. By
projecting the Hamiltonian at the desired test geometry $\mathbf{R}$ into this
many-body basis, we get

$\boldsymbol{\mathcal{H}}(\mathbf{R})\mathbf{X}(\mathbf{R})=\mathbf{E}(\mathbf{R})\,\boldsymbol{\mathcal{S}}\mathbf{X}(\mathbf{R}),$
(2)

with the eigenvectors, $\mathbf{X}(\mathbf{R})$, giving the amplitudes of the
training states defining the interpolated wave functions at the test geometry,
with inferred energy spectrum $\mathbf{E}(\mathbf{R})$. The electronic
Hamiltonian of the test geometry, $\boldsymbol{\mathcal{H}}(\mathbf{R})$, is
projected into a basis defined by the fixed probability amplitudes of the
many-body training states, written as

$\displaystyle\mathcal{H}_{ab}$
$\displaystyle=\sum_{ijkl}\sum_{\mathbf{n}\mathbf{n}^{\prime}}C_{\mathbf{n}}^{(a)*}C_{\mathbf{n}^{\prime}}^{(b)}\langle\mathbf{n}|{\hat{c}}^{\dagger}_{i}{\hat{c}}^{\dagger}_{j}{\hat{c}}_{l}{\hat{c}}_{k}|\mathbf{n}^{\prime}\rangle\,K_{ijkl}(\mathbf{R})$
$\displaystyle=\sum_{ijkl}\Gamma^{ijkl}_{ab}\,K_{ijkl}(\mathbf{R}),$ (3)

where $\mathbf{n}$ denotes the many-electron configurations in the SAO basis
of the test geometry, with
$|\mathbf{n}\rangle\equiv|n_{1},n_{2},\dots,n_{L}\rangle$, and with
$C_{\mathbf{n}}^{(a)}$ and $C_{\mathbf{n}}^{(b)}$ the fixed SAO probability
amplitudes of the training points at atomic geometries $a$ and $b$
respectively. This definition therefore implicitly transfers the probability
amplitudes between the Hilbert spaces of the training and test geometries. The
$K_{ijkl}(\mathbf{R})$ tensor is the two-electron reduced Hamiltonian defined
in the SAO basis of the test geometry (given explicitly in the Methods
section) with second-quantized Fermionic operators acting in this basis as
${\hat{c}}^{(\dagger)}$ shown.

Since the Hamiltonian is a sum of only two-electron interactions, the
contraction over the exponential many-electron configurations, $\mathbf{n}$,
is performed for all pairs of training states to give the transition two-body
density matrices $\boldsymbol{\Gamma}^{ijkl}$ defined in Eq. 3. Crucially,
since the training probability amplitudes are defined not to change with test
geometry, this contraction is only performed once on the training states, not
for each test geometry. The construction of the subspace Hamiltonian at a test
geometry therefore only requires the $\mathcal{O}[L^{4}]$ contraction of Eq. 3
– the same scaling as traditional (hybrid) density functional theory, with
only $K_{ijkl}(\mathbf{R})$ changing with test configuration. This scaling for
evaluation of the model at test geometries could also be further lowered with
factorizations exploiting the low-rank nature of the tensors in Eq. (3)
[doi:10.1021/acs.jctc.7b00605, 10.1063/1.5047207,
doi:10.1021/acs.jctc.8b00780, doi:10.1021/acs.jctc.3c00851]. The overlap
between the training states, $\boldsymbol{\mathcal{S}}$ of Eq. 2, can
similarly be precomputed during the training phase, as

$\mathcal{S}_{ab}=\sum_{\mathbf{n}}C_{\mathbf{n}}^{(a)*}C_{\mathbf{n}}^{(b)},$
(4)

due to the orthonormality of the SAO basis at all test geometries. This
ensures that the overlaps of the training states do not change with geometry,
despite the physical training wave functions changing with atomic
rearrangements as they transform between Hilbert spaces. In this way, the
exponential complexity of the many-electron states are completely avoided in
the inference of wave functions at new test points, by representing the
training states in the polynomially-compact tensors
$\boldsymbol{\Gamma}^{ijkl}$ and their overlaps.

The lowest energy eigenvector from the diagonalization of this Hamiltonian
(whose dimensionality scales only with the number of training points and is
independent of system size) defines the specific variationally optimal linear
combination of training probability amplitudes for the state, which can
subsequently be used to predict any electronic property at this test geometry.
Due to the variationality we have the desirable properties that the inferred
state at a point which coincides with a training geometry must necessarily be
exact, as well as the fact that each additional training point must
necessarily lower the inferred energy towards the exact electronic solution
across all possible test atomic configurations, assuming linear independence.
In this way, the method more closely resembles a reduced order method than a
machine learning model, yet both are useful viewpoints. Due to the fixed
training amplitudes across geometries, as well as the variational optimization
of the model, computing analytic atomic forces from the inferred state is also
straightforward via the Hellmann-Feynman theorem (see Methods), ensuring a
particular relevance of this acceleration in molecular dynamics applications.

Figure 1: Schematic overview of the proposed ‘Eigenvector Continuation’ scheme
for the high accuracy prediction across chemical space from few _ab initio_
data points. Utilizing training states in a Löwdin atomic orbital (SAO) basis,
states are predicted for new molecular geometries by diagonalization in the
subspace of the training states, allowing for a fast variational prediction of
highly accurate potential energy surfaces and force fields via the compact
one- and two-body transition density matrices between states. The scheme is
exemplified for symmetrically stretched chains of $6$ hydrogen atoms in a
minimal basis set (STO-6G) using one, two, and three training points.
$\boldsymbol{\mathcal{S}}$ and $\boldsymbol{\Gamma}^{ijkl}$ represent the
geometry-agnostic overlap and two-body transition density matrices between
training states. $\langle\mathbf{n}|\Psi^{(a)}\rangle$ shows the value of the
(exponentially many) wavefunction amplitudes in the SAO basis for the
$a^{\mathrm{th}}$ training state, and the test geometry shows the combination
of training states used to define its own SAO probability amplitudes.

This approach also builds on the perspective of ‘eigenvector continuation’
which was recently introduced in both nuclear physics and condensed matter
lattice models, where an eigenstate is analytically continued to different
parts of the phase diagram [PhysRevLett.121.032501, PhysRevC.101.041302,
duguet2023eigenvector, KONIG2020135814, PhysRevC.107.064316,
DRISCHLER2021136777, 10.1063/5.0141145]. Even more recently this was extended
to simple ab initio quantum chemistry applications, with a related scheme to
the one proposed [mejutozaera2023quantum]. However, a crucial difference was
the use of a non-orthogonal atomic basis, which necessitated evaluation of the
test point Hamiltonian directly from the many-body states. This retained the
exponential complexity of the many-body state for inference at each test
geometry, which is avoided here. The method of Mejuto-Zaera et. al. was
therefore presented instead as an approach for quantum computers, where many-
body unitary operations can be applied in polynomial complexity. In contrast,
the SAO basis for the interpolation formally breaks this requirement, ensuring
the approach is amenable to classical computation in the predictions at test
points with tractable mean-field computational cost.

A simple example of the scheme is shown in Fig. 1, for the symmetric stretch
of a chain of six hydrogen atoms, with up to three atomic displacements
considered in the training set. The compressed intermediate representation of
the overlaps and transition density matrices between the training states in
the SAO basis are shown, enabling variationally optimal predictions across the
whole potential energy surface as linear combinations of the many-body
training basis transformed between geometries. The predictions are found to
converge to near exactness for this system with only three training points,
with a guarantee of smoothness on adiabatic surfaces, and exactness at any
training point geometry.

This approach is invariant to translation and rigid body rotations of the
chemical system, provided a consistent ordering of the SAO representation is
maintained, which is straightforward to achieve. However, this is not trivial
for atomic permutations or point group symmetries which would change the SAO
ordering, and hence probability amplitudes of the state definition.
Furthermore, in contrast with building a force field based on local
descriptors, the inference requires the same dimensionality Hilbert space for
the electronic state, necessitating that the training and prediction points
are taken from the same sized system. Future work will look to relax this
constraint.

Figure 2: Mean relative energy error of the prediction for distorted ten-atom
hydrogen chains against the absolute displacement of each atom from the
equilibrium position. For each realization, a distorted chain was created by
moving each atom from their position in the equilibrium geometry by a fixed
displacement with a random direction. The comparison includes predictions from
Hartree-Fock (red, dotted), the Gaussian approximation potential framework
(green, dashed), as well as the variational continuation scheme from $5$
training states of the symmetrically stretched chain (blue, solid). Each data
point corresponds to the mean over 1000 randomly generated geometries. The
inset shows the mean squared force error obtained with the three methods,
where the shaded area denotes the range of the errors over the random
realizations. The training set of equidistant one-dimensional geometries
include the equilibrium length, with an interatomic distance of $\approx
1.79\,a_{0}$, as well as the $4$ symmetric stretches of the atoms where the
inter-atomic distance was increased and decreased by $0.5\,a_{0}$ and
$1\,a_{0}$.

While the proof-of-principle in Fig. 1 demonstrates excellent accuracy with
few training points, it is also interpolation within a simple one-dimensional
phase space of geometries. We now compare to a far larger phase space,
composed of averaging the errors in both the inferred energy and analytic
forces on the atoms over randomly oriented three-dimensional displacements of
each atom from a ten-atom linear hydrogen chain. This provides an
exponentially large phase space of distorted chain configurations to test,
where the radius of the displacements of each atom can be used to control the
magnitude of the geometric distortions from the parent linear chain. Only five
training points from the symmetric stretch of the equidistant linear chain are
used. We consider the increase in error as the magnitude of the displacements
are increased in Fig. 2, as the test configurations move further from these
training samples. We also compare these errors to a Gaussian approximated
potential (GAP); a widely used machine-learning approach based on Gaussian
process regression in a space of local descriptors from the superposition of
atomic potentials [PhysRevLett.104.136403, 10.1063/5.0160898]. This models a
force field directly from the same training energies, but results in a
materially larger error for the energy and forces over all displacements. We
note that five points would generally be a very small training set for GAP,
and that improved techniques to directly train on the forces of the training
data themselves or improved model definitions were not used
[PhysRevB.95.214302]. Nevertheless, a demonstration that inferring the wave
function amplitudes themselves can outperform traditional machine-learning
inference of the properties directly is noteworthy. Furthermore, we compare to
Hartree–Fock theory (HF), which neglects all correlated electron effects and
has the same computational scaling as the inference of the proposed
‘eigenvector continuation’ scheme. This is also significantly worse at small
distortions of the chain, though outperforms the largest distortions which are
far from the training geometries and deep in the extrapolation regime.

## Bridging timescales

While it is easy to envisage many applications of an interpolation scheme for
accurate correlated electronic structure, an obvious target is Born-
Oppenheimer molecular dynamics (MD) [marx2000ab]. In particular, the
variationality of the scheme allows for systematic and quantifiable
improvability to the exact solution of the electronic Schrödinger equation in
the inferred potential energy surface at each geometry, while retaining a
mean-field scaling with respect to the timescales which can be accessed. To
access larger systems and basis sizes we also turn to modern electronic
structure approximations for the evaluation of training states. In particular,
we use the density matrix renormalization group (DMRG) to obtain training
states with controllable accuracy to exactness [10.1063/1.4905329,
10.1063/5.0050902]. These DMRG calculations can either be performed directly
in the SAO basis or the state rotated into this basis after optimization, in
advance of computation of the required overlaps and transition density
matrices between the training states. As an alternative to DMRG, we are also
able to approximate the training states by restricting the space of
correlations to a low-energy complete active subspace (CAS) selected from the
low-energy orbitals of a mean-field calculation [10.1063/5.0042147].

Figure 3: Predicted frequency of the $a_{1}$ symmetric stretch of a water
molecule obtained from MD trajectories starting from a stretched initial
configuration, predicted with different numbers of training data geometries
for the continuation. We simulate the system in increasingly large 6-31G
(blue, dashed), cc-pVDZ (green, dotted) and cc-pVTZ (red, dash-dotted) basis
sets where the larger two bases use training data restricted to a CAS space of
$4$ electrons in $8$ Hartree–Fock orbitals. Horizontal lines give reference
values from MD results on a FCI surface in the 6-31G basis, and CAS
simulations in the cc-pVDZ and cc-pVTZ basis. The inset shows the oxygen-
hydrogen distance over the MD trajectory in the 6-31G basis, as obtained from
continuation with $N=6$ (blue, solid) and $N=2$ (red, dotted) training points,
as well as the reference trajectory from FCI (black, dashed).

We consider these approaches for constructing training states and the
subsequent MD of a water molecule in increasing basis sets in Fig. 3. In
particular, we consider convergence of the predicted vibrational frequency of
the $a_{1}$ symmetric stretching mode as the number of training points
increases. For the smallest basis, we find the full vibrational dynamics
converge with just three training points, where we can compare directly to
exact FCI calculations of the dynamics. As we increase the basis, FCI is
intractable and we restrict the training to a CAS of low-energy orbitals,
where the number of training points required grows modestly to seven and
thirteen training points in a cc-pVDZ and cc-pVTZ basis respectively. While we
compare to CAS-CI for these larger bases for a reference estimate, we note
that the inferred energies are necessarily variationally below these
approximate training states, as shown in the extended data.

It was found important to develop an active learning scheme for the selection
of appropriate atomic configurations to include in the training data for rapid
convergence. In Ref. [mejutozaera2023quantum], the energy variance was
motivated as an appropriate measure for the inclusion of data points, however
this is impractical in the current lower-scaling scheme as it would require
the evaluation of higher-body transition density matrices between training
states. Instead, we consider the addition of training points which will
maximize the improvement in the MD trajectories while respecting the
invariances in the model predictions. This is performed by selecting the point
on the trajectory where the Hamiltonian operator in the SAO basis
$K_{ijkl}(\mathbf{R})$ has changed most (in a least squares sense) compared to
the Hamiltonians used in the current training data set. Since the probability
amplitudes are uniquely defined by this Hamiltonian, it is a suitable measure
for the addition of new data points. Furthermore, due to the variationality of
the method, it is guaranteed that the potential energy with the enlarged
training data will be equal or lower to the previous predictions, across the
whole trajectory. This can therefore be used as a rigorous metric for the
systematic convergence of the potential energy surface for the MD, with more
details in the Methods section. We consider the potential energy surface over
the whole MD simulation fully converged when the maximum reduction in energy
for any point over the whole trajectory is less than $1\,mE_{\mathrm{h}}$ for
two consecutive increases in the data set size.

Figure 4: Molecular dynamics of a $30$ atom hydrogen chain in a STO-6G basis
from an initially symmetrically stretched geometry. The panels report the
distance between two of the hydrogen atoms (top panel: first and last atom in
the chain, bottom panel: first and second atom) against the simulation time,
showing that the first two hydrogen atoms form a stable vibrating dimer while
the chain lengthens. The trajectories were obtained from the eigenvector
continuation with $N=47$ (blue, solid), $N=10$ (green, dashed) and $N=2$ (red,
dotted) DMRG data points, together with the trajectories from DMRG (black,
dashed), Hartree-Fock (grey, dashed) and density functional theory with a PBE
exchange correlation (grey, dash-dotted) potential energy surfaces. Additional
snapshots depict the initial and final hydrogen chain arrangements obtained
from the converged eigenvector continuation ($N=47$).

A semi-infinite symmetric one-dimensional chain of hydrogen atoms has emerged
as a paradigmatic benchmark system of strongly correlated electronic structure
in recent years, as a platform towards larger ab initio and extended systems.
Almost all modern electronic structure methods have been applied to the system
with varying success, and it has motivated further developments in both theory
and understanding of its unexpectedly rich phase diagram [PhysRevX.7.031059,
PhysRevX.10.031058]. While the symmetric stretch of this system has been
considered extensively via single-point electronic structure, its full
dynamics at this level has not. In Fig. 4, we release the atoms to dynamically
move on a tightly-converged ground state surface of the DMRG-trained
continuation scheme, starting from a $\sim$10% symmetric stretching of thirty
atoms equally from the symmetric equilibrium structure. We find that along
with the vibrations of the bonds, the atoms rapidly dimerize and separate,
with the overall length of the chain increasing approximately linearly with
time. We are able to converge the dynamics of this dimerization and
dissociation exactly (albeit in a minimal basis) with only a small number of
DMRG calculations. We note that in comparison, DFT-based AI-BOMD significantly
underestimates the rate of dimerization of the chain, while Hartree–Fock
theory conversely results in a bond for the hydrogen dimers which is too
stiff, demonstrating the importance of an accurate treatment of the electronic
correlations in the dynamics.

## Realistic chemistry on the exact electronic surface

We consider the feasibility of converging faithful reactive chemistry on a
near-exact potential energy surface for the gas-phase dynamics of a Zundel
cation, comprising a water molecule and hydronium ion - a system whose
intricate potential energy surface poses a challenging test case for novel
numerical techniques, yet is particularly important in the understanding of
proton diffusion in aqueous solution [Huang2005AbIP,
McCoy2005FulldimensionalVC, KondatiNatarajan2015RepresentingTP,
Heindel2018BenchmarkES, 10.1063/1.4941091]. We start the simulation from a
high-energy configuration away from the traditional Grotthuss mechanism, where
the additional proton bridges the two ions, interrupting the hydrogen bond as
shown in Fig. 5. We rigorously converge the potential energy surface of the MD
with the continuation scheme based on DMRG training data, which required 84
training points to cover the $120\,fs$ of the simulation (see extended data
for evidence of this convergence with training data over the trajectory). Due
to the fact that an explicit representation of the electronic state is
retained over the trajectory, we also extract non-energetic electronic
properties of the system over time. We use this to consider the evolution of
the Mulliken charge as the electron density is redistributed around the system
in response to the atomic motion, beyond the physics considered in traditional
polarizable force fields.

Figure 5: Predicted Mulliken charges for the hydronium (top panel) and water
(bottom panel) sub-units from an MD simulation of the reaction from an
eigenvector continuation with $N=40$ (green, dashed) and $N=84$ (blue, solid)
DMRG data points. The system uses a 6-31G basis, with snapshots depicting the
evolution of the molecular geometry at four evenly spaced times. Reference
charges from simulations with DFT with a B3LYP exchange correlation (dash-
dotted) and HF (dotted) are included for comparison.

The positive charge is initially fairly evenly distributed amongst the sub-
units, with the anticipation being that the water would rotate to adopt a
lower-energy configuration. However, on the DMRG interpolated potential energy
surface we find that before this is able to occur, a (neutral) hydrogen is
ejected from the system leaving a bound hydronium and hydroxide radical
system. The charges on the sub-units of this reaction show the redistribution
of charge as the hydrogen oscillates a number of times before its eventual
ejection from the system. This behavior is not seen on the more approximate HF
or DFT electronic surfaces where the additional hydrogen remains bound over
timescales enabling the water to rotate its orientation, with the HF charge
distribution substantially in error even in the initial state.

While these results describe the converged dynamics of a single trajectory in
a relatively restricted 6-31G basis, more consideration of a statistical
ensemble of initial states, as well as potentially basis size, thermalization
at finite temperature, nuclear quantum and solvent effects may be needed
before predictions as to the nature of this physical reaction can be given
with confidence [10.1063/1.4941091]. However, the fact that qualitative
changes in dynamical behavior already result from the quality of the treatment
of electronic correlation effects in the determination of the potential energy
surface underlines the importance of a robust and systematically improvable
approach to compute this surface. The eigenvector continuation acceleration
allows computation of this surface with high-level quantum chemical methods,
and extends their scope to enable them to access timescales of atomic dynamics
with provable convergence.

## Outlook

We develop a practical approach for eigenvector continuation of many-body
electronic wave functions in ab initio settings. In contrast to the
traditional paradigm of machine-learning force fields from training energies,
this considers the interpolation of accurate wave functions across chemical
space, from which all electronic properties as well as atomic forces can be
efficiently computed via a variational ansatz, avoiding the exponential
complexity of the many-body states themselves. Using the scheme to converge
the potential surface for molecular dynamics, we find examples of
qualitatively different behavior to state-of-the-art techniques, demonstrating
the importance of systematically converging the electronic structure across
the timescales.

The acceleration scheme therefore holds huge potential to extend the scope of
modern highly accurate electronic structure to molecular dynamics
applications. However, the potential of reliable wave function interpolation
also goes beyond this, towards a consideration of non-adiabatic and beyond-
Born–Oppenheimer effects, efficient geometry optimization for ground,
transition states or conical intersections, as well as a general procedure for
vibrations and phonons, raising the possibility of the routine extraction of
thermodynamic variables from accurate quantum chemistry. The move from single-
point electronic internal energies to (thermo)dynamical quantities within
correlated electronic structure theory is a long saught-after ambition
[https://doi.org/10.1002/anie.201709943]. The use of wave function
interpolation with developments in solvers for the training data to extend
system sizes could bring this closer to reality.

## Methods

### Interpolating across _ab initio_ potential energy surfaces

At the core of the methodology lies the prediction of the ground state
electronic energy for given molecular arrangement of $N_{\mathrm{elec}}$
electrons based on few exemplary solutions of the electronic structure problem
at different molecular geometries. We define the ab initio electronic
Hamiltonian for a $3\times N_{\mathrm{nuc}}$ atomic configuration,
$\mathbf{R}$, in a discrete basis of electronic orbitals
$\\{\chi(\mathbf{r};\mathbf{R})\\}$ as [https://doi.org/10.1002/wcms.1123]

$\displaystyle\hat{H}(\mathbf{R})$
$\displaystyle=\sum_{ij}h^{(1)}_{ij}(\mathbf{R})\,\hat{c}^{\dagger}_{i}\hat{c}_{j}+\frac{1}{2}\sum_{ijkl}h^{(2)}_{ijkl}(\mathbf{R})\,\hat{c}^{\dagger}_{i}\hat{c}^{\dagger}_{j}\hat{c}_{l}\hat{c}_{k}$
$\displaystyle\qquad+E_{\mathrm{nuc}}({\mathbf{R}})$ (5)
$\displaystyle=\sum_{ijkl}K_{ijkl}(\mathbf{R})\,\hat{c}_{i}^{\dagger}\hat{c}_{j}^{\dagger}\hat{c}_{l}\hat{c}_{k}+E_{\mathrm{nuc}}({\mathbf{R}}),$
(6)

with Fermionic creation and annihilation operators, $\hat{c}^{\dagger}$ and
$\hat{c}$ acting on the orbitals. The one-electron terms,
$h^{(1)}_{ij}(\mathbf{R})$, are matrix elements of the electron-nuclear and
electronic kinetic operators, while the electron-electron repulsion integrals
are

$\displaystyle h^{(2)}_{ijkl}(\mathbf{R})$ $\displaystyle=\int\int
d\mathbf{r}_{1}d\mathbf{r}_{2}\,\chi_{i}^{*}(\mathbf{r}_{1};\mathbf{R})\,\chi_{j}^{*}(\mathbf{r}_{2};\mathbf{R})$
(7)
$\displaystyle\qquad\frac{1}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}\chi_{k}(\mathbf{r}_{1};\mathbf{R})\,\chi_{l}(\mathbf{r}_{2};\mathbf{R}),$
$\displaystyle=\langle ij|kl\rangle.$ (8)

A convenient reduced two-body Hamiltonian which subsumes the one-body into the
two-body term can be written as

$K_{ijkl}=\frac{1}{2}\langle
ij|kl\rangle+\frac{1}{2(N_{\mathrm{elec}}-1)}(\delta_{jl}h^{(1)}_{ik}+\delta_{ik}h^{(1)}_{jl}),$
(9)

with $E_{\mathrm{nuc}}(\mathbf{R})$ the classical nuclear-nuclear repulsion
term [PhysRevA.57.4219].

As described in the main text, the eigenvector continuation proceeds via the
definition of a symmetrically (Löwdin) orthonormalized atomic orbital basis
(SAO) [10.1063/1.1747632, https://doi.org/10.1002/qua.981]. This allows the
training wave functions to be transferred between the Hilbert spaces of
different geometries by fixing their many-body probability amplitudes in this
representation. These SAOs are defined with an orbital transformation of an
underlying non-orthogonal atom-centered ‘AO’ orbital basis set at each
geometry, $\\{\phi_{\alpha}(\mathbf{r};\mathbf{R})\\}$, as

$\chi_{i}(\mathbf{r};\mathbf{R})=\sum_{\alpha}[\mathbf{S}(\mathbf{R})]^{-1/2}_{\alpha
i}\,\phi_{\alpha}(\mathbf{r};\mathbf{R}).$ (10)

where $\mathbf{S}(\mathbf{R})$ is the atomic orbital overlap matrix

$S_{\alpha\beta}(\mathbf{R})=\int
d\mathbf{r}\,\phi^{\ast}_{\alpha}(\mathbf{r};\mathbf{R})\,\phi_{\beta}(\mathbf{r};\mathbf{R}).$
(11)

The continuation then proceeds according to the scheme outlined in the main
text, with the evaluation of the transition two-body density matrices
(t-2RDMs) and overlaps between the training points in their SAO
representations. Of particular importance for molecular dynamics is the
evaluation of analytic forces at each test geometry, which can be achieved via
application of the Hellmann-Feynman theorem, which is simple due to the fully
optimized variational nature of the interpolated states in the geometry-
independent basis [PhysRev.56.340, 10.1063/1.4927594]. This therefore only
required the derivatives of the electron integrals in the AO basis
[https://doi.org/10.1002/jcc.23981, https://doi.org/10.1002/wcms.1340,
10.1063/5.0006074], as well as the derivative of the transformation from the
atomic orbitals to the SAOs with respect to nuclear positions, which we
evaluate via first order perturbation theory [bamieh2022tutorial]. The
specifics of this evaluation is given in the Supplementary information.

#### Approximate training data: DMRG

Rather than relying on exact (FCI) training data, we also consider modern
numerically efficient approximations to the correlated electronic structure to
allow for access to larger systems, which are nevertheless systematically
improvable to the exact solution to the electronic Schrödinger equation for
training. These require not only the evaluation of accurate many-body wave
functions at the training geometries, but also the evaluation of the t-2RDMs
and overlaps between different training states.

Firstly, we consider the compression of the training wave functions in the
form of _Matrix Product States_ (MPS), optimized via the _density matrix
renormalization group_ (DMRG) algorithm [PhysRevLett.69.2863]. For this, we
used the spin-adapted implementation from the _block2_ library
[10.1063/1.3695642, 10.1063/5.0050902], working directly in the Fock space
defined by the SAOs. We optimize the training states with a schedule for
exponentially increasing bond dimension (a factor of $1.8$ per increase) and
decreasing noise in the MPS, a standard practice for stable ab initio DMRG
[10.1063/1.4905329], terminating when the difference of the energy upon fully
relaxing the state at an increased bond dimension is less than a specified
threshold, $\epsilon$. For all presented results, we set
$\epsilon=10^{-3}\,E_{\mathrm{h}}$ and start the MPS with an initial bond
dimension of $34$, giving training data confidently below the accepted
‘chemical accuracy’. For the reference data for the hydrogen chain evolution
(Fig. (4)), we set the tolerance to $\epsilon=10^{-5}\,E_{h}$, and starting
MPS bond dimension to $61$.

#### Approximate training data: CAS

In addition to the continuation from MPS training states optimized with DMRG,
we also present the use of _complete active space_ (CAS) solvers to access the
results of Fig. 3. These give an approximate ground state of the full
electronic structure problem according to

$|\Psi_{\mathrm{CAS}}\rangle=|0\rangle^{N_{\mathrm{vir}}}\otimes|\Psi_{\mathrm{AS}}\rangle\otimes|1\rangle^{N_{\mathrm{core}}},$
(12)

where $|\Psi_{\mathrm{AS}}\rangle$ represents the fully variationally
optimized state over all many-electron configurations within a chosen active
subspace of orbitals and electrons, while $|1\rangle^{N_{\mathrm{core}}}$
represents fully occupied orbitals spanning the remaining space of states that
are occupied in a mean-field (in this case Hartree–Fock) description of the
system, and $|0\rangle^{N_{\mathrm{vir}}}$ explicitly indicate that the
higher-energy virtual states are unoccupied. In this way, the electronic
fluctuations of a low-energy subspace are considered fully, with the choice of
active space in this work selected purely based on the mean-field orbital
energies about the chemical potential of the system.

While this state can be straightforwardly optimized within a ‘CAS-CI’ scheme
implemented in the _PySCF_ package [https://doi.org/10.1002/wcms.1340,
10.1063/5.0006074], we also require the evaluation of the overlap and the
t-2RDMs between training states in their SAO basis, while the state is defined
(and optimized) in a mean-field canonical basis. Therefore, it is necessary to
rotate these many-body states into their respective SAO bases before the
t-2RDMs and overlaps computed. We show this for the t-2RDM as

$\displaystyle\Gamma^{ijkl}_{ab}$
$\displaystyle=\langle\Psi^{(a)}_{\mathrm{CAS}}|\hat{U}^{\dagger}_{\mathbf{R}^{(a)}}\hat{c}^{\dagger}_{i}\hat{c}^{\dagger}_{j}\hat{c}_{k}\hat{c}_{l}\hat{U}_{\mathbf{R}^{(b)}}|\Psi^{(b)}_{\mathrm{CAS}}\rangle,$
(13)

where $|\Psi^{(a/b)}_{\mathrm{CAS}}\rangle$ denotes the CAS states at the
different training points and $\hat{U}_{\mathbf{R}^{(a/b)}}$ is the unitary
transformation from the state in the basis of its canonical orbitals to the
SAO basis for the corresponding training point. This is evaluated efficiently
as a double summation over the active space many-electron configurations
(including their core) of each training state

$\Gamma^{ijkl}_{ab}=\\\ \sum_{\mathbf{n},\mathbf{n}^{\prime}\in
AS}C^{(a)\ast}_{\mathbf{n}}C^{(b)}_{\mathbf{n}^{\prime}}\,\langle\mathbf{n}|\hat{U}^{\dagger}_{\mathbf{R}^{(a)}}\hat{c}^{\dagger}_{i}\hat{c}^{\dagger}_{j}\hat{c}_{l}\hat{c}_{k}\hat{U}_{\mathbf{R}^{(b)}}|\mathbf{n}^{\prime}\rangle,$
(14)

where $C^{(a/b)}_{\mathbf{n}/\mathbf{n}^{\prime}}$ are the CAS-CI probability
amplitudes of the active spaces. This single-particle unitary transformation
$\hat{U}^{\dagger}_{\mathbf{R}}$ can be formed as

$\hat{U}_{\mathbf{R},ix}=\sum_{\alpha,\beta}Z_{\alpha
i}\,S_{\alpha\beta}\,\tilde{Z}_{\beta x},$ (15)

where $Z_{\alpha i}$ is the transformation matrix from AO to SAO and
$\tilde{Z}_{\beta x}$ is the transformation matrix from AO to canonical
Hatree–Fock orbitals, while $S_{\alpha\beta}$ is the AO overlap matrix. All of
these quantities are dependent on the specific training geometry,
$\mathbf{R}$.

The inner products
$\langle\mathbf{n}|\hat{U}^{\dagger}_{\mathbf{R}^{(a)}}\hat{c}^{\dagger}_{i}\hat{c}^{\dagger}_{j}\hat{c}_{k}\hat{c}_{l}\hat{U}_{\mathbf{R}^{(b)}}|\mathbf{n}^{\prime}\rangle$
from Eq. (14) can be identified as a matrix element between two different non-
orthogonal Slater determinants [doi:10.1021/acs.jctc.5b01148]. The efficient
evaluation of such overlaps between different non-orthogonal Slater
determinants is discussed in Refs. [10.1063/5.0045442, 10.1063/5.0122094]. We
utilize the _libgnme_ package, together with its python interface _pygnme_ ,
to evaluate the overlaps and t-2RDMs between CAS states for the continuation
in the SAO basis. Due to the non-orthogonality of the different CAS spaces,
the double contraction of Eq. (14) results in a cost scaling quadratically in
the size of the active space, thus more expensive than the evaluation of
expectation values of a single point CAS state, however this cost could be
reduced in the future by rotating to an intermediate basis representing the
co-domain of the occupied CAS orbitals in a pair of CAS-CI training states.

We note here that the eigenvector continuation scheme based on approximate
solvers (especially CAS-CI) can result in inferred energies at training points
which are variationally lower than those of the corresponding training energy
at the same geometry. This is because it can mix contributions to the test
state from other training states, with the resulting linear combination of
approximate states lower in energy, resulting in an improved estimate for the
energy at a test geometry than its corresponding training energy. Fig. 6 of
the extended data section exemplifies this behaviour, indicating lower
predicted energies from the continuation of CAS states as compared to the
reference CAS energies along the trajectory of a vibrating water molecule.

#### Gaussian Approximation Potentials

We include comparison results obtained from the prediction of potential
energies via _Gaussian Approximation Potentials_ (GAP)
[PhysRevLett.104.136403] – a well-established framework for the prediction of
potential energy surfaces and force fields. The model is extracted by fitting
a data set of training geometries, $\\{\mathbf{R}^{(a)}\\}_{a=1}^{N}$, with
associated energies $\\{E^{(a)}\\}_{a=1}^{N}$ using a kernel model
[rasmussen2006gaussian] incorporating symmetries of atomic environments via
the _smooth overlap of atomic position_ (SOAP) descriptors
[doi:10.1021/acs.chemrev.1c00021]. We apply the GAP framework following
standard approaches from the literature [C6CP00415F, PhysRevB.95.214302,
10.1063/5.0160898], based on the implementation of the SOAP descriptors in the
_dscribe_ package [HIMANEN2020106949]. Additional details of this prediction
procedure can be found in the supplementary information.

### Eigenvector continuation for BOMD

The Born–Oppenheimer molecular dynamics trajectories were computed in vacuum
based on a microcanonical (NVE) ensemble using the Velocity-Verlet integration
implemented in _PySCF_ [https://doi.org/10.1002/wcms.1340, 10.1063/5.0006074,
10.1063/1.442716, marx2000ab]. This solves Newton’s equation for the position
of each nucleus,

$m_{i}\frac{\mathrm{d}^{2}\mathbf{R}_{i}}{\mathrm{d}t^{2}}=\mathbf{F}_{i}(\mathbf{R}),$
(16)

where $\mathbf{R}_{i}$ is the position of the nucleus indexed by $i$, with
mass $m_{i}$ and associated force component $\mathbf{F}_{i}$ according to the
analytic nuclear gradients derived for the eigenvector continuation in the
supplementary information. In all tests, we initialized the nuclei at rest,
and chose a fixed timestep of $\delta t=5\,a.u\approx 0.121\,fs$ for the
integration.

#### Active learning for data selection

For the molecular dynamics applications, we perform an active learning scheme
in which we identify and select suitable molecular configurations for training
the eigenvector continuation scheme on-the-fly. This scheme is based on
iteratively running the MD with a given training set, and selecting an
enlarged training dataset with a new molecular configuration from the sampled
trajectory. A correlated electronic structure calculation is performed at the
selected geometry which is then included in the training data for an improved
inferred potential energy surface and resulting MD trajectory in the next
step. Starting from just a single training state (the initial geometry) and
iteratively adding new configurations to the dataset in this way, the number
of costly electronic structure calculations can be minimized and the
trajectory can be systematically and rapidly converged, noting that adding
training geometries from the simulated trajectories guarantees an improved
prediction in each step.

To select the new training geometry, we develop a ‘distance’ heuristic for all
geometries along the trajectory, which can be used as a metric for the
addition of the data, and quantifies the suitability of the current training
data in describing the test state at each point. The point along the
trajectory with the largest measure is added to the training data set. Since
the (non-degenerate) ground states along the trajectory are uniquely defined
by the ab initio Hamiltonian at each geometry, we use the differences between
the Hamiltonian elements at the training points and all trajectory points in
defining this measure. Defining these elements in their respective SAO basis
of each geometry used for the inference also ensures that the invariances and
symmetries of the eigenvector continuation are also respected in this measure.
Specifically, we define this Hamiltonian distance between two geometries,
$D(\mathbf{R},\mathbf{R}^{\prime})$, as

$D(\mathbf{R},\mathbf{R}^{\prime})=\sum_{ij}|h^{(1)}_{ij}(\mathbf{R})-h^{(1)}_{ij}(\mathbf{R}^{\prime})|^{2}+\\\
\frac{1}{2}\sum_{ijkl}|h^{(2)}_{ijkl}(\mathbf{R})-h^{(2)}_{ijkl}(\mathbf{R}^{\prime})|^{2}.$
(17)

In addition to respecting the symmetries of the model, this ensures that two
geometries with similar Hamiltonians (and thus wave functions) are considered
similar, even though an evaluation of the Euclidean distance between these two
geometries might be large (e.g., for geometries from near a dissociated
limit). To add a new configuration, we evaluate
$D(\mathbf{R}(t),\mathbf{R}^{(a)})$ for all geometries $\mathbf{R}(t)$ from
the trajectory and each training geometry $\mathbf{R}^{(a)}$ already contained
in the training set. We then pick that configuration $\mathbf{R}(t_{add})$
from the trajectory for which the distance to the closest training
configuration is maximal, i.e., where

$\quad
t_{\mathrm{add}}=\operatorname*{arg\,max}_{t}\left(\operatorname*{min}_{a}(D(\mathbf{R}(t),\mathbf{R}^{(a)}))\right)$
(18)

To gauge the systematic convergence of the MD trajectory, we track the
variational lowering (and hence improvement) of the potential energy surface
as the dataset is enlarged. This is done by comparing the smaller and larger
PES along the same trajectory corresponding to the larger data set. Due to the
variationality of the method, it is guaranteed that the potential energy
inference with the larger dataset will be lower or equal to the predictions
with the smaller dataset, and we use the difference between the predicted
energies as a convergence measure. In our applications, we terminate the
simulation when the predicted energy with the enlarged data set stays within a
tolerance of $\epsilon=10^{-3}\,E_{\mathrm{h}}$ along the full MD trajectory
for two iterations in a row. Examples of this convergence are shown in the
extended data.

## Code and data availability

The code and inputs to fully reproduce the numerical experiments of this work
can be found at https://github.com/BoothGroup/evcont. Furthermore, we provide
the videos of the molecular dynamics trajectories presented in this work at
https://doi.org/10.5281/zenodo.10658975.

## Acknowledgments

We thank Huanchen Zhai for support with the _block2_ code, Hugh Burton for
support with the _libgnme_ code, Carlos Mejuto-Zaera for helpful insights
about an independently-developed related scheme, Oliver Backhouse for
technical help and Kemal Atalar and Lachlan Lindoy for additional feedback on
the manuscript. We gratefully acknowledge support from the Air Force Office of
Scientific Research under award number FA8655-22-1-7011 and the UK Materials
and Molecular Modelling Hub for computational resources, which is partially
funded by EPSRC (EP/T022213/1, EP/W032260/1 and EP/P020194/1). YR also
acknowledges the support of the Engineering and Physical Sciences Research
Council [grant EP/Y005090/1].

## Author contributions

YR and GHB jointly developed the methodology and wrote the manuscript. YR
implemented the approach and performed the numerical experiments.

## References

  *   * [1] Motta, M. _et al._ Towards the solution of the many-electron problem in real materials: Equation of state of the hydrogen chain with state-of-the-art many-body methods.  _Phys. Rev. X_ 7, 031059 (2017).
  * [2] Williams, K. T. _et al._ Direct comparison of many-body methods for realistic electronic hamiltonians.  _Phys. Rev. X_ 10, 011041 (2020).
  * [3] Eriksen, J. J. _et al._ The ground state electronic energy of benzene.  _J. Phys. Chem. Lett._ 11, 8922–8929 (2020).  PMID: 33022176.
  * [4] Olivares-Amaya, R. _et al._ The ab-initio density matrix renormalization group in practice.  _J. Chem. Phys._ 142, 034102 (2015).
  * [5] Zhai, H. & Chan, G. K.-L.  Low communication high performance ab initio density matrix renormalization group algorithms.  _J. Chem. Phys._ 154, 224116 (2021).
  * [6] Larsson, H. R., Zhai, H., Umrigar, C. J., Kin-Lic, G. & Chan.  The chromium dimer: Closing a chapter of quantum chemistry.  _J. Am. Chem. Soc._ 144, 15932 – 15937 (2022).  URL http://dx.doi.org/10.1021/jacs.2c06357.
  * [7] Cui, Z.-H., Zhai, H., Zhang, X. & Chan, G. K.-L.  Systematic electronic structure in the cuprate parent state from quantum many-body simulations.  _Science_ 377, 1192–1198 (2022).
  * [8] Li, Z., Guo, S., Sun, Q. & Chan, G. K.-L.  Electronic landscape of the p-cluster of nitrogenase as revealed through many-electron quantum wavefunction simulations.  _Nat. Chem._ 11, 1026–1033 (2019).
  * [9] Booth, G. H., Grüneis, A., Kresse, G. & Alavi, A.  Towards an exact description of electronic wavefunctions in real solids.  _Nature_ 493, 365–370 (2012).  URL https://doi.org/10.1038/nature11770.
  * [10] Eskridge, B., Krakauer, H. & Zhang, S.  Non-perturbative many-body treatment of molecular magnets.  _J. Chem. Phys._ 158, 234110 (2023).
  * [11] Ren, W., Fu, W., Wu, X. & Chen, J.  Towards the ground state of molecules via diffusion monte carlo on neural networks.  _Nat. Commun._ 14, 1860 (2023).
  * [12] Cuzzocrea, A., Moroni, S., Scemama, A. & Filippi, C.  Reference excitation energies of increasingly large molecules: A qmc study of cyanine dyes.  _J. Chem. Theory Comput._ 18, 1089–1095 (2022).  PMID: 35080893.
  * [13] Koch, H. & Dalgaard, E.  Linear superposition of optimized non-orthogonal slater determinants for singlet states.  _Chem. Phys. Lett._ 212, 193–200 (1993).  URL https://www.sciencedirect.com/science/article/pii/000926149387129Q.
  * [14] Li, Z., Li, J., Dattani, N. S., Umrigar, C. J. & Chan, G. K.-L.  The electronic complexity of the ground-state of the FeMo cofactor of nitrogenase as relevant to quantum simulations.  _J. Chem. Phys._ 150, 024302 (2019).
  * [15] Huron, B., Malrieu, J. P. & Rancurel, P.  Iterative perturbation calculations of ground and excited state energies from multiconfigurational zeroth‐order wavefunctions.  _J. Chem. Phys._ 58, 5745–5759 (2003).
  * [16] Choo, K., Mezzacapo, A. & Carleo, G.  Fermionic neural-network states for ab-initio electronic structure.  _Nat. Commun._ 11 (2020).
  * [17] Barrett, T. D., Malyshev, A. & Lvovsky, A. I.  Autoregressive neural-network wavefunctions for ab initio quantum chemistry (2021).  2109.12606.
  * [18] Zhao, T., Stokes, J. & Veerapaneni, S.  Scalable neural quantum states architecture for quantum chemistry (2022).  2208.05637.
  * [19] Wu, Y. _et al._ A real neural network state for quantum chemistry (2023).  2301.03755.
  * [20] Yang, P.-J., Sugiyama, M., Tsuda, K. & Yanai, T.  Artificial Neural Networks Applied as Molecular Wave Function Solvers.  _J. Chem. Theory Comput._ 16, 3513–3529 (2020).
  * [21] Rath, Y. & Booth, G. H.  Framework for efficient ab initio electronic structure with gaussian process states.  _Phys. Rev. B_ 107, 205119 (2023).
  * [22] Bortone, M., Rath, Y. & Booth, G. H.  Impact of conditional modelling for a universal autoregressive quantum state (2024).
  * [23] Pfau, D., Spencer, J. S., Matthews, A. G. D. G. & Foulkes, W. M. C.  Ab initio solution of the many-electron schrödinger equation with deep neural networks.  _Phys. Rev. Res._ 2, 033429 (2020).
  * [24] Hermann, J., Schätzle, Z. & Noé, F.  Deep-neural-network solution of the electronic schrödinger equation.  _Nat. Chem._ 1–7 (2019).  URL http://dx.doi.org/10.1038/s41557-020-0544-y.
  * [25] Gerard, L., Scherbela, M., Marquetand, P. & Grohs, P.  Oh, A. H., Agarwal, A., Belgrave, D. & Cho, K. (eds) _Gold-standard solutions to the schrödinger equation using deep learning: How much physics do we need?_ (eds Oh, A. H., Agarwal, A., Belgrave, D. & Cho, K.) _Advances in Neural Information Processing Systems_ (2022).  URL https://openreview.net/forum?id=nX-gReQ0OT.  2205.09438.
  * [26] Park, J. W., Al-Saadon, R., MacLeod, M. K., Shiozaki, T. & Vlaisavljevich, B.  Multireference electron correlation methods: Journeys along potential energy surfaces.  _Chem. Rev._ 120, 5878–5909 (2020).  PMID: 32239929.
  * [27] Del Ben, M., Schönherr, M., Hutter, J. & VandeVondele, J.  Bulk liquid water at ambient temperature and pressure from mp2 theory.  _J. Phys. Chem. Lett._ 4, 3753–3759 (2013).
  * [28] Luo, Y. & Sorella, S.  Ab initio molecular dynamics with quantum monte carlo.  _Front. Mater._ 2 (2015).
  * [29] Kapil, V., VandeVondele, J. & Ceriotti, M.  Accurate molecular dynamics and nuclear quantum effects at low cost by multiple steps in real and imaginary time: Using density functional theory to accelerate wavefunction methods.  _J. Chem. Phys._ 144, 054111 (2016).
  * [30] Marx, D. & Hutter, J.  Ab initio molecular dynamics: Theory and implementation.  _Modern methods and algorithms of quantum chemistry_ 1, 141 (2000).
  * [31] Cohen, A. J., Mori-Sánchez, P. & Yang, W.  Challenges for density functional theory.  _Chem. Rev._ 112 1, 289–320 (2012).  URL https://doi.org/10.1021/cr200107z.
  * [32] Kabylda, A., Vassilev-Galindo, V., Chmiela, S., Poltavsky, I. & Tkatchenko, A.  Efficient interatomic descriptors for accurate machine learning force fields of extended molecules.  _Nat. Commun._ 14 (2023).  URL https://doi.org/10.1038/s41467-023-39214-w.
  * [33] Fu, X. _et al._ Forces are not enough: Benchmark and critical evaluation for machine learning force fields with molecular simulations.  _ArXiv_ abs/2210.07237 (2022).  URL https://doi.org/10.48550/arXiv.2210.07237.
  * [34] Unke, O. T. _et al._ Machine learning force fields.  _Chem. Rev._ 121, 10142 – 10186 (2020).  URL https://doi.org/10.1021/acs.chemrev.0c01111.
  * [35] Gkeka, P. _et al._ Machine learning force fields and coarse-grained variables in molecular dynamics: application to materials and biological systems.  _J. Chem. Theory Comput._ (2020).  URL https://doi.org/10.1021/acs.jctc.0c00355.
  * [36] Musil, F. _et al._ Physics-inspired structural representations for molecules and materials.  _Chem. Rev._ 121, 9759–9815 (2021).  PMID: 34310133.
  * [37] Kulik, H. J. _et al._ Roadmap on machine learning in electronic structure.  _Electron. Struct._ 4, 023004 (2022).
  * [38] Chong, S. _et al._ Robustness of local predictions in atomistic machine learning models.  _J. Chem. Theory Comput._ 19, 8020–8031 (2023).  PMID: 37948446.
  * [39] Kovács, D. P. _et al._ Linear atomic cluster expansion force fields for organic molecules: Beyond rmse.  _J. Chem. Theory Comput._ 17, 7696–7711 (2021).  PMID: 34735161.
  * [40] Schütt, K. T., Gastegger, M., Tkatchenko, A., Müller, K.-R. & Maurer, R. J.  Unifying machine learning and quantum chemistry with a deep neural network for molecular wavefunctions.  _Nat. Commun._ 10 (2019).  URL https://doi.org/10.1038/s41467-019-12875-2.
  * [41] Fedik, N. _et al._ Extending machine learning beyond interatomic potentials for predicting molecular properties.  _Nat. Rev. Chem._ 6, 653 – 672 (2022).  URL https://doi.org/10.1038/s41570-022-00416-3.
  * [42] Westermayr, J., Gastegger, M. & Marquetand, P.  Combining schnet and sharc: The schnarc machine learning approach for excited-state dynamics.  _J. Phys. Chem. Lett._ 11, 3828–3834 (2020).  PMID: 32311258.
  * [43] Klawohn, S. _et al._ Gaussian approximation potentials: Theory, software implementation and application examples.  _J. Chem. Phys._ 159, 174108 (2023).
  * [44] Agmon, N.  The grotthuss mechanism.  _Chem. Phys. Lett._ 244, 456–462 (1995).  URL https://doi.org/10.1016/0009-2614(95)00905-J.
  * [45] Cukierman, S.  Et tu, grotthuss! and other unfinished stories.  _Biochim. Biophys. Acta_ 1757 8, 876–85 (2006).  URL https://doi.org/10.1016/j.bbabio.2005.12.001.
  * [46] Knowles, P. J. & Handy, N. C.  A determinant based full configuration interaction program.  _Comput. Phys. Commun._ 54, 75–83 (1989).  URL https://www.sciencedirect.com/science/article/pii/0010465589900337.
  * [47] Mejuto-Zaera, C. & Kemper, A. F.  Quantum eigenvector continuation for chemistry applications.  _Electron. Struct._ 5, 045007 (2023).
  * [48] Löwdin, P.  On the Non‐Orthogonality Problem Connected with the Use of Atomic Wave Functions in the Theory of Molecules and Crystals.  _J. Chem. Phys._ 18, 365–375 (2004).
  * [49] Mayer, I.  On löwdin’s method of symmetric orthogonalization*.  _Int. J. Quantum Chem._ 90, 63–65 (2002).
  * [50] Peng, B. & Kowalski, K.  Highly efficient and scalable compound decomposition of two-electron integral tensor and its application in coupled cluster calculations.  _J. Chem. Theory Comput._ 13, 4179–4192 (2017).  PMID: 28834428.
  * [51] Wei, H. & Neuscamman, E.  Reduced scaling Hilbert space variational Monte Carlo.  _J. Chem. Phys._ 149, 184106 (2018).
  * [52] Sabzevari, I. & Sharma, S.  Improved speed and scaling in orbital space variational monte carlo.  _J. Chem. Theory Comput._ 14, 6276–6286 (2018).  PMID: 30418769.
  * [53] Peng, L., Zhang, X. & Chan, G. K.-L.  Fermionic reduced density low-rank matrix completion, noise filtering, and measurement reduction in quantum simulations.  _J. Chem. Theory Comput._ 19, 9151–9160 (2023).  PMID: 38095484.
  * [54] Frame, D. _et al._ Eigenvector continuation with subspace learning.  _Phys. Rev. Lett._ 121, 032501 (2018).
  * [55] Demol, P. _et al._ Improved many-body expansions from eigenvector continuation.  _Phys. Rev. C_ 101, 041302 (2020).
  * [56] Duguet, T., Ekström, A., Furnstahl, R. J., König, S. & Lee, D.  Eigenvector continuation and projection-based emulators (2023).  2310.19419.
  * [57] König, S., Ekström, A., Hebeler, K., Lee, D. & Schwenk, A.  Eigenvector continuation as an efficient and accurate emulator for uncertainty quantification.  _Phys. Lett. B_ 810, 135814 (2020).  URL https://www.sciencedirect.com/science/article/pii/S0370269320306171.
  * [58] Yapa, N., Fossez, K. & König, S.  Eigenvector continuation for emulating and extrapolating two-body resonances.  _Phys. Rev. C_ 107, 064316 (2023).
  * [59] Drischler, C., Quinonez, M., Giuliani, P., Lovell, A. & Nunes, F.  Toward emulating nuclear reactions using eigenvector continuation.  _Phys. Lett. B_ 823, 136777 (2021).  URL https://www.sciencedirect.com/science/article/pii/S0370269321007176.
  * [60] Schrader, S. E. & Kvaal, S.  Accelerated coupled cluster calculations with Procrustes orbital interpolation.  _J. Chem. Phys._ 158, 114116 (2023).
  * [61] Bartók, A. P., Payne, M. C., Kondor, R. & Csányi, G.  Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons.  _Phys. Rev. Lett._ 104, 136403 (2010).
  * [62] Glielmo, A., Sollich, P. & De Vita, A.  Accurate interatomic force fields via machine learning with covariant kernels.  _Phys. Rev. B_ 95, 214302 (2017).
  * [63] Levine, B. G., Durden, A. S., Esch, M. P., Liang, F. & Shu, Y.  CAS without SCF—Why to use CASCI and where to get the orbitals.  _J. Chem. Phys._ 154, 090902 (2021).
  * [64] Motta, M. _et al._ Ground-state properties of the hydrogen chain: Dimerization, insulator-to-metal transition, and magnetic phases.  _Phys. Rev. X_ 10, 031058 (2020).
  * [65] Huang, X., Braams, B. J. & Bowman, J. M.  Ab initio potential energy and dipole moment surfaces for h5o2 +.  _J. Chem. Phys._ 122 4, 44308 (2005).  URL https://doi.org/10.1063/1.1834500.
  * [66] McCoy, A. B., Huang, X., Carter, S., Landeweer, M. Y. & Bowman, J. M.  Full-dimensional vibrational calculations for h5o2+ using an ab initio potential energy surface.  _J. Chem. Phys._ 122 6, 061101 (2005).  URL https://doi.org/10.1063/1.1857472.
  * [67] Natarajan, S. K., Morawietz, T. & Behler, J.  Representing the potential-energy surface of protonated water clusters by high-dimensional neural network potentials.  _Phys. Chem. Chem. Phys._ 17 13, 8356–71 (2015).  URL https://doi.org/10.1039/c4cp04751f.
  * [68] Heindel, J., Yu, Q., Bowman, J. M. & Xantheas, S. S.  Benchmark electronic structure calculations for h3o+(h2o) n, n = 0-5, clusters and tests of an existing 1,2,3-body potential energy surface with a new 4-body correction.  _J. Chem. Theory Comput._ 14 9, 4553–4566 (2018).  URL https://doi.org/10.1021/acs.jctc.8b00598.
  * [69] Grimme, S. & Schreiner, P. R.  Computational chemistry: The fate of current methods and future challenges.  _Angew. Chem. Int. Ed._ 57, 4170–4176 (2018).
  * [70] Jensen, F.  Atomic orbital basis sets.  _WIREs Comput. Mol. Sci._ 3, 273–295 (2013).
  * [71] Mazziotti, D. A.  Contracted schrödinger equation: Determining quantum energies and two-particle density matrices without wave functions.  _Phys. Rev. A_ 57, 4219–4234 (1998).
  * [72] Feynman, R. P.  Forces in molecules.  _Phys. Rev._ 56, 340–343 (1939).
  * [73] Thomas, R. E. _et al._ Analytic nuclear forces and molecular properties from full configuration interaction quantum Monte Carlo.  _J. Chem. Phys._ 143, 054108 (2015).
  * [74] Sun, Q.  Libcint: An efficient general integral library for gaussian basis functions.  _J. Comput. Chem._ 36, 1664–1671 (2015).
  * [75] Sun, Q. _et al._ Pyscf: the python-based simulations of chemistry framework.  _WIREs Comput. Mol. Sci._ 8, e1340 (2018).
  * [76] Sun, Q. _et al._ Recent developments in the PySCF program package.  _J. Chem. Phys._ 153, 024109 (2020).
  * [77] Bamieh, B.  A tutorial on matrix perturbation theory (using compact matrix notation) (2022).  2002.05001.
  * [78] White, S. R.  Density matrix formulation for quantum renormalization groups.  _Phys. Rev. Lett._ 69, 2863–2866 (1992).
  * [79] Sharma, S. & Chan, G. K.-L.  Spin-adapted density matrix renormalization group algorithms for quantum chemistry.  _J. Chem. Phys._ 136, 124121 (2012).
  * [80] Plasser, F. _et al._ Efficient and flexible computation of many-electron wave function overlaps.  _J. Chem. Theory Comput._ 12, 1207–1219 (2016).  PMID: 26854874.
  * [81] Burton, H. G. A.  Generalized nonorthogonal matrix elements: Unifying Wick’s theorem and the Slater–Condon rules.  _J. Chem. Phys._ 154, 144109 (2021).
  * [82] Burton, H. G. A.  Generalized nonorthogonal matrix elements. II: Extension to arbitrary excitations.  _J. Chem. Phys._ 157, 204109 (2022).
  * [83] Rasmussen, C. E., Williams, C. K. _et al._ _Gaussian processes for machine learning_ Vol. 1 (Springer, 2006).
  * [84] De, S., Bartók, A. P., Csányi, G. & Ceriotti, M.  Comparing molecules and solids across structural and alchemical space.  _Phys. Chem. Chem. Phys._ 18, 13754–13769 (2016).
  * [85] Himanen, L. _et al._ Dscribe: Library of descriptors for machine learning in materials science.  _Comput. Phys. Commun._ 247, 106949 (2020).  URL https://www.sciencedirect.com/science/article/pii/S0010465519303042.
  * [86] Swope, W. C., Andersen, H. C., Berens, P. H. & Wilson, K. R.  A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters.  _J. Chem. Phys._ 76, 637–649 (1982).

## Extended data

#### Accuracy of the potential energy surface along MD trajectories

In Fig. 6, we report the difference between the inferred and reference
energies at every point along the simulated molecular dynamics trajectory of
the water vibration model in the main text. The top panel shows the accuracy
of the potential energy surface with the 6-31G basis set, where we performed
the DMRG-based continuation scheme, compared to the exact potential energy
surface. It can be seen that increasing the number of training states rapidly
converges the full trajectory with $N=6$ training states achieving an accuracy
well below $10^{-4}\,E_{\mathrm{h}}$ across all points.

The variationality of the method guarantees that the predicted energies are
always an upper bound to the exact ground state energy at any geometry.
Nonetheless, when the continuation is based on approximate training wave
functions, the computed linear combination of training states may result in an
energy lower than the reference training energy at that geometry. This is
exemplified in the bottom panel of Fig. 6, detailing the energetic difference
between the prediction and CAS-CI energies used for the training data along
the trajectory in cc-pVDZ and cc-pVTZ basis sets. While the same active space
sizes were used for the training states as for the computation of the
reference energies, the inferred energies generally lie below the CAS-CI
reference energies. Although this improvement is small, mostly less than
$1\,mE_{\mathrm{h}}$, it is obtained for the majority geometries from the
converged trajectory, ultimately improving the accuracy beyond what is
obtained with the reference method.

Figure 6: Difference between the predicted energy and a reference method over
all geometries from the MD trajectory of a water molecule. Top panel: Error
compared to FCI of DMRG-based eigenvector continuation with increasing data
set in the 6-31G basis. Bottom panel: Difference between CASCI(4,8)-trained
eigenvector continuation and CASCI(4,8) energies at each geometry along the
trajectory for $N=9$ training points in a cc-pVDZ basis, and $N=14$ training
points in a cc-pVTZ basis. We note that variationality with respect to this
approximate training data is not expected.

Analogously to Fig. 6, Fig. 7 depicts the difference between the energy
predicted by the eigenvector continuation and the reference DMRG energy for
the converged MD trajectory of the $30$-atom hydrogen system of the main text.
It can be seen that, even though there is no prior knowledge about the
targetted trajectory, the on-the-fly dataset generation for the molecular
dynamics simulation achieves a potential energy surface resolved to roughly a
few millihartree along the whole trajectory.

Figure 7: Energy difference between reference DMRG energies and the DMRG-
trained eigenvector continuation for the trajectory of the $30$-atom hydrogen
chain (main results presented in Fig. 4), with $N=2$, $N=20$, and $N=47$
training data points.

#### Zundel cation energy surface convergence

Figure 8 provides additional insights into the convergence of the DMRG-trained
potential energy surface (PES) underpinning the Zundel cation dynamics. In
this case, the system and timescales are too large to obtain an exact
reference comparison along the length of the trajectory, and so we consider
other metrics to assess its accuracy. In the top panel, we show the final
$N=84$ inferred electronic energy along the trajectory, compared to
approximate Hartree–Fock (HF) and coupled cluster with single and double
excitations (CCSD). The obtained values improve significantly upon single
point energy estimates from Hartree-Fock, and are generally comparible with
CCSD over the whole simulation, noting that CCSD is not a variational method
and therefore can be above or below the exact energy. Nevertheless, we expect
CCSD to be generally an accurate method to describe this system over most of
the geometries visited. It is significant however that for some geometries
along the trajectory, CCSD could not converge with the default settings, while
the eigenvector continuation, which lacks any self-consistent or iterative
procedure (even mean-field) for each trajectory can be run in an entirely
robust and black-box fashion.

In the center panel of the figure we highlight the convergence of the PES as a
function of the number of training geometries used, showing the average
electronic energy over the final trajectory for different numbers of training
points (shifting the energies such that the $N=84$ point has a zero average
PES). While the variational character of the continuation guarantees a
monotonic decrease of this average potential energy as more training
configurations are added, various flat regions can be observed, indicating
training configurations not contributing to a substantial improvement to the
potential energy. This highlights the potential for further improvements to
the data selection scheme in order to identify the most appropriate training
configurations and further accelerate the convergence of the PES with training
points.

Finally, the bottom panel of Fig. 8 shows the energy at the $84$ training
configurations used to train the model in the final converged MD simulation.
At a number of these training geometries, the HF and/or CCSD were not able to
fully converge, and these points have been omitted. While optimizing
convergence parameters would likely have enabled convergence, this is
potentially a flag for stronger correlation effects at these points, where the
DMRG-trained continuation would be more valuable and HF and CCSD unsuitable.
Furthermore, there are clear indications of at least 13 training points where
the DMRG training data was not able to correctly converge, likely due to
getting stuck in local minima in the DMRG sweep algorithm that optimizes the
MPS.

This also highlights the lack of ‘black-box’ character for these highly
accurate wave function based electronic structure methods, which means that
much care is needed if they were to be applied to MD on their own. Likely the
convergence of the continuation with respect to training points would have
been improved had these been fully converged, however it is noteworthy that
this does not substantially impact upon the quality of the inference of the
eigenvector continuation at these training points, which is found to be
substantially more accurate (variationally lower) than their corresponding
training point at these points even where the optimization of the training
wave function failed. Furthermore, the surface for the MD trajectory
necessarily remains smooth and the gradients well-behaved. We note that the
nature of the active data selection means that the molecular geometries of
these training points are unlikely to actually feature in the final MD
trajectory, since each time the data set is enlarged, the PES changes and the
MD follows a different path. In the final MD simulation, the HF does converge
for all points along the trajectory, and CCSD with default optimization
parameters fails $\sim 4\%$ of the time.

We note however that the rationale for selection of additional data points
builds upon the assumption that data points add the correct training state to
the dataset and therefore significantly improve the potential energy surface
at that geometry. Adding data points for which an appropriately converged
solution to the ground state was not obtained (in this case with DMRG) can
thus negatively affect the performance of the data selection for subsequent
runs and the rate of convergence of the method to exactness. Improving the
convergence of the molecular dynamics simulation even in limits where the data
set comprises inaccurate solutions to the electronic structure problem remains
a subject of future research.

Figure 8: Characteristics of potential energy surfaces and single point
training energy estimates for the simulation of the water-hydronium reaction
from the eigenvector continuation (blue), Hartree-Fock (green), as well as
CCSD (red). The top panel (a) shows the potential energy for these methods as
a function of time over the $N=84$ converged MD trajectory. Panel (b) shows
the average energy over the final trajectory which systematically converges as
more training data points are included in the eigenvector continuation. Panel
(c) shows the single point energy estimates for each of the $N=84$ training
geometries. Blue crosses denote the energy of the DMRG-trained eigenvector
continuation at the training geometries, the gray line gives the single point
energy estimates from the individual DMRG calculations used for the training
points (noting the spikes corresponding to poor DMRG convergence), while red
dots and green crosses show the CCSD and HF values at the training geometries
which were able to be converged along the trajectory with default parameters.

## Supplementary information

#### Nuclear forces for eigenvector continuation

Based on the potential energy estimate $E(\mathbf{R})$ for a nuclear
configuration $\mathbf{R}$, the force can be defined as

$F(\mathbf{R})=-\frac{\partial
E(\mathbf{R})}{\partial\mathbf{R}}=-\left(\frac{\partial
E_{\mathrm{nuc}}(\mathbf{R})}{\partial\mathbf{R}}+\frac{\partial
E_{\mathrm{el}}(\mathbf{R})}{\partial\mathbf{R}}\right),$ (19)

which as shown can be split up into the electronic and nuclear contributions
to this force. Focusing on the electronic contribution from the eigenvector
continuation (and dropping the subscript for clarity), we apply the Hellmann-
Feynman theorem [PhysRev.56.340], within the restricted subspace spanned by
the training states, giving

$\frac{\partial
E(\mathbf{R})}{\partial\mathbf{R}}=\frac{\partial}{\partial\mathbf{R}}\left(\boldsymbol{\psi}^{\dagger}\,\boldsymbol{\mathcal{H}}\,\boldsymbol{\psi}\right)=\boldsymbol{\psi}^{\dagger}\,\frac{\partial\boldsymbol{\mathcal{H}}}{\partial\mathbf{R}}\,\boldsymbol{\psi},$
(20)

where $\boldsymbol{\psi}$ is the inferred ground-state eigenfunction from the
diagonalization of the Hamiltonian in the many-body basis of the fixed
training states (Eq. (2) of the main text). Substituting in the definition of
the Hamiltonian from Eq. (5), similar to the derivation in Ref.
[10.1063/1.4927594], we obtain

$\frac{\partial E(\mathbf{R})}{\partial\mathbf{R}}=\sum_{ijkl}\frac{\partial
K_{ijkl}(\mathbf{R})}{\partial\mathbf{R}}\,\Gamma^{ijkl}$ (21)

where $\Gamma^{ijkl}$ is the two-electron reduced density matrix of the
inferred state,

$\Gamma^{ijkl}=\sum_{a,b}^{N}\psi^{*}_{a}\,\Gamma^{ijkl}_{ab}\,\psi_{b}.$ (22)

The reduced Hamiltonian gradient, $\frac{\partial
K_{ijkl}(\mathbf{R})}{\partial\mathbf{R}}$, is obtained according to

$\frac{\partial
K_{ijkl}(\mathbf{R})}{\partial\mathbf{R}}=\frac{1}{2}\frac{\partial
h^{(2)}_{ijkl}(\mathbf{R})}{\partial\mathbf{R}}+\\\
\frac{1}{2(N_{\mathrm{elec}}-1)}\left(\delta_{jl}\frac{\partial
h^{(1)}_{ik}(\mathbf{R})}{\partial\mathbf{R}}+\delta_{ik}\frac{\partial
h^{(1)}_{jl}(\mathbf{R})}{\partial\mathbf{R}}\right)$ (23)

The evaluation of the the one- and two-electron integral gradients in the SAO
basis requires the derivative of the transformation $\mathbf{Z}(\mathbf{R})$
from the atomic orbitals to the SAO basis, as well as the gradients of the
integrals between atomic orbitals. For the one-electron integral derivative,
we obtain

$\frac{\partial h^{(1)}_{ij}}{\partial\mathbf{R}}=\sum_{\alpha\beta}Z_{\alpha
i}\,Z_{\beta
j}\,\frac{\partial\tilde{h}^{(1)}_{\alpha\beta}}{\partial\mathbf{R}}\\\
+\sum_{\alpha\beta}\frac{\partial Z_{\alpha i}}{\partial\mathbf{R}}\,Z_{\beta
j}\,\tilde{h}^{(1)}_{\alpha\beta}+\sum_{\alpha\beta}\frac{\partial Z_{\beta
j}}{\partial\mathbf{R}}\,Z_{\alpha i}\,\tilde{h}^{(1)}_{\alpha\beta},$ (24)

where greek indices refer to atomic orbitals, latin indices to SAO, and
$\tilde{h}^{(1)}_{\alpha\beta}$ denotes the one-electron integral in the
atomic orbital basis. The derivatives of the integrals over the atomic
orbitals, are readly accessible from standard libraries, and we utilize the
implementation in the _Libcint_ library [https://doi.org/10.1002/jcc.23981] as
bundled in the _PySCF_ package [https://doi.org/10.1002/wcms.1340,
10.1063/5.0006074].

The gradient of the SAO transformation, $\frac{\partial Z_{\alpha
i}}{\partial\mathbf{R}}$, follows directly from the definition of the Löwdin
orthogonalization, and application of the chain rule. This gives

$\frac{\partial Z_{\alpha
i}}{\partial\mathbf{R}}=\sum_{\beta\gamma}\frac{\partial
S_{\beta\gamma}}{\partial\mathbf{R}}\frac{\partial(\mathbf{S}^{-1/2})_{i\alpha}}{\partial
S_{\beta\gamma}}.$ (25)

In addition to the derivative of the overlap matrix between atomic orbitals,
this requires the evaluation of the derivative of the inverse square root
matrix elements with respect to the overlap elements,
$\frac{\partial(\mathbf{S}^{-1/2})_{i\alpha}}{\partial S_{\beta\gamma}}$. To
evaluate this derivative, we consider the spectral decomposition of
$\mathbf{S}$ according to

$\mathbf{S}=\mathbf{T}\,\mathbf{s}\,\mathbf{T}^{\dagger},$ (26)

from which we obtain

$\frac{\partial\mathbf{S}^{-1/2}}{\partial
S_{\beta\gamma}}=\frac{\partial\mathbf{T}}{\partial
S_{\beta\gamma}}\,\mathbf{s}^{-1/2}\,\mathbf{T}^{\dagger}+\\\
\mathbf{T}\,\frac{\partial\mathbf{s}^{-1/2}}{\partial
S_{\beta\gamma}}\,\mathbf{T}^{\dagger}+\mathbf{T}\,\mathbf{s}^{-1/2}\,\frac{\partial\mathbf{T}^{\dagger}}{\partial
S_{\beta\gamma}}.$ (27)

The eigenvector and eigenvalue derivatives can be evaluated via first order
perturbation theory of each element of $S_{\beta\gamma}$ and evaluating the
derivative at vanishing perturbation strength [bamieh2022tutorial]. The chain
rule then provides

$\frac{\partial s^{-1/2}_{i}}{\partial
S_{\beta\gamma}}=-\frac{1}{2}\,s^{-3/2}_{i}\,\frac{\partial s_{i}}{\partial
S_{\beta\gamma}}=-\frac{1}{2}\,s^{-3/2}_{i}\,T^{\ast}_{\beta i}\,T_{\gamma
i},$ (28)

and we obtain for the derivative of the eigenvector coefficients

$\frac{\partial T_{\alpha i}}{\partial S_{\beta\gamma}}=-\sum_{j\neq
i}\frac{T^{\ast}_{\beta j}\,T_{\gamma i}\,T_{\alpha j}}{s_{j}-s_{i}}.$ (29)

The sum in the above expression runs over all eigenstates except $i$.
Additional care needs to be taken if eigenvalue degeneracies are present in
Eq. (26), in which case we appeal to standard degenerate perturbation theory
[bamieh2022tutorial], ensuring we rotate the eigenvectors in each degenerate
subspace such that the applied perturbation in this degenerate subspace is
diagonal. In this case, Eq. (29) becomes valid by only summing over all
eigenstates which have a distinct eigenvalue from $s_{i}$.

#### Details of GAP model predictions

Within the application of the GAP framework, the potential energy for a test
geometry is predicted based on a kernel model according to

$E(\mathbf{R})=\sum_{a=1}^{N}w_{a}\,k(\mathbf{R},\mathbf{R}^{(a)}),$ (30)

where the weights $w_{a}$ may be interpreted as learned parameters of the
model, and the kernel function $k$ defines a similarity measure between
nuclear geometries.

Without fully optimizing the various design choices, we follow standard
approaches to construct a suitable kernel function within the GAP framework
based on the idea of a smooth overlap of atomic positions (SOAP) [C6CP00415F].
This representation defines a set of features for each atom in the system
based on its local atomic environment represented in a basis of spherical
harmonics together with radial basis functions automatically incorporating
important symmetries into the representation. We utilize the implementation of
the SOAP features from the _dscribe_ library [HIMANEN2020106949], implementing
the SOAP features as described in Ref. [C6CP00415F]. Without additional
hyperparameter tuning, we constructed the SOAP features in our tests by going
up to 20th order in the spherical harmonics and 10th order in the radial basis
functions, and chose a radial cutoff of $10\,\text{\AA}$.

Based on these SOAP features, we construct an averaged kernel function
[C6CP00415F], given as the average of scalar products between the SOAP feature
vectors over all pairs of atoms in the system. Denoting the vector of SOAP
features for the $i$-indexed atom as $\mathbf{p}(\mathbf{R}_{i})$, this
kernel, appropriately normalized, is thus given as

$k(\mathbf{R},\mathbf{\mathbf{R}}^{\prime})=\frac{\sum_{ij}\mathbf{p}(\mathbf{R}_{i})\cdot\mathbf{p}(\mathbf{R}^{\prime}_{j})}{\sqrt{\mathcal{N}(\mathbf{R})\times\mathcal{N}(\mathbf{R}^{\prime})}},$
(31)

with a normalization constant

$\mathcal{N}(\mathbf{R})=\sum_{ij}\mathbf{p}(\mathbf{R}_{i})\cdot\mathbf{p}(\mathbf{R}_{j}).$
(32)

Following standard Gaussian process regression techniques
[rasmussen2006gaussian], the weights of the model are obtained by inversion of
the kernel matrix according to

$\mathbf{w}=(\mathbf{K}+\sigma^{2}\mathbb{1})^{-1}\mathbf{E}.$ (33)

Here, $\mathbf{w}$ denotes the vector of weights for all data points,
$\mathbf{K}$ is the matrix of pairwise kernel values between the training
configurations, and $\mathbf{E}$ is a vector comprising the training energies.
For the results shown in Fig. 2 of the main text, we set the additional noise
parameter $\sigma^{2}$, effectively regularizing the fit, to a small fixed
value of $\sigma^{2}=10^{-15}$, without additional optimization.

While it is possible to train GAP models also on force data to improve the
prediction [PhysRevB.95.214302], here we simply estimate the force by
differentiating through the potential energy predictor according to Eq. (30).
The force prediction is therefore obtained as

$F(\mathbf{R})=-\sum_{a=1}^{N}w_{a}\frac{\partial
k(\mathbf{R},\mathbf{R}^{(a)})}{\partial\mathbf{R}},$ (34)

with the kernel function as defined in Eq. (31). It should be noted that the
SOAP features used for the kernel are atom-centered, and therefore the
derivative with respect to these basis functions needs to be taken into
account when evaluating the gradient, as implemented in the _dscribe_ library
[HIMANEN2020106949].

## References

  * [1] Bassam Bamieh.  A tutorial on matrix perturbation theory (using compact matrix notation), 2022.
  * [2] Sandip De, Albert P. Bartók, Gábor Csányi, and Michele Ceriotti.  Comparing molecules and solids across structural and alchemical space.  Phys. Chem. Chem. Phys., 18:13754–13769, 2016.
  * [3] R. P. Feynman.  Forces in molecules.  Phys. Rev., 56:340–343, Aug 1939.
  * [4] Aldo Glielmo, Peter Sollich, and Alessandro De Vita.  Accurate interatomic force fields via machine learning with covariant kernels.  Phys. Rev. B, 95:214302, Jun 2017.
  * [5] Lauri Himanen, Marc O.J. Jäger, Eiaki V. Morooka, Filippo Federici Canova, Yashasvi S. Ranawat, David Z. Gao, Patrick Rinke, and Adam S. Foster.  Dscribe: Library of descriptors for machine learning in materials science.  Comput. Phys. Commun., 247:106949, 2020.
  * [6] Carl Edward Rasmussen, Christopher KI Williams, et al.  Gaussian processes for machine learning, volume 1.  Springer, 2006.
  * [7] Qiming Sun.  Libcint: An efficient general integral library for gaussian basis functions.  J. Comput. Chem., 36(22):1664–1671, 2015.
  * [8] Qiming Sun, Timothy C. Berkelbach, Nick S. Blunt, George H. Booth, Sheng Guo, Zhendong Li, Junzi Liu, James D. McClain, Elvira R. Sayfutyarova, Sandeep Sharma, Sebastian Wouters, and Garnet Kin-Lic Chan.  Pyscf: the python-based simulations of chemistry framework.  WIREs Comput. Mol. Sci., 8(1):e1340, 2018.
  * [9] Qiming Sun, Xing Zhang, Samragni Banerjee, Peng Bao, Marc Barbry, Nick S. Blunt, Nikolay A. Bogdanov, George H. Booth, Jia Chen, Zhi-Hao Cui, Janus J. Eriksen, Yang Gao, Sheng Guo, Jan Hermann, Matthew R. Hermes, Kevin Koh, Peter Koval, Susi Lehtola, Zhendong Li, Junzi Liu, Narbe Mardirossian, James D. McClain, Mario Motta, Bastien Mussard, Hung Q. Pham, Artem Pulkin, Wirawan Purwanto, Paul J. Robinson, Enrico Ronca, Elvira R. Sayfutyarova, Maximilian Scheurer, Henry F. Schurkus, James E. T. Smith, Chong Sun, Shi-Ning Sun, Shiv Upadhyay, Lucas K. Wagner, Xiao Wang, Alec White, James Daniel Whitfield, Mark J. Williamson, Sebastian Wouters, Jun Yang, Jason M. Yu, Tianyu Zhu, Timothy C. Berkelbach, Sandeep Sharma, Alexander Yu. Sokolov, and Garnet Kin-Lic Chan.  Recent developments in the PySCF program package.  J. Chem. Phys., 153(2):024109, 07 2020.
  * [10] Robert E. Thomas, Daniel Opalka, Catherine Overy, Peter J. Knowles, Ali Alavi, and George H. Booth.  Analytic nuclear forces and molecular properties from full configuration interaction quantum Monte Carlo.  J. Chem. Phys., 143(5):054108, 08 2015.