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

- Department of Chemistry, Indian Institute of Technology Kanpur, Kanpur, India

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.

## 1. 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. (2016, 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.

## 2. Theory

### 2.1. 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 Laio et al. (2002) 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

having the shell *k* harmonically bound to the core atom $I\in {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}_{\text{b}}^{\prime}$ and ${E}_{\text{vdw}}^{\prime}$ 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 ($\left\{{q}_{I}^{c}\right\}$) and the shell ($\left\{{q}_{k}^{s}\right\}$) degrees of freedom with the electronic density $\rho (\overline{\text{r}})$, respectively. Electrostatic interactions are computed in the real space with the modified Coulomb kernel as in Laio et al. (2002).

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, ${\dot{\text{S}}}_{I}$ is the velocity of the center of mass of a core–shell pair (*I, k*), and ${\dot{\text{s}}}_{k}$ is the relative velocity of the shell *k* connected to a core atom *I*. Here,

${{M}}_{I}$ and ${\overline{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.

### 2.2. 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 ${\text{F}}_{I}^{(S)}$ 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 ${\text{f}}_{k}^{(\text{s})}$ is the force acting on the relative coordinate **s**_{k} of a core–shell pair (*I, k*). The thermostat variables for the relative motion are $\{{\eta}_{j}^{*}$}, having masses $\left\{{q}_{j}^{*}\right\}$, and their conjugate momenta are given by $\left\{{p}_{{\eta}_{j}^{*}}\right\}$. 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 ${\text{f}}_{k}^{(r)}$ 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).

## 3. 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.

## 4. Results and Discussion

### 4.1. 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.

**Figure 1**. The system of five water molecules used for validating the QM/p–MM implementation. The QM water molecule is highlighted by a circle. The core and the shell of O atoms are indicated by red spheres and green transparent spheres, respectively. H atoms are in white color.

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.

**Figure 2**. *NVE* MD simulation using QM/p–MM implementation for 1H_{2}O (QM) + 4H_{2}O (p–MM) system: **(A)** total energy (violet) and potential energy (green) **(B)** orbital kinetic energy, and **(C)** shell temperature. All the energies are in Hartree unit and temperature is in Kelvin. *E*_{Cons} is the drift in total energy per atom per ps.

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.

**Figure 3**. *NVT* MD simulation using QM/p–MM implementation for 1H_{2}O (QM) + 4H_{2}O (p–MM) system: **(A)** conserved energy (violet) and potential energy (green) **(B)** orbital kinetic energy **(C)** shell temperature *T*_{s} **(D)** distributions of physical temperature *T*_{phys} (red line with lower X–axis) and shell temperature *T*_{s} (blue line with upper X-axis). All the energies are in Hartree unit and the temperature is in Kelvin.

### 4.2. 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) level of theory for the QM part and p-MZHB potential for the MM part. We used a supercell of 8 × 8 × 8 (Si_{2048} O_{4096}) for QM/p–MM calculations. Optimized structure and lattice parameters using the p–MZHB MM potential were used here. Multiple QM/p–MM calculations were carried out with different QM sizes: 2T (Si_{2}O_{7}), 8T (Si_{8}O_{25}), 14T (Si_{14}O_{40}), and 26T (Si_{26}O_{67}), where T stands for SiO_{4} tetrahedral unit (see Figure 4).

**Figure 4**. QM subsystems used in the four QM/MM calculations of oxygen vacancy in α–Cristobalite: **(A)** 2T (Si_{2}O_{7}), **(B)** 8T (Si_{8}O_{25}), **(C)** 14T (Si_{14}O_{40}), and **(D)** 26T (Si_{26}O_{67}). In the vacancy calculations, the highlighted (by blue circle) O atom was removed from its lattice position. Note that MM atoms are not shown here. Atom colors: Si (yellow), O (red).

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).

**Table 2**. The average value of Si–O bond length (Å), O–Si–O, and Si–O–Si angles (°) of α–cristobalite computed from MD simulations in *NVT* ensemble at 300 K using p–MZHB, MZHB (Sahoo and Nair, 2015), QM/p–MM (14T), and QM potentials.

Next, the vacancy formation energy, (Δ*E*_{f}) in α–cristobalite silica was then computed, as

where *E*(O_{2}), *E*_{diss}(O_{2}), *E*(SiO_{2−x}), 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.16eV 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.

**Table 3**. Δ*E*_{f} computed from the single point energy calculations using “all–QM” with PBE density functional.

Δ*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 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.

**Figure 5**. Probability distribution of *r*_{SiSi} distance obtained from the MD simulations at 300 K using “all–QM” (3 × 3 × 3 supercell), QM/MM (14T) and QM/p–MM (14T) potentials. The vertical lines show the corresponding average values.

### 4.3. 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**).

**Figure 6. (A)** Mechanism of hydrogenation of ethene catalyzed by Rh_{3}H_{7}/Y. **(B)** Snapshots of the reactant(**1.2**), the intermediate(**1.3**) and the product(**1.4**). Atom colors: Rh–ochre, Si–yellow, Al–green, O–red, C–black, and H–white. **(C)** Free energy profiles computed from the QM/MM and the QM/p–MM metadynamics simulations are given.

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.

### 4.4. 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.

**Figure 7. (A)** Mechanism of the proton exchange reaction between a methane molecule and the H–ZSM–5 zeolite. **(B)** Snapshots of the reactant, transition state, and product states taken from QM/p–MM MD simulations. Here the atoms colors used are Si (yellow), Al (green), O (red), C (black), and H (white). **(C)** Free energy profiles constructed from TASS simulations are shown here. The free energy profiles were obtained by projecting the six–dimensional free energy surface on the CV, *d*_{H1−C2} − *d*_{C2−H3}, which is the difference in the distances H1–C2 and C2-H3.

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 (${\text{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}.

## 5. 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.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling Editor and reviewer TH declared their involvement as co-editors in the Research Topic, and confirm the absence of any other collaboration.

## 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).

## Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2018.00275/full#supplementary-material

## References

Arzumanov, S. S., Moroz, I. B., Freude, D., Haase, J., and Stepanov, A. G. (2014). Methane activation on in modified ZSM-5 zeolite. H/D hydrogen exchange of the alkane with Brønsted acid sites. *J. Phys. Chem. C* 118, 14427–14432. doi: 10.1021/jp5037316

Awasthi, S., and Nair, N. N. (2017). Exploring high dimensional free energy landscapes: temperature accelerated sliced sampling. *J. Chem. Phys.* 146:094108. doi: 10.1063/1.4977704

Bakowies, D., and Thiel, W. (1996). Hybrid models for combined quantum mechanical and molecular mechanical approaches. *J. Phys. Chem.* 100, 10580–10594.

Boulanger, E., and Thiel, W. (2012). Solvent boundary potentials for hybrid QM/MM computations using classical Drude oscillators: a fully polarizable model. *J. Chem. Theory Comput.* 8, 4527–4538. doi: 10.1021/ct300722e

Boulanger, E., and Thiel, W. (2014). Toward QM/MM simulation of enzymatic reactions with the Drude oscillator polarizable force field. *J. Chem. Theory Comput.* 10, 1795–1809. doi: 10.1021/ct401095k

Bučko, T., Benco, L., Hafner, J., and Ángyán, J. G. (2007). Proton exchange of small hydrocarbons over acidic chabazite: *Ab initio* study of entropic effects. *J. Catal.* 250, 171–183. doi: 10.1016/j.jcat.2007.05.025

Car, R., and Parrinello, M. (1985). Unified approach for molecular dynamics and density-functional theory. *Phys. Rev. Lett.* 55, 2471–2474.

Chu, Y., Xue, N., Xu, B., Ding, Q., Feng, Z., Zheng, A., et al. (2016). Mechanism of alkane H/D exchange over zeolite H-ZSM-5 at low temperature: a combined computational and experimental study. *Catal. Sci. Technol.* 6, 5350–5363.

CPMD (v 13.2). CPMD Program Package. *IBM Corp 1990-2018, MPI für Festkörperforschung Stuttgart 1997-2018*. Available online at: http://www.cpmd.org.

de Leeuw, N. H., and Parker, S. C. (1998). Molecular-dynamics simulation of MgO surfaces in liquid water using a shell-model potential for water. *Phys. Rev. B* 58, 13901–13908.

Dick, B. G., and Overhauser, A. W. (1958). Theory of the dielectric constants of alkali halide crystals. *Phys. Rev.* 112, 90–103.

Downs, R. T., and Palmer, D. C. (1994). The pressure behavior of alpha-cristobalite. *Am. Min.* 79, 9–14.

Gabrienko, A. A., Arzumanov, S. S., Toktarev, A. V., Freude, D., Haase, J., and Stepanov, A. G. (2011). Hydrogen H/D exchange and activation of C1n-C4 alkanes on Ga-modified zeolite BEA studied with 1H magic angle spinning nuclear magnetic resonance *in situ*. *J. Phys. Chem. C* 115, 13877–13886. doi: 10.1021/jp204398r

Gale, J. D. (1997). GULP: a computer program for the symmetry-adapted simulation of solids. *J. Chem. Soc. Faraday Trans.* 93, 629–637.

Geerke, D. P., Thiel, S., Thiel, W., and van Gunsteren, W. F. (2007). Combined QM/MM molecular dynamics study on a condensed-phase SN2 reaction at nitrogen: The effect of explicitly including solvent polarization. *J. Chem. Theory Comput.* 3, 1499–1509. doi: 10.1021/ct7000123

Grimme, S. (2006). Semiempirical GGA-type density functional constructed with a long-range dispersion correction. *J. Comput. Chem.* 27, 1787–1799. doi: 10.1002/jcc.20495

Illingworth, C. J., Gooding, S. R., Winn, P. J., Jones, G. A., Ferenczy, G. G., and Reynolds, C. A. (2006). Classical polarization in hybrid QM/MM methods. *J. Phys. Chem. A* 110, 6487–6497. doi: 10.1021/jp046944i

Jorgensen, J. D. (1978). Compression mechanisms in α–quartz structures-SiO_{2} and GeO_{2}. *J. Appl. Phys.* 49, 5473–5478.

Laio, A., and Parrinello, M. (2002). Escaping free-energy minima. *Proc. Natl. Acad. Sci. U.S.A.* 99, 12562–12566. doi: 10.1073/pnas.202427399

Laio, A., VandeVondele, J., and Rothlisberger, U. (2002). A Hamiltonian electrostatic coupling scheme for hybrid Car–Parrinello molecular dynamics simulations. *J. Chem. Phys.* 116, 6941–6947. doi: 10.1063/1.1462041

Lamoureux, G., and Roux, B. (2003). Modeling induced polarization with classical Drude oscillators: theory and molecular dynamics simulation algorithm. *J. Chem. Phys.* 119, 3025–3039. doi: 10.1063/1.1589749

Lev, B., Zhang, R., de la Lande, A., Salahub, D., and Noskov, S. Y. (2010). The QM-MM interface for CHARMM–deMon. *J. Comput. Chem.* 31, 1015–1023. doi: 10.1002/jcc.21387

Liang, A. J., and Gates, B. C. (2008). Time-resolved structural characterization of formation and break-up of rhodium clusters supported in highly dealuminated Y zeolite. *J. Phys. Chem. C* 112, 18039–18049. doi: 10.1021/jp805917g

Lide, D. R. (2005). *Crc Handbook. CRC Handbook of Chemistry and Physics, Internet Version 2005.* Boca Raton, FL: CRC Press. Available online at: http://www.hbcpnetbase.com

Lindan, P. J. D., and Gillan, M. J. (1993). Shell-model molecular dynamics simulation of superionic conduction in CaF_{2}. *J. Phys. Condens. Matter* 5:1019.

Loco, D., Lagardère, L., Caprasecca, S., Lipparini, F., Mennucci, B., and Piquemal, J.-P. (2017). Hybrid QM/MM molecular dynamics with AMOEBA polarizable embedding. *J. Chem. Theory Comput.* 13, 4025–4033. doi: 10.1021/acs.jctc.7b00572

Loco, D., Polack, É., Caprasecca, S., Lagardére, L., Lipparini, F., Piquemal, J. P., et al. (2016). A QM/MM approach using the AMOEBA polarizable embedding: from ground state energies to electronic excitations. *J. Chem. Theory Comput.* 12, 3654–3661. doi: 10.1021/acs.jctc.6b00385

Lu, Z., and Zhang, Y. (2008). Interfacing *ab initio* quantum mechanical method with classical Drude osillator polarizable model for molecular dynamics simulation of chemical reactions. *J. Chem. Theory Comput.* 4, 1237–1248. doi: 10.1021/ct800116e

Martyna, G. J., Klein, M. L., and Tuckerman, M. (1992). Nosé–hoover chains: the canonical ensemble via continuous dynamics. *J. Chem. Phys.* 97, 2635–2643.

Marx, D., and Hutter, J. (2009). *Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods.* Cambridge: Cambridge University Press.

Marx, D., Tuckerman, M. E., and Martyna, G. J. (1999). Quantum dynamics via adiabatic *ab initio* centroid molecular dynamics. *Comp. Phys. Commun.* 118, 166–184.

Mitchell, P. J., and Fincham, D. (1993). Shell model simulations by adiabatic dynamics. *J. Phys.: Condens. Matter* 5:1031.

Nasluzov, V. A., Ivanova, E. A., Shor, A. M., Vayssilov, G. N., Birkenheuer, U., and Rösch, N. (2003). Elastic polarizable environment cluster embedding approach for covalent oxides and zeolites based on a density functional method. *J. Phys. Chem. B* 107, 2228–2241. doi: 10.1021/jp026742r

Niklasson, A. M., Steneteg, P., Odell, A., Bock, N., Challacombe, M., Tymczak, C. J., et al. (2009). Extended Lagrangian Born–Oppenheimer molecular dynamics with dissipation. *J. Chem. Phys.* 130:214109. doi: 10.1063/1.3148075

Perdew, J. P., Burke, K., and Ernzerhof, M. (1996). Generalized gradient approximation made simple. *Phys. Rev. Lett.* 77, 3865–3868.

Plimpton, S. (1995). Fast parallel algorithms for short-range molecular dynamics. *J. Comput. Phys.* 117:1–19.

Riahi, S., and Rowley, C. N. (2014). The CHARMM-Turbomole interface for efficient and accurate QM/MM molecular dynamics, free energies, and excited state properties. *J. Comput. Chem.* 35, 2076–2086. doi: 10.1002/jcc.23716

Rowley, C. N., and Roux, B. (2012). The solvation structure of Na+ and K+ in liquid water determined from high level *ab initio* molecular dynamics simulations. *J. Chem. Theory Comput.* 8, 3526–3535.

Sahoo, S. K., and Nair, N. N. (2015). A potential with low point charges for pure siliceous zeolites. *J. Comput. Chem.* 36, 1562–1567. doi: 10.1021/ct300091w

Sahoo, S. K., and Nair, N. N. (2016). CPMD/GULP QM/MM interface for modeling periodic solids: implementation and its application in the study of Y–zeolite supported Rh_{n} clusters. *J. Comput. Chem.* 37, 1657–1667. doi: 10.1002/jcc.23968

Sangster, M., and Dixon, M. (1976). Interionic potentials in alkali-halides and their use in simulations of molten-salts. *Adv. Phys.* 25, 247–342.

Schlick, T. (2010). *Molecular Modeling and Simulation: An Interdisciplinary Guide, 2nd Edn.* New York, NY: Springer.

Sprik, M. (1991). Computer simulation of the dynamics of induced polarization fluctuations in water. *J. Phys. Chem.* 95, 2283–2291.

Sulimov, V. B., Sushko, P. V., Edwards, A. H., Shluger, A. L., and Stoneham, A. M. (2002). Asymmetry and long-range character of lattice deformation by neutral oxygen vacancy in α-quartz. *Phys. Rev. B* 66:024108. doi: 10.1103/PhysRevB.66.024108

Todorov, I. T., Smith, W., Trachenko, K., and Dove, M. T. (2006). DL_POLY_3: new dimensions in molecular dynamics simulations via massive parallelism. *J. Mater. Chem.* 16, 1911–1918. doi: 10.1039/b517931a

Truitt, M. J., Toporek, S. S., Rovira-Truitt, R., and White, J. L. (2006). Alkane C–H bond activation in zeolites: evidence for direct protium exchange. *J. Am. Chem. Soc.* 128, 1847–1852. doi: 10.1021/ja0558802

Tuma, C., and Sauer, J. (2015). Quantum chemical *ab initio* prediction of proton exchange barriers between CH_{4} and different H-zeolites. *J. Chem. Phys.* 143:102810. doi: 10.1063/1.4923086

Vollmer, J. M., and Truong, T. N. (2000). Mechanisms of hydrogen exchange of methane with h-zeolite y: an *ab initio* embedded cluster study. *J. Phys. Chem. B* 104, 6308–6312. doi: 10.1021/jp0008445

Wilson, M., and Madden, P. A. (1993). Polarization effects in ionic systems from first principles. *J. Phys. Condens. Matter* 5:2687.

Woodcock, H. L., Hodošček, M., Gilbert, A. T. B., Gill, P. M. W., Schaefer, H. F., and Brooks, B. R. (2007). Interfacing Q-Chem and CHARMM to perform QM/MM reaction path calculations. *J. Comput. Chem.* 28, 1485–1502. doi: 10.1002/jcc.20587

Keywords: QMMM simulations, MD, POLARIZED MM, catalysis, CPMD-GULP

Citation: Sahoo SK and Nair NN (2018) Interfacing the Core-Shell or the Drude Polarizable Force Field With Car–Parrinello Molecular Dynamics for QM/MM Simulations. *Front. Chem*. 6:275. doi: 10.3389/fchem.2018.00275

Received: 22 February 2018; Accepted: 18 June 2018;

Published: 10 July 2018.

Edited by:

Sam P. De Visser, University of Manchester, United KingdomReviewed by:

Thomas S. Hofer, University of Innsbruck, AustriaHugo Gattuso, University of Liége, Belgium

Copyright © 2018 Sahoo and Nair. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Nisanth N. Nair, nnair@iitk.ac.in

^{†}Present Address: Sudhir K. Sahoo, Department Chemie, Universität Paderborn, Paderborn, Germany