MM Boundary Artifacts and Correction in

Article Cite This: J. Chem. Theory Comput. 2019, 15, 3917−3928

pubs.acs.org/JCTC

Quantitative Analysis of QM/MM Boundary Artifacts and Correction in Adaptive QM/MM Simulations Hiroshi C. Watanabe*,†,‡ and Qiang Cui§ †

Quantum Computing Center, Keio University, 3-14-1 Hiyoshi, Kohoku-ku, Yokohama 223-8522, Japan Japan Science and Technology Agency, PRESTO, 4-18, Honcho, Kawaguchi, Saitama 332-0012, Japan § Departments of Chemistry, Physics, and Biomedical Engineering, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, United States Downloaded via NOTTINGHAM TRENT UNIV on August 6, 2019 at 12:52:44 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

S Supporting Information *

ABSTRACT: A quantum chemical treatment of solvation effects using the standard quantum mechanical/molecular mechanical (QM/MM) molecular dynamics simulations is challenging due to QM and MM solvent exchange near the QM solute. To this end, various adaptive QM/MM methods have been proposed; free solvent exchanges are allowed via flexible switching of their identities between QM and MM systems depending on their distances from the QM solute. However, temporal and spatial discontinuities remain in the standard implementations of adaptive QM/MM approaches and continue to hamper stable and accurate dynamics simulations. We previously demonstrated that the size-consistent multipartitioning (SCMP) method achieves temporal continuity while, to some extent, avoiding spatial discontinuities. In the present study, we demonstrate that residual spatial discontinuities may lead to severe artifacts under certain conditions. Through quantitative analyses, we show that all multiscale including QM/MM methods might experience these problems, which so far have not been investigated in depth. To alleviate these artifacts, we propose a correction scheme in the framework of the SCMP approach and demonstrate its effectiveness using bulk water simulations.

■

INTRODUCTION Quantum chemical effects of solvation have a significant impact on various chemical processes in solution and biomolecules. To capture these effects in molecular simulations, quantum mechanical (QM) models that explicitly consider the electronic structure of the solution are required. However, due to the high computational cost of QM calculations, it is typically not feasible to treat the entire solution system using QM models. Accordingly, hybrid computational models that combine QM and molecular mechanical (MM) potentials have been widely applied, in which the region of central interest is treated by QM, while the rest is described by MM.1−4 In general, there are two possible choices for modeling the solution system with respect to QM/MM partitionings: (i) only the solute molecule of interest is treated by QM; (ii) the solute molecule and the nearby solvent molecules are treated by QM (Figure 1). Intuitively, the second case appears to be more accurate because it can capture quantum chemical interactions between the solute and nearby solvent. However, there remains an important technical issue for molecular dynamics (MD) simulations with the standard QM/MM implementations, in which the partitioning of the system into QM and MM regions remains fixed during the MD simulation. Solvent molecules treated by QM, which are usually in the vicinity of the QM solute molecule at the beginning of the simulation, diffuse away in the course of the MD simulations and are replaced by MM solvent, leading to a less accurate treatment of solute−solvent interactions. Therefore, simulation frameworks that maintain the QM solvent © 2019 American Chemical Society

Figure 1. Two possible QM/MM partitionings for a solution simulation. QM and MM molecules are represented van der Waals and ball/stick models, respectively. (left) The solute molecule of interest is modeled by QM. (right) Both the solute molecule of interest and the nearby solvent molecules are modeled by QM.

near the QM solute, while preserving the validity of the simulation, are required. A naive treatment of the solvent exchange problem is to update the identity of a solvent as QM or MM as soon as it enters or leaves the region near the solute during the MD simulation, which is termed abrupt partitioning Received: February 21, 2019 Published: May 16, 2019 3917

DOI: 10.1021/acs.jctc.9b00180 J. Chem. Theory Comput. 2019, 15, 3917−3928

Article

Journal of Chemical Theory and Computation

M number of possible QM/MM partitionings is ∑m ( M m ), where M and m denote the total number of solvent molecules in the system and the number of QM solvent molecules in a partitioning, respectively. Due to the limited computational resources, it is not feasible to treat all the partitionings. Thus, only representative partitionings are considered in practice, and partitionings whose QM region is fragmented are replaced by those whose QM region is compact during the MD simulation (partitioning update). However, abrupt updates lead to temporal discontinuities in Veff. Thus, the weight function σ(n) is designed to suppress discontinuities in Veff such that the partitionings included and excluded during the update have a weight of zero. The continuity of the potential can be ensured via the conservation of the Hamiltonian, which can be achieved through the SAP, PAP, or SCMP methods. However, note that the PAP method is not practical because it has to count all possible partitionings with respect to the solvent molecules in the transition region. Therefore, in practice, the PAP is replaced by a modified PAP (mPAP) method, where the potential energy expansion of PAP is truncated at the second or third order to limit the number of partitionings considered, albeit the potential energy is no longer continuous in such case.17 The third aspect of adaptive QM/MM methods is further divided into two issues: (i) the formulation of the Hamiltonian and (ii) the artifacts from the QM/MM spatial boundary. Although the first issue has been discussed extensively in the literature, the second issue has rarely been discussed because it is assumed that such artifacts can be largely suppressed via smooth switching of the QM vs MM identity of a solvent. In this study, we demonstrate that, while this assumption is partially correct, the residual artifacts are not negligible. Moreover, we discuss how these two issues are rooted in the same cause and are closely intertwined. To better explain these points, it is useful to briefly review the formulation of potential-oriented and force-oriented approaches. In potential-oriented approaches, the Hamiltonian H is written as

update. However, this abrupt partitioning update leads to spatial and temporal discontinuities. The spatial discontinuity is caused by the artificial boundary between the QM and MM regions in the solvation structure, as reflected by the solvent−solute radial distribution functions (RDFs). This occurs because, although the QM and MM methods describe the same molecular species, the properties of the solvent at the two levels are not identical. The temporal discontinuity is caused by the abrupt change in the potential energy before and after the QM/MM partitioning update, which leads to monotonic temperature increase (40− 900 K/ps in a previous report5) and instability of the MD simulation. To deal with such problems, several research groups have proposed revised QM/MM techniques over the past two decades, which can be broadly classified into “constrained QM/ MM” and “adaptive QM/MM”. Most constrained QM/MM techniques are based on a single QM/MM partitioning scheme, such as FIRES,6 BEST,7 and BCC.8 In these schemes, the boundary between the QM and MM systems is regarded as closed because the bias potentials are introduced to prevent solvent exchange between the two regions and the QM vs MM identities of the solvent molecules remain fixed during the simulations. Owing to the constraining forces, the resulting dynamics are no longer realistic; instead, these methods focus on reproducing equilibrium properties. For instance, they cannot be used to treat reaction or diffusion dynamics, but equilibrium properties such as the solvation structure can be evaluated. In contrast, adaptive QM/MM methods are open-boundary approaches and allow exchange of solvent molecules between the QM and MM regions; therefore, they aim at computing not only equilibrium properties but also realistic dynamics. To this end, the molecular properties of solvents have to be smoothly switched between QM and MM models depending on their distance from the QM solute. Most adaptive QM/MM methods are based on multipartitioning schemes, whereby various QM/ MM partitionings that have different combinations of QM solvents are defined. In these methods, the effective potential energy, Veff, is defined as eff

V (q) =

∑σ n

(n)

H=

(n)

(q)V (q)

∑ j

(1)

pj 2 2mj

+ V eff (2)

where mj and pj represent the mass and the moment, respectively, of the jth particle. Then, the effective force, Feff, used to propagate the MD trajectory is,

where V(n) and σ(n) are the QM/MM potential energy and the normalized weight function, respectively, of the nth partitioning and q indicates the atomic coordinates. The differences between the various adaptive QM/MM methods lie in the definition of the partitionings and the specific form of the weight function. To date, different adaptive QM/MM methods have been proposed, each having respective advantages and limitations.9−16 To compare these methods, there are three important aspects worth considering: (1) the existence of a Hamiltonian, (2) the differentiability of Veff, and (3) the physical validity of the resulting dynamics and ensemble. The first aspect can be paraphrased as whether the equations of motion are properly derived from the Hamiltonian. The dynamics based on the Hamiltonian ought to be symplectic, which ensures the time reversibility and stability of the MD simulation. Note that, as described below, this aspect does not depend on whether the Hamiltonian is potential-oriented or force-oriented. To our knowledge, the ONIOM-XS,9 PAP,10 SAP,10 DAS,11 and SCMP15 methods satisfy this requirement. The second aspect, that is, the continuity of the potential, is mainly related to the weight function, σ(n), in eq 1. In general, the

F jeff = −

∂V eff ∂qj

N

N

=∑ σ (n)(q)f j(n) − n

= Fj +

∑ ∂σ n

Fjtrans

(n)

∂q

V (n)(q) (3)

(n)

where f is the force in the nth QM/MM partitioning, given by the negative derivative of V(n) with respect to q(j). The second term, Ftrans, sometimes referred to as the transition force, is not physically realistic because it arises from the derivative of the artificial weight function σ(n). Because the transition force can lead to unrealistic dynamics and solvation structures,6 in potential-oriented approaches it is important to minimize Ftrans. Note that the contribution of Ftrans to Feff is significant in magnitude in multisize partitioning approaches, such as PAP, SAP, and DAS, in which the partitionings have different 3918


Article


properties,23,24 which are important challenges in popular DFT methods.25−27 In a recent investigation of the sensitivity of the bulk water properties on the DFTB3 parameters, Goyal and coworkers employed an inverse Monte Carlo approach to adjust the O−H repulsive potential and reproduce the radial distribution function of bulk water under ambient conditions. They observed that the resulting 3OBw model, which features a small perturbation (on the order of kBT) in the O−H repulsive potential, was able to considerably improve the structural properties of bulk water; some dynamical properties (diffusion and orientational relaxation) and the infrared spectrum of water were also described reasonably well. However, it was noted that the 3OBw model reduced the degree of oversolvation of 3OB at the cost of considerable errors in the energetic properties, such as the enthalpy of vaporization, suggesting that a transferable DFTB model requires improvements of the electronic components of DFTB3; the implication of this is discussed below, where we compare DFTB3/MM simulations with 3OB and 3OBw models. In the next section, we also demonstrate that artifacts arise from the QM/MM system boundary through a quantitative analysis. Further, we show that, although previous studies assumed that the QM/MM boundary effect is effectively removed in adaptive QM/MM methods, considerable artifacts can still remain. To suppress the bookkeeping term beyond the size-consistency requirements, further improvements in existing methods are required. We propose a simple scheme to reduce such errors and improve the performance of adaptive QM/MM method-based simulations.

numbers of QM solvent molecules. Because the potential energy of a QM model is usually much more negative than that of an MM model, the QM/MM energy V(n) becomes increasingly negative as the number of QM solvent molecules increases. Because Feff should work to minimize Veff, the solvation structure can be distorted such that partitionings with larger QM regions have a larger value of σ(n). In contrast, because all partitionings have the same number of QM solvent molecules in the SCMP method, the potential energy differences among partitionings are much smaller than those in multisize partitioning approaches and, correspondingly, the artifacts are smaller. In force-oriented approaches, the Hamiltonian H is written as H=

pj 2

∑

2mj

j

−

∑∫ n

+

∑ σ(n)(q)V (n)(q) n (n)

∂σ V (n)(q) dq ∂q

=K (p) + V eff (q) + B(q)

(4)

where K(p) is the kinetic energy and B(q) is a bookkeeping term, originally proposed by Ensing et al.18 for atomistic/coarsegrained study and introduced to adaptive QM/MM models by Bulo et al.11 Because the derivative of the bookkeeping term cancels out the transition force, the effective force is N

F jeff =

∑ σ(n)(q)f j(n) (q) n

■

(5)

SIZE-CONSISTENT MULTIPARTITIONING QM/MM METHOD Suppose a system is partitioned into two parts: a QM region consisting of one solute and m solvent molecules and an MM region containing M solvent molecules. While the number of possible QM/MM partitionings for an MD snapshot is very large, the SCMP method considers only a limited number of partitionings in which the QM region has the same number of solvent molecules in different combinations. Then, the force f and the potential energy V are evaluated via an ordinary QM/ MM scheme for the respective partitioning, and the effective force Feff used for coordinate update in the MD simulation is evaluated as a weighted sum of forces f following eq 5. Because the partitionings in the disordered QM region are replaced by those in the ordered QM region in the SCMP method (Scheme 1), the weight functions σ(n) for the nth partitioning have to satisfy the two following conditions for the temporal continuity: (i) σ(n) = 0 when the QM region of the nth partitioning is disordered; (ii) σ(n) = 0 when the QM region of the nth partitioning is perfectly ordered (Scheme 1). To this end, we introduce the fade-out function, O(n), for the nth partitioning:

Unlike potential-oriented approaches, force-oriented approaches do not exhibit transition forces due to the weight function and yield more realistic dynamics. The bookkeeping term can decrease monotonically over the course of the MD simulation. Correspondingly, the internal energy, which consists of K and Veff, increases monotonically while the Hamiltonian is conserved. As a result, the trajectory does not yield an accurate ensemble regardless of whether the temperature is held constant with thermostat. To avoid this, the bookkeeping term should be kept constant. It is worth noting that, in potential-oriented approaches, the bookkeeping term is defined as an integral of the transition force, Ftrans. Because the SCMP method yields a smaller transition force, the drift of the bookkeeping term in SCMP should be also smaller in magnitude than that of multisize partitioning methods. Indeed, previous studies demonstrated that temperature drifts caused by the bookkeeping term obtained in SCMP simulations were less than 0.06 K/ps.15,19 This is in contrast to multisize methods, which caused temperature drifts of more than 100 K/ps,5 although these may partially arise from potential discontinuities.5 However, in the Results section, we demonstrate that, even with the SCMP method, the bookkeeping term can exhibit a large drift depending on the molecular model employed; this implies that other factors can cause significant drifts in the bookkeeping term. Throughout this study, we employ the approximate density functional tight binding method, DFTB3,20−22 as the QM model in the SCMP method. DFTB3 is an approximation of density functional theory (DFT) via expansion of the DFT energy up to the third order.21 Although DFTB3 is an attractive choice for QM−MD simulations owing to its computational efficiency, it causes well-documented errors in the description of hydrogenbonding interactions involving water,21 especially for bulk

m (n) OQM =

∏ λj(n)(q) j

(6)

where λ(n) is a progress function of the jth QM solvent j molecules, ranging from 0 to 1 depending on the distance between the QM solute and the jth solvent molecule. In addition, the fade-in function, I(n), is defined as (n) (n) IQM = 1 − OQM

3919

(7) DOI: 10.1021/acs.jctc.9b00180 J. Chem. Theory Comput. 2019, 15, 3917−3928

Article

Journal of Chemical Theory and Computation Scheme 1. Weight Function, σ, and Partitioning Update in the SCMP Methoda

a

When the QM region in a partitioning is completely ordered or disordered, σ will become 0. A partitioning whose QM region is disordered is replaced by one whose QM region is ordered.

We also introduce the fade-out and fade-in functions for the MM solvent molecules as M (n) OMM =

∏ Λ(jn)(q)

(n) (n) IMM = 1 − OMM

j

(8)

Λ(n) j

where is a progress function of the jth MM solvent molecule. The specific forms of the progress functions can be found in previous works.15,28 The weight function σ for the nth partitioning is σ (n) =

■

(n) (n) (n) (n) OQM IQMOMMIMM N

(k) (k) (k) (k) ∑k OQM IQMOMMIMM

(9)

RESULTS The Bookkeeping Term. Figure 2a shows that the Hamiltonian is numerically conserved in the SCMP simulations with DFTB3/3OBw, which is consistent with previous studies that used DFTB3/3OB.15,19 This highlights the continuity of the potential function in the SCMP method regardless of the DFTB parameters. On the other hand, the bookkeeping terms in Figure 2b show monotonic drifts over the course of the MD simulations with both DFTB parameters. However, note that the magnitude of the drift depends significantly on the DFTB parameters, which are −8.5 kJ/(mol·ps) for the 3OB and −39.2 kJ/(mol·ps) for 3OBw.29 Because the Hamiltonian is conserved in SCMP, the drift of the bookkeeping term is coupled to the increase in the internal energy, which consists of kinetic and effective potential energies. Indeed, the simulation with the 3OBw parameters shows a temperature drift of 0.18 K/ps under microcanonical (NVE) conditions, whereas the 3OB parameters result in a much smaller temperature drift of 0.04 K/ps, in agreement with previous studies.15,19,30 As expected, this increase is consistently observed also in the effective potential (data not shown). We note that the continuous drifts of the bookkeeping term imply systematic differences in the potential energy between ordered and disordered partitionings. Suppose that the nth partitioning is perfectly ordered, satisfying σ(n) = 0 at the MD time point t1. Over the course of the MD simulation, the partitioning is bound to become disordered due to solvent diffusion, leading to σ(n) > 0. Eventually, it will be completely disordered at t = t2, and σ(n) will become 0 again. The contribution of the nth partitioning to the bookkeeping term for the time period between t1 and t2 is

Figure 2. Time evolution of Hamiltonian (a), temperature (b), and the bookkeeping term (c) over the course of the MD simulations. The black and red lines represent SCMP simulations using DFTB3/3OB and DFTB3/3OBw as the QM model, respectively.

ΔB(n) = −

∫q

q2

1

∂σ (n) (n) V (q) dq ∂q

(10)

where q1 = {q(t1)} and q2 = {q(t2)}. Therefore, if ΔV(n) = V(n)(q2) − V(n)(q1) ≠ 0, then ΔB(n) ≠ 0. Because the bookkeeping term is the sum of ΔB(n) over all partitionings, a monotonic decrease indicates that the QM/MM potential energies statistically satisfy ⟨ΔV(n)⟩ < 0. Furthermore, in the SCMP simulations, because the disordered partitionings are updated with an ordered partitioning, a nonzero ΔB accumulates continuously, resulting in the constant drift of the 3920


Article


partitioning has statistically higher potential energy than a disordered one implies that there exist driving forces that promote solvent diffusion. In other words, on average, a systematic set of repulsive forces act on the QM solvent molecules from the QM center or attractive forces act on the MM solvent. Because a homogeneous system modeled with a single molecular model does not have such driving forces, they are clearly artifacts produced by the spatial discontinuity associated with the QM/MM partitioning. Potentials of the Mean Force. To gain a deeper understanding of these artifacts, we calculate the potentials of the mean forces (PMFs) acting on the oxygen and hydrogen atoms in the QM and MM water molecules, as shown in Figure 4. These PMFs are evaluated with f(n) (eq 5) via conventional QM/MM calculations in a postprocessing step following the SCMP simulations. Because the PMF represents an averaged influence, including direct and indirect effects, on the solvent from the QM solute oxygen, it should attenuate and reach a plateau as the distance increases. However, this is not the case

bookkeeping term. Thus, the present results suggest that, for an arbitrary time period, the averaged potential energy difference ⟨ΔV(n)⟩ for 3OBw is greater in magnitude than that for 3OB. Figure 3 shows the potential energy curves of the water dimers as a function of the distance between the oxygen atoms. A

Figure 3. Potential energy curves for a water dimer in vacuo. The black, red, green, and blue lines represent homodimers modeled with the CAM-B3LYP with 6-311++G(3fd, 2pd) basis set, SPC-Fw, DFTB3/ 3OB, and DFTB/3OBw models, respectively. The violet line represents a heterodimer in which hydrogen binding donors and acceptors were modeled with DFTB3/3OBw and SPC-Fw, respectively.

reference calculation using the CAM-B3LYP functional with 6311++G(3df, 3pd) basis set yields a binding energy of approximately 23.1 kJ/mol at a distance of 2.88 Å from the oxygen atom. The homodimer modeled by the SPC-Fw31 overestimates the minimum depth as 28.8 kJ/mol at 2.73 Å, whereas DFTB3/3OB and DFTB3/3OBw show weaker hydrogen bonds, with minima of approximately 19.3 and 16.8 kJ/mol, respectively. In addition, the heterodimer pairs of QM and MM waters yield the same depth of 23.4 kJ/mol regardless of whether 3OB or 3OBw is used, where the QM and MM models are proton acceptor and donor, respectively. Note that the potential energy curve is consistent with an interchange between donor and acceptor in a modest distance range (data not shown). Therefore, there exist three types of hydrogen bonds in the present SCMP simulations, namely, QM−QM, QM−MM, and MM−MM, and the energy differences between these three types of hydrogen bonds are more distinct for DFTB3/3OBw than for DFTB3/3OB. Suppose that a partitioning divides a homogeneous system into QM and MM regions modeled by methods A and B, respectively. If models A and B are identical, the QM/ MM potential energy V(n) is independent of the shape and size of the QM region. In contrast, as the energy difference between A and B becomes large, the potential energy becomes more dependent on the shape of the QM region. Thus, considering the hydrogen bonding energies shown in Figure 3, the QM/MM potential energy V(n) obtained with the 3OBw parameters varies depending on the shape of the QM region more than that obtained with the 3OB parameters, implying that ΔV3OBw < ΔV3OB < 0. Therefore, the hydrogen bonding energy shown in Figure 3 is in accordance with the observation, shown in Figure 2b, that the 3OBw parameters result in a larger drift in the bookkeeping term. Furthermore, the fact that an ordered

Figure 4. Potential of the mean force of oxygen (top) and hydrogen (bottom) versus distance from the QM center. The obtained PMFs for the QM and MM solvent molecules are represented by dashed and dotted lines, respectively. The resulting PMF curves obtained using the 3OB and 3OBw parameters are colored in red and black, respectively. The effective PMFs of oxygen, represented by solid lines, are estimated by combining the QM and MM PMFs using the QM profile (see Figure 5). 3921


Article

Journal of Chemical Theory and Computation for the computed PMFs of the QM solvent molecules. On the contrary, the PMF slopes become steeper at longer distances for both oxygen and hydrogen atoms. It is clear that these unexpected characteristics arise from the artificial transition between the QM and MM regions. In addition, both oxygen and hydrogen PMFs exhibit a distinct minimum at approximately 6− 7 Å. This implies that, on average, the QM solvent molecules, are subject to repulsive forces from the QM center toward the minimum around the QM/MM boundary, which features stronger water−water interactions (see Figure 3). In contrast to the QM solvent, the MM solvent showed a flat PMF at the long-distance limit. Instead, the short-range behavior seemed to be distorted. As shown in Figure 4, the PMFs of MM oxygen and hydrogen have a maximum at the positions that correspond to the minima of the QM oxygen and hydrogen, respectively, although their positions are slightly shifted to the shorter side. Note that both PMFs drop abruptly from the peak to the shorter side, which indicates a significant attraction to the QM solute. In particular, the hydrogen can freely approach the QM center (QM solute oxygen), whereas the oxygen presumably has a minimum at approximately 4.5 Å (Supporting Information). This indicates that the MM solvent molecules are more stable in the vicinity of the QM center than the QM solvent molecules, thus driving the MM solvent molecules toward the QM region. Overall, the resulting PMFs showed that the artifacts on the QM and MM solvent molecules emerge as driving forces that promote solvent exchange. In other words, the obtained PMFs imply that ⟨ΔV(n)⟩ < 0, which is in accordance with the observed monotonic decrease of the bookkeeping term shown in Figure 2b. Note that these artifacts are not limited to SCMP simulations but are a general consequence of the QM/MM setup because the resulting PMFs are evaluated from conventional QM/MM calculations using f(n) in eq 5 and, therefore, are free from SCMP force compositions. Thus, we expect that most QM/MM calculations are subject to these artifacts, although they presumably emerge more distinctively in dynamics simulations. Moreover, note that these artifacts are alleviated to some extent in the SCMP method. To capture this, the QM profile ω that we previously proposed15 is a useful metric and is defined as

Figure 5. QM profiles versus distance from the QM center obtained using SCMP (for the oxygen of the QM solute water). The red and black lines represent the results obtained with the 3OB and 3OBw parameters, respectively.

effective PMFs for 3OB and 3OBw in Figure 4 indicate that QM characteristics do not contribute to the solvation dynamics around ω = 0. Likewise, the resulting dynamics are free from the artifacts seen in the PMFs of the MM solvents around ω = 1. Note that the effective PMF of oxygen, shown in Figure 4a, decreases by approximately 1.5 kJ/mol in the range from 6.0 to 8.0 Å, indicating that the solvent is subject to repulsive forces from the QM solute. This leads to an underestimation of the RDF, as discussed below, which indicates that QM/MM boundary artifacts cannot be completely removed in the SCMP method. Because the SCMP simulations with the DFTB3/3OB parameters result in a smaller drift of the bookkeeping term, it is insightful to compare the resulting PMF curves with those obtained using the 3OBw parameters. While the hydrogen PMFs appear to be similar in shape to those of the 3OBw simulations, the oxygen PMFs show distinct differences; for example, in the transition region, the QM PMFs increase as the distance increases while the MM PMFs decrease. Although these behaviors are not realistic, the resulting curve appears to be symmetric around 6 Å in the transition region as well as in the QM profile, and thus, these artifacts are canceled out in the effective force ∑σ(n)f(n), presumably yielding flat PMFs in the analyzed range. It should be noted that we employed different transition ranges described by parameters s and t (details can be found in the literature9,15) for the 3OB and 3OBw simulations because they yield different densities around the QM regions, and the proper transition range for efficient simulations differs accordingly, as discussed in the literature.28 To verify the influence of the transition range on the PMF, we have carried out additional SCMP simulations using different values of s and t (Supporting Information) and confirmed that, although the transition range can affect the maximum and minimum positions of the PMFs to some extent, it does not cause changes as drastic as those seen for different DFTB parameters. Thus, we conclude that the significant PMF differences shown in Figure 4 mainly arise from differences in the QM models rather than from the transition parameters and that the 3OB parameters yield better results in QM/MM simulations because of the stronger QM− QM hydrogen bonds and error cancellations.

N

ωj =

∑ σ(n)δj(n) n

(11)

where the subscript j represents the jth nearest solvent from the QM center. l o1 if j ∈ S(n) δj(n) = o m o o 0 otherwise n

(12)

and S(n) is a subset of the QM solvent in the nth partitioning. The QM profile indicates the extent to which the respective solvent molecules behave quantum mechanically and its values range from 0 to 1, where 0 and 1 indicate that the corresponding solvent molecule behaves as an MM and a QM solvent molecule, respectively. As shown in Figure 5, the QM profile dampens smoothly to zero as the distance to the QM center increases from 5.0 to 8.5 Å. Note that one can roughly estimate an effective solvent PMF using the QM profile ω(j) as the mixing ratio of the PMFs of the QM and MM solvents. Although the exact PMF can be obtained with the radial distribution function (RDF), the use of the QM profile provides useful insights into the artifacts via the solvent characteristics in the transition region. The obtained 3922


Article

Journal of Chemical Theory and Computation Radial Distribution Function (RDF). The RDF provides an alternative way to analyze the effective PMF in SCMP simulations and evaluate the error suppression. First, Figure 6

⟨F eff (r )⟩ = kBT

d ln g (r ) dr

(13)

we can evaluate the deviation of the PMFs based on the SCMP effective force Feff from experimentally determined ones (Supporting Information). The PMFs of both oxygen and hydrogen atoms exhibited deviations of at most 1.0 kJ/mol in a distance range where the QM profile ω is 0 < ω < 1. Artifact Correction. In this subsection, we propose a correction to alleviate the artifacts arising from the QM/MM boundary, namely, the residual spatial discontinuity. To this end, we introduce additional pairwise potentials, U(n)(q), which describe interactions between the QM solute and solvent molecules. Then, we replace the respective QM/MM potential V(n) in eq 4 by V′(n) = V(n)(q) + U(n)(q). Accordingly, the Hamiltonian is written as H=

∑ j

−

pj 2 2mj

∑∫ n

+

∑ σ(n)(q)V ′ (n)(q) n (n)

∂σ V ′ (n)(q) dq ∂q

(14)

(n)

because the derivative of U in the bookkeeping term is canceled out. The corrective force Fcorr acting on the ith solvent i molecule is F jcorr = −∑ σ (n) n

∂U (n) =∑ σ (n)ℑ(jn) ∂qj n

(15)

(n)

Note that U is arbitrary except for three requirements: (A) correction of the solvation structure, (B) suppression of the bookkeeping term drift, and (C) continuity in the region where it applies. In addition, note that, if necessary, U (n) can work simultaneously on the QM and MM solvents. Moreover, the corrective potential U(n) is allowed to have different forms, namely, UQM and UMM, for the QM and MM solvents, respectively. Both UQM and UMM are preferably monotonic, and the difference between the corrective potential as solvent exchange progresses should satisfy ΔU(n) = U(n)(q2) − U(n)(q1) < 0 to suppress the drift of the bookkeeping term. In other words, F(n) acting on the QM solvent should be attractive toward the QM center, whereas F(n) acting on the MM solvent should be repulsive. Because of requirement C, UQM should be continuous where r < tQM, whereas UMM should be continuous where r > sMM, with tQM and sMM being transition parameters that define the progress functions λQM and λMM, respectively (details can be found in the literature15,28). Here, the underestimation of the RDFs when using the 3OBw parameters implies that the artifacts from the QM/MM boundary lead to an overall repulsive force. Thus, to satisfy requirements A and B, we chose F(n) to be an attractive force acting on the QM solvent molecules. Note that the progress function λ, which ranges between 0 and 1, satisfies the monotonic and continuous requirements stated above. In the present study, we define UQM as

Figure 6. RDFs of oxygen (top) and hydrogen (bottom) versus distance from the oxygen of the QM solute (a water molecule in the center). The black and red lines represent the results for the SCMP simulations using DFTB3/3OBw and DFTB3/3OB, respectively. The full QM simulation using DFTB3/3OBw is represented in magenta. The green line represents the experimental results.32

shows that the first peak of the SCMP simulation with the 3OB parameters is higher than the experimental RDF, indicating oversolvation, and that the position of the second peak is shifted toward the long-distance side; this is consistent with previous full QM simulation.19,28,30 Next, as reported in a previous study,32 the RDF obtained from the full QM/MD simulation with DFTB3/3OBw shows good agreement with experimental results,29 although the heights of the first peaks of both oxygen and hydrogen RDFs are slightly overestimated. Nevertheless, the use of the 3OBw parameters in the SCMP approach resulted in an overall underestimation in the short-distance range compared with the full QM simulation, although the peak positions are in good agreement. The sparse QM center seems to be in accordance with the estimated effective PMF, as shown in Figure 4, which implies that artificial forces statistically promote solvent exchange. Because the RDFs are linked to the PMF as

(n) UQM (q) = U ·λ(q)

(16)

where U is a constant parameter that corresponds to the height of the additional potential barrier that the QM solvent molecules have to overcome during solvent exchange. Although the progress functions λQM and λMM are not unique, in the present study they are described by three spline curves, as proposed in a 3923


Article

Journal of Chemical Theory and Computation previous study.28 Likewise, λ in eq 16 consists of the three spline curves defined in the range sMM < r < tQM (details can be found in the Supporting Information). Benchmark Simulations. We evaluate the time evolution of the bookkeeping term with various values of U in eq 16 by applying the same value to oxygen and hydrogen. First, we confirm that the Hamiltonian is still numerically conserved, as shown in Figure 7a, verifying the continuity of U. Next, Figure 7b shows that the drift of the bookkeeping term depends largely on the value of U. As U becomes larger, the drift becomes smaller. When U is 2.5 kJ/mol, the bookkeeping term remains almost zero during the MD simulation, although it fluctuates around the average value, which is even smaller than that obtained with the 3OB parameters (Figure 2). This result indicates that, with a properly chosen U, the drift of the bookkeeping term can be drastically suppressed. Correspondingly, the temperature is kept almost constant, as shown in Figure 7c. Therefore, the resulting trajectories represent an accurate ensemble under microcanonical conditions. Finally, we compare the RDFs computed with the correction term with the original SCMP results. As shown in Figure 8, the correction leads to an increase in the height of the peaks in both Oc-O and Oc-H RDFs, whereas the peak positions show no distinct differences. As expected, the value of U dictates the increase in the RDF peak heights. It is also worth noting that the best U value for the suppression of the drift in the bookkeeping term, namely, U = 2.5 kJ/mol, also leads to the best agreement in the RDFs with those obtained from the full QM calculations. Therefore, the results confirm that the artifacts in the QM/MM potential V(n) are manifested in both distortion of the solvation structure and drift in the bookkeeping term. On the other hand, the first minimum of both oxygen and hydrogen RDFs are still underestimated, which suggests that the present correction cannot remove all artifacts from the QM/MM boundary.

■

DISCUSSION The proposed correction approach is similar to the HAMBC, which was proposed by Boereboom et al.33 using SAP and implemented into mPAP.34 However, the basic concepts and their practical implementation are different. In the HAMBC method, an additional potential was included in the potentialoriented Hamiltonian of a multisize partitioning method as H=

pj 2

∑ j

−

2mj ij

+

∑ σ(n)(q)V (n)(q) yz

n

∑ jjjjjσ(n)(q) ∑ δj(n)C zzzzz n

j k

z {

j

(17)

where δ(n) j equals unity when the jth solvent is QM in the nth partitioning and C is a constant parameter, which is usually chosen depending on the energy differences in the solvent between the QM and MM models in the gas phase. The effective forces can be described as f jeff =

∑ σ(n)f j(n)

yz − VMM(q)zzzz z { n

−

n

(n) i

jj (n) jjV − ∂qj jj k

∑ ∂σ

Figure 7. Time evolution of the Hamiltonian (a), temperature (b), and the bookkeeping term (c) over the course of the SCMP simulations with a correction potential (eq 16). The green, red, and black lines represent the results for the SCMP simulations using U = 1.0, 2.0, and 2.5 kJ/mol, respectively. The bookkeeping term and temperature obtained with the 3OBw parameters and U = 0.0 kJ/mol are also presented (blue line) as a reference.

∑ δk(n)C k

where VMM is a potential energy fully described by the MM model. Note that eqs 17 and 18 contain a sum of the constant value C, which is thus proportional to the QM size in a partitioning. Therefore, the HAMBC approach is based on the

(18) 3924


Article


models. After comparing the results for different DFTB parameters, it was found that the original DFTB3/3OB parameters can be used to carry out reliable SCMP simulations for at least several hundred picoseconds, as previously demonstrated.19,30 On the other hand, the 3OBw parameters lead to more severe structural discontinuities in SCMP simulations, although they reproduce more accurate solvation structures when employed in full QM simulations. Therefore, the use of other QM models that have a more balanced hydrogen bonding energy compared to MM models, likely results in more stable and accurate QM/MM simulations. This conclusion is consistent with the recent work of Jiang et al.,35 who proposed a tailored MM modeling approach by tuning the MM parameters to reproduce the energetics of a QM model and demonstrated the excellent agreement between their RDF results with experimental data. Although their tailored MM approach was based on the single-partitioning method FIRES, we expect that it is also effective for adaptive QM/MM methods. However, tailoring MM models seems to require a tremendous amount of effort, in particular for inhomogeneous systems. As an alternative, including a correction potential to balance interactions across the QM/MM boundary as we have shown here is more straightforward. We note that hydrogen atoms in most MM water models, such as SPC-Fw, do not have LennardJones parameters; including explicit Lennard-Jones interactions for hydrogens may further improve the consistency across the QM/MM interface.

■

CONCLUDING REMARKS In conclusion, we conducted a quantitative analysis of the spatial discontinuities arising at the QM/MM boundary, which distort the solvation structure as shown by the RDFs. Although these artifacts are more critical to dynamics simulations than to static ones, the present results indicated that all QM/MM simulations are potentially subject to these artifacts, which has not been discussed in depth in previous works. In addition, the present study also demonstrated that, although the obtained dynamics are free from the transition forces present in adaptive QM/MM methods, the bookkeeping term can drift as a result of these boundary artifacts, which can be linked to temperature drifts in the resulting Hamiltonian dynamics. Thus, the artifacts not only destabilize the simulation but also prevent the methods from yielding an accurate statistical ensemble. These artifacts presumably originate from energy imbalances between the QM and MM solvent models, in particular with respect to hydrogen bonding interactions between water molecules. Therefore, these artifacts are more evident when using the DFTB3/3OBw parameters in combination with SPC-Fw than when using the DFTB3/3OB parameters because of the large energy differences between the DFTB3/3OBw and MM water interactions. Furthermore, these artifacts can cause large errors in the RDF in DFTB3/3OBw/MM simulations. To alleviate the boundary artifacts, we proposed a correction potential that can balance the QM and MM interactions across the QM/MM boundary better, leading to the suppression of artificial solvent exchanges; this is similar in spirit to the adjustment of van der Waals potentials between QM and MM components to balance QM−QM and MM−MM interactions.36,37 The correction was observed to significantly suppress the drift of the bookkeeping term and improve the RDFs in DFTB3/3OBw/MM simulations. However, it should be noted that the simple correction adopted here does not remove all spatial discontinuities. Moreover, although there is increasing

Figure 8. RDFs of oxygen (top) and hydrogen (bottom) obtained from the SCMP simulations with a correction potential. The blue, red, green, and black lines represent the SCMP simulation results using U = 0.0, 1.0, 2.0, and 2.5 kJ/mol, respectively. The results of the full QM simulation with DFTB3/3OBw are also presented (magenta line).

same concept as the SCMP method for minimizing the transition force via the size-consistency of the partitionings. However, the SCMP method is advantageous because it ensures size-consistency without any parametrization. In other words, the SCMP method can remove most of the artifacts caused by the derivative of σ(n). Therefore, the SCMP method can attribute the drift of the bookkeeping term purely to the QM/MM boundary artifacts in V(n), which allows detailed analyses and performing a simple correction. In addition, although in this work we apply the correction potential only to the QM solvent, it is also possible to apply different corrections to the QM and MM solvents and to oxygen and hydrogen, respectively. However, this also means that the choice of correction potentials is not unique. Further, the form of a suitable correction function can vary depending on the system. Because we found the proper conditions by changing the value of U to suppress the drift of the bookkeeping term, the present approach has the disadvantage of requiring preliminary calculations for accurate production runs. The present study also demonstrates a clear guideline for better spatial continuity because the energy differences between the QM and MM models must be minimized. The DFTB models (DFTB3/3OB and 3OBw) used in the present study underestimate the hydrogen bonding energy more than the MM 3925



■

demand for multiscale simulations in practical applications, the proper form of the correction potential is not obvious for anisotropic systems, such as gas−liquid interfaces and biomolecular surfaces or binding pockets.38−41 Considerable efforts have been made to improve approximate QM methods to achieve a good trade-off with the computational efficiency.29,42−45 Moreover, significant efforts have been made to develop polarizable MM models,46−49 including in QM/MM simulations.50 In addition to these efforts, the present study further highlights that, as emphasized in previous studies,51,52 it is crucial to balance the interaction between QM and MM components at the QM−QM, QM−MM, and MM−MM levels.

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.9b00180.

■

COMPUTATIONAL DETAILS We implemented the SCMP method in a local version of the GROMACS 5.0.7 package.53−55 In all the SCMP simulations, a total of 80 QM/MM partitionings were considered. The system consisted of 2048 water molecules in a periodic cubic box with a side length of 40.28 Å, so that an experimental density of 0.997 g/cm3 was realized for bulk water at 300 K. Each partitioning contained one QM solute water molecule and 32 solvent water molecules. The transition parameters sQM, tQM, sMM, and tMM used for the SCMP method are listed in Table 1. In the

tMM(sQM)

tQM

3OBw 3OB

4.5 3.8

6.6 6.0

8.8 8.8

Details of progress functions and the transition parameters, details of evaluation of the effective PMF, QM profiles based on the solvent index, cumulative coordinate number of solvent molecules, mean forces versus distance from the QM center, oxygen and hydrogen PMFs shown with larger range of potential energy and with different transition parameters, and the used parameter set (Table S1). (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Hiroshi C. Watanabe: 0000-0003-4379-8633 Qiang Cui: 0000-0001-6214-5211 Author Contributions

H.C.W. designed this work and conducted all the simulations and analyses. The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

Table 1. Transition Parameters (Å) Used in the SCMP Simulations sMM

ASSOCIATED CONTENT

S Supporting Information *

■

DFTB parameter set

Article

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Yohich Suzuki and Yutaka Shikano for the fruitful discussions. This research was supported by JSPS KAKENHI Grant Number JP17K15101 and JST PRESTO Grant Number JPMJPR17GC.

partitioning updating protocol, we allowed updated partitionings to have a degree of order of 75% for efficiency, as discussed in previous work.28 We employed the SPC-Fw water model31 for the MM model in the SCMP method, as well as the LennardJones terms of the QM/MM interactions. For the QM part, we employed DFTB320,21 implemented in GROMACS, as reported previously.56 For the DFTB3 parameters, we used the standard 3OB22 and the 3OBw sets.29 The electrostatic interactions in the MM−MM and QM−MM models were calculated with the particle-mesh Ewald method57 with a switching function in the range from 8.5 to 9.0 Å. Likewise, the electrostatic potentials on the QM atoms induced by the charges of the MM atoms were obtained with the smooth particle-mesh Ewald method with a switching function for electrostatic interactions (electrostatic embedding), whereas van der Waals interactions were damped to zero with a switching function in the range between 8.5 and 9.0 Å, as detailed in the literature.56 All MD simulations were conducted for 100 ps with a 0.5 fs time step under NVE conditions with initial velocities generated with a Maxwell−Boltzmann distribution at 300 K. The PMFs were computed from QM/MM calculations with f(n) in eq 5 evaluated using 5000 atomic coordinates sampled through three independent 100 ps SCMP simulation trajectories. Note that 80 QM/MM calculations based on different partitionings were carried out for each atomic coordinate set. For comparison, we conducted a full QM MD simulation with DFTB3 for a bulk water system consisting of 2048 water molecules under periodic boundary conditions over 1 ps. With the DFTB+ program,58 the water molecules were simulated in a cubic box with a side length of 41.48 Å under NVE conditions.

■

ABBREVIATIONS SCMP, size-consistent multipartitioning; QM/MM, quantum mechanical/molecular mechanical; MD, molecular dynamics; DFTB, density functional tight binding; DFT, density functional theory; PMF, potential of the mean force; RDF, radial distribution function

■

REFERENCES

(1) Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions Dielectric, Electrostatic and Steric Stabilization of Carbonium-Ion in Reaction of Lysozyme. J. Mol. Biol. 1976, 103 (2), 227−249. (2) Singh, U. C.; Kollman, P. A. A Combined Abinitio QuantumMechanical and Molecular Mechanical Method for Carrying out Simulations on Complex Molecular-Systems - Applications to the Ch3cl + Cl- Exchange-Reaction and Gas-Phase Protonation of Polyethers. J. Comput. Chem. 1986, 7 (6), 718−730. (3) Field, M. J.; Bash, P. A.; Karplus, M. A Combined QuantumMechanical and Molecular Mechanical Potential for MolecularDynamics Simulations. J. Comput. Chem. 1990, 11 (6), 700−733. (4) Gao, J. L.; Ma, S. H.; Major, D. T.; Nam, K.; Pu, J. Z.; Truhlar, D. G. Mechanisms and free energies of enzymatic reactions. Chem. Rev. 2006, 106 (8), 3188−3209. (5) Bulo, R. E.; Michel, C.; Fleurat-Lessard, P.; Sautet, P. Multiscale Modeling of Chemistry in Water: Are We There Yet? J. Chem. Theory Comput. 2013, 9 (12), 5567−5577. (6) Rowley, C. N.; Roux, B. The Solvation Structure of Na+ and K+ in Liquid Water Determined from High Level ab Initio Molecular Dynamics Simulations. J. Chem. Theory Comput. 2012, 8 (10), 3526− 3535. 3926


Article


using high quality meta-GGA functionals. Chem. Sci. 2017, 8 (5), 3554−3565. (28) Watanabe, H. C. Improvement of Performance, Stability and Continuity by Modified Size-Consistent Multipartitioning Quantum Mechanical/Molecular Mechanical Method. Molecules 2018, 23 (8), 1882. (29) Goyal, P.; Qian, H. J.; Irle, S.; Lu, X. Y.; Roston, D.; Mori, T.; Elstner, M.; Cui, Q. Molecular Simulation of Water and Hydration Effects in Different Environments: Challenges and Developments for DFTB Based Models. J. Phys. Chem. B 2014, 118 (38), 11007−11027. (30) Watanabe, H. C.; Banno, M.; Sakurai, M. An adaptive quantum mechanics/molecular mechanics method for the infrared spectrum of water: incorporation of the quantum effect between solute and solvent. Phys. Chem. Chem. Phys. 2016, 18 (10), 7318−7333. (31) Wu, Y. J.; Tepper, H. L.; Voth, G. A. Flexible simple point-charge water model with improved liquid-state properties. J. Chem. Phys. 2006, 124 (2), 024503. (32) Soper, A. K. The radial distribution functions of water as derived from radiation total scattering experiments: is there anything we can say for sure? ISRN Phys. Chem. 2013, 2013, 279463. (33) Boereboom, J. M.; Potestio, R.; Donadio, D.; Bulo, R. E. Toward Hamiltonian Adaptive QM/MM: Accurate Solvent Structures Using Many-Body Potentials. J. Chem. Theory Comput. 2016, 12 (8), 3441− 3448. (34) Duster, A. W.; Wang, C. H.; Lin, H. Adaptive QM/MM for Molecular Dynamics Simulations: 5. On the Energy-Conserved Permuted Adaptive-Partitioning Schemes. Molecules 2018, 23 (9), 2170. (35) Jiang, T.; Simko, S.; Bulo, R. E. Accurate Quantum Mechanics/ Molecular Mechanics Simulation of Aqueous Solutions with Tailored Molecular Mechanics Models. J. Chem. Theory Comput. 2018, 14 (8), 3943−3954. (36) Freindorf, M.; Gao, J. L. Optimization of the Lennard-Jones parameters for a combined ab initio quantum mechanical and molecular mechanical potential using the 3-21G basis set. J. Comput. Chem. 1996, 17, 386−395. (37) Riccardi, D.; Li, G.; Cui, Q. The importance of van der Waals interactions in QM/MM simulations. J. Phys. Chem. B 2004, 108, 6467−6478. (38) Pavlov, E.; Taiji, M.; Scukins, A.; Markesteijn, A.; Karabasov, S.; Nerukh, D. Visualising and controlling the flow in biomolecular systems at and between multiple scales: from atoms to hydrodynamics at different locations in time and space. Faraday Discuss. 2014, 169, 285− 302. (39) Watanabe, H. C.; Welke, K.; Sindhikara, D. J.; Hegemann, P.; Elstner, M. Towards an understanding of channelrhodopsin function: simulations lead to novel insights of the channel mechanism. J. Mol. Biol. 2013, 425 (10), 1795−1814. (40) Goyal, P.; Yang, S.; Cui, Q. Microscopic basis for kinetic gating in cytochrome c oxidase: insights from QM/MM analysis. Chem. Sci. 2015, 6 (1), 826−841. (41) Duster, A. W.; Garza, C. M.; Aydintug, B. O.; Negussie, M. B.; Lin, H. Adaptive Partitioning QM/MM for Molecular Dynamics Simulations: 6. Proton Transport through a Biological Channel. J. Chem. Theory Comput. 2019, 15 (2), 892−905. (42) Gaus, M.; Cui, Q.; Elstner, M. Density functional tight binding: application to organic and biological molecules. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4 (1), 49−61. (43) Dral, P. O.; Wu, X.; Sporkel, L.; Koslowski, A.; Thiel, W. Semiempirical Quantum-Chemical Orthogonalization-Corrected Methods: Benchmarks for Ground-State Properties. J. Chem. Theory Comput. 2016, 12, 1097−1120. (44) Dral, P. O.; von Lilienfeld, O. A.; Thiel, W. Machine Learning of Parameters for Accurate Semiempirical Quantum Chemical Calculations. J. Chem. Theory Comput. 2015, 11, 2120−2125. (45) Cui, Q. Quantum mechanical methods in biochemistry and biophysics. J. Chem. Phys. 2016, 145, 140901.

(7) Shiga, M.; Masia, M. Boundary based on exchange symmetry theory for multilevel simulations. I. Basic theory. J. Chem. Phys. 2013, 139 (4), 044120. (8) Takahashi, H.; Kambe, H.; Morita, A. A simple and effective solution to the constrained QM/MM simulations. J. Chem. Phys. 2018, 148 (13), 134119. (9) Kerdcharoen, T.; Morokuma, K. ONIOM-XS: an extension of the ONIOM method for molecular simulation in condensed phase. Chem. Phys. Lett. 2002, 355 (3−4), 257−262. (10) Heyden, A.; Lin, H.; Truhlar, D. G. Adaptive partitioning in combined quantum mechanical and molecular mechanical calculations of potential energy functions for multiscale simulations. J. Phys. Chem. B 2007, 111 (9), 2231−2241. (11) Bulo, R. E.; Ensing, B.; Sikkema, J.; Visscher, L. Toward a Practical Method for Adaptive QM/MM Simulations. J. Chem. Theory Comput. 2009, 5 (9), 2212−2221. (12) Bernstein, N.; Varnai, C.; Solt, I.; Winfield, S. A.; Payne, M. C.; Simon, I.; Fuxreiter, M.; Csanyi, G. QM/MM simulation of liquid water with an adaptive quantum region. Phys. Chem. Chem. Phys. 2012, 14 (2), 646−656. (13) Takenaka, N.; Kitamura, Y.; Koyano, Y.; Nagaoka, M. The number-adaptive multiscale QM/MM molecular dynamics simulation: Application to liquid water. Chem. Phys. Lett. 2012, 524, 56−61. (14) Waller, M. P.; Kumbhar, S.; Yang, J. A Density-Based Adaptive Quantum Mechanical/Molecular Mechanical Method. ChemPhysChem 2014, 15 (15), 3218−3225. (15) Watanabe, H. C.; Kubar, T.; Elstner, M. Size-Consistent Multipartitioning QM/MM: A Stable and Efficient Adaptive QM/ MM Method. J. Chem. Theory Comput. 2014, 10 (10), 4242−4252. (16) Field, M. J. An Algorithm for Adaptive QC/MM Simulations. J. Chem. Theory Comput. 2017, 13 (5), 2342−2351. (17) Pezeshki, S.; Lin, H. Adaptive-Partitioning QM/MM for Molecular Dynamics Simulations: 4. Proton Hopping in Bulk Water. J. Chem. Theory Comput. 2015, 11 (6), 2398−2411. (18) Ensing, B.; Nielsen, S. O.; Moore, P. B.; Klein, M. L.; Parrinello, M. Energy conservation in adaptive hybrid atomistic/coarse-grain molecular dynamics. J. Chem. Theory Comput. 2007, 3 (3), 1100−1105. (19) Watanabe, H. C.; Kubillus, M.; Kubar, T.; Stach, R.; Mizaikoff, B.; Ishikita, H. Cation solvation with quantum chemical effects modeled by a size-consistent multi-partitioning quantum mechanics/molecular mechanics method. Phys. Chem. Chem. Phys. 2017, 19 (27), 17985− 17997. (20) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-consistent-charge densityfunctional tight-binding method for simulations of complex materials properties. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58 (11), 7260−7268. (21) Gaus, M.; Cui, Q. A.; Elstner, M. DFTB3: Extension of the SelfConsistent-Charge Density-Functional Tight-Binding Method (SCCDFTB). J. Chem. Theory Comput. 2011, 7 (4), 931−948. (22) Gaus, M.; Goez, A.; Elstner, M. Parametrization and Benchmark of DFTB3 for Organic Molecules. J. Chem. Theory Comput. 2013, 9 (1), 338−354. (23) Choi, T. H.; Liang, R.; Maupin, C. M.; Voth, G. A. Application of the SCC-DFTB Method to Hydroxide Water Clusters and Aqueous Hydroxide Solutions. J. Phys. Chem. B 2013, 117 (17), 5165−5179. (24) Goyal, P.; Elstner, M.; Cui, Q. Application of the SCC-DFTB Method to Neutral and Protonated Water Clusters and Bulk Water. J. Phys. Chem. B 2011, 115 (20), 6790−6805. (25) Fernandez-Serra, M. V.; Artacho, E. Network equilibration and first-principles liquid water. J. Chem. Phys. 2004, 121 (22), 11136− 11144. (26) Kuo, I. F. W.; Mundy, C. J.; McGrath, M. J.; Siepmann, J. I.; VandeVondele, J.; Sprik, M.; Hutter, J.; Chen, B.; Klein, M. L.; Mohamed, F.; Krack, M.; Parrinello, M. Liquid water from first principles: Investigation of different sampling approaches. J. Phys. Chem. B 2004, 108 (34), 12990−12998. (27) Ruiz Pestana, L.; Mardirossian, N.; Head-Gordon, M.; HeadGordon, T. Ab initio molecular dynamics simulations of liquid water 3927


Article

Journal of Chemical Theory and Computation (46) Ishiyama, T.; Morita, A. Analysis of anisotropic local field in sum frequency generation spectroscopy with the charge response kernel water model. J. Chem. Phys. 2009, 131 (24), 244714. (47) Babin, V.; Medders, G. R.; Paesani, F. Toward a Universal Water Model: First Principles Simulations from the Dimer to the Liquid Phase. J. Phys. Chem. Lett. 2012, 3 (24), 3765−3769. (48) Medders, G. R.; Paesani, F. Infrared and Raman Spectroscopy of Liquid Water through ″First-Principles″ Many-Body Molecular Dynamics. J. Chem. Theory Comput. 2015, 11 (3), 1145−1154. (49) Laury, M. L.; Wang, L. P.; Pande, V. S.; Head-Gordon, T.; Ponder, J. W. Revised Parameters for the AMOEBA Polarizable Atomic Multipole Water Model. J. Phys. Chem. B 2015, 119 (29), 9423−9437. (50) Kratz, E. G.; Walker, A. R.; Lagardere, L.; Lipparini, F.; Piquemal, J.-P.; Cisneros, G. A. LICHEM: A QM/MM program for simulations with multipolar and polarizable force fields. J. Comput. Chem. 2016, 37, 1019−1029. (51) Geerke, D. P.; Thiel, S.; Thiel, W.; van Gunsteren, W. F. QMMM interactions in simulations of liquid water using combined semiempirical/classical Hamiltonians. Phys. Chem. Chem. Phys. 2008, 10 (2), 297−302. (52) Lu, X. Y.; Fang, D.; Ito, S.; Okamoto, Y.; Ovchinnikov, V.; Cui, Q. QM/MM free energy simulations: recent progress and challenges. Mol. Simul. 2016, 42 (13), 1056−1078. (53) Bjelkmar, P.; Larsson, P.; Cuendet, M. A.; Hess, B.; Lindahl, E. Implementation of the CHARMM Force Field in GROMACS: Analysis of Protein Stability Effects from Correction Maps, Virtual Interaction Sites, and Water Models. J. Chem. Theory Comput. 2010, 6 (2), 459− 466. (54) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1, 19−25. (55) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008, 4 (3), 435−447. (56) Kubar, T.; Welke, K.; Groenhof, G. New QM/MM implementation of the DFTB3 method in the gromacs package. J. Comput. Chem. 2015, 36 (26), 1978−1989. (57) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald - an N.Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98 (12), 10089−10092. (58) Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB+, a sparse matrix-based implementation of the DFTB method. J. Phys. Chem. A 2007, 111 (26), 5678−5684.

3928


MM Boundary Artifacts and Correction in

Recommend Documents