{'fraction_non_alphanumeric': 0.02521563589649477, 'fraction_numerical': 0.00880895577170123, 'mean_word_length': 4.809381663113006, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reactive molecular dynamics at constant pressure via non-reactive force fields: extending the Empirical Valence Bond method to the isothermal-isobaric ensemble

The Empirical Valence Bond (EVB) method offers a suitable framework to obtain reactive potentials through the coupling of non-reactive force fields. However, most of the implemented functional forms for the coupling terms depend on complex spatial coordinates, which precludes the computation of the stress tensor for condensed phase systems and prevents the possibility to carry out EVB molecular dynamics in the isothermal-isobaric (NPT) ensemble. In this work, we make use of coupling terms that depend on the energy gaps, defined as the energy differences between the participating non-reactive force fields, and derive an expression for the EVB stress tensor suitable for computations. Implementation of this new methodology is tested for a model of a single reactive malonaldehyde solvated in non-reactive water. Computed densities and classical probability distributions in the NPT ensemble reveals a negligible role of the reactive potential in the limit of low concentrated solutions, thus corroborating the validity of standard approximations customarily adopted for EVB simulations.


Introduction
Molecular dynamics (MD) simulations offer a powerful computational tool to derive atomistic insight of complex phenomena from organic chemistry and biochemistry to heterogeneous catalysis. [1][2][3] The interatomic interactions in classical MD simulations are based on force field (FF) descriptors 1 which allows for very fast computation of the interactions and access to simulate very large systems. Commonly, these FF descriptors have simple functional forms, with parameters either fitted to experimental data or derived from quantum mechanical calculations. 4,5 In most of the available FF libraries, functional forms and fitted parameters remain unchanged during the course of the MD simulation. In reactive processes, however, the interactions inevitably change due to the breaking and/or formation of chemical species.
Thus, standard FFs are not suitable to simulate chemical reactions and they are referred to as non-reactive.
An alternative to simulate chemical reactions with MD is given by Reactive FFs (RFFs), [6][7][8][9][10][11] that are designed to model interatomic interactions of multiple states representing different chemical species. The task of designing RFFs, however, is very challenging 2 and, despite the enormous progress over the last years, 7,19 a general parameterization is not yet available.
The Empirical Valence Bond (EVB) method [20][21][22][23][24][25] offers, instead, a simple general framework to model reactive processes through the coupling of multiple non-reactive FFs, where each FF corresponds to a different chemical state for the system. In this method, a suitable EVB matrix is built using the computed energies of the involved chemical states as well as appropriate coupling terms. Matrix diagonalization at each time step allows computation of reactive energy landscapes that account for the change in chemistry when sampling conformations between the participating, chemically different, states. 1 A force field is a mathematical construction to model the interactions between atoms without having to compute the electronic Schrödinger equation. This construction reduces dramatically the computational effort to obtain energies and forces (and stress tensors for extended systems), allowing the sampling of the phase space up to nano-seconds, depending on the system and the computational resources available. 2 Indeed, designing RFFs requires a high level of expertise to tackle a multi-dimensional problem, 11 where the modelled interactions are often expressed by complicated functional forms with many strongly coupled parameters that are optimized via the use of sophisticated tools. [12][13][14][15][16][17][18] In contrast to RFFs, the advantage of the EVB method lies in the large availability of standard non-reactive FFs libraries, which has offered an appealing strategy for computational implementation and development over the past four decades. [26][27][28][29] Moreover, despite the tedious initial task to calibrate the coupling terms against reference data, research has demonstrated that these couplings are invariant to the surrounding electrostatics, making it possible to simulate the same reactive unit in different environments. 30 This convenient feature of the EVB method has widely increased its recognition as a powerful tool within the computational chemistry community. 24 For condensed phase systems, reported MD simulations with the EVB method (herein, MD-EVB simulations) are conducted either in the microcanonical (NVE) or the canonical (NVT) ensemble. However, MD-EVB simulations at constant pressure and temperature, i.e. using the isothermal-isobaric (NPT) ensemble have not been addressed in the literature. In fact, the standard protocol for EVB simulations in condensed phase is to first consider only one of the possible chemical states of the system (preferably the state with lowest free energy) together with the surrounding, non-reactive, environment and carry out a standard NPT simulation (without EVB) at the target pressure and temperature. 31 The converged volume is then fixed and the MD-EVB simulation is performed using the NVT ensemble. This procedure appears to be a sensible strategy for very large, homogeneous soft-matter systems.
However, for smaller systems or higher concentration of solutes, the validity of this standard protocol to approximate real experimental conditions at constant pressure and temperature has never been corroborated to date.
In this work, we demonstrate that the use of the standard formulation to compute the stress tensor cannot be directly applied to derive the components of the EVB stress tensor. We argue that this limitation explains the absence of MD-EVB simulations in the NPT ensemble.
In contrast to using complex spatial variables to fit the coupling terms of the EVB matrix, we propose to make use of energy gaps, 32 defined as the energy difference between the nonreactive FFs. With this choice, we derive an expression for the EVB stress tensor suitable for computational implementation, not only offering a solution to an overlooked limitation of EVB but also extending the applicability of MD-EVB simulations to NPT ensembles for the first time. The computational implementation of this new formalism is tested using a model of a solvated reactive malonaldehyde molecule in water. MD-EVB simulations at 300 K and 1 atm are used to quantify the role of the reactive potential in the computed density and classical probability distributions of the energy gaps obtained from the sampling of the configurational space. Results also allow to evaluate the validity of the standard protocol for MD-EVB simulations, while the derived method offers an opportunity to explore new strategies for future implementation and development of the EVB method.
The structure of this paper is as follows. The fundamentals of the EVB method developed over the years are presented in a convenient notation in section 2. In section 3, we discuss the limitation of the standard formulation to calculate the EVB stress tensor, and propose a new alternative method. A brief overview of the computational implementation is given in section 4. Section 5 discusses general aspects of the coupling terms within the framework of the present paper. Details of the model and MD computational setting are provided in section 6, which is followed by section 7 with the results and discussion. Concluding remarks are finally addressed in section 8.

The EVB method
In this section we present the fundamentals of the EVB formalism in a convenient notation.
Let us assume an atomic system composed of N p particles with positions described by the set of vectors {R}. The non-reactive force field (FF) for the chemical state (m) is described by the configurational energy E   can be generally written as a sum of different terms as follows 33 are the interactions representing core-shell polarization, tethered particles, chemical bonds, valence angles, dihedrals, inversion angles, three-body, four-body, Tersoff, metallic, van der Waals and coulombic contributions, respectively. Following Eq. (1), the forces can be expressed using a similar decomposition. In the current notation, we shall use indexes m and k for the chemical states (and FFs), I and J for atoms and Greek letters for Cartesian coordinates. Indexes in parenthesis are used to emphasize the particular chemical state.
The purpose of the EVB method is to couple N F non-reactive force fields to obtain a reactive potential. These FFs are coupled through the HamiltonianĤ EVB with a matrix representation H EVB ∈ R N F ×N F that has the following components where each diagonal element corresponds to the configurational energy E whereas the off-diagonal terms C mk are the couplings between states m and k. For convenience in the notation, we shall omit hereinafter the dependence on the set of coordinates {R} for the particles. Even though there are different possible choices for the coupling terms, in the above definition we have set C mk to depend on −ǫ km , where ǫ mk is commonly referred to as energy gap and defines a possible reaction co-ordinate for the process. 22,24,32,34 Since the H EVB matrix is Hermitian by construction and the C mk terms are real, the condition of C mk = C km must be imposed to the off-diagonal elements. Diagonalization of H EVB leads to N F possible eigenvalues {λ 1 , ..., λ N F } with The EVB energy, E EVB , is defined as the lowest eigenvalue with the corresponding normalized EVB eigenvector and Since the eigenvector Ψ EVB is real and normalized we have from which we can interpret |Ψ (k) EVB 2 as the fraction of the chemical state (k) being part of the EVB state. The eigenvector Ψ EVB can also be represented as a column vector EVB is the element of the k-row. Thus, Eq. (6) is expressed as a matrix multiplication whereα corresponds to each of the orthonormal Cartesian vectors and From Eq. (2) the matrix components of the operator ∂Ĥ EVB ∂ R J α are given as follows where C ′ mk = dC mk (ǫ mk ) The above equations define the standard EVB force field (EVB-FF). Even though the EVB formalism was first developed to compute molecular systems, EVB is also applicable to extended systems, customarily modelled using the supercell approximation and periodic boundary conditions (PBCs). Nevertheless, MD-EVB simulations have only been conducted for the NVE and NVT ensembles, to the best of our knowledge, as there is no evidence of a previously reported method to compute the EVB stress tensor. In the next section, we discuss the intricacies related to computing the stress tensor using the standard formulation and propose a new method that allows extending the applicability of MD-EVB to NPT ensembles the first time.

The EVB stress tensor
The key requirement for a NPT simulation with the EVB method is to being able to compute the EVB stress tensor σ EVB . Similarly to the energy and forces, the configurational stress tensor for the force field m, σ c(m) , can be decomposed in a general expression equivalent to Eq. (1), where each contribution is computed separately using well-known functional forms. 33,36 For bonded interactions, for example, the αβ contribution to the stress tensor from particle J due to the bonded interactions with the surrounding particles, σ bond(m) J,αβ , is given by where R JI,α is the α component of the vector separation R JI = R J − R I between particles I and J, and f bond(m) IJ the bond force over particle J from its bonded interaction with particle I. In Eq. (13) the sum runs over all particles I interacting with particle J via bonds.
Analogously, we could in principle propose an expression for the αβ component of the EVB stress tensor resulting from the EVB bonded forces, f EVB IJ,β , as follows, In the present case of bonded interactions, the evaluation of Eq. (14) requires of each individual EVB-bonded force over particle J from interaction with particles I, given by f EVB
Nevertheless, the EVB force given in Eq. (10) represents the total force, F EVB J , resulting from the interaction of particle J with all the neighboring particles, which generally include other type of interactions apart from bonding interactions. As far as we can discern, each individual contribution to the force f EVB IJ cannot be computed from the EVB formalism presented in last section and, consequently, the evaluation of the stress tensor via Eq. (14) is not possible. The same reasoning applies to other type of interactions. This limitation precludes the computation of the stress tensor within the EVB formalism via standard formulae and, consequently, MD simulations using the NPT ensemble. Surprisingly, this inherent limitation of the EVB method has not been previously discussed in the literature, to the best of our knowledge.
To circumvent this problem, we propose to use the well-known relation between the configurational energy and the configurational stress tensor 37 where h is the set of lattice vectors of the supercell with volume V =det(h). Multiplying to the left by h νβ and summing over β we obtain the inverse relation to Eq.
which can be used to define the EVB stress tensor Similar to the definition of the EVB force, we evaluate ∂E EVB /∂h αγ using the Eq. (6) and the Hellman-Feynman theorem 35 The matrix components of the operator ∂Ĥ EVB ∂ h αβ follow from the definition of the EVB matrix (2) and the use of relation (15) ∂H mk Finally, the EVB stress tensor of Eq. (17) can be expressed as a matrix multiplication These expressions provide an alternative to compute the stress tensor σ EVB from the configurational stress tensors of each non-reactive FF, σ αγ . It is important to note that this new scheme to compute σ EVB can only be derived if one uses functional forms for C mk that depend on the energy differences ǫ mk , for which one can evaluate ∂E (m) c ∂h αβ − ∂E (m) c ∂h αβ and use relation (15) with the computed configurational stress tensor for each chemical state. In contrast, if the choice was to use coupling terms that do not depend on ǫ mk but other degrees of freedom such as spatial coordinates (see refs. 38-44), we cannot discern a clear logic to derive an expression for σ EVB , which might explain the fact there is no evidence of any previous reported method to compute the stress tensor using EVB.
So far we have presented an alternative to compute the stress tensor σ EVB αβ but have not discussed the total virial V EVB . Similarly to the stress tensor, the inability to compute individual contributions of the EVB force prevents the evaluation of the virial using the standard formulation, 33 and the usual decomposition of the virial depending of the type of interaction under consideration. Within the presented formalism, we compute the virial V EVB from σ EVB αβ as follows In contrast to the EVB energy, it is possible to decomposed the virial into different type of interactions, as we discuss in section S2 of the Supporting Information. The total stress tensor, σ T , is given by the following general expression where σ kin , σ RB and σ bc are the contributions to the stress tensor from the kinetic energy, rigid bodies (RB) and bond constraints (bc), respectively. The EVB method only accounts for the configurational interactions, as described. The kinetic stress tensor is computed as usual from the instantaneous velocities of the particles. 33 where F I B is the total force over particle I of rigid body B and d I B the vector distance from atom I B to the center of mass of the rigid body B. In the above expression, index B runs over all the rigid bodies. Each rigid body is composed of η B particles. Since, by definition, the topology of rigid bodies remain unaltered during the simulation, the use of RBs within in the present framework is meaningful only to model the environment interacting reactive EVB site. A common example is the use of rigidly constrained water molecules to model a solution.
Contributions to the stress tensor from bond constraints, σ bc αβ , are obtained using the SHAKE/RATTLE algorithm 45,46 during the course of the simulation. This algorithm is independent of the EVB formalism, and corrects for the dynamics of the constrained particles. Finally, frozen particles do not contributed to the stress tensor and are not considered in the formalism. It is important to note that the topology defined via the setting of RBs, frozen atoms and bond constraints must be the consistent for all the coupled FFs, as they impose well defined conditions for the dynamics. For example, if a group of atoms form a rigid body, they must remain a rigid body independently of chemical state under consideration.

Overview of the computational implementation
The EVB method described in section 2 and its extension for the computation of the stress tensor (section 3) were implemented within the DL_POLY_4 code. 47,48 In the standard format, DL_POLY_4 reads the initial coordinates, velocities and forces from the CONFIG file. Each particle is labelled according to its specification in the FIELD file, which contains the information of the FF type and parameters for the interactions between the particles.
Settings for the MD simulation are specified in the CONTROL file. Initially, the code was modified to allow i) reading multiple (N F ) CONFIG and FIELD files, ii) allocating arrays of dimension N F for the relevant quantities, iii) checking consistency of specification between all force fields and initial coordinates (including any possible constraint such as rigid bodies), iv)

Coupling terms
The quality of EVB method depends on the choice for the coupling terms C mk , particularly to The implementation of such procedure necessarily requires the use of force-fields i) consistent with the level of theory that is used to compute the explicit electronic problem for the reaction and ii) accurate enough far from the reference geometry for which they were fitted.
Ultimately, meeting these requirements is a non-trivial challenge, generally impossible in many cases, particularly for large systems. In addition, previous research claimed that for several reactions the resulting EVB energy leads to large errors along the MEP, especially in the transition region where, artificial minima are created in the worst cases. To overcome these limitations, a combination of Gaussian functions were proposed to model the coupling terms, 43 thus offering a promising route for future calibration and development of EVB potentials.

Model study and MD settings
To test the implementation of EVB method and its extension to the NPT ensemble, we  Using the initial arrangement of atoms, the system was initially computed in the NVT ensemble at 300K using only FF 1 for MA and a Nose-Hoover 61,62 thermostat with a relaxation time of 0.5 ps. Equilibration was conducted for 5 ps, scaling the system temperature every 5 fs and resampling instantaneous system momenta distribution every 9 fs. Production MD followed for 30 ps. The last snapshot with positions, velocities and forces served as the starting point of a NPT simulation at 1 atm, this time using a Nose-Hoover thermostat and barostat with relaxation times of 3.0 and 1.0 ps for the thermostat and barostat, respectively.
Equilibration was conducted for 2 ps while allowing for a variation of 10% in the system density. This was followed by a 200 ps of MD production run. The average supercell dimension was used for all the EVB-NVT simulations and as starting point for the EVB-NPT runs. All Even though this choice represents a departure from a more realistic chemistry, lower energy barriers are more convenient to MD, as less computational time is needed to switch between conformations, thus allowing a better sampling of the configurational space.
On the other hand, energy barriers computed with the EVB method depend on the term C 12 . Here, we have assumed C 12 to be a constant. Note that the choice of a constant for the coupling term, even trivial, complies with the functional form requirement for the coupling terms to compute the EVB stress tensor. Figure 4 shows that the computed EVB barrier for MA decreases as the value of C 12 increases. For the adopted OPLS-2005-FFs and to the purpose of comparison with the DFT energy barriers of Fig. 3, we have only considered values of C 12 in the range between 46.5 and 50.0 kcal/mol. These results constitute an example of how coupling terms can be used to calibrate EVB potentials against a reference value, in this case obtained from DFT. The advantage of the EVB formalism lies in the assumption that the coupling terms calibrated for reactive molecules in gas phase do not change significantly when transferring the reactive system from one phase to the another. 24 This approximation has been rigorously validated via Constrained-DFT calculation. 30 Moreover, the use of constants for the coupling terms is the most common choice in the execution of EVB simulations for solvated reactive sites, as a constant can be finely adjusted until a calculated property (usually free energy) agrees with the experimental value. 24,31 Here, we are not interested in comparing with experiments but evaluating how the results are affected Clearly, the reactive EVB potential is different from any of the individual non-reactive FFs, particularly in the TS region. Thus, it is natural to argue as to which extent this reactive EVB potential affects the stress tensor and the converged volume (and density), and how results compare with the density resulting from a NPT simulation using one of the involved FF, as in the standard protocol. If the energy barrier is sufficiently large, the system will only sample the vicinity of one of the possible configurations. Even though the system might occasionally swap conformation, the TS region will be hardly sampled. Thus, for cases where the conformations are chemically equivalent and the barrier is large enough, the standard protocol appears to be a sensible approximation. In contrast, if the barrier is sufficiently low, the TS region will be better sampled during the course of a MD simulation and the average potential will depart from any of the individual FFs.
Following a NPT simulation of the solvated MA using only one of the two FFs at 1 atm and 300 K, the computed density of the system is predicted to be in the region between the To the best of our knowledge, such probability distributions have not been reported before within the framework of EVB. Moreover, these computed distributions support previous research that suggest the convenience of ǫ as an alternative reaction coordinate. 22,24,32,34 In fact, in contrast to the MEP, ǫ represents a coordinate that implicitly accounts not only for the reactive site but also includes the effect of the surrounding solvent.
The symmetric distributions obtained for the NVT and NPT ensembles indicate a good level of agreement. Nevertheless, it is convenient to quantify this agreement for a better comparison. To this purpose, we first estimate the width for both peaks, which is of the order of 320 kcal/mol. This range for ǫ 12 was used to define the TS domain, indicated by the region between the dashed lines of Fig. 6, located at ǫ T S 12 = ±160 kcal/mol. This TS region is also adopted to be independent of the reactive model for MA. We define the probability for the the system to be in the TS region as follows where P(ǫ 12 ) is the normalized distribution obtained from the histograms, as shown in Fig.   6. Clearly, the computed TS probability depends on the choice for the extension of the TS region, which is completely arbitrary. However, Eq. (23) provides a method to quantify the relevance of the TS region and compare the probability distributions for different ensembles. Figure 7 shows the computed values for the TS Probabilities for the different reactive models of MA. As expected, the classical probability for the system to sample the TS region increases as the barrier reduces. Additionally, the computed probability distributions do not depend on the assumed ensemble for the present model.
Based on these results, we conclude that the standard protocol to compute single reactive sites in water solution with EVB using the NVT ensemble is already a remarkably good approximation to the self-consistent EVB-NPT simulation, at least for concentrations of reactive sites lower than 9.19 × 10 −2 molality. In the present case study of solvated MA, the flexible-reactive MA is composed of nine atoms and contributes with 27 degrees of freedom  It would be interesting to investigate if the dynamics will differ on increasing the concentration of the reactive MA. As discussed in section 6, due to restriction on the size of the In comparison to the standard protocol of converging the system volume using only one of the non-reactive FFs, results from EVB-NPT simulations for different reactive models of MA (i.e. different energy barriers) demonstrate a negligible effect of the EVB potential in the computed density of the solution.
In addition, we performed EVB-NVT simulations using the converged volume from the standard protocol. To compare the sampling of the configurational space with respect to EVB-NPT simulations, we use the energy gap as a variable to compute probability distributions of the reactive system. Detailed analysis of these distributions also demonstrate a negligible difference between both ensembles. We attribute these findings to the relative low concentration for the model of MA in water, where the non-reactive dynamics of the rigid water molecules dominates over the reactive dynamics of the single MA. Therefore, future EVB simulations of solvated reactive molecules with higher concentrations would be beneficial to quantify the role of the reactive dynamics on the whole system. To this purpose, developments to extend the EVB method are required to accommodate multiple reactive sites at the same time.

Settings for DFT simulations
Geometry optimization and NEB calculations for the two conformations of MA were conducted using the library DL_FIND 89  where the index type runs over all type of possible interactions (bonds, angles, coulombic, etc). In contrast, C mk (ǫ mk ) cannot be decomposed in terms of the individual contributions U (m) type . Consequently, matrix H EVB cannot be decomposed for each type of interaction. One might consider the particular case of constant coupling terms C mk (ǫ mk ) = C mk , ∀m, k = 1, · · · , N F , with m = k to check if a separation into individual terms is possible. For the sake of simplicity, let us consider the case of two FFs with C 12 = C 21 . Without loss of generality, we can write C 12 as a sum of a set of constants

S2-On the EVB stress tensor and virial decomposition
The purpose of this section is to demonstrate that the EVB stress tensor σ EVB αβ and virial V EVB can be decomposed in different components according to the type of the interaction.
We note that σ