Metadynamics simulations with Bohmian-style bias potential

Here, we present a parametrization of the metadynamics simulations for reactions involving breaking the chemical bonds along a single collective variable coordinate. The parameterization is based on the similarity between the bias potential in metadynamics and the quantum potential in the de Broglie – Bohm formalism. We derive the method and test it on two prototypical reaction types: proton transfer and breaking of the cyclohexene cycle (reversed Diels – Alder reaction).


| INTRODUCTION
Metadynamics (MTD) is a powerful molecular-dynamics-based simulation technique that is known to provide the chemical reaction pathways without prior knowledge of the transition states the system goes through. [1][2][3][4] MTD is an example of so-called enhanced sampling techniques, 5,6 among which there are such classical schemes as local elevation 7 or conformational flooding. 8 To include the quantum nature of the nuclear motions, 9 the MTD can be combined with the path-integral molecular dynamics. 10,11 In recent years, MTD was also successfully extended to automatic conformational search. 12 The idea of the MTD is the propagation of molecular dynamics with an effective potential given as [1][2][3][4] where V is the potential energy surface (PES) of the nuclei, and the V MTD is the so-called bias potential. The latter marks the previous configurations of the system by increasing their energies, making it less probable for the system to go back, thus forcing the molecule to explore new regions of the PES.
The bias potential is usually formulated using collective variables (CV). These are just functions of the internal coordinates, parametrizing the expected reaction coordinate of the system. [2][3][4] The classical numerical representation of the bias potential is a sum of Gaussian functions of the user-defined widths and heights. These basic components are added with some periodicity during the MTD simulation, and the centers of these Gaussians are placed at the current values of the CVs. Therefore, at the end of the simulation, ÀV MTD represents a coarse and averaged snapshot of the sampled part of the PES. [1][2][3][4] The interpretation of this bias potential can also be made through the probability distributions. 13 In this article, we propose an alternative view of the MTD bias potential. We use its similarity to the quantum potential from the de Broglie-Bohm formalism of quantum mechanics. 14 We discuss the application of this reformulated MTD simulation to a specific class of reactions that involve the breaking of chemical bonds. In particular, we demonstrate the applicability of the proposed approach using the two prototypic reactions: proton transfer in the malonaldehyde 15,16 and the inverse Diels-Alder reaction [17][18][19] of the cyclohexene. In the end, we provide conclusions and perspectives on the current approach.

| Bias potential shape
To find the building block for the bias potential, let us consider the onedimensional quantum potential from the De Broglie-Bohm formalism 14,20 : Abbreviations: MTD, metadynamics; CV, collective variable; PES, potential energy surface.
where q is the CV, μ is the effective mass of the molecular motion along the CV, and A is the wavefunction's (ψ) amplitude (A ¼ jψj). Let us assume that at the current value of the CV, q ¼ Q, we have a Gaussian wavepacket centered at the position q ¼ Q 20 : where σ is the standard deviation (SD) of the CV q from the mean position Q, and N ¼ ffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffiffiffi ffi 2π p σ p À1 is the normalization factor. Substitution of the Equation (2) into the Equation (3) gives us the quantum potential of the shape We want our bias potential for the MTD to be non-negative. Therefore, we can take only the part of the Equation (3) corresponding to V quant ðqÞ ≥ 0 (ðq À QÞ 2 ≤ 2σ 2 ), and set the rest of the potential to zero, thus having the following building blocks for the bias potential: where notation ðqjQÞ denotes that the bias component's argument is q and this function is centered at the position q ¼ Q, and E b is the bias energy: The total bias potential at a given point of the MTD simulation will thus be where N is the total number of marked CV values Q n .

| Collective variable definition
To model the bond breaking in the molecules, the most native coordinate is the distance between the atoms: , where i and j are the atomic indexes, and r i , r j denote atomic position vectors. However, since there might be cases when more than two bonds are broken during a reaction, such as the double proton transfer or the breaking of the cycles, we might need more than one distance in the CV definition. Therefore, we will take the following definition of the CV: where index k enumerates the distances in the definition of CV, jc k j > 0 are the arbitrary coefficients, and R k ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðr ak À r bk Þ 2 q is the distance between atoms number a k and b k .
To define the bias energy (Equation 5), we require a reduced mass of the motion. In the case of the q ¼ R ij , it is defined simply as the reduced mass of the two atoms forming the bond, that is, as where m i and m j are the masses of the atoms number i and j, respectively. 21 We will consider a simplified approach based on the classical equipartition theorem 22,23 to define the reduced mass μ q for the CV coordinate given by Equation (7).
Let us assume that the following equalities hold: T is the absolute temperature, and _ x ¼ dx=dt is the velocity of the value x. With the definition of the CV (Equation 7), we can also write: Assuming that the motions along different bonds are uncorrelated which gives the expression for the μ q : This expression provides a unique reduced mass value for any CV of the type given in Equation (7). In the standard MTD simulation, the user defines the CV, the width of the Gaussian bias potential, and the height of this potential. However, in this formulation, based on the quantum potential, the user needs to provide a definition of the CV and the width of the bias potential σ to define the MTD simulation because the reduced mass μ q and the bias energy are computed automatically with Equations (8) and (5), respectively.

| Bias force
The last issue to the MTD simulations with the bias potential given in Equation (6) is the evaluation of the forces. The bias force acting on the atom i with coordinate r i is given with the equation 24 However, some terms of this sum can be zero because points Q n can be out of the action range from the current value q (those that fulfill the inequality ðq À Q n Þ 2 > 2σ 2 , see Equation 4). Thus, we can order the Q n by their proximity to q, that is, in the increasing order of jq À Q n j, and take only thoseÑ ≤ N values that correspond to nonzero contributions to the potential (i.e., ðq À Q n Þ 2 ≤ 2σ 2 for 1 ≤ n ≤Ñ, and ðq À Q n Þ 2 > 2σ 2 for n >Ñ). The derivatives of the components for 1 ≤ n ≤Ñ are easy to compute: which upon substitution to Equation (9) will give where F is the force magnitude and ∂q ∂r i is the direction vector of the force for atom number i. The latter can be computed in a straightforward fashion using Equation (7) as The vector ∂R k ∂r i ¼ 0 if atom i is not a part of the distance R k between atoms a k and b k (i.e., i ≠ a k and i ≠ b k ). In the other case, when, say, i ¼ a k and j ¼ b k , and thus R k ¼

| NUMERICAL DEMONSTRATION
We implemented the proposed MTD routine into the BBBMTD script of the PyRAMD software. 25,26 To demonstrate the algorithm's performance, we chose two prototypical reactions: proton transfer between two oxygens in the enol form of the malonaldehyde and decomposition of the cyclohexene into the 1,3-butadiene and ethylene, which is the inverse of the Diels-Alder reaction. The schematics of these reactions are given in the insets (A) of Figures 1 and 2, respectively.
To perform quantum-chemical calculations, including the computation of the gradients of the PES in the MTD, we used the ORCA 5 software. 27 The level of theory in all the calculations was BLYP-D3BJ/6-31G. [28][29][30][31] First, we optimized the molecular geometries of the malonaldehyde and cyclohexane with subsequent calculations of the harmonic frequencies at the chosen level of theory. From these equilibrium geometries, we started the MTD simulations, generating initial velocities of the atoms from Maxwell-Boltzmann distribution at the temperature of 300 K, 22 and then this temperature was kept by a Berendsen thermostat 32 with a relaxation time of 100 fs. The time step in every MTD simulation was 1 fs.
The MTD simulations for malonaldehyde had the following settings. The CV was where rðO À HÞ and rðOÁÁÁHÞ denote the distances between oxygens and hydrogen, corresponding to the OÀH covalent and OÁÁÁH hydrogen bonds in the equilibrium geometry, respectively (see Figure 1).
The reduced mass of this CV was thus μ q ¼ 1:9 a.m.u. The width of the bias potential σ was 0.2 Å, which subsequently gave the bias energy E b (Equation 5) of 111 cm À1 . The bias potential was updated with the current value of CV every 30 fs. We accumulated 16 MTD trajectories in total, each 2 ps in duration.
In the case of cyclohexene, the CV was where r 1 ðC À CÞ and r 2 ðC À CÞ the distances corresponding to the two CÀC bonds broken in the reaction (see the inset (A) of Figure 2). This definition gave the μ q of 3.0 a.m.u. Nine MTD trajectories with the bias potential update time of 10 fs and width σ ¼ 0:1 Å, resulting in bias energies E b ¼ 281 cm À1 , with total simulation duration of 5 ps. To prevent the dissociation of the ethylene and 1,3-butadiene, an artificial constraining potential was applied. It was defined as where k B is the Boltzmann's constant, T ¼ 300 K is the temperature of the MTD simulation, r i is the distance of the ith atom to the center of mass of the system, and fðrÞ is defined as  15,16 The estimated barrier height was found to be 590 AE 60 cm À1 . We also computed the electronic activation energy (ΔE ≠ ) at the chosen level of theory for the same reaction as a reference value for comparison. We used a nudged elastic band with the transition state optimization afterward (NEB-TS) to find the transition state molecular geometry. 35 The resulting value of ΔE ≠ ¼ 623 cm À1 is within the uncertainty of the estimated value from the MTD simulations.
In the case of cyclohexene, the estimated activation energy from the fine MTD simulations was 19,300 AE 500 cm À1 . For comparison, the ΔE ≠ of the inverse Diels-Alder reaction was also computed. The

| CONCLUSIONS
Here, we presented an alternative approach to constructing the meta- Open Access funding enabled and organized by Projekt DEAL.

CONFLICT OF INTEREST
The authors have no conflicts of interest to declare.

DATA AVAILABILITY STATEMENT
The data that supports the findings of this study are available in Data S1.