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.

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,

$\begin{array}{l}{{L}}_{\text{CP}/\text{QMMM}}(\text{R},\dot{\text{R}},\varphi ,\dot{\varphi })={\displaystyle \sum _I}\frac{1}{2}{M}_I{\dot{\text{R}}}_I^2+{\displaystyle \sum _i}\frac{1}{2}\mu \langle \dot{{\varphi }_i}|\dot{{\varphi }_i}\rangle \\ -E_{\text{KS}}(\text{R},\varphi )-E_{\text{MM}}(\text{R})-E_{\text{QM}/\text{MM}}(\text{R},\varphi )\\ +{\displaystyle \sum _{i,j}}{\Lambda }_{ij}(\langle {\varphi }_i\mid {\varphi }_j\rangle -{\delta }_{ij}),& (1)\end{array}$

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,

$\begin{array}{l}{{L}}_{\text{CP}/\text{Shell}}(\text{R},\dot{\text{R}},\text{r},\dot{\text{r}},\varphi ,\dot{\varphi })={\displaystyle \sum _I}\frac{1}{2}{M}_I{\dot{\text{R}}}_I^2+{\displaystyle \sum _i}\frac{1}{2}\mu \langle \dot{{\varphi }_i}|\dot{{\varphi }_i}\rangle +{\displaystyle \sum _{k\in {P}}^{n_{\text{s}}}}\frac{1}{2}{m}_k{\dot{\text{r}}}_k^2\\ -E_{\text{KS}}(\text{R},\varphi )-E_{\text{MM}}(\text{R},\text{r})-E_{\text{QM}/\text{MM}}(\text{R},\text{r},\varphi )\\ +{\displaystyle \sum _{i,j}}{\Lambda }_{ij}(\langle {\varphi }_i\mid {\varphi }_j\rangle -{\delta }_{ij}),& (2)\end{array}$

where, n_s is the number of shells, and m_k is the fictitious mass of a shell k. Also,

$\begin{array}{ll}{E}_{\text{MM}}(\text{R},\text{r})=E_{\text{b}}(\text{R},\text{r})+E_{\text{nb}}(\text{R},\text{r})+{\displaystyle \sum _k^{n_{\text{s}}}}\frac{1}{2}{\kappa }_k{s}_k^2,& (3)\end{array}$

where

$s_k=\left|{\text{s}}_k\right|=|{\text{R}}_I-{\text{r}}_k|,$

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,

$\begin{array}{l}{E}_{\text{QM}/\text{MM}}(\text{R},\text{r})=E_{\text{b}}^{\prime }(\text{R},\text{r})+E_{\text{vdw}}^{\prime }(\text{R},\text{r})-{\displaystyle \sum _I}{q}_I^{\text{c}}{\displaystyle \int }d\overline{\text{r}}\rho (\overline{\text{r}})\frac{1}{|{\text{R}}_I-\overline{\text{r}}|}\\ -{\displaystyle \sum _k}{q}_k^{\text{s}}{\displaystyle \int }d\overline{\text{r}}\rho (\overline{\text{r}})\frac{1}{|{\text{r}}_k-\overline{\text{r}}|},& (4)\end{array}$

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

$\begin{array}{ll}{\displaystyle \underset{\left\{\varphi ,\text{r}\right\}}{min}}{\displaystyle \underset{T_{\varphi },T_{\text{s}}\to 0}{lim}}{{L}}_{\text{CP}/\text{Shell}}(\text{R},\dot{\text{R}},\text{r},\dot{r},\varphi ,\dot{\varphi })\to {L}(\text{R},\dot{\text{R}}).& (5)\end{array}$

The physical temperature (T_phys) is defined as,

$\begin{array}{ll}{T}_{\text{phys}}=\frac{1}{N_{\text{f}}{k}_{\text{B}}}{\displaystyle \sum _I}{{M}}_I{\dot{\text{S}}}_I^2,& (6)\end{array}$

while the shell temperature is defined as,

$\begin{array}{ll}{T}_{\text{s}}=\frac{1}{3n_{\text{s}}{k}_{\text{B}}}{\displaystyle \sum _k^{n_{\text{s}}}}{\overline{m}}_k{\dot{\text{s}}}_k^2.& (7)\end{array}$

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,

${\text{S}}_I=\frac{1}{{{M}}_I}(M_I{\text{R}}_I+m_k{\text{r}}_k),$

${{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:

${{M}}_I=M_I+m_k$

${\overline{m}}_k=\frac{M_I{m}_k}{M_I+m_k}$

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,

$\begin{array}{l}{{M}}_I{\ddot{\text{S}}}_I={\text{F}}_I^{(\text{S})}-{{M}}_I{\dot{\text{S}}}_I\frac{p_{{\eta }_1}}{q_1}\\ {ṗ}_{{\eta }_1}={\displaystyle \sum _I}{{M}}_I{\dot{\text{S}}}_I^2-N_{\text{f}}{k}_{\text{B}}{T}_{\text{phys}}-\frac{p_{{\eta }_2}}{q_2}{p}_{{\eta }_1}\\ {ṗ}_{{\eta }_j}=\frac{p_{{\eta }_{j-1}}^2}{q_{j-1}}-k_{\text{B}}{T}_{\text{phys}}-\frac{p_{{\eta }_{j+1}}}{q_{j+1}}{p}_{{\eta }_j},j=2,\cdots ,n_{\text{c}}-1,\\ {ṗ}_{{\eta }_j}=\frac{p_{{\eta }_{j-1}}^2}{q_{j-1}}-k_{\text{B}}{T}_{\text{phys}},j=n_{\text{c}},\\ {\dot{\eta }}_j=\frac{p_{{\eta }_j}}{q_j},j=1,\cdots ,n_{\text{c}}.\end{array}$

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,

$\begin{array}{l}{\overline{m}}_k{\ddot{\text{s}}}_k={\text{f}}_k^{(s)}-{\overline{m}}_k{\dot{\text{s}}}_k\frac{p_{{\eta }_1^{*}}}{q_1^{*}}\\ {ṗ}_{{\eta }_1^{*}}={\displaystyle \sum _k}{\overline{m}}_k{\dot{\text{s}}}_k^2-3n_{\text{s}}{k}_{\text{B}}{T}_{\text{s}}-\frac{p_{{\eta }_2^{*}}}{q_2^{*}}{p}_{{\eta }_1^{*}}\\ {ṗ}_{{\eta }_j^{*}}=\frac{p_{{\eta }_{j-1}^{*}}^2}{q_{j-1}^{*}}-k_{\text{B}}{T}_{\text{s}}-\frac{p_{{\eta }_{j+1}^{*}}}{q_{j+1}^{*}}{p}_{{\eta }_j^{*}},j=2,\cdots ,n_{\text{c}}^{*}-1,\\ {ṗ}_{{\eta }_j^{*}}=\frac{p_{{\eta }_{j-1}^{*}}^2}{q_{j-1}^{*}}-k_{\text{B}}{T}_{\text{s}},j=n_{\text{c}}^{*},\\ {\dot{\eta }}_j^{*}=\frac{p_{{\eta }_j}}{q_j^{*}},j=1,2,\cdots ,n_{\text{c}}^{*}.\end{array}$

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,

$\begin{array}{l}{\text{F}}_I^{(S)}={\text{F}}_I+{\text{f}}_k^{(r)},\\ {\text{f}}_I^{(s)}=\frac{1}{{{M}}_I}\left(m_k{\text{F}}_I-M_I{\text{f}}_k^{(r)}\right).\end{array}$

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.

TABLE 1

Table 1. The p–MZHB force field parameters for siliceous zeolites.

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₂O(QM)+4H₂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

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⁻⁷ a.u. atom⁻¹ ps⁻¹, 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

Figure 2. NVE MD simulation using QM/p–MM implementation for 1H₂O (QM) + 4H₂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⁻⁷ a.u. atom⁻¹ ps⁻¹, 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

Figure 3. NVT MD simulation using QM/p–MM implementation for 1H₂O (QM) + 4H₂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₂₀₄₈ O₄₀₉₆) 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₂O₇), 8T (Si₈O₂₅), 14T (Si₁₄O₄₀), and 26T (Si₂₆O₆₇), where T stands for SiO₄ tetrahedral unit (see Figure 4).

FIGURE 4

Figure 4. QM subsystems used in the four QM/MM calculations of oxygen vacancy in α–Cristobalite: (A) 2T (Si₂O₇), (B) 8T (Si₈O₂₅), (C) 14T (Si₁₄O₄₀), and (D) 26T (Si₂₆O₆₇). 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

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

$\begin{array}{ll}\Delta E_{\text{f}}=\frac{x}{2}\left[E({\text{O}}_2)+E_{\text{diss}}({\text{O}}_2)\right]+E({\text{SiO}}_{2-x})-E(\text{Si}{\text{O}}_2)& (8)\end{array}$

where E(O₂), E_diss(O₂), E(SiO_2−x), and, E(SiO₂) are the energies of O₂ molecule (in the triplet electronic ground state), the dissociation energy of O₂ molecule, the energy of bulk silica with oxygen vacancy, and the energy of pure bulk silica, respectively. E(O₂) was computed from “all–QM” calculations whereas E(SiO_2−x) and E(SiO₂) were computed either by “all–QM” or QM/p–MM calculations. E_diss(O₂) = 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

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

TABLE 4

Table 4. ΔE_f computed from single point energy calculations using QM/p–MM potential.

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

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₃H₇ 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₁₅₃₆O₃₀₇₂) 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₃H₇ 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

Figure 6. (A) Mechanism of hydrogenation of ethene catalyzed by Rh₃H₇/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⁻¹ 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⁻¹ 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⁻¹), 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):

$\begin{array}{ll}\text{H}-\text{Z}S\text{M}+{\text{H}}^{\prime }-{\text{CH}}_3\to {\text{H}}^{\prime }-\text{Z}S\text{M}+\text{H}-{\text{CH}}_3.& (9)\end{array}$

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

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⁻¹, 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⁻¹.

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | Google Scholar

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.

Google Scholar

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.

Google Scholar

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

Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Jorgensen, J. D. (1978). Compression mechanisms in α–quartz structures-SiO₂ and GeO₂. J. Appl. Phys. 49, 5473–5478.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Leach, A. R. (2010). Molecular Modelling: Principle and Applications, 2nd Edn. Harlow: Pearson.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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.

Google Scholar

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

Google Scholar

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.

Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | Google Scholar

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

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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.

PubMed Abstract | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Google Scholar

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

Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, X., and Blowers, P. (2005). A computational study of alkane hydrogen-exchange reactions on zeolites. J. Mol. Catal. A Chem. 242, 18–25. doi: 10.1016/j.molcata.2005.07.029

CrossRef Full Text | Google Scholar