Interfacing the Core-Shell or the Drude Polarizable Force Field With Car–Parrinello Molecular Dynamics for QM/MM Simulations

We report a quantum mechanics/polarizable–molecular mechanics (QM/p–MM) potential based molecular dynamics (MD) technique where the core–shell (or the Drude) type polarizable MM force field is interfaced with the plane-wave density functional theory based QM force field which allows Car–Parrinello MD for the QM subsystem. In the QM/p-MM Lagrangian proposed here, the shell (or the Drude) MM variables are treated as extended degrees of freedom along with the Kohn–Sham (KS) orbitals describing the QM wavefunction. The shell and the KS orbital degrees of freedom are then adiabatically decoupled from the nuclear degrees of freedom. In this respect, we also present here the Nosé–Hoover Chain thermostat implementation for the dynamical subsystems. Our approach is then used to investigate the effect of MM polarization on the QM/MM results. Especially, the consequence of MM polarization on reaction free energy barriers, defect formation energy, and structural and dynamical properties are investigated. A low point charge polarizable potential (p–MZHB) for pure siliceous systems is also reported here.


INTRODUCTION
Hybrid quantum mechanical/molecular mechanical (QM/MM) calculations offer a powerful way to bridge the length scales in a chemically complex system where a small region of the system of interest is treated by QM techniques, while the rest is described by computationally cheap MM force-fields. Widely used MM force fields employ a fixed point charge model for accounting the electrostatic interactions between MM atoms. The QM/MM implementations with fixed charge MM models enable polarization of QM charge density due to MM electrostatic potential. However, such approaches cannot take into account the polarization of MM atoms due to the QM electrostatic potential. Inclusion of polarization of MM atoms in QM/MM simulations demands usage of polarizable MM force fields, i.e., QM/polarized-MM (QM/p-MM) methods. It was reported that inclusion of polarization of MM atoms has significant effects on various properties (Bakowies and Thiel, 1996;Illingworth et al., 2006;Geerke et al., 2007;Lu and Zhang, 2008;Boulanger and Thiel, 2014), for instance, free energy barriers of chemical reactions are affected by about 10% (Lu and Zhang, 2008;Boulanger and Thiel, 2014).
The shell model (Dick and Overhauser, 1958) (or the Drude oscillator model) is widely used to describe polarization of MM atoms. Molecular dynamics (MD) simulations with the shell model based MM force fields can be carried out in two ways. In the conventional scheme, the position of the shells are minimized (Sangster and Dixon, 1976;Lindan and Gillan, 1993) at every MD step while keeping the core coordinates fixed. In an alternative scheme, the shells are treated as extended degrees of freedom and these are propagated classically to avoid minimization of their positions (Sprik, 1991;Mitchell and Fincham, 1993;Wilson and Madden, 1993) in the spirit of the Car-Parrinello MD method (Car and Parrinello, 1985). Often, the shell variables are assigned a mass smaller compared to that of the nuclear masses. Due to this reason, a smaller time step than used in a conventional MD is required for this approach. In practice, the shell temperature is kept close to 0 K, and most importantly, dynamics of shell degrees of freedom is made adiabatically decoupled from the rest of the system. Different QM/MM MD schemes have been proposed to interfere the shell model with ab initio methods (Sulimov et al., 2002;Nasluzov et al., 2003;Woodcock et al., 2007;Geerke et al., 2007;Lu and Zhang, 2008;Lev et al., 2010;Boulanger and Thiel, 2012;Rowley and Roux, 2012;Boulanger and Thiel, 2014;Riahi and Rowley, 2014). Ideally, wavefunction of the QM subsystem and positions of shells are minimized at every MD step. In the approach by Lu and Zhang (2008) the positions of the shells were either minimized or updated only once in every MD step. The shell variables are treated as extended degrees of freedom in the QM/MM scheme proposed by Boulanger and Thiel (2012). Similar method was also reported by Rowley and Roux (2012). Recently Loco et al. (2016Loco et al. ( , 2017 have presented a QM/MM coupling method using the AMOEBA polarizable force fields. Here, SCF (self consistent field) calculations were carried out for the induced dipoles, while either SCF or the extended Lagrangian variant of Born-Oppenheimer MD (Niklasson et al., 2009) technique was employed for the wavefunction update. Incorporating the shell model in the extended Lagrangian scheme for Car-Parrinello MD within a QM/MM implementation is, however, not straightforward, and has not been attempted before, to the best of our knowledge.
In this paper we present an extended Lagrangian scheme to carry out Car-Parrinello MD for the QM subsystem which is coupled to a polarizable shell model based MM force field. First, we discuss the theory and the technical details of our method. Next, the implementation is validated by taking a system of water cluster composed of five water molecules. A new polarizable MM potential for silica with low point charges is then developed. Using our implementation and this new force-field for silica, we study three problems: (a) Oxygen vacancy in α-cristobalite silica; (b) Hydrogenation of ethene catalyzed by Rh cluster supported in Y-zeolite; (c) Proton exchange between methane and H-ZSM-5 zeolite.

Formulation of the Extended Lagrangian QM/p-MM Method
The Lagrangian for the conventional QM/MM Car-Parrinello MD simulation is, where {M I } and {µ i } are the masses of the ionic and the orbital degrees of freedom, respectively, and {R I } and {φ i } are the nuclear coordinates (in Cartesian) and the Kohn-Sham orbitals, respectively. Here, E KS , E MM , and E QM/MM are the energy of the QM subsystem, the energy of the MM subsystem and the energy due to QM-MM non-bonding interactions, respectively. The last term in the Lagrangian invokes orthonormality constraints during the classical evolution of the Kohn-Sham orbitals (Marx and Hutter, 2009). For details of this implementation, see  and Sahoo and Nair (2016).
In the case of QM/p-MM implementation, MM polarization is included by augmenting the degrees of freedom by the shells {r k }, which are connected to a selected set of polarizable atoms P. We propose the QM/p-MM Lagrangian, where, n s is the number of shells, and m k is the fictitious mass of a shell k. Also, where s k = |s k | = |R I − r k | , having the shell k harmonically bound to the core atom I ∈ P.
Here E b refers to the sum of all the bonding terms in the MM potential, which is conventionally defined over the shell variables.
The total non-bonding interaction energy, E nb , is the sum of the dispersive and the electrostatic interactions within the MM subsystem. The dispersive interactions are defined over the shell degrees of freedom, while the electrostatic interactions span over the cores and the shells degrees of freedom. Further, where E ′ b and E ′ vdw are the energy contributions due to the bonding and the dispersive interactions between the QM and the MM atoms, respectively. The last two terms in the above equation account for the electrostatic interaction between the point charges of the core ({q c I }) and the shell ({q s k }) degrees of freedom with the electronic density ρ(r), respectively. Electrostatic interactions are computed in the real space with the modified Coulomb kernel as in .
For the success of this approach it is crucial that the Lagrangian in Equation (2) leads to dynamics close to that on the Born-Oppenheimer surface. This is taken care by starting the MD simulations with the optimized {φ} and {r} , while maintaining the temperatures of the orbitals (T φ ) and the shells (T s ) degrees of freedom close to zero, considering The physical temperature (T phys ) is defined as, while the shell temperature is defined as, In the above equations,Ṡ I is the velocity of the center of mass of a core-shell pair (I, k), andṡ k is the relative velocity of the shell k connected to a core atom I. Here, M I and m k are the total mass of a core-shell pair (I, k), and the reduced mass of a shell k, respectively: In the above, k B is the Boltzmann constant, N f is the total nuclear degrees of freedom and n s is the total number of shell variables, respectively. We have implemented this method in the CPMD/GULP QM/MM interface program, as developed in Sahoo and Nair (2016), where the plane wave density functional theory (DFT) based CPMD (CPMD, 132) code is interfaced with the MM based GULP (Gale, 1997) program. At this stage, the following points are noted: 1. Both the orbital and the shell degrees of freedom have to be adiabatically separated from the nuclear degrees of freedom. While the nuclei could be kept hot, both the orbitals and the shells variables should remain cold throughout the dynamics. 2. As per Equation (5), it is not critical to have an adiabatic separation between shell and orbital degrees of freedom, provided that both sets of variables remain cold.
Accordingly, we have strategized the application of our implementation. The time step of integration and the masses of both shell and orbital degrees of freedom have to be chosen such that adiabatic separation between the nuclear subsystem and the subsystem containing shells and orbitals is maintained. We choose m k = µ (i.e., the masses of the shell and the orbital degrees of freedom are taken to be the same), which allows us to choose the same time step for integrating the equations of motion for all the subsystems.

Implementation of Nosé-Hoover Chain Thermostat for Shell Dynamics
For obtaining stable dynamics and to achieve a canonical ensemble, it is crucial to couple the dynamical subsystems with thermostats. We have implemented three separate sets of Nosé-Hoover Chain (NHC) thermostats (Martyna et al., 1992). The system temperature is maintained at T phys using one set of NHC thermostats whereas the shell and the orbital variables are maintained close to 0 K using two separate thermostats. We coupled the nuclear and the shell NHC thermostats to the center of mass motion and the relative motion of the core-shell pairs, and the corresponding equations of motion are given by, Here F (S) I is the force acting on the center of mass coordinates S I , {p η i } and {q i } are the momentum and the masses of the thermostat variables η. Number of thermostat variables, n c , is chosen to be more than one.
In order to thermostat the shell dynamics, we write the shell equations of motion in relative coordinates as, Here f (s) k is the force acting on the relative coordinate s k of a coreshell pair (I, k). The thermostat variables for the relative motion are {η * j }, having masses {q * j }, and their conjugate momenta are given by {p η * j }. In our practical implementation, we transform the Cartesian coordinates of a core-shell pair to the corresponding relative and the center of mass coordinates at every MD steps. Transformation of forces for a core-shell pair (I, k) is achieved by, Here, F I and f (r) k are the Cartesian forces on an atom (or core) I and on a shell k connected to core I, respectively. A similar coordinate transformation from cartesian to normal mode and vice-versa was reported by Marx et al. (1999) and can be also found in the work of Lamoureux and Roux (2003).

LOW POINT CHARGE POLARIZABLE FORCE FIELD FOR SILICA
Herein, the previously reported low point charge potential (MZHB) for silica (Sahoo and Nair, 2015) is further extended to include polarization of O atoms using the core shell model. This new potential is termed as p-MZHB hereafter. We have re-parameterized k θ of O-Si-O and Si-O-Si angles and optimized the core-shell coupling parameters {κ k }. During the parameterization, point charges of the core and the shell of O atoms were allowed to vary, while their sum was fixed to −0.35 e. The re-parameterization was done using the GULP (Gale, 1997) program to reproduce the experimental structure of α-Quartz (Jorgensen, 1978). The final set of parameters are given in Table 1. The detailed validation of the p-MZHB potential is discussed in the Supporting Information.

Validation of the Implementation in CPMD/GULP Interface Program
Total energy conservation during MD runs is tested to validate the method and the implementation. For the benchmarking purpose, we used 1H 2 O(QM)+4H 2 O(p-MM) system (Figure 1), and carried out NVE MD runs for 36 ps using DFT/PBE (Perdew et al., 1996) to describe the QM subsystem and polarizable de Leeuw- Parker (de Leeuw and Parker, 1998) potential to treat the MM subsystem; see Supporting Information for other technical details. In this MM potential, shells are added on the oxygen atoms of the water molecules. For the cross Lennard-Jones parameters, the Lorentz-Berthelot (Leach, 2010;Schlick, 2010) combination rule is applied. The total energy, the orbital kinetic energy and the shell temperature of the system were monitored; see Figure 2. The drift in total energy is only of the order of 10 −7 a.u. atom −1 ps −1 , indicating that the total energy conservation is fairly good. The plot of kinetic energy of orbitals shows only a small long-time drift, while the temperature of the shell variables are maintained at low temperature. Such small drifts in kinetic energy could be controlled by connecting the orbital degrees of freedom with a thermostat, as demonstrated below.
As next, we repeated the same simulation, but in the NVT ensemble. Here, three thermostats were added on the nuclear, the orbitals and the shells degrees of freedom, and these were thermostated to 300 K (T phys , 0.0007 Hartree and 1 K, respectively). Since the NHC thermostat posseses a conserved quantity, drift in this conserved quantity allows us to verify our implementation further. The drift in the conserved quantity is only of the order of 10 −7 a.u. atom −1 ps −1 , confirming the correctness of our implementation (see Figure 3A). The orbital kinetic energy and the shell temperature (see Figure 3B-D) plots clearly show that the dynamics of the extended variables is stable and is well thermostated.

Benchmark Studies Using α-Cristobalite Silica
To further benchmark the performance of the developed method, we carried out MD simulation of α-cristobalite silica in the NVT ensemble at 300 K. Here we use the DFT/PBE (Perdew et al., 1996) Figure 4).
The structure (inner QM atoms only) obtained from QM/p-MM simulation was compared with QM data ("all-QM") and MM data ("all-MM") data; see Table 2. Difference between "all-QM" and QM/p-MM data for Si-O bond length is only 0.01 Å, O-Si-O angle is only 0.2 • and Si-O-Si angle is only 0.8 • . However, structures near the boundary are deviating more from the "all-QM" data; see Supporting Information. This is expected due to the boundary effects in QM/MM calculations (as also seen in Sahoo and Nair, 2016).
Next, the vacancy formation energy, ( E f ) in α-cristobalite silica was then computed, as , and, E(SiO 2 ) are the energies of O 2 molecule (in the triplet electronic ground state), the dissociation energy of O 2 molecule, the energy of bulk silica with oxygen vacancy, and the energy of pure bulk silica, respectively. E(O 2 ) was computed from "all-QM" calculations whereas E(SiO 2−x ) and E(SiO 2 ) were computed either by "all-QM" or QM/p-MM calculations. E diss (O 2 ) = 5.16 eV was taken from the available experimental data (Lide, 2005). E f values computed from "all-QM" calculations with varying supercell sizes are listed in Table 3. The converged value of E f (w.r.t supercell size) is 8.73 eV. E f values computed from QM/p-MM calculations listed in Table 3, show that E f is nearly converged to 8.77 eV which is in excellent agreement with the "all-QM" data. However, it may  These results are also compared with experimental data (Downs and Palmer, 1994). For details see text. a Using p-MZHB MM potentials. b Using PBE density functional. c Using PBE density functional and p-MZHB MM force-field. be noted that the values predicted by QM/p-MM calculations are sensitive to the size of the QM subsystem. The E f value of 9.43 eV using 2T site is a very poor estimate compared with the "all-QM" data. Clearly, the computed value converges close to "all-QM" data with increase in the size of the QM subsystem (see Table 4).
The bond length distribution of the Si-Si bond (r SiSi ) at the defect site is then computed from the NVT trajectory using "all-QM" (using 3×3×3 supercell), QM/MM, and QM/p-MM potentials; see Figure 5. The average r SiSi value using "all-QM" potential is 2.38 Å. Using the non-polarized MZHB QM/MM simulations, we obtained a slightly increased value of 2.41 Å, while this was 2.40 Å in the polarized case. Overall, Figure 5 shows that the mean and the standard deviation of r SiSi from QM/p-MM MD agree better to the "all-QM" MD results than that from the QM/MM MD using the non-polarizable MM force-field.

Application: Hydrogenation of Ethene Catalyzed by Rh Clusters Supported in Y-Zeolite
In order to investigate the effect of polarization of MM atoms on the free energy barriers, we revisit the study of hydrogenation of ethene catalyzed by Rh clusters supported in Y-zeolites, as in our previous study (Sahoo and Nair, 2016). Based on the experimental results from the Gates group (Liang and Gates, 2008) and the previous study (Sahoo and Nair, 2016), we choose the hydrogenated cluster Rh 3 H 7 in Y-zeolite as the catalyst model. In order to model Y-zeolite, a supercell of size 2 × 2 × 2 of pure siliceous Y-zeolite (Si 1536 O 3072 ) was taken. The metal cluster, its ligands, and the metal-coordinated T25 site of the Y-zeolite were treated in the QM region (at the level of DFT/PBE+D2; Perdew et al., 1996;Grimme, 2006), while the rest of the zeolite was treated using MZHB or p-MZHB MM potentials.
The free energy profile for ethene hydrogenation was computed employing metadynamics (Laio and Parrinello, 2002) using the QM/MM and QM/p-MM implementations; see Figure 6 for the mechanism of the reaction studied here. For the details of the metadynamics simulation setup, see Supporting Information. In the metadynamics simulations using QM/MM and QM/p-MM methods, we observed the same reaction mechanism. Here, one of the hydrogen atoms from Rh 3 H 7 cluster moved to one of the C atoms of ethene forming an intermediate 1.3, which further reacted with a hydrogen atom on the Rh cluster to form ethane (1.4).
Interestingly, we observed differences of the order of 1 kcal mol −1 only in the free energy barriers computed from QM/p-MM and QM/MM computations; see Figure 6C. The main difference is in the stability of the ethyl intermediate 1.3 and on the free energy barrier for 1.3→1.4. However, the difference in free energy is only ∼ 1 kcal mol −1 which is close to the error associated with the metadynamics simulation. Thus, we conclude that for this specific reaction, the effect of MM polarization in the free energy estimates and on the reaction mechanism are negligible.

Application: Proton Exchange Between Methane and H-ZSM-5 Zeolite
Methane activation is one of the important steps in various industrially relevant processes such as conversion of methane to higher hydrocarbons and methanol. In this process, the C-H bond, which is quite inert as evident from the high bond formation energy (> 100 kcal mol −1 ), is activated for further functionalization. Of great importance, protonated zeolites, particularly H-ZSM-5, is known to activate the methane C-H bonds (Arzumanov et al., 2014;Chu et al., 2016). Here we look at the following chemical reaction ( Figure 7A): One of the acidic protons in the H-ZSM-5 zeolite is exchanged with one of the protons in the methane molecule during this reaction. The mechanism for C-H activation in zeolites is previously studied in the literature (Vollmer and Truong, 2000;Zheng and Blowers, 2005;Truitt et al., 2006;Bučko et al., 2007;Gabrienko et al., 2011;Arzumanov et al., 2014;Tuma and Sauer, 2015;Chu et al., 2016). Two types of mechanisms are proposed for this reaction: (a) direct and (b) bimolecular. In the direct mechanism, hydrogen is exchanged through a carbonium ion intermediate, i.e., via the formation of the penta-coordinated carbocation. In the bimolecular mechanism, the alkane molecule dissociates one of H atoms to the zeolite framework, forming an alkoxyl intermediate (Truitt et al., 2006), which subsequently undergoes reactions with other alkane molecules. Chu et al. (2016) have reported from both experimental and theoretical studies that higher alkanes react through bimolecular mechanism whereas the lower alkanes react via the direct mechanistic route. Our interest here is the direct mechanism as it involves the formation of a charged reactive intermediate and we anticipate some effect of MM polarization in the free energy estimates.
We could successfully simulate the proton exchange reaction using the Temperature Accelerated Sliced Sampling (TASS) method (Awasthi and Nair, 2017). This method allowed us to explore a high-dimensional free energy landscape composed of five collective variables (CVs). More details about the CVs used here and other technical details of the TASS simulation are given in the Supporting Information. The computed free energy profiles (projected along one of the crucial CVs) are given in Figure 7C. In TASS+QM/MM and TASS+QM/p-MM simulations, the hydrogen exchange reaction was found to proceed through the formation of a carbonium ion (CH + 5 ; see structure 2.2 in Figure 7B). The free energy barrier ( F ‡ ) for the reaction computed from QM/MM and QM/p-MM simulations are 39.5 and 36.0 kcal mol −1 , respectively. As expected, we are observing significant difference in the free energy barriers when polarized MM force-field is used. It is also noted in passing that the free energy barrier computed here are close to the potential energy barriers computed for similar reactions in zeolites in Vollmer and Truong (2000), Zheng and Blowers (2005), Bučko et al. (2007), and Tuma and Sauer (2015), which are in the range 29-38 kcal mol −1 .

CONCLUSIONS
An extended Lagrangian based implementation of QM/p-MM method that allows to perform conventional Car-Parrinello MD for the QM subsytem is presented here. In particular, we have discussed a combined scheme where the extended Lagrangian dynamics of the shell (or the Drude) variables are performed together with the Car-Parrinello dynamics of the KS orbitals. Inclusion of polarization does not increase the computational cost of the QM/MM Car-Parrinello MD simulations within our approach, mainly because we use the same time step as that of the conventional Car-Parrinello MD.
We find that, invoking polarization of MM atoms in QM/MM calculations only marginally improve the predictions of equilibrium structure and dynamics of non-charged systems.
On the other hand, when the transition state is charged, free energy barriers are considerably affected on the inclusion of MM polarization. We believe that the methods and the strategies developed here for a QM/polarized-MM implementation will be useful to study more complex problems in catalysis, reactions in solid-liquid interfaces, crystallization etc.
Other than these, we also report here a new low point charge polarizable potential for silica, which is suitable for performing QM/MM calculations. The potential is shown to be performing well and is able to reproduce bulk structures of various silica polymorphs. The developed MM potential has simple and commonly used potential functions with fewer parameters, making it easy to use in various simulation packages such as GULP, DL_POLY (Todorov et al., 2006) and LAMMPS (Plimpton, 1995).

AUTHOR CONTRIBUTIONS
NN has planned the work and SS have executed the implementation and computations. Both authors have analyzed the results and have written the paper.

ACKNOWLEDGMENTS
Authors acknowledge the HPC facility, IIT Kanpur. SS thanks CSIR New Delhi and IIT Kanpur for the Ph.D. fellowship. Authors are thankful to Dr. Shalini Awasthi (IIT Kanpur) for her help in setting up the metadynamics and the TASS simulations. Authors are also grateful for the valuable inputs from Prof. Alessandro Laio (SISSA, Italy).