Reducing the cost of neural network potential generation for reactive
  molecular systems
Although machine-learning potentials have recently had substantial impact on
molecular simulations, the construction of a robust training set can still
become a limiting factor, especially due to the requirement of a reference ab
initio simulation that covers all the relevant geometries of the system.
Recognizing that this can be prohibitive for certain systems, we develop the
method of transition tube sampling that mitigates the computational cost of
training set and model generation. In this approach, we generate classical or
quantum thermal geometries around a transition path describing a conformational
change or a chemical reaction using only a sparse set of local normal mode
expansions along this path and select from these geometries by an active
learning protocol. This yields a training set with geometries that characterize
the whole transition without the need for a costly reference trajectory. The
performance of the method is evaluated on different molecular systems with the
complexity of the potential energy landscape increasing from a single minimum
to a double proton-transfer reaction with high barriers. Our results show that
the method leads to training sets that give rise to models applicable in
classical and path integral simulations alike that are on par with those based
directly on ab initio calculations while providing the computational speed-up
we have come to expect from machine-learning potentials.
Krystof Brezina, Hubert Beck, Ondrej Marsalek
# Reducing the cost of neural network potential generation for reactive molecular systems

###### Abstract

Although machine-learning potentials have recently had substantial impact on molecular simulations, the construction of a robust training set can still become a limiting factor, especially due to the requirement of a reference _ab initio_ simulation that covers all the relevant geometries of the system. Recognizing that this can be prohibitive for certain systems, we develop the method of transition tube sampling that mitigates the computational cost of training set and model generation. In this approach, we generate classical or quantum thermal geometries around a transition path describing a conformational change or a chemical reaction using only a sparse set of local normal mode expansions along this path and select from these geometries by an active learning protocol. This yields a training set with geometries that characterize the whole transition without the need for a costly reference trajectory. The performance of the method is evaluated on different molecular systems with the complexity of the potential energy landscape increasing from a single minimum to a double proton-transfer reaction with high barriers. Our results show that the method leads to training sets that give rise to models applicable in classical and path integral simulations alike that are on par with those based directly on _ab initio_ calculations while providing the computational speed-up we have come to expect from machine-learning potentials.

## I Introduction

Owing to the detailed atomistic insight into the structure and dynamics of molecular systems and materials, the relevance of computer simulations of molecular dynamics (MD) in current research is undeniable. MD simulations represent a valuable analytic and predictive tool in multiple fields of both basic and applied research including physical chemistry, material science, or drug design.[1; 2; 3; 4; 5] They also provide a way to explain and corroborate experimental data that might be difficult to interpret otherwise. For many systems of interest, MD simulations can be routinely performed under the Born-Oppenheimer approximation in the electronic ground state, with the nuclei being treated either classically or quantum-mechanically within the imaginary-time path integral formalism. This makes the choice of the potential energy surface (PES) a key decision that determines the accuracy of the resulting simulation. Out of the available options, _ab initio_ molecular dynamics[6] (AIMD) is a state-of-the-art methodology that relies on a full, on-the-fly quantum electronic structure calculation [7; 8] at every step of the simulation to evaluate the potential energy and forces. This is most commonly performed at the level of density functional theory[9; 10; 8] (DFT), which provides correlated electronic energies at a computational cost accessible in practice, but for smaller systems, the use of correlated wave function methods is feasible as well.[11; 12; 13] In any case, the computational cost of AIMD simulations is typically large -- especially so for advanced hybrid DFT functionals in the condensed phase -- and can easily become prohibitive in the light of the ever-growing demand for larger time and length scales of the relevant simulations.

This issue can be mitigated by the recent development of the so-called machine learning potentials (MLPs).[14; 15] These use various machine learning approaches[16; 17; 18; 19; 20] to faithfully approximate the desired _ab initio_ PES by training on a reference data set consisting of a relatively modest number of _ab initio_ geometries and their corresponding energies and, optionally, forces. As such, they indeed combine the best of the two worlds: they are able to maintain the accuracy of the parent _ab initio_ method, but they also circumvent the need for explicit electronic structure calculations at each step of the MD simulation. Thus, they evaluate the potential energy and forces at a significantly reduced computational cost.[21] One particular flavor of MLPs of major practical importance is represented by neural network potentials (NNPs), which rely on artificial neural networks combined with a set of appropriate atomic descriptors to accurately represent the molecular geometry-to-energy relationship, including all its symmetries.[16; 22] NNPs have repeatedly proved their worth in modeling a plethora of various molecular
systems ranging from liquids and solutions to interfaces and solids.[20; 23; 24; 25; 26] Our recent work focusing on NNPs[27; 26] shows that rather than using a single NNP to represent the PES, it is advantageous to build a model as a committee[28] of NNPs (C-NNP) that comprises a small number of NNPs, each trained individually to a subset of the main training set. The advantage is twofold: first, the energy prediction obtained as the committee average is known to be a better approximation of the _ab initio_ energy than the estimates of the individual members.[29; 30] Second, the committee disagreement,[31] represented by the standard deviation of the individual member estimates of energies or forces, serves as a valuable indicator of prediction reliability and can be used to monitor and optionally ensure the stability of the simulation.[27] Crucially, this disagreement can be used as the key ingredient of the active learning process called query by committee[32] (QbC) that systematically builds the training set in a data-driven way.[31; 27; 33]

An accurate and stable NNP can only be obtained on a foundation of robust, high-quality training data. This is typically based on a reference AIMD trajectory, from which geometries are selected for the training set, together with the corresponding energies and forces. However, the trajectory is highly correlated in time and thus most of the expensive AIMD data does not contribute useful information for the training of the model. This selection has been approached in different ways from random sampling and manual selection to more data-driven procedures,[25; 34; 35; 36; 37; 38; 39; 40; 41] with QbC being a particularly efficient method. QbC considers a set of candidate structures, in this case the whole AIMD trajectory, and iteratively builds up the training set. It starts by training a C-NNP on a very small set of initial configurations and using its disagreement to screen the candidate configurations for those with the most uncertain prediction. A small number of these configurations are then added to the training set, a new C-NNP is trained, and the process is iteratively repeated until some convergence criteria are met. In comparison to random selection, this approach is known to generate more compact training sets that give rise to robust models of similar accuracy.[35; 37] Even though the initial AIMD trajectory is typically the most expensive part of the procedure, numerous successful MLPs have been generated on top of reasonably short AIMD simulations.[26]

However, for many purposes, this process involving AIMD is still too expensive to be practical. For instance, the requirements on a high-level electronic structure method can raise the computational demands above a reasonable threshold. One might also be interested in a system that features rare events, such as chemical or conformational changes, which will happen quickly and occur infrequently or not at all in a direct AIMD simulation. In turn, these crucial configurations are underrepresented in the set of candidates and enhanced sampling simulations would be required in order to construct a robust training set, which typically raises the computational cost further by one or more orders of magnitude.

In case such a situation occurs, one needs to adhere to an approximate method of candidate generation that relieves some of the computational expenses while maintaining the quality of the resulting candidate set. For simple systems with a single potential energy minimum, the solution is fairly straightforward. In this case, one can benefit from random sampling of displacements in the directions of a fixed set of normal modes to obtain a set of distorted configurations. This approach, sometimes called normal mode sampling (NMS) in the literature,[18] avoids the cost of a full AIMD simulation by replacing it with a more manageable combination of a Hessian matrix evaluation and a number of single-point electronic structure calculations for the generated geometries. The sampling of the known normal mode distribution itself yields uncorrelated samples by definition and requires no _ab initio_ calculations, therefore its cost is negligible. Various versions of this approach were successfully used to generate structures for the training of MLPs. Using a scaled uniform random sampling of the normal modes, the method was used to obtain auxiliary structures used in model validation[42] and with approximate thermal distortions in NNP training set generation around configurational minima[18] as well as to construct an NNP model for a gas-phase ammonia molecule.[40] Clearly, the utility of NMS is limited when the harmonic approximation becomes insufficient. This can be the case if individual modes are strongly anharmonic or coupled, or if the system features conformational changes or reactions, where multiple local minima come into play.

In this work, we propose transition tube sampling (TTS), a robust and general approach to the generation of training sets and models that are able to accurately describe processes that feature transitions over potential energy barriers, which includes both conformational flexibility and chemical reactivity. We achieve this by generating thermally distorted candidate geometries along a reaction pathway with the help of multiple normal mode expansions and screening these candidates using QbC. The role of the minimum geometry in NMS is taken by the minimum energy path (MEP) that describes the course of the reaction through configuration space. Local harmonic expansions are performed in a small number of relevant configurations along the MEP and physically relevant candidate configurations are generated with uniform distribution along the MEP and with classical or quantum thermal weights in all perpendicular directions based on one of the sets of normal modes. An arbitrary number of these candidate configurations can be generated at a negligible computational cost and submitted to the QbC process, which selects the most important ones to have _ab initio_ calculations performed and to be included in the training set. This results in compact and robust training sets and models that maintain consistent accuracy along the reaction path, making them suitable for MD simulations of the reactive process, including enhanced sampling simulations, while no AIMD trajectories
are required as part of this process. We test this method on three different molecules in the gas phase to illustrate its capabilities.

The rest of the paper is organized as follows. In Section II, we begin by formalizing thermal NMS, which samples the exact classical or quantum canonical distribution under the harmonic approximation. With the obtained framework, we then proceed to introduce the MEP into the picture and describe the technical details of TTS. In Section III, we describe how we used TTS to create C-NNP models, the simulations performed with these models, and other related computational details. In Section IV, we apply this approach to three different gas-phase systems of increasing complexity represented by the molecules of benzene, malonaldehyde and 2,5-diaminonenzoquinone-1,4-diimine (DABQDI) and discuss its successes and possible pitfalls. Section V concludes the paper and offers outlooks concerning the generalization and the limitations of the method beyond the gas phase.

## II Theory

In this section, we discuss the theoretical basis of the TTS method. In this approach, we rely on the harmonic approximation and the vibrational normal mode formalism to obtain _ab initio_ training data for the construction of C-NNPs for reactive systems without the need to run expensive sampling simulations, such as full AIMD. First, we present the simple key idea behind NMS which relies on the harmonic approximation to describe the underlying PES and thus is expected to work well for systems that are close to harmonic around a single given minimum geometry at the temperature of interest. Clearly, this does not yield a flexible and general method, since the harmonic approximation is readily challenged by many realistic systems, notably those that exhibit more pronounced configurational flexibility or chemical reactivity. Therefore, we propose a more general approach to sampling candidate geometries based on NMS which is applicable even to systems described by multiple minima separated by barriers. This is achieved by using the harmonic expansion locally along an MEP in a way that eventually generates a balanced training set.

### NMS for thermal sampling around minimum geometries

To open the discussion of the theory behind TTS, we first turn our attention to the simple case represented by a PES with a single minimum geometry \(\mathbf{R}_{0}\) on which the nuclear motion is described by classical mechanics. Assuming a reasonable extent of validity of the harmonic approximation to capture the thermally accessible potential energy landscape, the classical thermal probability density \(\rho_{\mathrm{c}}\) at temperature \(T\) is approximated by

\[\rho_{\mathrm{c}}(\Omega_{1},\dots,\Omega_{N_{\mathrm{int}}})\propto\prod_{i=1 }^{N_{\mathrm{int}}}\exp\left(-\frac{1}{2}\beta\omega_{i}^{2}\Omega_{i}^{2} \right). \tag{1}\]

In this expression, \(\omega_{i}\) and \(\Omega_{i}\) denote, respectively, the natural frequency and the normal coordinate corresponding to the \(i\)-th normalized vibrational normal mode vector \(\mathbf{\Omega}_{i}\), and \(\beta\) is the inverse temperature equal to \(1/k_{\mathrm{B}}T\) (with \(k_{\mathrm{B}}\) representing the Boltzmann constant). \(N_{\mathrm{int}}\) is the total number of internal degrees of freedom of the species, typically \(3N-6\) for \(N\) atoms. Hence, in the harmonic approximation, the thermal density is described as a multivariate, yet uncoupled normal distribution where each \(i\)-th orthogonal degree of freedom has the standard deviation of \(\sigma_{i}=1/\sqrt{\beta\omega_{i}^{2}}\).

As such, it is straightforward to generate completely uncorrelated thermal geometries \(\mathbf{R}\) by distorting the minimum geometry \(\mathbf{R}_{0}\) independently in the direction of each of the normal mode vectors. The appropriate magnitude of the distortions is given by a randomly generated value of the corresponding normal coordinate \(\Omega_{i}\) from the distribution in Equation 1. The instrumental prescription for this procedure is the inverse coordinate transformation from normal modes back to Cartesian coordinates

\[\mathbf{R}=\mathbf{R}_{0}+\mathbb{M}^{-\frac{1}{2}}\sum_{i=1}^{N_{\mathrm{int }}}\Omega_{i}\mathbf{\Omega}_{i}, \tag{2}\]

where \(\mathbb{M}\) represents the diagonal mass matrix. Thus, by drawing samples of normal coordinates and transforming them, we obtain correctly distributed thermal samples in Cartesian coordinates.

We can now perform thermal NMS by sampling from this auxiliary harmonic ensemble as a source of candidate geometries to be potentially included in the training set of an MLP. The auxiliary ensemble is thus never used directly and no expectation values are calculated over it. It only needs to provide good coverage of the thermally accessible region of the PES, which will be the case as long as the harmonic approximation is reasonably accurate for the system of interest. Specifically, we construct a training set in our active learning procedure by generating a large number of these NMS candidate geometries and screening them in a QbC process using a C-NNP model. In each QbC iteration, electronic structure calculations are performed only for a small number of selected structures to obtain their potential energies and possibly forces, which then comprise the final training set once the process converges. The computational cost is thus determined primarily by the geometry optimization procedure, the Hessian calculation, the C-NNP prediction required for screening, and the electronic structure calculations for the selected geometries. The cost of the sample generation is negligible. This approach is substantially less computationally demanding when compared to the more conventional approach of sampling the candidate geometries for QbC from an AIMD trajectory, which requires a
large number of electronic structure calculations for very similar geometries that do not contribute diversity to the training set. In contrast to that, NMS generates fully decorrelated geometries by construction, and electronic structure calculations are only needed for the relatively small number of the most important geometries selected by the subsequent QbC process.

So far, we have focused on the situation where NMS is used to sample a classical distribution on the studied PES. However, since the harmonic approximation describes a molecule as a set of independent one-dimensional harmonic oscillators, we can readily generalize the above classical case to a quantum one as the analytic solution of the quantum harmonic oscillator is known. Specifically, it is straightforward to show (see Section S1 of the Supporting Information) that the canonical thermal density of a quantum harmonic oscillator at a given temperature is Gaussian just as its classical counterpart, but broader. This broadening is encoded in the quantum effective inverse temperature [43]

\[\beta^{*}(\beta,\omega)=\frac{2}{\hbar\omega}\tanh\left(\frac{\beta\hbar \omega}{2}\right) \tag{3}\]

at which a classical harmonic oscillator would have the same thermal width as a quantum harmonic oscillator at a reference inverse temperature \(\beta\). Since \(\beta^{*}\) is by definition a frequency-dependent quantity, one cannot describe the whole molecule by a single quantum effective temperature, but instead has to assign one to each individual mode. In turn, the quantum thermal density is given by

\[\rho_{\rm q}(\Omega_{1},\ldots,\Omega_{N_{\rm int}})\propto\prod_{i=1}^{N_{\rm int }}\exp\left[-\frac{1}{2}\beta^{*}(\beta,\omega_{i})\omega_{i}^{2}\Omega_{i}^{ 2}\right]. \tag{4}\]

This simple modification allows one to generate an auxiliary quantum ensemble at practically the same cost as the classical one that would otherwise need to be approached from a significantly more demanding perspective, perhaps based on sampling techniques using the imaginary time path integral formalism.

### Transition tube sampling

Up to this point, we have relied on the ability of the harmonic expansion around a single minimum to approximate the real PES so that the generated samples cover sufficiently all the relevant regions for the purpose of generating an MLP. Arguably, this is a reasonable requirement for most stable molecules with a single minimum geometry where the onset of the anharmonic region connected to the dissociation of the molecule is not thermally accessible. On the other hand, it is a stringent requirement for molecular systems which display conformational changes or chemical reactivity and, therefore, are represented by multiple PES minima connected by MEPs: features not captured by a single harmonic expansion. However, in such cases, it is desirable for the resulting model to be able to describe the potential energy landscape not only around local potential energy minima, but also in the transition regions. This is vital in the case of low-\(k_{\rm B}T\) barriers, where spontaneous transitions occur during direct MD. Nonetheless, it cannot be omitted even in the case of high-\(k_{\rm B}T\) barriers where an enhanced-sampling simulation would be required to cross the barrier. Even if the transition does not actually occur, the presence of the transition state may introduce substantial anharmonicity within the original PES basin that an eventual MLP should learn. However, in the case of barrier transitions it is not desirable to attempt to build the C-NNP model starting from a candidate set representing the true thermal ensemble, even if we could obtain it, since this would lead to a possibly detrimental under-representation of the high-energy configurations close to the transition state in the resulting

Figure 1: A schematic depiction of the TTS approach proposed for reactive systems. Panel A: An illustrative MEP winding through a model two-dimensional configuration space. Panel B: Control points are selected and their local normal modes (gray arrows) are calculated. Here, the control points are taken as the two endpoint minima on the MEP (purple and brown dots) and the transition state (yellow–green). Panel C: A much denser set of uniformly distributed reference geometries is generated along the MEP. Panel D: Each of the reference geometries is distorted using the local modes of their assigned control point (as detailed in panel E) to become a candidate geometry. Panel E: A detailed view showing how each reference geometry is assigned to a set of local modes. A set of reference geometries on the MEP is assigned to each control point following Equation 5. At the MEP edges, standard Gaussian thermal NMS is performed outside of the reaction coordinate (decaying purple and brown tails of the total density).
candidate set and, in turn, to poor performance of the resulting model in the transition regions.

Therefore, we propose the TTS method: a generalization of NMS for systems with transitions that employs local normal modes along an MEP to sample uniformly along the path and with proper thermal weights in all perpendicular directions. This leads to an auxiliary harmonic ensemble that differs significantly from the true thermal one but enables construction of MLPs with uniform accuracy along the whole transition. The TTS method naturally reduces to thermal NMS as described above for single-minimum systems in the zero MEP length limit. The process, illustrated in Figure 1, starts by finding the MEP \(\mathbf{R}(\xi)\) on the given PES (panel A). Here, \(\xi\) is a dimensionless reaction coordinate along the MEP curve through configuration space normalized to the interval from 0 to 1. In the following, we shall assume that the MEP is available as a continuous, differentiable function of the parameter \(\xi\). In practice, this can be achieved by spline fitting of the discretized representation of the MEP originating from, for instance, a nudged elastic band calculation.[44] Note that by definition, the MEP is a minimum of the PES in all directions perpendicular to it. Once the relevant MEP is known, we continue by selecting a sparse set of control points \(\mathbf{R}_{c},c=1,\ldots,N_{\mathrm{p}}\) along the MEP at positions \(\xi_{c}\) for which the Hessian matrices are calculated and diagonalized to give the set of local normal mode vectors \(\mathbf{\Omega}_{c,i}\) and their corresponding frequencies. For instance, this can be the two end-point minima and the transition state between them (Figure 1, panel B), although there is no constraint on how densely one might select the control points along the MEP other than the limiting computational expense of the Hessian matrix calculation. Formally, the expansion of the PES along the MEP becomes exact under the harmonic approximation in the limit of a large number of control points \(N_{\mathrm{p}}\). Since we want to achieve uniform sampling along the MEP, we now proceed to the generation of reference geometries on the MEP that do have this property. Specifically, to each control point \(\mathbf{R}_{c}\) we first assign a probability distribution \(p_{c}(\xi)\) defined as

\[p_{c}(\xi)=\begin{cases}\sin^{2}\left[\frac{\pi}{2|\xi_{c}-\xi_{c-1}|}(\xi-\xi_ {c-1})\right],\xi_{c-1}\leq\xi\leq\xi_{c}\\ \cos^{2}\left[\frac{\pi}{2|\xi_{c}-\xi_{c+1}|}(\xi-\xi_{c})\right],\xi_{c}\leq \xi\leq\xi_{c+1}\\ 0\text{ otherwise}\end{cases} \tag{5}\]

for all \(\xi\in[0,1]\) so that the identity

\[p(\xi)=\sum_{c=1}^{N_{\mathrm{p}}}p_{c}(\xi)=1 \tag{6}\]

holds over the whole length of the MEP (Figure 1, panel E over the range of \(\xi\)). Note that the choice of squares of harmonic functions is only one out of many possibilities, and any other pair of complementary functions that sum up to unity would work in this case. Next, we generate an arbitrary number of reference geometries \(\mathbf{R}_{0}(\xi)\) at a chosen linear density by drawing random values of \(\xi\) from the above distributions and passing them to the continuous prescription of the MEP, all while keeping track of the parent \(c\)-th control point (Figure 1, panel C). Analogously to the distortion of the minimum geometry in the single-minimum case through Equation 2, we distort each of these reference geometries using the normal modes and frequencies of its parent control point using

\[\mathbf{R}=\mathbf{R}_{0}(\xi)+\mathbb{M}^{-\frac{1}{2}}\sum_{i}{}^{{}^{\prime }}\Omega_{c,i}\left[1-\mathbb{P}(\xi)\right]\mathbf{\Omega}_{c,i}, \tag{7}\]

where the normal coordinate values \(\Omega_{c,i}\) are sampled thermally according to Equation 1 or 4 (Figure 1, panel D); the prime indicates that modes with imaginary frequencies are omitted from the sum. The matrix \(\mathbb{P}(\xi)\) is the projector on the tangent direction at point \(\xi\) which can be constructed analytically from \(\mathrm{d}\mathbf{R}(\xi)/\mathrm{d}\xi\). This is used to obtain distortions strictly perpendicular to the MEP and thus to correct for the approximate validity of the normal mode expansion calculated at \(\xi_{c}\) for all the displaced geometries. However, the use of the decaying probability distributions (Equation 5) favors the use of the local modes close to their origin. Through this procedure, one obtains a set of candidate geometries distributed inside a tube around the MEP the width of which is given thermally. At this point, this tube still has open ends cut sharply by the planes defined by normal vectors equal to the MEP tangent vector at the endpoints of the MEP. Since these points are (usually) also well-defined minima on the PES, the presence of these sharp edges is easily sanitized by appending the usual thermal NMS samples at these minima, although only adding the configurations away from the MEP (Figure 1, decaying gray dashed lines). In other words, just one half of the multivariate Gaussian is appended to the tube that does not overlap with it. In our TTS implementation, we ensure that the uniform density of the sampling along the MEP and the one at the peak of the half-Gaussian are seamlessly matched (as described in Section S1 of the Supporting Information).

Using the described sampling approach leads to an auxiliary ensemble of candidates that does not correspond to the true thermal ensemble, but contains a balanced selection of geometries distributed uniformly along the MEP with classical or quantum thermal displacements around it. Just like with plain NMS, we submit these samples as candidates to the QbC procedure, where in each iteration a large number of them is screened and a small number of those is selected to be included in the training set. _Ab initio_ calculations are only required for these selected geometries. This enables the building of diverse training sets in which all representative structures that might be encountered in a future simulation are contained so that the resulting model is, in fact, able to accurately describe the PES along the whole MEP, even in regions that have negligible thermal populations. Similar to NMS, the computational cost of TTS is de
termined primarily by the MEP optimization procedure, the Hessian calculation, the C-NNP prediction required for screening, and the electronic structure calculations for the selected geometries. In general, this can be expected to be substantially less computationally demanding than executing direct, or even enhanced-sampling, classical or path integral AIMD simulations and sampling from their trajectories.

## III Computational details

### Ab initio electronic structure

Two different levels of electronic structure theory were used in the simulations presented in this work. In both cases, we used the implementation provided by the CP2K software package [45] with its Quickstep DFT module.[46; 47] We described the electronic structure of the benzene molecule in the gas phase at the self-consistent charge density-functional tight binding [48] (SCC-DFTB) level with third-order diagonal corrections. The system was enclosed in a 10 A wide cubic box with open boundary conditions. For malonaldehyde and DABQDI systems, we used the revPBE0-D3 hybrid density functional [49; 50; 51; 52] combined with the TVZ2P Gaussian basis set [53; 46; 54] to represent the molecular orbitals and a plane wave basis with a 600 Ry cutoff to represent the density. The core electrons of the heavy atoms were represented using Goedecker-Tetter-Hutter pseudopotentials.[55] In addition, we used the auxiliary density matrix method [56] with the cpFIT3 fitting basis set for the DABQDI molecule. Both systems using hybrid DFT were centered in a 15 A wide cubic box with open boundary conditions and the wavelet Poisson equation solver.

### C-NNP model generation

Throughout all of our investigations, we used committees consisting of eight different Behler-Parrinello neural network potentials,[16] where for each of them, a different initialization of weights and a different 90 % subset of the full training data set was used to ensure a diverse committee. The models consisted of two hidden layers of 20 nodes each and were trained using the multistream [57] adaptive extended Kalman [58; 59] filter with 32 streams. The input features were a standard set of atom-centered symmetry functions.[26] The training of the model was done using the n2p2 package [57] and the selection of training structures by QbC was done with a development version of our AML package.[60] For benzene, 20 structures were randomly sampled initially and in each of the 40 QbC iterations, 10 new structures with the highest committee force disagreement were added to the data set, for a total of 420 training geometries. The final NNPs were trained for 2000 epochs. For the first generation of malonaldehyde and DABQDI models, 20 initial structures were sampled randomly and then the 15 structures with the highest disagreement were selected in each of the 40 QbC iterations, for a total of 620 training geometries. For malonaldehyde, where additional generations of models were required (as detailed in Section IV), the original training set was supplemented by additional structures QbC-sampled from an MD trajectory which was produced using the previous C-NNP model. Here, 15 structures were added in each iteration until the force committee disagreements of the selected structures and the remaining candidate structures were similar. All reference calculations were done using CP2K and the electronic structure settings described above.

### Geometry optimization and vibrational analysis

The optimization of the minimum reference geometries for benzene and DABQDI was executed natively in the CP2K software. It was performed using the BFGS optimizer [61] combined with threshold criteria of 0.07 eV A\({}^{-1}\) for the maximum change in force components, 0.009 A for the change in atomic positions and 0.13 eV for the change in total energy. For the malonaldehyde molecule, we employed the Atomic Simulation Environment (ASE) [62] together with CP2K and performed the optimization using the FIRE optimizer [63] while specifying only a force criterion of 0.01 eV A\({}^{-1}\). Additional constrained optimizations in the case of DABQDI needed for the relaxed PES scan were performed using the constraint functionality provided by ASE together with its FIRE optimizer. The Hessian matrix evaluation on the optimized structures was performed in each case using CP2K and a Cartesian atomic displacement of 0.0005 A.

### Nudged elastic band calculations

The relevant MEPs needed for the TTS procedure were obtained through the climbing-image [64] nudged elastic band [44] (CI-NEB) optimization procedure as implemented in CP2K. The initial band geometries in this work consisted of 15 replicas of the molecule in question including the two fixed, pre-optimized endpoints, and were obtained through linear interpolation. The spring constant of the harmonic links between the neighboring replicas was kept constant at the value of 4.86 eV A\({}^{-2}\). We used a force convergence criterion of 0.007 eV A\({}^{-1}\) and the minimization of the band energy was performed using a DIIS optimizer.

### MD simulations

All MD simulations involving both _ab initio_ as well as C-NNP potentials [30] were run using the CP2K package.
The simulations with classical representation of the nuclei were propagated at a temperature of 300 K using a time step of 0.5 fs to numerically integrate the Langevin equation with the friction coefficient \(\gamma\) of 0.02 fs\({}^{-1}\) to achieve canonical sampling. The path integral simulations that include nuclear quantum effects were performed using imaginary-time ring polymers consisting of 64 replicas using the RPMD propagator. The canonical distribution at 300 K was sampled using the path integral Langevin equation thermostat [65] with the time constant for the centroid motion of 200 fs while the integration time step was kept at 0.25 fs.

### Umbrella sampling

The initial conditions for each umbrella sampling window were extracted from a steered MD trajectory, which was performed in the CP2K v2022.1 software package combined with the PLUMED plugin [66, 67, 68]. In this case, the value of the proton-sharing coordinate \(\delta_{1}\) (as detailed in Section IV) was biased from \(-1.2\,\mathrm{\SIUnitSymbolAngstrom}\) to \(1.2\,\mathrm{\SIUnitSymbolAngstrom}\) during a 10 ps long simulation using a moving harmonic restraint with the force constant \(\kappa\) of 500.0 kJ mol\({}^{-1}\) A\({}^{-2}\). The simulation was performed classically with an integration time step of 0.5 fs in the canonical ensemble at 300 K using a local CVR thermostat [69] with a time constant of 50 fs.

30 equidistant umbrella sampling windows separated by 0.08 A were set up from the above steered MD simulation. Individually in each window, the value of \(\delta_{1}\) was biased by a static harmonic restraint of 500.0 kJ mol\({}^{-1}\) A\({}^{-2}\) and simulated for 50 ps using the same setup as for the steered MD simulation above. The overlap of the corresponding histograms of \(\delta_{1}\) values observed in each simulation window is shown in Section S2 of the Supporting Information. The value of \(\delta_{2}\) was kept unbiased in each simulation window but was monitored for use in the following analysis. The biased configurations were reweighed to the unbiased ensemble using a Python implementation of the multistate Bennet acceptance ratio [70, 71] (MBAR) procedure to obtain both a one-dimensional free energy profile for the proton-transfer along \(\delta_{1}\) as well as a two-dimensional free energy surface showing the dependence on both proton-sharing coordinates. This was done by determining the thermal weight associated with each configuration in the biased simulations and using these to obtain the probability distribution in the \(\delta_{1}\), \(\delta_{2}\) subspace, and from that the corresponding free energy surface.

## IV Results and discussion

To showcase the performance of the TTS procedure in the creation of models for realistic potentials, we select three different gas-phase molecules with increasing complexity of their PES. We begin with benzene, which represents a single-minimum system with a close-to-harmonic potential at room temperature and thus allows us to illustrate the simple thermal NMS procedure. This is followed by a study of the end form of 1,3-propanediol (malonaldehyde), which exhibits reactivity by sharing the acidic proton between the two oxygen atoms spontaneously at ambient conditions [72]. Finally, we focus on a more involved proton-sharing system represented by 2,5-diaminobenzzoquinone-1,4-diimine (DABQDI), which has two proton-sharing sites [73]. Spontaneous proton transfer is hindered by a barrier thermally insurmountable at room temperature, and an enhanced sampling simulation is necessary to determine the free energy profile.

### Benzene

To lead off the discussion of the ability of TTS to seed a training set for the creation of C-NNP models for realistic systems in the gas phase, we focus on the benzene molecule. It represents an ideal example to illustrate the basic idea of thermal NMS using a single normal mode expansion at an optimal geometry, since it features a single configurational minimum and the surrounding PES exhibits almost no anharmonic effects.

To prepare the ground for comparison with the relevant C-NNP data, we initially performed one 250 ps AIMD simulation of gas-phase benzene at 300 K at the DFTB level using a classical representation of the atomic nuclei as well as a 100 ps PIMD simulation using 64 replicas to approximate the imaginary time path. Two C-NNP models were then based on candidate sets obtained from a thermal NMS of gas-phase benzene using a Hessian matrix calculated at the same DFTB level of theory as the (PI)-AIMD simulations at 300 K for the classical model and with the appropriate effective temperatures at 300 K for the quantum one. The resulting models were evaluated on test sets consisting of 1000 structures sampled from the two AIMD trajectories. Both models performed very well with an energy root mean square error (RMSE) of 1.66 meV and 5.90 meV for the model constructed for the use without and with path integral structures, respectively. The RMSE for a single force component was \(14.9\,\mathrm{meVA}^{-1}\) and \(30.4\,\mathrm{meVA}^{-1}\). Subsequently, the models were used to obtain new 500 ps long MD and 100 ps PIMD simulations at 300 K.

The comparison of the C-NNP models to the corresponding (PI)-AIMD trajectories in terms of molecular geometry properties is summarized in Figure 2. In general, we can see the expected broadening of probability distributions due to nuclear quantum effects when we compare the left and right columns of Figure 2. In both the classical and quantum case, we observe a perfect match between the _ab initio_ (green shading) and C-NNP distributions (blue dashed lines) in C-C bond lengths (panels A and B), C-H bond lengths (panels D and E), C-C-C angles (panels E and F), and C-C-C dihe
drals (panels G and H). The two types of covalent bonds have expected distributions; the mean of the C-C-C angle is located at 120\({}^{\circ}\) which shows the average hexagonal arrangement of the aromatic ring subject to planarity, which is, in turn, demonstrated by the (signed) C-C-C dihedral angle peaking at 0\({}^{\circ}\) as expected. This level of agreement suggests that the final models used for production MD represent excellent approximations to the original DFTB PES. Additionally, we show the distributions of the NMS structures (orange dotted lines) alongside the anharmonic distributions. These exhibit significant overlap with both the (PI)-AIMD and C-NNP data. This suggests that the harmonic approximation to the original ensemble is relatively good and confirms the assumed high degree of harmonicity of the 300 K gas-phase benzene PES, even in the quantum case. However, note that the match of the NMS data with the (PI)-AIMD data is not nearly as perfect as that of the C-NNP data and certain deviations are, in fact, present. As discussed in Section II, these are to be expected since the NMS ensemble is only auxiliary and its goal is to provide sufficient coverage of the accessible PES which ultimately leads to an accurate C-NNP model. Using thermal NMS, we were thus able to construct a C-NNP that accurately describes the original PES of benzene based on a single Hessian evaluation and 420 single-point electronic structure calculations.

### Malonaldehyde

The end form of malonaldehyde is a simple organic molecule that has been used in MD simulations [72] to illustrate a simple intramolecular proton-transfer reaction in which the proton is moved from the end oxygen to the aldehyde oxygen with a simultaneous electronic rearrangement which altogether causes the reactant and the product to become symmetrically mirrored, chemically identical structures. As such, malonaldehyde is a convenient molecule to demonstrate the ability of TTS to describe a simple reaction.

With the aim to describe the proton-sharing process accurately, we decided to model the original _ab initio_ PES at the hybrid DFT level using the revPBE0-D3 functional. Thanks to the ability to use quantum normal coordinate distributions in the TTS method, we produced both classical and quantum models and classical and PIMD trajectories for malonaldehyde to test and showcase this functionality. However, we focus mostly on the classical case in the main text and discuss the complementary quantum results, for which qualitatively similar conclusions arise, in the Supporting information, Section S2. The starting point of the TTS procedure is the proton-sharing MEP, which was discretized into 15 replicas and optimized using the CI-NEB procedure and the revPBE0-D3 density functional. For illustration purposes, we decided to use the full MEP with both the reactant and product explicitly represented: this is strictly speaking not necessary since the symmetry of the reaction allows the use of only one non-trivial half of the MEP for TTS. Out of the optimized full-length MEP, three control points were selected in the two minima (reactant and product) and in the transition state. For the visualization of the multidimensional configuration data, we choose the reduction into a 2D space of two geometric parameters: the proton-sharing coordinate \(\delta(\mathbf{R})=|\mathbf{R}_{\mathrm{O}}-\mathbf{R}_{\mathrm{H}}|-|\mathbf{ R}_{\mathrm{O^{\prime}}}-\mathbf{R}_{\mathrm{H}}|\) and the oxygen-oxygen \(d_{\mathrm{OO^{\prime}}}(\mathbf{R})=|\mathbf{R}_{\mathrm{O}}-\mathbf{R}_{ \mathrm{O^{\prime}}}|\) where O and O' denote the two oxygen atoms that share the proton H. The parameters are illustrated in the snapshot on the left of Figure 3. The obtained optimized MEP and the selected control

Figure 2: Thermal geometry properties of benzene in the gas phase at 300 K from classical MD (left column) and path integral MD (right column) compared between simulations using the reference DFTB potential, the harmonic TTS ensemble, and simulations using a C-NNP model building on the thermal NMS geometries. Panels A and B show the distribution of C–C bond lengths, panels C and D the distribution of C–H bond lengths, panels E and F the distribution of C–C–C angles, and, finally, panels G and H the distribution of the C–C–C dihedral angles.
points in this representation are shown in the top left panel of Figure 3 in orange. The TTS classical candidate structures were then generated using the procedure outlined in Section II at the temperature of 300 K, linear sampling density along the MEP of \(1\times 10^{5}\,\mathrm{\SIUnitSymbolAngstrom}^{-1}\) and matched-density sampling at the minima. The distribution of the obtained configurations is shown in the top left panel of Figure 3 as blue contours. The same distribution colored by the assignment of each candidate to the individual control points (corresponding to the situation shown in panel D of Figure 1) is shown in Section S2 of the Supporting information. Note that this particular presentation of the data does not do justice to the uniformity of the sample distribution along the MEP as the regions around the minima seem to be more populated than the transition state. This is an effect of the deformation of the configuration space by the projection on the selected subspace; the samples are distributed uniformly in the full dimension. After passing the resulting set of candidates through the QbC selection and training a C-NNP model on the obtained training set, the model was used to run a direct 250 ps MD simulation of gas-phase malonaldehyde at 300 K. A subset of the obtained configurations is shown in the form of a scatter plot in the top right panel of Figure 3 colored by the decadic logarithm of the norm of the force committee disagreement on carbon atoms in the usual \(\delta\) and \(d_{\mathrm{OO}^{\prime}}\) representation. Additionally, a 1D free-energy profile obtained by a Boltzmann inversion of the probability density of configurations along the \(\delta\)-axis is shown in the bottom panel of Figure 3; the height of the barrier is approximately 120 meV which corresponds to roughly 5 \(k_{\mathrm{B}}T\) at 300 K. This accounts for the expected low, but existing population surrounding the transition state at \(\delta=0\,\mathrm{\SIUnitSymbolAngstrom}\). Alongside the free energy profile, we show the corresponding average potential energy as a function of \(\delta\).

An important observation can be made from the presented data. By comparison of the two distributions in the top panels of Figure 3, it is clear that the TTS distribution populates a smaller volume of the configuration space than the data obtained from the MD simulation. In this particular case, it means that TTS does not directly provide good enough coverage and the resulting model is undertrained in the missing, yet thermally accessible regions. Specifically, the C-NNP model performs poorly in the large \(d_{\mathrm{OO}^{\prime}}\) tails of the shown distribution, as quantified by the larger force disagreement values. On the contrary, in regions around the proton-sharing MEP, the coverage is good and the force disagreement remains small, despite the tiny thermal population. Regardless of the elevated model uncertainty for some configurations, these MD simulations remain stable. The observed increased disagreement can be interpreted in the following way: going in the opposite direction from the minima as the proton-sharing MEP, the true anharmonic PES has a potential wall that grows slower than the harmonic wall captured by the TTS distribution and, therefore, the thermal coverage of the TTS configurations cannot reach far enough. In principle, this behavior is caused by either the strongly anharmonic character of the chemical bonds leading to bond dissociation or, more likely in this case, the presence of another reactive process leading to a new transition state. In the following two paragraphs, we present two possible solutions to the issue. The first one relies on an active-learning-based iterative improvement of the model, which has the advantage of requiring no knowledge of the origin of the anharmonicity, but is tedious to perform since several intermediate simulations and model generations need to be created. Meanwhile, the other solution relies on the chemical intuition of the user to identify the reactive nature of the issue with the aim to extend the initial TTS, which leads to a fully capable C-NNP model straight away.

The QbC process can be used to fill in an already existing training data set that has gaps, perhaps due to an incomplete TTS candidate set in the QbC selection for the initial model. We illustrate this process in Figure 4, where the left panel shows the same data as the top right panel of Figure 3 as a starting point. Regions of configuration space not covered well in the training

Figure 3: The TTS sampling of the malonaldehyde proton transfer and the MD simulation with the resulting C-NNP model. The top left panel shows the relevant MEP (orange) with the three selected control points in the two configurational minima and the transition state highlighted and the distribution of the TTS geometries (blue). The top right panel shows a scatter plot of a subset of geometries (selected with a stride of 37.5 fs) obtained during a 250 ps MD simulation using the C-NNP model built on top of the TTS ensemble. Each point is colored by the norm of the force committee disagreement on the carbon atoms and the mean of the quantity is shown in the box. Note the high force disagreement in the high \(\delta\) tails of the distribution. The bottom panel shows the Boltzmann-inverted free energy profile (red) and the corresponding binned average potential energy of the system (black, aligned to zero) along the proton transfer reaction as observed in the MD simulation.
set of this generation 1 model can be easily identified by the high committee disagreement, as can be seen in the tails of the distribution. Hence, one can start a new QbC using the existing training data set augmented by structures from an MD simulation performed with the initial C-NNP. Depending on the size of the gaps in the initial training data set, only a few iterations of QbC are typically necessary. However, adding these structures to the training data set can lead to substantial changes in the previously inaccurate regions of the potential energy surface, resulting in an MD simulation that again reaches new regions of the configuration space where the shape of the potential energy surface is yet unknown to the model and the committee disagreement is high. This can be seen in the generation 2 model shown in the middle panel of Figure 4. Therefore, multiple repetitions of the MD-QbC cycle might be necessary until a highly accurate model that exhibits low and uniform disagreements over the sampled data is reached, as is the case for the generation 3 model in the right panel of Figure 4. As such, the approach could become practically cumbersome when the regions of high disagreement coincide with regions of high free energy and long MD simulations are needed to uncover these structures, but nonetheless represents a general solution to the anharmonicity problem. Overall, repeating the cycle of sampling MD configurations with a given generation of a C-NNP model followed by training a new generation on a training set enhanced by high-disagreement QbC-selected structures from the previous MD simulations leads to a force disagreement that is lower in the problematic PES regions and, therefore, more uniform overall. In addition, we observe a decrease of the mean of the force disagreement of the sampled configurations due to the fact that the size of the training set increases in each generation. Specifically, 620 structures were used for the original model in the left panel of Figure 4, 800 structures for the model in the middle panel, and 950 structures for the last model in the right panel. This approach can be beneficial in systems where it is difficult to identify the origin of the anharmonicity of the original PES, but is rather demanding from the point of view of both computational requirements and user involvement due to the need for the semi-supervised iterative procedure.

However, in the case of malonaldehyde, the general active learning iterative procedure to improve the model might be excessive. The possible reasons for the softer-than-harmonic wall due to conformational flexibility are few and the particular direction against which the first generation of the model is pushing can be easily identified with the _s-cis_ and _s-trans_ torsional isomerism

which is mediated by rotation around the C-C single bond in the propane backbone. The optimized MEP corresponding to this torsion displays a perfect continuation in the correct direction when projected into the \(\delta\) and \(d_{\text{OO}^{\prime}}\) subspace, as shown in the top left panel of Figure 5 in orange. Although this pair of descriptors is not appropriate for the whole torsion MEP, which entails a more complicated motion, it is accurate enough at small deviations from the equilibrium geometry. To include structures along this MEP into the initial (first generation) proton-sharing TTS, two control points were chosen in the minimum (shared by the two MEPs) and in the new transition state (not shown in Figure 5, as it is around \(d_{\text{OO}^{\prime}}=3.8\,\text{\AA}\)). We do not need to use the _s-trans_ minimum at all, as we are not interested in including the

Figure 4: Evolution of the force disagreement of the carbon atoms through multiple instances of QbC. The left panel shows a subset of configurations originating from an MD simulation using a model trained on data selected directly from the TTS candidate set (identical data as in the top right panel of Figure 3 is shown). The force disagreement (depicted using the color scale) in the vicinity of the two configurational minima and along the proton-sharing reaction is adequately low, however, for structures with a high absolute \(\delta\) it is more than one order of magnitude higher. The central panel shows configurations and disagreements obtained from an MD simulation using a model trained on the initial training data augmented by QbC-selected high-disagreement configurations from the data in the left panel. In turn, the right panel shows data obtained by improving the model by the new data sampled in the simulation shown in the central panel. Most notably, structures at the tails of the populated configuration space are substantially improved. The mean force disagreement over all configurations in each data set is shown in the framed box in each panel. The \(\delta\) and \(d_{\text{OO}^{\prime}}\) coordinates are illustrated in the snapshot to the left of the panels.
transition itself, only the shape of the PES on the side of the global minimum. A new TTS was performed between these control points with the same parameters as the initial one and the resulting distribution of the combined sets of configurations is shown in blue in the top left panel of Figure 5. Running a 250 ps long MD simulation with a new C-NNP model trained on the QoC-selected training set from this combined candidate set leads to the distribution shown in the top right panel in Figure 5. Clearly, the distribution reaches all the expected thermally accessible regions, the force disagreement is evened out across the configurations, and the high-disagreement tails are no longer present. The mean value of the disagreement is comparable to that of the first-generation model, as these two models are based on training sets of the same size. The bottom panel of Figure 5 shows the potential energy and the free energy profile along \(\delta\), where the latter is shown for both classical and path-integral data. In all three curves, one can observe a softening of the barrier in the high \(|\delta|\) regions in comparison to the data shown in Figure 3 resulting from the present C-NNP being aware of the anharmonic nature of the PES in these regions. The transition-state free energy of 43 meV for the quantum simulation is \(\sim\)4 times lower than in the classical case, which is a manifestation of proton tunneling through the barrier. This is qualitatively consistent with existing literature but quantitatively different from the results reported for the BLYP functional, where the classical free energy barrier is somewhat higher and the quantum effect is substantially smaller, decreasing the barrier by a factor of only \(\sim\)2.[72]

For further insights, a test set independent of the training set data was created by generating 500 structures using TTS and sampling 500 structures from the classical MD trajectory shown in Figure 5, and evaluating their energies and forces with the revPBE0-D3 reference method. The generation 3 C-NNP of the iterative approach as well as the extended TTS C-NNP perform well with an energy root mean square error (RMSE) of 1.80 meV and 3.44 meV and a force component RMSE of 18.3 meV\(\text{\AA}^{-1}\) and 24.4 meV\(\text{\AA}^{-1}\), respectively. The slightly elevated RMSE of the extended TTS approach is due to the broader coverage of the training set. It includes a range of geometries along the C-C single bond torsion, even in regions that are not populated during MD, leading to a less dense coverage of the rest of the configuration space. The validation errors of the intermediate models of the iterative approach and the distribution of errors within configuration space are discussed in more detail in the Supporting Information, Section S2.

Both of the approaches above thus yield highly accurate models for the description of the proton-sharing reaction in malonaldehyde for classical and quantum nuclei. The advantage of one over the other therefore depends mostly on the specific situation in which the need for any of them should arise: if the reaction coordinate of the complementary transition can be identified, then the latter approach using the extended TTS reaches the desired result with higher efficiency. Note that this approach can also be used when multiple different transitions are to be included in a single model.

## Dabqdi

The most complex reactive system used to demonstrate the performance of the TTS method is represented by the DABQDI molecule. This nitrogenated benzoquinone derivative can exchange two protons between the neighboring amine and iming groups

\[\text{H}_{2}\text{N}\]

\[\text{NH}_{2}\]

\[\text{NH}_{2}\]
where each proton-sharing coordinate describes a single proton-sharing site, was obtained at the revPBE0-D3 level of electronic structure theory through a relaxed scan of the molecular potential energies while applying appropriate constraints and is shown in the left panel of Figure 6. The shown data was aligned so that the global minimum of the PES corresponds to the zero-energy level. The typical structure of the PES featuring four distinct configurational minima and four transition states corresponds to a sequential double proton transfer at the level of an MEP. Here, one proton is first fully exchanged to reach an intermediate state located at a higher potential energy and only then the other proton follows. The height of the potential energy barrier for this sequential process of roughly 0.8 eV indicates that its thermal rate should be negligible. The alternative concerted proton transfer path that is seen in other species including carboxylic acid dimers [74] is classically disallowed in this case by a tall (\(>\)1.2 eV) potential barrier in the middle of the presented PES which represents a second-order saddle point and, as such, no MEP can go through it. The symmetry of the DABQDI PES allows us to explicitly address only a single proton transfer: unlike in the previous example, we exploit this feature here for the C-NNP model generation. The relevant non-trivial MEP connecting the two chemically distinct minima was obtained using the CI-NEB optimization at the revPBE0-D3 level of theory and is shown in white in the left panel of Figure 6. From there, the straightforward TTS candidate generation was performed using the three usual control points in the two minima and in the transition state at 300 K with the linear sampling density of \(1\times 10^{3}\,\mathrm{\SIUnitSymbolAngstrom}^{-1}\). A subset of the candidate geometries is shown in the left panel of Figure 6 as green points. The obtained C-NNP model was used to recreate the 2D proton transfer PES which is shown in the middle panel of Figure 6 with the energies aligned in the same way as in the DFT case. Qualitatively, the C-NNP model captures all the features of the original _ab initio_ PES including the position of the minima, the transition state, and the barriers, as well as the potential energy values. Note that the good agreement in the representation of the central barrier in spite of the lack of corresponding geometries in the TTS candidates is due to the successful extrapolation by the model. The quantitative difference between the _ab initio_ and C-NNP potential energy landscape is captured in the right panel of Figure 6 and shows that the potential energies are reproduced as closely as approximately \(\pm 0.05\,\mathrm{eV}\). This difference could be decreased further, if desired, by changing the architecture of the C-NNP and by increasing the size of the training set beyond the current 620 structures.

Since DABQDI features barriers that are not practically accessible by direct MD, it serves as a useful example to illustrate the power of the TTS-based C-NNP model to perform an enhanced sampling calculation to correctly estimate the free energy profile of the double proton transfer at 300 K. This was obtained using an umbrella sampling simulation in the coordinate \(\delta_{1}\) with the C-NNP PES (as described in Section III) followed by a multistate Bennet acceptance ratio (MBAR) reweighting of the biased configurations. The obtained 1D free energy profile in \(\delta_{1}\) is shown in blue in the top panel of Figure 7. The transition state is located at roughly 0.8 eV above the global minimum. Comparing this with

Figure 6: Comparison between the reference _ab initio_ revPBE0-D3 and the C-NNP proton-sharing PESs of the DABQDI molecule. The left panel shows the projection is the reference DFT PES into the \(\delta_{1}\), \(\delta_{2}\) subspace using a color scale as well as individual isoenergetic contours. Furthermore, the minimal non-trivial MEP which describes a single proton transfer is depicted in white with the selected control points highlighted. A sparse subset of the TTS geometries sampled around the MEP at 300 K is shown in green. The middle panel shows the corresponding PES projection calculated with the resulting C-NNP. Finally, the right panel shows the difference between the two PESs and quantifies the maximum error of the C-NNP fit in the classically thermally accessible regions on the order of \(10^{-3}\) to \(10^{-2}\) eV. The two \(\delta\) coordinates are illustrated in the snapshot of the DABQDI molecule to the left of the panels.
the value of the corresponding potential energy suggests that the entropic contribution in the gas phase system is small and that the population at the barrier is clearly negligible at 300 K. To validate the obtained free energy profile, we perform a reweighing of a subset of the obtained configurations in each umbrella window to the original DFT ensemble. This is achieved by additionally multiplying each MBAR-obtained unbiased weight by the corresponding factor \(e^{-\beta\Delta E}\) where the energy difference \(\Delta E\) is the difference between the C-NNP and DFT potential energy for each configuration. For this purpose, we used a total of 3000 configurations obtained by selecting 100 geometries evenly spaced in time from each of the 30 umbrella sampling windows. The resulting profile, which is a good approximation to the full-DFT free energy profile, is displayed as the orange dashed line in Figure 7 and shows very good correspondence with the profile obtained using the C-NNP model alone. This procedure thus at the same time validates the C-NNP model and provides DFT data for a fraction of the cost of the hypothetical purely _ab initio_ enhanced sampling simulation. Monitoring the values of the collective variable \(\delta_{2}\) along the umbrella sampling simulation and using the thermal weights obtained from the MBAR treatment of the biased simulations also allows recovering the 2D free energy surface in \(\delta_{1}\) and \(\delta_{2}\). Its symmetrized version is shown in the bottom panel of Figure 7. It is worth noting that the computational acceleration of the C-NNP umbrella simulation in comparison to the naive execution with the original DFT method is substantial. To illustrate the computational savings, we can compare the times required for one MD step with the implementations in CP2K used in this work. With the hybrid functional, one MD step takes 272 s on a single core or 17 s on 32 cores (a full node) of our EPYC-based cluster. With the C-NNP, one step takes 0.006 s and does not scale meaningfully to more cores due to the small system size. This yields a speedup of \(\sim\)45000\(\times\) on identical resources or \(\sim\)2800\(\times\) with more resources given to the DFT calculation. Obviously, the specific numbers will depend on the details of the electronic structure setup and the MLP architecture used, as well as the specific implementations and hardware used, but this behavior of our particular setup should provide a general idea.

## V Conclusions

In this work, we have introduced the TTS method to sample thermal geometries around MEPs that describe barrier-crossing transitions in molecular systems. The goal of the method is to provide a physically meaningful set of candidate structures for the creation of MLPs without the need to run computationally demanding _ab initio_ simulations. In our case specifically, we submit these geometries to QbC and construct a C-NNP model, but the same candidates could be used for other types of models as well. The execution of the TTS protocol as a whole entails a relatively modest computational cost with respect to the original _ab initio_ method that is given by the MEP optimization, several Hessian evaluations, and a small number of single-point _ab initio_ calculations for the QbC-selected geometries. In terms of application to realistic systems, the TTS method yields highly accurate C-NNP models in all studied cases. This was achieved either by using the generated candidate set directly, or by letting the resulting C-NNP model undergo additional active learning generations to compensate for a pronounced anharmonic effect as seen in the case of malonaldehyde. As such, the performance of TTS in the presented test systems demonstrates its robustness and efficiency and suggests applicability in most gas-phase systems, including highly anharmonic cases.

A noteworthy feature of the TTS method is its ability to provide thermal geometries sampled from the quantum thermal distribution at essentially the computational cost of the classical case. As such, models that are appropriate for use in path integral simulations are made readily available without the need to run expensive PI-AIMD simulations at all. However, it is important to recognize

Figure 7: Umbrella sampling simulation of the single proton transfer in the DABQDI molecule along the \(\delta_{1}\) collective variable using the C-NNP potential. The top panel shows the obtained free energy profile in blue. For validation purposes, we also show the DFT free energy profile (orange, dashed) obtained by reweighting the C-NNP configurations as described in the text. Note that no umbrella sampling using the DFT potential was performed to obtain the DFT free energy profile. The bottom panel shows the full 2D free energy surface obtained by weighing the distribution in the two proton-sharing coordinates by the thermal Boltzmann factors extracted from the biased simulation and symmetrizing the resulting histogram.
that although the present formulation of TTS can address quantum behavior, it has limitations in this regard that derive from the fundamentally classical nature of the MEP. Nuclear quantum effects, in particular quantum tunneling through the potential barrier, can cause the configuration-space probability density of the system to deviate from the transition tube around the MEP in a way that renders the coverage by TTS samples insufficient.

To account for this, the above formulation of TTS can be straightforwardly generalized from sampling around classical MEPs to ring-polymer instantons, [75] which represent the paths of optimal tunneling. While this modification requires essentially no adaptation of the TTS theory and implementation itself, one can anticipate an elevated computational cost due to the required instanton optimization at the explicit _ab initio_ level. The approach will find applications beyond the gas phase, in systems where vibrational normal modes are a meaningful concept, such as crystals and surface-bound molecules. Since our research anticipates the need to address such systems in the near future, we expect TTS to be a valuable tool in the creation of accurate, yet computationally accessible potentials that will enable the accurate description of these more complex systems at unprecedented sizes and simulation time scales.

###### Acknowledgements.

 K.B. and H.B. acknowledge funding from grant schemes at Charles University, reg. n. CZ.02.2.69/0.0/0.0/19_073/0016935, and from the IMPRS for Quantum Dynamics and Control. O.M. acknowledges support from the Czech Science Foundation, project No. 21-27987S.

## References

* Section C **104**, 142-164 (2008).
* [2] Marco De Vivo, Matteo Masetti, Giovanni Bottegoni, and Andrea Cavalli, "Role of Molecular Dynamics and Related Methods in Drug Discovery," Journal of Medicinal Chemistry **59**, 4035-4061 (2016).
* [3] Scott A. Hollingsworth and Ron O. Dror, "Molecular Dynamics Simulation for All," Neuron **99**, 1129-1143 (2018).
* [4] Volker L. Deringer, Miguel A. Caro, and Gabor Csanyi, "Machine Learning Interatomic Potentials as Emerging Tools for Materials Science," Advanced Materials **1902765**, 1902765 (2019).
* [5] Nan Yao, Xiang Chen, Zhong Fu, and Qiang Zhang, "Applying Classical, Ab Initio, and Machine-Learning Molecular Dynamics Simulations to the Liquid Electrolyte for Rechargeable Batteries," Chemical Reviews **122**, 10970-11021 (2022).
* [6] Dominik Marx and Jurg Hutter, _Ab initio molecular dynamics: basic theory and advanced methods_ (Cambridge University Press, New York, 2009).
* [7] Attila Szabo and Neil S. Ostlund, _Modern quantum chemistry: introduction to advanced electronic structure theory_ (Dover Publications, Inc., Mineola, N.Y., 1996).
* [8] Robert G. Parr and Weitao Yang, _Density-functional theory of atoms and molecules_ (Oxford University Press, New York, 1989).
* [9] P. Hohenberg and W. Kohn, "Inhomogeneous Electron Gas," Physical Review **136**, B864-B871 (1964).
* [10] W. Kohn and L. J. Sham, "Self-Consistent Equations Including Exchange and Correlation Effects," Physical Review **140**, A1133-A1138 (1965).
* [11] Ondrej Marsalek, Pei Yang Chen, Romain Dupuis, Magali Benoit, Merlin Meheut, Zlatko Bacic, and Mark E. Tuckerman, "Efficient calculation of free energy differences associated with isotopic substitution using path-integral molecular dynamics," Journal of Chemical Theory and Computation **10**, 1440-1453 (2014).
* [12] Thomas Spura, Hossam Elgabarty, and Thomas D. Kuhne, ""On-the-fly" coupled cluster path-integral molecular dynamics: impact of nuclear quantum effects on the protonated water dimer," Physical Chemistry Chemical Physics **17**, 14355-14359 (2015).
* [13] V. Kapil, J. VandeVondele, and M. Ceriotti, "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," The Journal of Chemical Physics **144**, 054111 (2016), arXiv:1512.00176.
* [14] Jorg Behler, "Neural network potential-energy surfaces in chemistry: a tool for large-scale simulations," Physical Chemistry Chemical Physics **13**, 17930-17955 (2011).
* [15] Pascal Friederich, Florian Hase, Jonny Proppe, and Alan Aspuru-Guzik, "Machine-learned potentials for next-generation matter simulations," (2021).
* [16] Jorg Behler and Michele Parrinello, "Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces," Physical Review Letters **98**, 146401 (2007).
* [17] Albert P. Bartok, Mike C. Payne, Risi Kondor, and Gabor Csanyi, "Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons," Physical Review Letters **104** (2010), 10.1103/PhysRevLett.104.136403, arXiv:0910.1019.
* [18] J. S. Smith, O. Isayev, and A. E. Roitberg, "ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost," Chemical Science **8**, 3192-3203 (2017), arXiv:1610.08935.
* [19] K. T. Schutt, P. J. Kindermans, H. E. Sauceda, S. Chmiela, A. Tkatchenko, and K. R. Muller, "SchNet: A continuous-filter convolutional neural network for modeling quantum interactions," Advances in Neural Information Processing Systems **2017-Decem**, 992-1002 (2017), arXiv:1706.08566.
* [20] Simon Batzner, Albert Musaelian, Lixin Sun, Mario Geiger, Jonathan P. Mailoa, Mordechai Kornbluth, Nicola Molinari, Tse E. Smid, and Boris Kozinsky, "E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials," Nature Communications **13** (2022), 10.1038/s41467-022-29939-5, arXiv:2101.03164.
* [21] Pavlo O. Dral, "Quantum Chemistry in the Age of Machine Learning," (2020).
* [22] Jorg Behler, "Atom-centered symmetry functions for constructing high-dimensional neural network potentials," Journal of Chemical Physics **134** (2011), 10.1063/1.3553717.
* [23] Suresh Kondati Natarajan and Jorg Behler, "Neural network molecular dynamics simulations of solid-liquid interfaces: Water at low-index copper surfaces," Physical Chemistry Chemical Physics **18**, 28704-28752 (2016).
* A deep learning architecture for molecules and materials," Journal of Chemical Physics **148** (2018), 10.1063/1.5019779, arXiv:1712.06113.
* [25] Ganesh Sivaram, Anand Narayanan Krishnamoorthy, Matthias Baur, Christian Holm, Marius Stan, Gabor Csanyi, Chris Benmore, and Alvaro Vazquez-Mayagoitia, "Machine-learned interatomic potentials by active learning: amorphous and
liquid hafnium dioxide," npj Computational Materials **6** (2020), 10.1038/s41524-020-00367-7.
* Schran _et al._ [2021]Christoph Schran, Fabian L. Thiemann, Patrick Rowe, Erich A. Muller, Ondrej Marsalek, and Angelos Michaelides, "Machine learning potentials for complex aqueous systems made simple," Proceedings of the National Academy of Sciences of the United States of America **118**, e211007118 (2021), arXiv:2106.00048.
* Schran _et al._ [2020]Christoph Schran, Krystof Brezina, and Ondrej Marsalek, "Committee neural network potentials control generalization errors and enable active learning," The Journal of Chemical Physics **153**, 104105 (2020), arXiv:2006.01541.
* Hansen and Salamon [1990]Lars Kai Hansen and Peter Salamon, "Neural Network Ensembles," IEEE Transactions on Pattern Analysis and Machine Intelligence **12**, 993-1001 (1990).
* Gastegger _et al._ [2017]Michael Gastegger, Jorg Behler, and Philipp Marquetand, "Machine learning molecular dynamics for the simulation of infrared spectra," Chemical Science **8**, 6924-6935 (2017).
* Schran _et al._ [2020]Christoph Schran, Jorg Behler, and Dominik Marx, "Automated Fitting of Neural Network Potentials at Coupled Cluster Accuracy: Protonated Water Clusters as Testing Ground," Journal of Chemical Theory and Computation **16**, 88-99 (2020), arXiv:1908.08734.
* Krogh and Vedelsby [1995]Anders Krogh and Jesper Vedelsby, "Neural Network Ensembles, Cross Validation, and Active Learning," Advances in Neural Information Processing Systems **7**, 231-238 (1995).
* COLT '92_ (ACM Press, New York, New York, USA, 1992) pp. 287-294.
* Chen _et al._ [2020]Lei Chen, Ivan Sukuba, Michael Probst, and Alexander Kaiser, "Iterative training set refinement enables reactive molecular dynamics via machine learned forces," RSC Advances **10**, 4293-4299 (2020).
* Morawietz _et al._ [2016]Tobias Morawietz, Andreas Singraber, Christoph Dellago, and Jorg Behler, "How van der Waals interactions determine the unique properties of water," Proceedings of the National Academy of Sciences **113**, 8368-8373 (2016).
* Smith _et al._ [2018]Justin S. Smith, Ben Nebgen, Nicholas Lubbers, Olexandr Isayev, and Adrian E. Roitberg, "Less is more: Sampling chemical space with active learning," Journal of Chemical Physics **148** (2018), 10.1063/1.5023802, arXiv:1801.09319.
* Musl _et al._ [2018]Felix Musl, Sandip De, Jack Yang, Joshua E. Campbell, Graeme M. Day, and Michele Ceriotti, "Machine learning for the structure-energy-property landscapes of molecular crystals," Chemical Science **9**, 1289-1300 (2018).
* Vandermause _et al._ [2020]Jonathan Vandermause, Steven B. Torrisi, Simon Batzner, Yu Xie, Lixin Sun, Alexie M. Kolpak, and Boris Kozinsky, "On-the-fly active learning of interpretable Bayesian force fields for atomistic rare events," npj Computational Materials **6** (2020), 10.1038/s41524-020-0283-z, arXiv:1904.02042.
* Wang _et al._ [2020]Wujie Wang, Tzhusing Yang, William H. Harris, and Rafael Gomez-Bombarelli, "Active learning and neural network potentials accelerate molecular screening of ether-based soluic ionic liquids," Chemical Communications **56**, 8920-8923 (2020).
* Miksch _et al._ [2021]April M. Miksch, Tobias Morawietz, Johannes Kastner, Alexandre Urban, and Nongnuch Artiti, "Strategies for the construction of machine-learning potentials for accurate and efficient atomic-scale simulations," (2021), arXiv:2101.10468.
* Schwabbe-Koda _et al._ [2021]Daniel Schwabbe-Koda, Aik Rui Tan, and Rafael Gomez-Bombarelli, "Differentiable sampling of molecular geometries with uncertainty-based adversarial attacks," Nature Communications **12**, 1-12 (2021), arXiv:2101.11588.
* Vandenhaute _et al._ [2023]Sander Vandenhaute, Maarten Cools-Ceuppens, Simon DeKeyser, Toon Verstraelen, and Veronique Van Speybroeck, "Machine learning potentials for metal-organic frameworks using an incremental learning approach," npj Computational Materials **9** (2023), 10.1038/s41524-023-00969-x.
* Rupp _et al._ [2015]Matthias Rupp, Raghunathan Ramakrishnan, and O. Anatole Von Lilienfeld, "Machine Learning for Quantum Mechanical Properties of Atoms in Molecules," Journal of Physical Chemistry Letters **6**, 3309-3313 (2015), arXiv:1505.0035.
* Ceriotti _et al._ [2009]Michele Ceriotti, Giovanni Bussi, and Michele Parrinello, "Nuclear Quantum Effects in Solids Using a Colored-Noise Thermostat," Physical Review Letters **103**, 030603 (2009).
* Jonsson _et al._ [1998]Hannes Jonsson, Greg Mills, and Karsten W. Jacobsen, "Nudged elastic band method for finding minimum energy paths of transitions,", 385-404 (1998).
* Hutter _et al._ [2014]Jurg Hutter, Marcella Iannuzzi, Florian Schiffmann, and Joost Vandevondele, "Cp2k: Atomistic simulations of condensed matter systems," Wiley Interdisciplinary Reviews: Computational Molecular Science **4**, 15-25 (2014).
* Vandevondele _et al._ [2005]Joost Vandevondele, Matthias Krack, Fawzi Mohamed, Michele Parrinello, Thomas Chassaing, and Jurg Hutter, "Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach," Computer Physics Communications **167**, 103-128 (2005).
* Quickstep: Efficient and accurate electronic structure calculations," The Journal of Chemical Physics **152**, 194103 (2020), arXiv:2003.03868.
* Porezag [1995]D. Porezag, Th Frauenheim, Th Kohler, G. Seifert, and R. Kaschner, "Construction of tight-binding-like potentials on the basis of density-functional theory: Application to carbon," Physical Review B **51**, 12947 (1995).
* Perdew _et al._ [1996]John P. Perdew, Kieron Burke, and Matthias Ernzerhof, "Generalized gradient approximation made simple," Physical Review Letters **77**, 3865-3868 (1996).
* Zhang and Yang [1998]Yingkai Zhang and Weitao Yang, "Comment on "generalized gradient approximation made simple"," Physical Review Letters **80**, 890 (1998).
* Adamo and Barone [1999]Carlo Adamo and Vincenzo Barone, "Toward reliable density functional methods without adjustable parameters: The PBE0 model," Journal of Chemical Physics **110**, 6158-6170 (1999).
* Goerigk and Grimme [2011]Lars Goerigk and Stefan Grimme, "A thorough benchmark of density functional methods for general main group thermochemistry, kinetics, and noncovalent interactions," Physical Chemistry Chemical Physics **13**, 6670-6688 (2011).
* Guidon _et al._ [2008]Manuel Guidon, Florian Schiffmann, Jurg Hutter, and Joost Vandevondele, "Ab initio molecular dynamics using hybrid density functionals," Journal of Chemical Physics **128**, 1-15 (2008).
* Guidon _et al._ [2009]Manuel Guidon, Jurg Hutter, and Joost Vandevondele, "Robust periodic Hartree-Fock exchange for large-scale simulations using Gaussian basis sets," Journal of Chemical Theory and Computation **5**, 3010-3021 (2009).
* Condensed Matter and Materials Physics **54**, 1703-1710 (1996), arXiv:9512004 [mtrl-th].
* Guidon _et al._ [2010]Manuel Guidon, Jurg Hutter, and Joost Vandevondele, "Auxiliary density matrix methods for Hartree-Fock exchange calculations," Journal of Chemical Theory and Computation **6**, 2348-2364 (2010).
* Singraber _et al._ [2019]Andresa Singraber, Tobias Morawietz, Jorg Behler, and Christoph Dellago, "Parallel Multistream Training of High-Dimensional Neural Network Potentials," Journal of Chemical Theory and Computation **15**, 3075-3092 (2019).
* Shah _et al._ [1992]Samir Shah, Francesco Palmieri, and Michael Datum, "Optimal filtering algorithms for fast learning in feedforward neural networks," Neural Networks **5**, 779-787 (1992).
* Blank and Brown [2015]Thomas B. Blank and Steven D. Brown, "Adaptive, global, ex
tended Kalman filters for training feedforward neural networks," Journal of Chemometrics **8**, 391-407 (1994).
* [60] "AML public GitHub repository. [https://github.com/marselekgroup/aml](https://github.com/marselekgroup/aml), Accessed: 2023-03-20,".
* [61] C. G. Broyden, "The Convergence of a Class of Double-rank Minimization Algorithms 1. General Considerations," IMA Journal of Applied Mathematics **6**, 76-90 (1970).
* A Python library for working with atoms," Journal of Physics Condensed Matter **29**, 273002 (2017).
* [63] Erik Bitzek, Pekka Koskinen, Franz Gahler, Michael Moseler, and Peter Gumbsch, "Structural relaxation made simple," Physical Review Letters **97**, 170201 (2006).
* [64] Graeme Henkelman, Blas P. Uberuaga, and Hannes Jonsson, "Climbing image nudged elastic band method for finding saddle points and minimum energy paths," Journal of Chemical Physics **113**, 9901-9904 (2000).
* [65] Michele Ceriotti, Michele Parrinello, Thomas E Markland, and David E Manolopoulos, "Efficient stochastic thermostatting of path integral molecular dynamics." The Journal of chemical physics **133**, 124104 (2010).
* [66] Massimiliano Bonomi, Davide Branduardi, Giovanni Bussi, Carlo Camilloni, Davide Provasi, Paolo Raiteri, Davide Donadio, Fabrizio Marinelli, Pabio Pietrucci, Ricardo A. Broglia, and Michele Parrinello, "PLUMED: A portable plugin for free-energy calculations with molecular dynamics," Computer Physics Communications **180**, 1961-1972 (2009), arXiv:0902.0874.
* [67] Massimiliano Bonomi, Giovanni Bussi, Carlo Camilloni, Gareth A. Tribello, Pavel Banas, Alessandro Barducci, Mattia Bernetti, Peter G. Bolhuis, Sando Bottaro, Davide Branduardi, Riccardo Capelli, Paolo Carloni, Michele Ceriotti, Andrea Cesari, Haochuan Chen, Wei Chen, Francesco Colizzi, Sandip De, Marco De La Pierre, Davide Donadio, Viktor Drobot, Bernd Ensing, Andrew L. Ferguson, Marta Filizola, James S. Fraser, Hao-hao Fu, Piero Gasparotto, Francesco Luigi Gervasio, Federico Giberti, Alejandro Gil-Ley, Toni Giorgino, Gabriella T. Heller, Glen M. Hocky, Marcella Iannuzzi, Michele Invernizzi, Kim E. Jelfs, Alexander Jussupp, Espery Kirilin, Alessandro Laio, Vittorio Limongelli, Kresten Lindorff-Larsen, Thomas Lohr, Fabrizio Marinelli, Layla Martin-Samos, Matteo Masetti, Ralf Meyer, Angelos Michaelides, Carla Molteni, Tetsuya Morishita, Marco Nava, Cristina Passoni, Elena Papaleo, Michele Parrinello, Jim Pfaendtner, Pablo Piaggi, Giovanni Matic Piccini, Adriana Pietropaolo, Fabio Pietrucci, Silvio Pipolo, Davide Provasi, David Quigley, Paolo Raiteri, Stefano Ranioio, Jakabv Rydzewski, Matteo Salvalaglio, Gabriele Cesare Sosso, Vojtech Spiwok, Jiri Spomer, David W.B. Swenson, Pratyush Tiwary, Omar Valsson, Michele Vendruscolo, Gregory A. Voth, and Andrew White, "Promoting transparency and reproducibility in enhanced molecular simulations," Nature Methods **16**, 670-673 (2019).
* [68] Gareth A. Tribello, Massimiliano Bonomi, Davide Branduardi, Carlo Camilloni, and Giovanni Bussi, "PLUMED 2: New feathers for an old bird," Computer Physics Communications **185**, 604-613 (2014), arXiv:1310.0980.
* [69] Giovanni Bussi, Davide Donadio, and Michele Parrinello, "Canonical sampling through velocity rescaling," Journal of Chemical Physics **126**, 014101 (2007), arXiv:0803.4060.
* [70] Michael R. Shirts and John D. Chodera, "Statistically optimal analysis of samples from multiple equilibrium states," Journal of Chemical Physics **129**, 124105 (2008), arXiv:0801.1426.
* [71] "WHAM, [https://github.com/apallath/wham](https://github.com/apallath/wham), Accessed: 2023-03-20,".
* [72] Mark E Tuckerman and Dominik Marx, "Heavy-Atom Skeleton Quantization and Proton Tunneling in "Intermediate-Barrier" Hydrogen Bonds," Physical Review Letters **86**, 4946-4949 (2001).
* [73] Ales Cahlfk, Jack Hellerstedt, Jesus I. Mendieta-Moreno, Martin Svec, Vijai M. Santhini, Simon Pascal, Diego Soler-Polo, Sigurdur I. Erlingsson, Karel Vyborny, Ping Mutombo, Ondrej Marsalek, Olivier Siri, and Pavel Jelinek, "Significance of Nuclear Quantum Effects in Hydrogen Bonded Molecular Chains," ACS Nano **15**, 10357-10365 (2021), arXiv:2007.14657.
* [74] Zorka Smedarchina, Willem Siebrand, Antonio Fernandez-Ramos, and Ruben Meana-Paneda, "Mechanisms of double proton transfer. Theory and applications," Zeitschrift fur Physikalische Chemie **222**, 1291-1309 (2008).
* [75] Johannes Kastner, "Theory and simulation of atom tunneling in chemical reactions," WIREs Computational Molecular Science **4**, 158-168 (2014).
**Supporting information for: Reducing the cost of neural network potential generation for reactive molecular systems**

Krystof Brezina, Hubert Beck, and Ondrej Marsalek

[0.5cm] Charles University, Faculty of Mathematics and Physics, Ke Karlovu 3, 121 16 Prague 2, Czech Republic

[0.5cm]

(Dated: 29 March 2023)

## S1 Additional theory details

The following section contains additional details on the theory discussed in the main text. First, the derivation of the effective quantum temperature of the quantum harmonic oscillator is presented. This is followed by the discussion of the density-matching condition needed in the TTS procedure for appending thermal NMS of the minima away from MEPs.

### Quantum effective temperature

To derive the expression for the quantum effective temperature, we start from the standard Hamiltonian for a 1D harmonic oscillator in a mass-weighted coordinate \(\Omega\) with the natural frequency \(\omega\)

\[\widehat{H}=-\frac{\hbar^{2}}{2}\frac{\mathrm{d}^{2}}{\mathrm{d}\Omega^{2}}+ \frac{1}{2}\omega^{2}\Omega^{2}.\] (S1)

The \(n\)-th solution for the corresponding time-independent Schrodinger equation (\(n=0,1,2,\ldots,\infty\)) is in the well-known form based on the physicist's Hermite polynomials \(H_{n}\)

\[\psi_{n}(\Omega)=\left(\frac{1}{2^{n}n!\sqrt{\pi}}\right)^{\frac{1}{2}}\left( \frac{\omega}{\hbar}\right)^{\frac{1}{4}}e^{-\frac{\omega\Omega^{2}}{2\hbar}} H_{n}\left(\sqrt{\frac{\omega}{\hbar}}\Omega\right),\] (S2)

which can be simplified to

\[\psi_{n}(y)=\left(\frac{\alpha}{2^{n}n!\sqrt{\pi}}\right)^{\frac{1}{2}}e^{- \frac{y^{2}}{2}}H_{n}(y)\] (S3)

with the substitutions \(\alpha=\sqrt{\frac{\omega}{\hbar}}\) and \(y=\alpha\Omega\). The corresponding energy levels are

\[E_{n}=\hbar\omega\left(n+\frac{1}{2}\right).\] (S4)

Ultimately, we are interested in the thermal density of the harmonic oscillator \(\rho(y)\), which is formally obtained as the diagonal elements of the full density matrix \(\rho(y,y^{\prime})\). A formula for the density matrix in terms of an infinite sum over all states can be obtained as

\[\begin{split}\rho(y,y^{\prime})&=\bra{y^{\prime}} \ket{e^{-\beta\widehat{H}}}{y}\\ &=\sum_{n,m}\bra{y^{\prime}}\psi_{n}\bra{\psi_{n}}e^{-\beta \widehat{H}}\ket{\psi_{m}}\bra{\psi_{m}}y\\ &=\sum_{n}\psi_{n}^{*}(y^{\prime})\psi_{n}(y)e^{-\beta E_{n}}\\ &=e^{-\frac{\beta\hbar\omega}{2}}\frac{\alpha}{\sqrt{\pi}}\sum_ {n}\frac{1}{2^{n}n!}e^{-\beta\hbar\omega n}.\\ &\cdot e^{-\frac{y^{2}+y^{\prime 2}}{2}}H_{n}(y)H_{n}(y^{\prime}), \end{split}\] (S5)

which can be closed using the so-called Mehler kernel [11]

\[\sum_{n}\frac{\chi^{n}}{2^{n}n!}H_{n}(x)H_{n}(y)=\frac{1}{\sqrt{1-\chi^{2}}}e ^{-\frac{\chi^{2}(x^{2}+y^{2})-2\chi x}{1-\chi^{2}}}:\] (S6)

a mathematical identity for a parameter \(\chi\) that allows summing over the above product of Hermite polynomials. With this, one obtains

\[\begin{split}\rho(y,y^{\prime})&=e^{-\frac{\beta \hbar\omega}{2}}\frac{\alpha}{\sqrt{\pi}}e^{-\frac{y^{2}+y^{\prime 2}}{2}}\frac{1}{ \sqrt{1-e^{-2\beta\hbar\omega}}}\\ &\cdot e^{-\frac{-2\beta\hbar\omega(y^{2}+y^{\prime 2})-2\epsilon- \beta\hbar\omega\,yy^{\prime}}{1-e^{-2\beta\hbar\omega}}}.\end{split}\] (S7)

The thermal density is then simply obtained by setting \(y=y^{\prime}\), which gives

\[\begin{split}\rho(y)&=e^{-\frac{\beta\hbar\omega}{2 }}\frac{\alpha}{\sqrt{\pi}}\frac{1}{\sqrt{1-e^{-2\beta\hbar\omega}}}\\ &\cdot e^{-y^{2}\left[1+\frac{2e^{-\beta\hbar\omega}\left(e^{- \beta\hbar\omega}-i\right)}{1-e^{-2\beta\hbar\omega}}\right]}\end{split}\] (S8)

after trivial rearrangements. As discussed in the main text, the thermal density is Gaussian with an \(\omega\)-dependent normalization factor (first line of Equation S8) and a non-trivial scaling of the exponent. The latter will be rearranged to extract the quantum effective temperature from Equation 3 of the main text. Substituting \(t=e^{-\beta\hbar\omega}\), we get

\[\begin{split}& 1+\frac{2e^{-\beta\hbar\omega}\left(e^{-\beta \hbar\omega}-1\right)}{1-e^{-2\beta\hbar\omega}}=1+\frac{2t(t-1)}{1-t^{2}}\\ &=\frac{t^{2}-2t+1}{1-t^{2}}=\frac{1-t}{1+t}=\frac{1-e^{-\beta \hbar\omega}}{1+e^{-\beta\hbar\omega}}\\ &=\tanh\left(\frac{\beta\hbar\omega}{2}\right),\end{split}\] (S9)
which gives

\[\rho(\Omega)\propto e^{-\frac{\omega}{\hbar}\tanh\left(\frac{\beta\hbar\omega}{2} \right)\Omega^{2}}\] (S10)

for the thermal density in the original coordinates. Comparing the hyperbolic tangent factor with the form of a classical distribution (\(\rho_{\mathrm{cl}}(y)\propto e^{-\beta V}=e^{-\frac{1}{2}\beta\omega^{2}\Omega^{ 2}}\)) of the harmonic oscillator at a new inverse temperature \(\beta^{*}\), we immediately get

\[\begin{split}&\frac{\omega}{\hbar}\tanh\left(\frac{\beta\hbar \omega}{2}\right)\equiv\frac{1}{2}\beta^{*}\omega^{2}\\ &\beta^{*}(\beta,\omega)=\frac{2}{\hbar\omega}\tanh\left(\frac{ \beta\hbar\omega}{2}\right),\end{split}\] (S11)

as stated in the main text. Since the hyperbolic tangent approaches unity as its argument grows, it is always lower-valued than a linear function and thus always \(\beta^{*}<\beta\). In turn, this means that the quantum harmonic oscillator is always at an effective temperature higher than the reference classical temperature and its density is thus always wider. Moreover, note that the classical limit \(\beta^{*}\rightarrow\beta\) is reached as \(\beta\to 0\) or \(\omega\to 0\): at a given temperature, the higher frequency modes will have a more pronounced quantum effect and, in turn, at a given frequency, the quantum effect will be more pronounced at a lower temperature.

### Density matching

An obvious technical problem with appending the thermal NMS at the minima at the MEP endpoints is the matching of the sampling density so that it is a smooth continuation of the linear MEP sampling (as indicated in Figure 1 of the main text). For practical purposes, one can achieve sufficient smoothness by separately choosing a linear density for the MEP sampling and a fixed number of points that are sampled thermally around the minima by a trial-and-error approach. However, this does not lead to an analytically smooth transition. That can be achieved by integrating the thermal Gaussian distributions at the minima over all dimensions but the one defined by the tangent vector to the MEP \(\mathbf{\tau}\) at the minimum and finding a scaling factor that brings the density at the maximum of the Gaussian equal to the linear density \(\rho_{0}\) of sampling along the MEP, yielding an analytic continuation. Rather inconveniently from a mathematical perspective, the MEP tangent \(\mathbf{\tau}\) is, in general, not identifiable with any of the normal modes \(\mathbf{\Omega}_{i}\), which diagonalize the Gaussian distribution. Therefore, to perform the integration over everything but the \(\tau\)-direction, one needs to set up a new orthogonal basis built around this vector in which, however, the Gaussian is no longer diagonal. Still, the result of such integration will be a 1D Gaussian and its width will be given by the effective frequency \(\widetilde{\omega}\)

\[\rho(\tau)=\int\rho(\tau,\tau_{2},\ldots,\tau_{3N})\ \mathrm{d}\tau_{2} \ldots\mathrm{d}\tau_{3N}\propto e^{-\frac{1}{2}\beta\widetilde{\omega}^{2} \tau^{2}}.\] (S12)

It can be shown (with the full proof being outside of the scope of the present work), that this effective frequency can be calculated from the following compact formula

\[\widetilde{\omega}^{2}=\frac{\det(\mathbb{H})}{\det_{\tau,\tau}(\mathbb{H})},\] (S13)

where \(\mathbb{H}\) is the Hessian matrix of the potential at the minimum and the symbol \(\det_{\tau,\tau}\) denotes the minor of the Hessian matrix taken over the row and the column corresponding to the \(\tau\)-dimension. Note that, unlike the full determinant, the minor is a basis-dependent operation and, therefore, the Hessian matrix must be transformed to the \(\tau\)-oriented basis prior to the calculation of \(\widetilde{\omega}\). The effective standard deviation of the 1D distribution is \(\widetilde{\sigma}=1/\sqrt{\beta\widetilde{\omega}^{2}}\). The aim at this point is to find a scaling factor \(\alpha\) which ensures that the density at the configurational minimum \(\alpha\rho(0)\) matches \(\rho_{0}\). However, one cannot simply equate the two values directly, since the \(\rho(\tau)\) comes from a mass-weighted coordinate system (note that the \(\mathbf{\tau}\)-based basis was obtained as a simple rotation of the normal mode basis, which itself is naturally mass-weighted) and \(\rho_{0}\) is understood as a density along the non-weighted Cartesian MEP curve. Therefore, we need to calculate the effective mass \(\mu\) of the \(\tau\)-direction, which is the following weighted average

\[\mu=\frac{\mathbf{\tau}^{T}\mathbb{M}\mathbf{\tau}}{\mathbf{\tau}^{T}\mathbf{\tau}},\] (S14)

where \(\mathbb{M}\) is the diagonal mass matrix. Then, we obtain the analytic expression

\[\alpha\rho(0)=\frac{\alpha}{\widetilde{\sigma}\sqrt{2\pi}}=\sqrt{\frac{1}{\mu} }\rho_{0}\] (S15)

which immediately gives the crucial density-matching condition

\[\alpha=\widetilde{\sigma}\sqrt{\frac{2\pi}{\mu}}\rho_{0}.\] (S16)

The scaling parameter \(\alpha\) thus captures how many times the original normalized density \(\rho(\tau)\) needs to be increased in order to match at its peak the sampling density at the MEP. In practice, where the sampling of configurations is discrete and finite, one needs to sample the nearest-integer-to-\(\alpha\) points to match the sampling density at the MEP.

## S2 Additional results

### TTS of malonaldehyde proton transfer

The top panel of Figure S1 shows the distribution of the TTS geometries of malonaldehyde obtained at 300 K
using the classical harmonic distributions in the mode directions perpendicular to the proton-sharing MEP. Three control points were selected for the TTS: the two symmetric minima at the edges of the MEP and the transition state at \(\delta=0\,\)A. Each candidate geometry is colored by the assignment to its control point following the methodology described in Section II of the main text and illustrated in Figure 1 therein. The bottom panel of Figure S1 then shows the same situation but using the quantum formulation of TTS in which \(\beta^{*}\) is calculated for each degree of freedom. The use of the quantum densities leads to the characteristic broadening of the distribution of the geometries.

### Path integral MD for malonaldehyde

The process described in the main paper for classical MD simulation can be done analogously for PIMD. TTS was used to generate structures along the proton-sharing MEP at an effective temperature of 300 K and 620 structures were selected during one QbC run for the initial model. Similar to the classical model, force disagreement during PIMD was low around the two free energy minima and along the MEP, but increased at the tails of the distribution, as can be seen in the left panel of Figure S2 To improve the MLP in these areas of the configuration space, two strategies were pursued. An iterative process, in which PIMD simulations and QbC processes with the structures of the PIMD trajectory as candidates are alternated until force disagreement is low along the whole PIMD trajectory. Or a single QbC run, in which TTS is used to create candidate structures along the MEP of the proton-transfer reaction as well as the MEP of the torsion of the C-C single bond of the propane back bone up until the transition state to include the strongly anharmonic vibrations in this direction of the configuration space. The distribution of the force disagreement for the second-generation MLP of the iterative process can be seen in the middle panel of Figure S2. Unlike the classical case in the main paper, a single additional PIMD-QbC sequence, which added 435 structures to the training data set, was sufficient to obtain an MLP capable of accurate PIMD simulations. In the right panel of Figure S2, the results from the extended TTS MLP can be seen. As for the classical case, the resulting model gives a low committee disagreement along the whole PIMD trajectory.

### Validation of malonaldehyde MLPs

In the main text of the paper, the energy and force RMSEs on an independent test set for the final models were reported. The complete set of validation errors can be seen in table S1. While the later generations of the iterative approach and the extended TTS model performed well, the initial model performed poorly. However, the errors averaged over the whole test set tell only a part of the story, because they are not evenly distributed for all MLPs. As shown in Figure S3, the first generation model displays low force errors for most structures with a low \(d_{\mathrm{O-O}}\), which accounts for structures along the proton-sharing transition and the minima, but exceedingly large errors of up to 603 meV (energy) and 3908 \(\mathrm{meV\AA}^{-1}\) (forces) for \(d_{\mathrm{O-O}}>2.8\). Due to the anharmonicity of the system, no structures from this region of the configuration space are included in the TTS along the proton sharing coordinate and hence such structures are absent from the training data set. Therefore, the model is not suitable for accurate MD simulations. The same effect can be observed in the second generation model,
however less severe. The final models of both approaches, on the other hand, exhibit an even distribution of force RMSEs across the complete span of oxygen-oxygen distances.

### Overlap of umbrella windows in DABQDI

Figure S4 shows that the histograms of \(\delta_{1}\) values from the individual windows of the umbrella sampling simulation of the proton transfer in the DABQDI molecule exhibit sufficient overlap. Therefore, the combination of window spacing and harmonic restraint stiffness selected to obtain the results shown in Figure 7 of the main text represents an appropriate choice.
