The outcome of a Molecular Dynamics (MD) simulation depends on the adopted Force Field (FF) to describe the interactions between the different constituents of the system under study. Such interactions are often modelled by simple equations with parameters fitted to experimental data or quantum mechanical calculations. In most of the classical FFs available, functional forms and fitted parameters remain unchanged during the MD simulation. In reactive process, however, the nature of the interactions is expected to change due to the formation of new species. Thus, standard FFs fail to describe chemical reactions. Here we shall refer to them as non-reactive FFs.
Reactive FFs (R-FF) offer an alternative to simulate chemical reactions with MD. Figure 1 shows a simplistic 1D picture of the potential energy landscape of a system with two possible and chemically different states (State 1 and State 2) and a Transition State (TS). Assuming that the black curve represents the “exact” energy landscape between State 1 and State 2, a R-FF is fitted to model the interactions in the whole domain (dashed blue).
Nevertheless, designing R-FFs for MD simulations is very challenging, since it necessarily entails to tackle a multi-dimensional problem, where the modelled interactions require of highly complicated functional forms with many strongly coupled parameters that are optimised via a difficult search. Still, R-FFs have evolved considerably in the last years and led to computation of systems up to billions of atoms. Unfortunately, a general parametrization for the commonly used R-FFs is not available yet and, instead, parameters are tuned to specific chemical systems and environments.
Figure 1. Schematic 1D picture of a reactive process between State 1 and State 2 with the “exact” energy landscape (solid black). Left: A R-FF is fitted to model the interactions in the whole domain (dashed blue). Right: E1 is a non-reactive FF fitted that models the interactions in the region around the State 1 (dashed red), whereas E2 is a non-reactive FF fitted that models the interactions in the region around the State 2 (dashed green). The EVB method aims to find the best coupling between these two non-reactive FFs to approximate the energy landscape in the whole domain (EEVB, dotted violet).
In order to avoid the complications of building reliable R-FFs, we have chosen a different strategy for DL-POLY and decided to implement the Empirical Valence Bond (EVB) method, inspired in the seminal work of A. Warshel [1]. In the EVB method, simulation of reactive processes is conducted via the coupling of non-reactive FFs, as schematically shown in the 1D-view of Fig. 1 (right), one non-reactive FF for the reactant E1 (dashed red) and one non-reactive FF for the product E2 (dashed green) [2]. As discussed above, if we used the FFs built for the State 1 (E1), State 2 will never be sampled, since the values of E1 in the region of State 2 will be exceedingly large. Analogously, State 1 will never be sampled if we used E2 as the FF to describe the interactions. Instead, the EVB method defines a Hamiltonian (HEVB) whose matrix representation (Eqn. 1) has each of the non-reactive FFs as diagonal components, whereas the off-diagonal terms are given by the coupling C12 between the FFs in the reaction:
where all terms obviously depend on the set of atomic coordinates R. Following diagonalization of HEVB we obtain two possible eigenvalues, λ±,
We define the EVB energy (EEVB) as the lowest value with the corresponding EVB eigenvector ψEVB
Consequently, the EVB force over the atom “I” is given by
Equations (3) and (4) define the EVB force field (EVB-FF). The main advantage of the EVB method with respect to the use of R-FF is that non-reactive FFs are already available in the literature. Nevertheless, no-reactive FFs must be sufficiently accurate for configurations sufficiently far from the reference geometry for which they were fitted. In addition, the quality of the EVB method to model the interactions in the intermediate region in between states (in the region of the Transition State) depends on the choice of the coupling term C12. The strategy is to fit the parameters of C12 using a suitable coordinate that connects State 1 and 2 such that the difference between the reference energy (black curve of Fig. 1) and EEVB along this path is minimised. Following the work of Hartke et. al [3] a convenient coordinate is the zero-temperature minimum energy path.
As a preliminary example of the implementation of the EVB method in DL-POLY we show the case of a single malonaldehyde in gas phase, as schematically depicted in Figure 2. This molecule exhibits two isomeric conformations (Fig. 2a-2c) with a transition state where the hydrogen atom is equally shared between the two oxygen atoms (Fig. 2b). If we used only a FF fitted for the conformation where the H is bound to the oxygen at the left, the configuration with the H bound to the oxygen at the right is never sampled (Fig. 3a). In contrast, when we couple the FFs fitted for conformations of Fig 2a. and Fig. 2c via the EVB method, both configurations are sampled, as long as the total energy of the molecule is sufficiently high to overcome the energy barrier in between both conformations (Fig. 3b).
Figure 2. Schematic representation of the isomerization of malonaldehyde.
Figure 3. MD simulation of an isolated malonaldehyde molecule. Left: the conformation with the H atom bound to the oxygen in the right (as shown in Fig 2c) is not sampled when we only use the FF fitted for the reference geometry with the H atom bound to the left oxygen (Fig. 2a). Right: isomerization is obtained when coupling the two FFs for the two different conformations via the EVB method.
[1] A. Warshel, J. Am. Chem. Soc., 1980, 102, 6218.
[2] The EVB method can be easily generalised to N-coupled non-reactive FFs.
[3] B.Hartke and S. Grimme, Phys.Chem.Chem.Phys., 2015, 17, 16715.