MM Boundary Artifacts and Correction in

Subscriber access provided by UNIV AUTONOMA DE COAHUILA UADEC

Dynamics

Quantitative analysis of QM/MM boundary artifacts and correction in adaptive QM/MM simulations Hiroshi C Watanabe, and Qiang Cui J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.9b00180 • Publication Date (Web): 16 May 2019 Downloaded from http://pubs.acs.org on May 17, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Quantitative analysis of QM/MM boundary artifacts and correction in adaptive QM/MM simulations Hiroshi C. Watanabe1,2*, Qiang Cui3 1 Quantum

Computing Center, Keio University, 3-14-1 Hiyoshi, Kohoku-ku, Yokohama, 2238522 Japan

2Japan

Science and Technology Agency, PRESTO, 4-18, Honcho, Kawaguchi, Saitama 3320012, Japan

3Departments

of Chemistry, Physics, Biomedical engineering, Boston University, 590 Commonwealth Avenue, Boston, MA 02215

KEYWORDS. QM/MM, adaptive QM/MM method, dynamics, solvation, condensed phase, size-consistent multi-partitioning

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 54

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.


2

Page 3 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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


3


Page 4 of 54

interactions. Therefore, simulation frameworks that maintain the QM solvent 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 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 BCC8. 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


4

Page 5 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 multi-partitioning schemes, whereby various QM/MM partitionings that have different combinations of QM solvents are defined. In these methods, the effective potential energy 𝑉eff is defined as 𝑉eff(𝑞) = ∑𝑛𝜎(𝑛)(𝑞)𝑉(𝑛)(𝑞),

(1)

where 𝑉(𝑛) and 𝜎(𝑛) are the QM/MM potential energy and the normalized weight function of the n-th partitioning, respectively, 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.


5


Page 6 of 54

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 𝑉eff, 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, i.e., the continuity of the potential, is mainly related to the weight function 𝜎(𝑛) in Eq. 1. In general, the number of possible QM/MM

( )

𝑀 𝑀 partitionings is ∑𝑚 𝑚 , where M and m denotes 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


6

Page 7 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(partitioning update). However, abrupt updates lead to temporal discontinuities in 𝑉eff. Thus, the weight function 𝜎(𝑛) is designed to suppress discontinuities in 𝑉eff such that the partitionings included/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 partitioning 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.


7


Page 8 of 54

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 𝐻=

𝑝2𝑗

∑2𝑚 + 𝑉 𝑗

𝑗

eff.

(2)

where, mj and pj represent the mass and the moment of the m-th particle, respectively. Then, the effective force 𝐹eff used to propagate the MD trajectory is, 𝐹eff =

―∂𝑉eff ∂𝑞 𝑁

=

∑𝜎

(𝑛)

𝑁

(𝑛)

(𝑞)𝑓

𝑛

―

∑ 𝑛

∂𝜎(𝑛) (𝑛) 𝑉 (𝑞) ∂𝑞

(3)

= 𝐹 + 𝐹trans, where 𝑓(𝑛) is the force in the n-th QM/MM partitioning, given by the negative derivative of 𝑉(𝑛). The second term, 𝐹trans, sometimes referred to as the transition force, is not physically realistic because it arises from the derivative of the artificial weight function 𝜎(𝑛). Because the transition force can lead to unrealistic dynamics and solvation structures,6 in potential-oriented approaches it is important to minimize 𝐹trans. Note that the contribution of 𝐹trans to 𝐹eff is significant in magnitude in multi-scale partitioning approaches, such as PAP, SAP, and DAS, in which the partitionings have different numbers of QM solvent molecules. Because the potential energy of a QM model is usually much more negative than that of


8

Page 9 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


an MM model, the QM/MM energy 𝑉(𝑛) becomes increasingly negative as the number of QM solvent molecules increases. Because 𝐹eff should work to minimize 𝑉eff, the solvation structure can be distorted such that partitionings with larger QM regions have a larger value of 𝜎(𝑛). 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 multi-scale partitioning approaches and, correspondingly, the artifacts are smaller. In force-oriented approaches, the Hamiltonian H is written as 𝐻=

𝑝2𝑗

∑2𝑚 ∑ 𝑗

+

𝜎(𝑛)(𝑞)𝑉(𝑛)(𝑞) ―

𝑛

∑∫ 𝑛

∂𝜎(𝑛) (𝑛) 𝑉 (𝑞)𝑑𝑞 ∂𝑞

(4)

= 𝐾(𝑝) + 𝑉eff(𝑞) + 𝐵(𝑞), where K(p) is the kinetic energy and B(q) is a bookkeeping term, originally proposed by Ensing et al.18 for atomistic/coarse-grained 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 𝑁

𝐹eff =

∑𝜎

(𝑛)

(𝑞)𝑓(𝑛)(𝑞)

(5)

𝑛

Unlike potential-oriented approaches, force-oriented approaches do not exhibit artifacts due to the weight function and yield more realistic dynamics. The bookkeeping term can decrease monotonically over the course of the MD


9


Page 10 of 54

simulation. Correspondingly, the internal energy, which consists of K and 𝑉eff, 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. 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 𝐹trans. 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 multi-scale 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 with multi-scale 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


10

Page 11 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


errors in the description of hydrogen-bonding interactions involving water,21 especially for bulk 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 co-workers employed an inverse Monte Carlo approach to adjust the O-H repulsive potential and reproduce the radial distribution function of bulk water at ambient condition. 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 systems 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


11


Page 12 of 54

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.

2. Size-consistent multi-partitioning 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 whose 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 𝐹eff 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 with the ordered QM region in the SCMP method (Scheme 1), the weight functions 𝜎(𝑛) for the n-th partitioning have to satisfy the two following conditions for the temporal continuity: (i) 𝜎(𝑛) = 0 when the QM region of the nth partitioning is disordered; (ii) 𝜎(𝑛) = 0 when the QM region of the n-th partitioning is perfectly ordered. (Scheme 1) To this end, we introduce the fadeout function 𝑂(𝑛) for the n-th partitioning:


12

Page 13 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


𝑚

𝑂(𝑛) QM

=

∏λ

(6)

(𝑛) 𝑗 (𝑞)

𝑗

where λ(𝑛) is a progress function of the j-th QM solvent molecules, ranging from 𝑗 0 to 1 depending on the distance between the QM solute and the j-th solvent molecule. In addition, the fade-in function 𝐼(𝑛) is defined as (7)

(𝑛) 𝐼(𝑛) QM = 1 ― 𝑂QM

We also introduce the fade-out and fade-in functions for the MM solvent molecules as 𝑀

𝑂(𝑛) MM

=

∏Λ

(𝑛) 𝑗 (𝑞)

(8)

𝑗

(𝑛) 𝐼(𝑛) MM = 1 ― 𝑂MM

where Λ(𝑛) is a progress function of the j-th MM solvent molecule. The specific 𝑗 forms of the progress functions can be found in previous works.15, 28　The weight function 𝜎 for the n-th partitioning is (𝑛)

𝜎

=

(𝑛) (𝑛) (𝑛) 𝑂(𝑛) QM𝐼QM𝑂MM𝐼MM 𝑁

(𝑘) (𝑘) (𝑘) ∑𝑘 𝑂(𝑘) QM𝐼QM𝑂MM𝐼MM

(9)

3. RESULTS The Bookkeeping term: Figure 2a shows that the Hamiltonian is numerically conserved in the SCMP simulations with DFTB3/3OBw, which is consistent with


13


Page 14 of 54

previous studies that used DFTB3/3OB15, 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 n-th partitioning is perfectly ordered, satisfying 𝜎(𝑛) = 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 𝜎(𝑛) > 0. Eventually, it will be completely disordered at t = 𝑡2 and 𝜎(𝑛) will become 0 again. The contribution of the n-th partitioning to the bookkeeping term for the time period between 𝑡1 and 𝑡2 is


14

Page 15 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


∆𝐵

(𝑛)

=―

∫

𝑞2∂𝜎(𝑛) 𝑞1

∂𝑞

𝑉(𝑛)(𝑞)𝑑𝑞

(10)

where 𝑞1 = {𝑞(𝑡1)} and 𝑞2 = {𝑞(𝑡2)}. Therefore, if ∆𝑉(𝑛) = 𝑉(𝑛)(𝑞2)–𝑉(𝑛)(𝑞1) ≠ 0, then ∆𝐵(𝑛) ≠ 0. Because the bookkeeping term is the sum of ∆𝐵(𝑛) over all partitionings, a monotonic decrease indicates that the QM/MM potential energies statistically satisfy ⟨∆𝑉(𝑛)⟩ < 0. Furthermore, in the SCMP simulations, because the disordered partitionings are updated with an ordered partitioning, a non-zero ∆𝐵 accumulates continuously, resulting in the constant drift of the bookkeeping term. Thus, the present results suggest that, for an arbitrary time period, the averaged potential energy difference ⟨∆𝑉(𝑛)⟩ 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 reference calculation using the CAMB3LYP functional with 6-311++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


15


Page 16 of 54

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 𝑉(𝑛) 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 𝑉(𝑛) obtained with the 3OBw parameters varies depending on the shape of the QM region more than that obtained with the 3OB parameters, implying that ∆𝑉3OBw < ∆𝑉3OB < 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 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


16

Page 17 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


center and/or attractive forces act on the MM solvent. Because a homogenous 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 𝑓(𝑛) (Eq. 5) via conventional QM/MM calculations in a post-processing 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 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).


17


Page 18 of 54

In contrast to the QM solvent, the MM solvent showed a flat PMF at the longdistance 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 unboundedly 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 ⟨∆𝑉(𝑛)⟩ < 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 𝑓(𝑛) 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


18

Page 19 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 𝑁

𝜔𝑗 = ∑𝑛 𝜎(𝑛)𝛿(𝑛) 𝑗 ,

(11)

where the subscript j represents the j-th nearest solvent from the QM center.

{

1 if 𝑗 ∈ 𝑆(𝑛) 𝛿(𝑛) 𝑗 = 0 otherwise

(12)

and 𝑆(𝑛) is a subset of the QM solvent in the n-th 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 a QM and an MM 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 𝜔(𝑗) 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 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


19


Page 20 of 54

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 artifacts cannot be completely removed in the SCMP method. Because the SCMP simulations with the DFTB3/3OB parameters results in a smaller drift of the bookkeeping term, it is insightful to compare the resulting PMF curves with those obtained using the 3OBw parameters. Whilst 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 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 ∑𝜎(𝑛)𝑓(𝑛), 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


20

Page 21 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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. Radial distribution function (RDF): The RDF provides an alternative way for analyzing the effective PMF in SCMP simulations and evaluating the error suppression. First, Figure 6 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 longdistance 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


21


Page 22 of 54

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 𝑑

〈𝐹eff(𝑟)〉 = 𝑘𝐵𝑇𝑑𝑟ln𝑔(𝑟) ,

(13)

we can evaluate the deviation of the PMFs based on the SCMP effective force 𝐹eff 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 𝑈(𝑛)(𝑞), which describe interactions between the QM solute and solvent molecules. Then, we replace the respective QM/MM potential 𝑉(𝑛) in Eq. 4 by 𝑉′(𝑛) = 𝑉(𝑛)(𝑞) + 𝑈(𝑛)(𝑞). Accordingly, the Hamiltonian is written as 𝑝2𝑗

𝐻 = ∑𝑗2𝑚 + ∑𝑛𝜎(𝑛)(𝑞)𝑉′(𝑛)(𝑞) ― ∑𝑛∫

∂𝜎(𝑛) ∂𝑞

𝑉′(𝑛)(𝑞)𝑑𝑞,

(14)

because the derivative of 𝑈(𝑛) in the bookkeeping term is canceled out. The corrective force 𝐹corr acting on the i-th solvent molecule is 𝑖 𝐹corr = ―

∑𝜎

(𝑛)

𝑛

=

∑𝜎

∂𝑈(𝑛) ∂𝑞

(15)

(𝑛) (𝑛)

𝔉

𝑛


22

Page 23 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Note that 𝑈(𝑛) 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, 𝑈(𝑛) can work simultaneously on the QM and MM solvents. Moreover, the corrective potential 𝑈(𝑛) is allowed to have different forms, namely 𝑈QM and 𝑈MM, for the QM and MM solvents, respectively. Both 𝑈QM and 𝑈MM are preferably monotonic, and the difference between the corrective potential as solvent exchange progresses should satisfy ∆ 𝑈(𝑛) = 𝑈(𝑛)(𝑞2)–𝑈(𝑛)(𝑞1) < 0 to suppress the drift of the bookkeeping term. In other words, 𝐹(𝑛) acting on the QM solvent should be attractive toward the QM center, whereas 𝐹(𝑛) acting on the MM solvent should be repulsive. Because of requirement (C), 𝑈QM should be continuous where 𝑟 < 𝑡QM, whereas 𝑈MM should be continuous where 𝑟 > 𝑠MM, with 𝑡QM and 𝑠MM 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 𝐹(𝑛) to be an attractive force acting on the QM solvent molecules. Note that the progress function 𝜆, which ranges between 0 and 1, satisfies the


23


Page 24 of 54

monotonic and continuous requirements stated above. In the present study, we define 𝑈QM as 𝑈(𝑛) QM(𝑞) = 𝑈 ∙ 𝜆(𝑞),

(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 previous study.28 Likewise, 𝜆 in Eq. 16 consists of the three spline curves defined in the range 𝑠MM < 𝑟 < 𝑡QM (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


24

Page 25 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 𝑉(𝑛) 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 multi-scale partitioning method as


25


𝐻=

𝑝2𝑗

∑2𝑚 + ∑𝜎

(𝑛)

𝑗

(𝑞)𝑉(𝑛)(𝑞) ―

𝑛

∑ (𝜎

(𝑛)

∑𝛿

)

(𝑛) 𝑗 𝐶

(𝑞)

𝑛

Page 26 of 54

𝑗

,

(17)

where 𝛿(𝑛) equals unity when the j-th solvent is QM in the n-th 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 𝐹eff =

∑𝜎

(𝑛) (𝑛)

𝑉

𝑛

―

∑ 𝑛

∂𝜎(𝑛) (𝑛) 𝑉 ― ∂𝑞

(

∑

𝛿(𝑛) 𝑗 𝐶

)

(18)

― 𝑉MM(𝑞) ,

𝑗

where 𝑉MM 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 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 parameterization. In other words, the SCMP method can remove all the artifacts caused by the derivative of 𝜎(𝑛). Therefore, the SCMP method can attribute the drift of the bookkeeping term purely to the QM/MM boundary artifacts in 𝑉(𝑛), 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


26

Page 27 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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


27


Page 28 of 54

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 Lennard–Jones 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


28

Page 29 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 demand for multi-scale 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/binding pockets.38-41 Considerable efforts have been made to improve approximate QM methods to achieve a good trade-off with the computational


29


efficiency.29,

42-45

Page 30 of 54

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 studies51-52, it is crucial to balance the interaction between QM and MM components at the QM–QM, QM–MM, and MM–MM levels.

Computational details We implemented the SCMP method in a local version of the GROMACS 5.0.7 package53-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 𝑠QM, 𝑡QM, 𝑠MM, and 𝑡QM used for the SCMP method are listed in Table 1. In the 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 Lennard– Jones terms of the QM/MM interactions. For the QM part, we employed DFTB32021

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


30

Page 31 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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

Supporting Information


31


Page 32 of 54

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

AUTHOR INFORMATION Corresponding Author * Hiroshi C. Watanabe, email: [email protected]

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.

ACKNOWLEDGMENT


32

Page 33 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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

ABBREVIATIONS SCMP,

size-consistent

multi-partitioning;

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.

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.


33


Page 34 of 54

(left) One 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.


34

Page 35 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. The 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.


35


Page 36 of 54

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 6311++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.


36

Page 37 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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


37


Page 38 of 54

lines, are estimated by combining the QM and MM PMFs using the QM profile (see Figure 5).

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.


38

Page 39 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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


39


Page 40 of 54

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


40

Page 41 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


bookkeeping term and temperature obtained with the 3OBw parameters and U = 0.0 kJ/mol are also presented (blue line) as a reference.

Figure 7. 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,


41


Page 42 of 54

respectively. The results of the full QM simulation with DFTB3/3OBw are also presented (magenta line).

Scheme 1. The weight function 𝜎 and partitioning update in the SCMP method. 𝜎 will become 0, when the QM region in a partitioning is completely ordered or disordered. A partitioning whose QM region is disordered is replaced by one whose QM region is ordered.


42

Page 43 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1. Transition parameters (Å) used in the SCMP simulations. DFTB parameter set

𝑠𝑀𝑀

𝑡𝑀𝑀(𝑠𝑄𝑀)

𝑡𝑄𝑀

3OBw

4.5

6.6

8.8

3OB

3.8

6.0

8.8

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 Quantum-Mechanical

and Molecular Mechanical Method for Carrying out Simulations on Complex Molecular-Systems - Applications to the Ch3cl + Cl- Exchange-Reaction and GasPhase Protonation of Polyethers. J. Comput. Chem. 1986, 7 (6), 718-730.


43


3.

Page 44 of 54

Field, M. J.; Bash, P. A.; Karplus, M., A Combined Quantum-Mechanical

and Molecular Mechanical Potential for Molecular-Dynamics 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. 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.


44

Page 45 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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), 22122221. 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), 32183225. 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.


45


16.

Page 46 of 54

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 density-functional tight-binding method for simulations of complex materials properties. Phys. Rev. B 1998, 58 (11), 7260-7268. 21.

Gaus, M.; Cui, Q. A.; Elstner, M., DFTB3: Extension of the Self-Consistent-

Charge Density-Functional Tight-Binding Method (SCC-DFTB). J. Chem. Theory Comput. 2011, 7 (4), 931-948.


46

Page 47 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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.

Pestana, L. R.; Mardirossian, N.; Head-Gordon, M.; Head-Gordon, T., Ab

initio molecular dynamics simulations of liquid water using high quality metaGGA functionals. Chem. Sci. 2017, 8 (5), 3554-3565.


47


28.

Page 48 of 54

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


48

Page 49 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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


49


39.

Page 50 of 54

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), 826841. 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.; O. A. von Lilienfeld; W. Thiel, Machine Learning of Parameters

for Accurate Semiempirical Quantum Chemical Calculations. J. Chem. Theory Comput. 2015, 11, 2120-2125.


50

Page 51 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


45.

Cui, Q., Quantum mechanical methods in biochemistry and biophysics.

J. Chem. Phys. 2016, 145, 140901. 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


51.

Page 52 of 54

Geerke, D. P.; Thiel, S.; Thiel, W.; van Gunsteren, W. F., QM-MM

interactions

in

simulations

of

liquid

water

using

combined

semi-

empirical/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.; Páll, S.; Smith, J. C.; Hess, B.;

Lindahl, E., GROMACS: High performance molecular simulations through multilevel 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), 19781989.


52

Page 53 of 54 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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


53



Page 54 of 54

MM Boundary Artifacts and Correction in

Recommend Documents