MM Simulations - Journal of Chemical

Article pubs.acs.org/JCTC

An Algorithm for Adaptive QC/MM Simulations Martin J. Field* Dynamo Team/DYNAMOP Group, UMR5075, Université Grenoble I, CEA, CNRS, Institut de Biologie Structurale, 71 Avenue des Martyrs, CS 10090, 38044 Grenoble Cedex 9, France ABSTRACT: An algorithm is proposed for the simulation of molecular systems with hybrid quantum chemical (QC) and molecular mechanical (MM) potentials that permits the adaptive partitioning of the atoms in the system between QC and MM regions. In contrast to existing methods, the algorithm requires only a single QC calculation of the QC/ MM system per energy calculation and yet has consistent energy and forces, which makes it suitable for geometry optimizations and molecular dynamics calculations within the microcanonical ensemble, in addition to other types of simulation. This article describes the algorithm and its implementation, presents some simple test cases with both semiempirical and density functional theory QC/MM potentials, and discusses perspectives for future work. overview of the literature.5−8 Nevertheless, a number of pertinent adaptive algorithms will be highlighted. One of the earliest of these, and probably the one that has been most used for applications, is the hot-spot family of methods developed by Rode and co-workers.9,10 These methods are quite efficient and ensure a smooth transition of atoms in the buffer region between QC and MM regions by defining their forces to be the weighted superposition of forces arising from the pure QC and MM potentials, respectively. Unfortunately, the formulation of the methods is such that they suffer from a number of deficiencies, including the absence of an energy consistent with the forces, which makes them unsuitable for many types of calculation. A contrast to the hot-spot methods is provided by the permuted adaptive partitioning (PAP) algorithm of Heyden et al.11 This approach, which is outlined in detail in section 2, is exact in the sense that it has consistent energy and forces, but it has the huge disadvantage that its cost scales as 2N, where N is the number of distinct mobile groups in the buffer region. To overcome this, in the same paper, Heyden and co-workers developed a more efficient version of the PAP algorithm sorted adaptive partitioningwith a cost that scales as ∼O(N). A related algorithmdistance adaptive partitioningwith similar scaling was also later suggested by Bulo et al.12 Although only a few of the available QC/MM adaptive algorithms have been mentioned above, all of which the author is aware are either significantly more costly than a QC/MM calculation with fixed partitioning or introduce approximations that limit their usefulness to a restricted class of problem or simulation type. The aim of this work is to make an initial attempt to propose an algorithm that eliminates these problems

1. INTRODUCTION Hybrid potentials that combine quantum chemical (QC) and molecular mechanical (MM) methods have become standard tools in molecular simulation.1−4 Most simulations with hybrid potentials have been done with a fixed partitioning of the atoms in the system between QC and MM regions, which is chosen beforehand. There are, however, significant classes of application in which a fixed partitioning is problematic. Notable examples are solution or enzyme reactions that involve solvent. A proper treatment of these cases requires the inclusion of one or more solvent molecules in the QC region, but a specific choice is rendered difficult by the mobility of the solvent molecules, which can diffuse in and out of range of the solute or active site. A number of approaches that permit the adaptive or dynamic partitioning of atoms between QC and MM regions in a hybrid potential simulation have been proposed. Most employ a framework that is illustrated schematically in Figure 1, in which a system is divided into a central QC region, an outer MM region, and an intermediate transition or buffer region where atoms transit between QC and MM representations. The boundaries between the different regions are typically defined using a geometrical criterion, such as the distance from a given atom or group of atoms in the QC region. Note that it is possible to mix static and dynamic partitioning, depending upon the application. Thus, for example, in a solute/solvent system, it would be usual to statically assign the solute to the QC region but allow the solvent to be able to move freely between regions. By contrast, in the study of an enzyme reaction, it would be more common to statically partition the protein and substrate between QC and MM regions and permit only the more mobile solvent molecules to evolve adaptively. There is not sufficient space here to provide a comprehensive review of published adaptive QC/MM methods. Readers are, instead, referred to recent reviews that provide an up-to-date © 2017 American Chemical Society

Received: January 31, 2017 Published: April 6, 2017 2342

DOI: 10.1021/acs.jctc.7b00099 J. Chem. Theory Comput. 2017, 13, 2342−2351

Article

Journal of Chemical Theory and Computation

multiplied by the product of the MM weights of the MM fragments within the partition. That is, wP =

∏ I ∈ QCP

λI

∏ J ∈ MMP

(1 − λJ ) (1)

The total PAP energy of the system, EPAP, is then given by the weighted sum of the energies of each partition, EP, noting that the weights are automatically normalized due to their construction. Thus, one has 2N

E PAP =

∑ wPEP P=1

(2)

If the partition energies and weights are continuous and wellbehaved, then quantities derived from them, such as the atomic forces, should also be so. 2.2. The Scaled Interaction Single Partition Adaptive Algorithm. The PAP method deals with the problem of the transition of fragments between QC and MM regions by generating a smooth combination of all possible QC/MM partitionings involving the buffer fragments. This results in an accurate algorithm but one that is impractically expensive, except for very small numbers of fragments in the transition region. Is there a way around this? In place of weighting distinct QC/MM partitionings, each with their own well-defined wave function, might it not be feasible instead to have a single partitioning, with an “average” wave function, that somehow encapsulates the mixed QC and MM nature of the buffer region fragments? There are undoubtedly several ways to achieve this, but the method discussed here employs scaling to appropriately weight the interactions involving the QC and MM representations of the buffer region fragments. This is most easily explained schematically, and so readers are referred to Figure 2, which

Figure 1. A schematic of the partitioning of the atoms in a system in an adaptive QC/MM simulation. The QC region is defined with respect to the central glycine zwitterion. All molecules of water (shown in gold) within rc Å of the α-carbon of the glycine are also defined as being in the QC region. Water molecules (cyan lines) further than rb Å from the glycine are in the MM region. The intermediate water molecules (in silver) are in the buffer or transition region.

and that allows adaptive calculations with similar cost, precision, and utility to those with fixed partitioning.

2. METHODS This section first describes the PAP algorithm,11 as this is used as the reference adaptive QC/MM method in the paper, before going on to outline the new adaptive algorithm. It finishes with some implementation details. 2.1. The PAP Algorithm. In what follows, we suppose that the system is partitioned as shown in Figure 1 and that only a subset of the atoms in the system can transform between QC and MM representations. In addition, we limit transformations to entities that we term fragments and do not allow atoms within a multiatom fragment to transform independently. Typical examples of fragments would be solvent molecules or ions in solution. Each time an energy is calculated in the PAP method, a weight is determined for each fragment, λI, that gives its fraction of QC character. Its corresponding fraction of MM character is simply 1 − λI. The weighting function is designed so that it is unity for fragments within the QC region, is zero within the MM region, and tapers smoothly between these values in the buffer. Now suppose that there are N fragments within the buffer region. It is not difficult to show that it is possible to generate 2N distinct QC/MM partitionings using these fragments, starting from one in which all fragments are treated as MM and ending with one in which they are all treated as QC. The weight of one of these partitions, wP, is given by the product of the QC weights of the QC fragments within the partition

Figure 2. A schematic of the scalings used for the different fragment energy contributions. Red, MM; orange, self and vacuum; green, QC/ MM; blue, QC. The subscripts I and J refer to different fragments in the buffer region.

gives the scalings for the case of pairwise interactions. As in the PAP method, each fragment is assigned a weighting, λI, that gives its fraction of QC character. This weighting goes from unity when the fragment is in the QC region to zero when it is in the MM region. Within the buffer, a fragment will have both QC and MM representations. These are shown in orange for fragment I in Figure 2. All interactions between the QC representation of fragment I and the atoms in the QC and MM regions are scaled 2343


Article

Journal of Chemical Theory and Computation by the fragment’s QC weighting, λI. Likewise, all interactions between its MM representation and the QC and MM region atoms are scaled by its MM weighting, 1 − λI. The scalings for interactions between two fragments in the buffer region follow the same pattern, except that now the scaling is a product of the QC or MM weight factors of each fragment. This is illustrated in Figure 2 for the interactions between fragments I and J. Counting up the scalings for each of the different types of interaction involving a fragment shows that they each add up to one as they should. The remaining question concerns what to do about intrafragment interactions. In principle, these should be scaled by λI and by 1 − λI for the fragment’s QC and MM representations, respectively. This can indeed be done directly for the MM representation as the most commonly used MM potentials have well-defined atomic contributions. For the QC representation, however, it poses a problem as scaling would severely compromise the internal electronic structure of the fragment arising from the QC calculation. Instead, therefore, the intrafragment interactions are left unscaled, and a set of corrections, shown on the right of Figure 2, is applied for the representations of each buffer fragment to avoid overcounting. In this work, these corrections are determined by isolating each fragment in vacuum with the same geometry that it has in the buffer and calculating its energy using its MM and QC representations. These energies are then scaled as shown in the figure and added into the total energy for the QC/MM system. Alternative correction techniques can be envisaged. For example, one method that was initially considered was to gradually mutate the QC representation of a fragment to that of a known reference state as the MM boundary was approached. This technique would merit investigation in the futureit is more complicated to implement but has the advantage that it dispenses completely with the additional fragment vacuum QC calculations. We denote the above scheme the scaled interaction single partition adaptive, or SISPA, algorithm. To summarize, it requires a single QC/MM calculation on the complete system, accompanied by a number of smaller and much cheaper QC and MM calculations of the fragments isolated from the buffer region. In the QC/MM calculation, various scalings are applied to the QC/QC, QC/MM, and MM/MM interactions, and the QC part of the energy evaluation requires the inclusion of all QC core and buffer fragment atoms as QC atoms. What assumptions are involved in the formulation of the algorithm? The principal one concerns the behavior of the QC representation of a fragment as it traverses the buffer region. The method relies on the fact that the energy and wave function (or electron density) of the fragment arising from the QC/MM’s QC calculation will tend to that of the isolated fragment as its interactions with the remaining atoms in the system are decoupled on approaching the boundary of the MM region. Alternatively, the method assumes that the wave function of the QC/MM system will tend to a Hartree product of the wave functions of the fragment and the remaining QC atoms. This assumption can be checked. For the tests performed in this article, it appears to be valid, and we have not seen any notable increase in difficulty in solving the QC equations for the system with scaled interactions. However, this might not always be so. In these cases, it might simply be a technical difficulty in obtaining a converged wave function or electron density, which could be tackled by using an improved initial

guess for the wave function, a better solver, or a constraint technique that ensures convergence to the desired solution. On the other hand, it might also indicate an inappropriate choice of buffer region where the separability assumption itself breaks down. This could be due, for example, to a long-range quantum mechanical effect that significantly perturbs the fragment wave function, even at the boundary of the MM region. 2.3. Implementation. The PAP and SISPA algorithms were implemented independently of each other in different developmental versions of the pDynamo molecular modeling program.13 The setup for simulations was common for both algorithms. As is standard for pDynamo QC/MM calculations, a full MM specification of the system is initially needed. This is then followed by the user identifying which atoms in the system are to be in the QC core (at least one), those that are always to be in the MM region (if any), and the remaining atoms that are allowed to transform adaptively. The latter are divided into chemically distinct user-defined fragmentsfor example, solvent molecules or counterionsthat transform as single entities. The setup is completed by specifying the geometry of the adaptive buffer region, which is currently limited to be a spherical shell around the weighted geometric center of the atoms in the QC core. The implementation of the PAP method itself was straightforward and was done using a pDynamo script, with no intervention required in the existing code. As a result, the PAP method can be employed with all of pDynamo’s existing QC/MM potentials. The algorithm works by replacing the existing energy and gradient function for a system with one which determines the number, N, of fragments in the buffer region and then sequentially generates and calculates the weighted sum of the energy and gradients of each of the 2N possible QC/MM partitions. No optimization of the method was attempted, with QC/MM calculations being performed independently for all PAP partitions, no matter what the PAP weight from eq 1. By contrast, to avoid unreasonable computational cost, it was essential at each energy evaluation to check the number of buffer fragments and abort the calculation if a preset maximum was exceeded. Implementation of the SISPA algorithm was also straightforward, but much more invasive than that of the PAP method as it required a selective scaling of the various interactions in the energy and gradient calculation. This was done for pDynamo’s semiempirical, density functional theory (DFT), and Hartree− Fock (HF) QC methods and for its MM potentials. For completeness, the scaling rules we used for this initial version of the algorithm were as follows: 1. All interactions involving nonfragment atoms are calculated as normal. 2. All interactions between atoms of the same fragment within the same representation are calculated as normal. 3. There are no interactions between atoms of the same fragment with different representations (i.e., the QC and MM representations of a fragment do not see each other). 4. For an n-body interaction involving fragments, the scaling is a product of n factors, using 1 for the scaling of nonfragment atoms, λI for a QC interaction, and 1 − λI for an MM interaction. 5. When an n-body interaction involves atoms from the same fragment, the scaling factors for those atoms are contracted. This means that terms of the form λkI , where k < n, are replaced by λI. Note that when all atoms in the interaction are in the same fragment (k = n), the scaling is 1 (from rule 2). 2344


Article


3. RESULTS This section outlines the systems and potentials that were employed for testing, before going on to describe the results of these tests. 3.1. Test Systems. The PAP and SISPA algorithms were tested on a number of molecular systems, including a variety of chlorohydrocarbon clusters and liquids, and aqueous solvated systems. The results of these tests, however, were broadly similar, and so in what follows attention is restricted to four pure water systemsa dimer, a cluster consisting of 22 molecules, a sphere of 1000 molecules with a radius of 19.277 Å, and a cubic box of 1000 molecules with side 31.074 Å. For the cubic system, periodic boundary conditions were employed, whereas for the sphere, piecewise half-harmonic restraining potentials were added between a central molecule and all the others in the system to prevent evaporation. These potentials had a force constant of 2000 kJ mol−1 Å−2 and were applied whenever the distance between molecule centers exceeded 18.777 (= 19.277−0.5) Å. In all tests, the flexible CHARMM TIP3P model with Lennard-Jones parameters on hydrogen was used for the MM potential,14,15 together with a variety of semiempirical, DFT, and HF QC methods. For the semiempirical approaches, we started with the full range of MNDO-like methods that pDynamo provides,13 most notably AM1,16 but subsequently implemented a number of others. These included AM1-FS1,17 which is like AM1 but includes dispersion and hydrogenbonding corrections, PM3-PIF,18 which gives good results for liquid water simulations, 19 and PMO1, 20,21 which has dispersion and an improved description of polarizability. For the DFT potentials, we employed two of pDynamo’s standard combinations,13 namely, (i) the LDA functional, a STO-3G orbital basis, and a DeMon fitting basis and (ii) the BP86 functional, a 3-21G orbital basis, and a DeMon fitting basis. In addition, we also developed an MM pseudo-QC method, which was useful for testing the implementation of the algorithms for example, for verifying that MM/MM PAP and SISPA simulations gave the same results as those with pure MM potentials. Finally, long-range MM/MM and QC/MM interactions were determined either without truncation, where possible, or using pDynamo’s standard atom-based forceswitching technique.13,22 For the vacuum and sphere systems, comparisons were made between both approaches to see how the results were affected. For the fixed-partitioning QC/MM calculations, one water was selected as being in the QC region. For the adaptive calculations, one water was assigned to the QC core, whereas the remaining waters were treated as fragments, which could migrate in and out of the QC core and buffer regions as necessary. All core water-fragment and interfragment distances were determined between the molecules’ weighted geometric centers, using atomic numbers as weights. 3.2. Basic Tests. Both the PAP and SISPA algorithms underwent extensive testing. This showed that they possessed consistent energies and forces, which was confirmed in a number of ways, including finite-difference calculation of the gradients from energies, energy minimization, and molecular dynamics simulations in the microcanonical ensemble using a velocity Verlet algorithm. To illustrate the latter, the rootmean-square (RMS) deviations in the total energy from 25 ps NVE simulations of the water box at 300 K with a pure MM potential were 7.4, 1.9, and 0.5 kJ mol−1 using timesteps of 1.0,

As in the PAP method, the buffer region fragments were identified at the beginning of each energy calculation and the scalings of each of the fragment (and nonfragment) atoms determined. The scaling rules outlined above were then applied directly to the different types of interaction appearing in the QC/MM energy expression as they were calculated. In the case of the QC electronic interactions, the scalings were applied to the one- and two-electron integrals over basis functions and, in the integration of the DFT exchange-correlation potential, to basis function products. We note that pDynamo uses atomcentered basis functions for all its QC methods, so each basis function takes on the scaling of the fragment to which its parent atom belongs. Explicit examples of scaling, in which all atoms involved in the interaction are from different fragments, are (1−λI)(1−λJ) for the MM/MM electrostatic and Lennard-Jones energies; λI(1 − λJ) for the QC/MM Lennard-Jones energy; λIλJ for the QC/ QC nuclear energy and the overlap and kinetic energy oneelectron integrals; λIλJλK for the three-center QC/QC electron−nuclear one-electron integrals; and λI λJ (1−λK) for the three-center QC/MM one-electron integrals involving basis functions from fragments I and J, and an MM charge from fragment K. Apart from scaling of the interactions, the remainder of SISPA’s QC/MM energy calculation proceeds in a standard fashion, with the solution of the QC equations to obtain the electronic energy and wave function and the computation of the system’s energy derivatives. Finally, the determination of the energy is completed by computing the weighted QC and MM fragment correction energies in vacuum that appear in Figure 2. Weighting functions of various forms were implemented and tested for both the PAP and SISPA algorithms. The tests in this paper were carried out with a fifth-order function that has continuous first and second derivatives at its boundaries: x = (r − rc)/(rb − rc) λ (x ) = 1

rc > r

= 1 − 10x 3 + 15x 4 − 6x 5 rb > r > rc = 0

r > rb

(3)

Here, r is the distance between the centers of the QC region and fragment and rc and rb are the distances that define the positions of the boundaries between the different regions (with rb > rc). To facilitate testing, it was found easier to accumulate separately terms of the form ∂E/∂λ and ∂E/∂(1 − λ) in the interaction energy code. These could then be postprocessed using the exact form of the weighting function and added into the full system coordinate derivatives. As a final remark, we note that all adaptive QC/MM methods that seek to have a continuous energy function must determine the energies of a fragment in its different representations with respect to a consistent zero of energy. This ensures, for example, that the fragments in the system do not partition themselves artificially due to one representation having a lower, but arbitrary, value of the energy than another. In this work, for both the PAP and SISPA algorithms, we always calculate a fragment’s energy in a given QC or MM representation with respect to its equivalent geometryoptimized energy in vacuum. 2345


Article


Figure 3. Potential energy curves for dissociation of a water dimer using the AM1 (top), LDA (middle), and BP86 (bottom) QC methods. The SISPA and PAP calculations were performed with the appropriate QC method for the panel, and with the buffer boundaries (inner/outer) that are specified in the legends. Also shown in each panel are the curves with fixed QC/MM partitioning and with the pure QC and MM potentials.

configuration was chosen for the comparison PAP calculations because this was the largest that was feasible computationally for all the test systems employed in this paper. It is evident from the figures that the rather simple individual QC and MM potentials used for the tests show a large variability in how they treat the energetics of the systems. As a consequence, the resulting adaptive QC/MM curves also display significant variability that depends upon the buffer geometry, and so upon how the separate potentials are mixed. Nevertheless, all the curves are quite smooth, and it is noteworthy that the PAP and SISPA algorithms with the 3−4 Å buffer geometry have very similar behavior. We also analyzed the electronic densities, notably the atomic charges, arising from the QC and QC/MM calculations (data not shown). These show the expected behavior for both the PAP and SISPA methods, in that the values of the charges on the atoms vary smoothly between those typical of the QC method in the QC core and the constant values of the MM potential in the MM environment. For the PAP algorithm, the

0.5, and 0.25 fs, respectively. The equivalent values for these RMS deviations were 7.1, 1.8, and 0.5 kJ mol−1 for a fixedpartitioning AM1/MM simulation; 7.6, 1.8, and 0.4 kJ mol−1 for a PAP AM1/MM simulation with a buffer region between 3 and 4 Å; and 7.0, 1.8, and 0.4 kJ mol−1 for a SISPA AM1/MM simulation with the same buffer geometry. In all cases, plots of the total energy showed no discernible drift. Examples of potential energy curves for two dissociation processes using the SISPA and PAP algorithms are shown in Figures 3 and 4. These show, respectively, the curves that are generated when the water dimer’s oxygen−oxygen distance is increased, and when the 22-molecule water cluster is expanded by displacing simultaneously and equally the enveloping 21 waters away from the central water in the cluster. Results are given for three different QC methods (AM1, LDA, and BP86), and with three different SISPA buffer configurations. Also shown in each panel are the results of PAP calculations with a 3−4 Å buffer, with a fixed partitioning QC/MM scheme, and with the pure QC and MM potentials. The 3−4 Å buffer 2346


Article


Figure 4. Potential energy curves for expansion of the 22-molecule water cluster using the AM1 (top), LDA (middle), and BP86 (bottom) QC methods. The expansion was performed by simultaneously and equally displacing 21 molecules away from a central molecule along the oxygen− oxygen direction vectors. See the legend to Figure 3 for further details.

semiempirical AM1 potential was used as the QC method in all cases. The data in Figure 3 indicate that the AM1/AM1 and AM1/ MM water−water interactions are less strong than the MM/ MM one. This can be seen by comparing the pure MM and fixed partitioning AM1/MM results in Figure 5, which show that the oxygen population in the first solvation shell in the AM1/MM simulation is slightly lower than in the MM case. It is even more evident in the adaptive SISPA and PAP simulations, as in these the QC core and buffer regions become almost devoid of water molecules. Clearly, when the molecules have the freedom to choose their own representation, it is energetically more favorable for them to be pure MM, than to be either AM1 or a mixture of both. As a final remark about Figure 5, we note that the SISPA and PAP results are very similar. Their costs though were very dissimilar, with the PAP simulation being significantly more expensive (∼60). Analysis of the fragment population in the buffer showed that it averaged about 5, although it could go up

charges of the fragments in the buffer tend toward those of the MM potential as the MM environment is approached because the weights of partitions with these fragments as MM are higher. The same occurs for the SISPA algorithm but for a different reason. In this case, the charges of the atoms in a fragment close to the MM environment tend to those of the isolated fragment in vacuum due to a decoupling of the fragment from its environment by the scaling factors. These charges are then removed by the charges coming from the QC fragment correction, leaving only those of the MM potential. 3.3. Condensed Phase Tests. Condensed phase tests were carried out with the water sphere and water box systems. Radial distribution functions (RDFs) and neighbor functions (NFs) resulting from some initial water box molecular dynamics simulations of 100 ps duration at 300 K with a Langevin thermostat and a 1 fs time step are plotted in Figure 5. These calculations were done with the SISPA and PAP algorithms with a 3−4 Å buffer geometry, with a fixed partitioning QC/ MM approach, and with the pure MM potential. The 2347


Article


Figure 5. Oxygen−oxygen radial distribution (top) and integrated neighbor (bottom) functions from periodic water box simulations. The calculations were performed between the oxygen of the water in the QC core and all other oxygens in the system. The SISPA and PAP calculations were done with the AM1 QC method with a buffer region between 3 and 4 Å, shown as the dotted vertical lines. Also plotted are the fixed partitioning AM1/MM and pure MM results.

As can be seen, the QC/MM SISPA distributions show the same behavior as in Figure 5, with the buffer region almost empty of water molecules. Despite this, there is a gradual improvement in the NFs on going from AM1 to PM3-PIF, which just modifies the QC/QC nuclear interactions, and then to PMO1, which adds dispersion and has an improved electronic treatment of polarization. At long range, PMO1/ MM has an NF that agrees well with the MM results, although it significantly overestimates the water density at short range. As a mimic of these results, Figure 6 also shows the curves obtained from two SISPA simulations with MM/MM potentials. In both cases, the TIP3P MM potential is used for the QC and MM regions, but in the MMa/MM simulation, the QC/MM electrostatic (charge−charge) interactions are scaled by 0.5, whereas in the second MMb/MM simulation, the baseline energy of the QC representation of each fragment is penalized by 10 kJ mol−1. The first modification leads to a depopulation of the buffer region, whereas the second causes a

to as much as 11, equivalent to 32 and 2048 QC/MM partitions, respectively. These problems with AM1, and equivalent semiempirical approaches, are well-known. For example, pure QC simulations of liquid water with AM119 gave a heat of vaporization of 32.1 kJ mol−1, compared to the 43.7 kJ mol−1 obtained with the TIP3P MM potential.14 As an attempt to alleviate them, we tested a number of more recent MNDO-like semiempirical methods with improved intermolecular interactions. The results for two of these, PM3-PIF18 and PMO1,20,21 are shown in Figure 6, along with results with the AM1 method. All three of these simulations were performed with the SISPA algorithm, and a larger 4−6 Å buffer geometry. The latter is more realistic than the 3−4 Å buffer geometry of Figure 5, as it moves the QC−MM transition zone away from the first solvation shell of the QC core’s water molecule. Also shown are the results of two SISPA simulations with MM/MM potentials (to be described below), and the standard MM curves. 2348


Article


Figure 6. Oxygen−oxygen radial distribution (top) and integrated neighbor (bottom) functions from periodic water box simulations. The calculations were performed between the oxygen of the water in the QC core and all other oxygens in the system. SISPA calculations were done with the AM1, PM3-PIF, and PMO1 semiempirical QC methods, with a buffer region between 4 and 6 Å, shown as the dotted vertical lines. Also plotted are the curves arising from two SISPA MM/MM simulations (see the text for details), and from a standard MM simulation.

marked reduction in the number of QC fragments in both the QC core and buffer regions. Insights into the origins of these results can be had by analyzing a fragment’s energies and forces during a simulation, although in what follows we shall concentrate on the energies as the force analysis gives much the same picture. Decomposition of the total QC/MM energy into fragment contributions is straightforward with the semiempirical QC and MM potentials that are employed here. The only three-center term is the MM angle, but this is intrafragment and so can be assigned directly to that fragment’s energy. The remaining terms are either one-center or two-center, with the latter being divided equally between the participating atoms. Figure 7 shows the averages of various fragment energy components as a function of the distance from the central water in the QC core. For clarity, the contributions from the latter are omitted. The values were obtained by calculating, with the appropriate QC/MM potential, the fragment energies of

structures selected at 100 fs intervals from a water box simulation with a pure MM potential of 200 ps duration (making 2001 structures in all). The fragment contributions were then histogrammed with respect to distance from the center of the QC core and the different values averaged within each histogram bin. We note that we chose configurations from an MM simulation, so as to ensure a proper water distribution. The top panel of Figure 7 shows the energy components arising from a SISPA calculation with a 4−6 Å buffer geometry and an MM/MM potential. We have verified that the total energy in this plot is identical to that obtained with a pure MM potential, and we note that its value is basically constant at distances greater than about 3 Å, indicative of a uniform water distribution. The bottom panel of the figure gives the QC (both intra- and interfragment interactions) and QC/MM electrostatic energy components for three semiempirical QC/MM potentials, in addition to those for the MM/MM potential. We only plot the QC and QC/MM electrostatic terms as the 2349


Article


Figure 7. Average fragment energy components as a function of distance from the center of the QC core region. The averages were obtained using 2001 structures from a water box simulation with a pure MM potential and were calculated with the SISPA algorithm with a 4−6 Å buffer geometry (shown as the dotted vertical lines). Top panel: the values arising from a MM/MM potential. The MM and QC fragment energies comprise the intrafragment and fragment correction energies, whereas the MM/MM and QC/QC interaction energies are the sums of both the interfragment electrostatic and Lennard-Jones energies. Bottom panel: the QC and QC/MM electrostatic interaction energies for four different QC/MM potentials. Only these two contributions are shown as the MM fragment, QC/MM Lennard-Jones, and MM/MM interaction terms are the same as those in the top panel. The MM/MM QC/MM term is identical to the QC/MM electrostatic term of the top panel, whereas the MM/MM QC term is equal to the sum of the QC fragment and QC/QC interaction terms of the top panel.

remaining terms are pure MM and, so, are identical to those of the top panel. Inspection of the curves in the bottom panel provides evidence of the behavior of the semiempirical methods. The AM1 QC fragment energy is always positive (i.e., destabilizing), and its QC/MM electrostatic energy is only about two-thirds that of the MM/MM potential. Modification of the QC/QC nuclear terms in the PM3-PIF potential improves somewhat the latter’s QC results (although they are still largely destabilizing) but leaves the QC/MM energies essentially unchanged. By contrast, the enhanced treatment of polarization in PMO1 means that its QC/MM electrostatic energy mirrors very closely that of the MM/MM potential, whereas the inclusion of both polarization and dispersion makes PMO1’s QC fragment energy mostly stabilizing, except in the crucial region toward the center of the buffer. This difference between the QC values of the PMO1 and MM/MM potentials is sufficient to cause the distortion in the water distribution observed in Figure 6. Clearly these results are unsatisfactory, and it would be unwise to use these QM/MM potentials in combination with SISPA in most applications. We have made substantial attempts to modify the separate QC and MM potentials and their coupling to overcome the problems that are observed. Indeed, in individual cases, we can obtain improved results by, for example, scaling various energy contributions coming from an energy calculation (similar to the MM/MM curves in Figure 6). Unfortunately, though, we have been unable to find a physically

convincing and reliable general strategy. As a consequence, we prefer to leave further discussion of these points to a future paper, in which we will employ alternative and more flexible semiempirical, DFT, and ab initio QC methods and MM potentials, and consider additional approaches, such as specific parametrization procedures.

4. CONCLUSIONS In this paper, the SISPA algorithm for adaptive QC/MM simulations has been proposed. Compared to previous algorithms, it has the advantage of having a formulation that gives consistent energy and forces, while only requiring a single QC/MM energy calculation. This makes it not significantly different in accuracy, cost, and utility to QC/MM calculations with fixed partitioning. The method does, however, have disadvantages. First, its implementation requires a nonnegligible, although straightforward, modification of the energy-calculating code of the QC/MM program, and second, it is still necessary to do a QC calculation that encompasses all the atoms in the QC core and buffer regions of the QC/MM system. This might be too expensive in some situations (for example, when the buffer region is very large), and so alternative algorithms, such as that of Watanbe et al.,23 might be preferable. The results of initial tests of the SISPA algorithm, described in section 3, are promising, but clearly there is significant scope for future developmental work. It is apparent that the extra 2350


Article


(11) 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, 2231−2241. (12) Bulo, R. E.; Ensing, B.; Sikkema, J.; Visscher, L. Toward a Practical Method for Adaptive QM/MM Simulations. J. Chem. Theory Comput. 2009, 5, 2212−2221. (13) Field, M. J. The pDynamo Program for Molecular Simulations using Hybrid Quantum Chemical and Molecular Mechanical Potentials. J. Chem. Theory Comput. 2008, 4, 1151−1161. (14) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Water. J. Chem. Phys. 1983, 79, 926−935. (15) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W., III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (16) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. Development and Use of Quantum Mechanical Molecular Models. 76. AM1: a New General Purpose Quantum Mechanical Molecular Model. J. Am. Chem. Soc. 1985, 107, 3902−3909. (17) Foster, M. E.; Sohlberg, K. A New Empirical Correction to the AM1Method for Macromolecular Complexes. J. Chem. Theory Comput. 2010, 6, 2153−66. (18) Bernal-Uruchurtu, M. I.; Martins-Costa, M. T. C.; Millot, C.; Ruiz-López, M. F. Improving Description of Hydrogen Bonds at the Semiempirical Level: Water−Water Interactions as Test Case. J. Comput. Chem. 2000, 21, 572−581. (19) Monard, G.; Bernal-Uruchurtu, M. I.; van der Vaart, A.; Merz, K. M., Jr; Ruiz-López, M. F. Simulation of Liquid Water using Semiempirical Hamiltonians and the Divide and Conquer Approach. J. Phys. Chem. A 2005, 109, 3425−3432. (20) Zhang, P.; Fiedler, L.; Leverentz, H. R.; Truhlar, D. G.; Gao, J. Polarized Molecular Orbital Model Chemistry. 2. The PMO Method. J. Chem. Theory Comput. 2011, 7, 857−867. (21) Zhang, P.; Fiedler, L.; Leverentz, H. R.; Truhlar, D. G.; Gao, J. Erratum: Polarized Molecular Orbital Model Chemistry. 2. The PMO Method. J. Chem. Theory Comput. 2012, 8, 2983−2983. (22) Field, M. J.; Albe, M.; Bret, C.; Proust-De Martin, F.; Thomas, A. The Dynamo Library for Molecular Simulations using Hybrid Quantum Mechanical and Molecular Mechanical Potentials. J. Comput. Chem. 2000, 21, 1088−1100. (23) 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, 4242−4252. (24) Pezeshki, S.; Lin, H. Adaptive-partitioning Redistributed Charge and Dipole Schemes for QM/MM Dynamics Simulations: On-the-fly Relocation of Boundaries that Pass Through Covalent Bonds. J. Chem. Theory Comput. 2011, 7, 3625−3634.

freedom permitted by adaptive algorithms in allowing the transformation of parts of a system between QC and MM representations places extra demands on the quality of the QC/ MM coupling, and on ensuring the compatibility of the separate QC and MM potentials. Likewise, for the SISPA algorithm itself, it would be worthwhile to investigate alternative fragment correction and scaling procedures, to extend the method to be able to handle covalently bound fragments,24 and to extend it to other types of multiscale simulation. Finally, the availability of robust, fast, and accurate adaptive QC/MM algorithms should facilitate the simulation of a range of molecular problems that hitherto have been difficult to study effectively.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: martin.fi[email protected]. ORCID

Martin J. Field: 0000-0001-8674-7997 Notes

The author declares no competing financial interest.

■

ACKNOWLEDGMENTS The author would like to thank Anirban Bhattacharjee, Ramon Crehuet, Thomas Hofer, and Melchor Sanchez for some helpful discussions about adaptive QC/MM methods.

■

REFERENCES

(1) Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme. J. Mol. Biol. 1976, 103, 227−249. (2) Singh, U. C.; Kollman, P. A. A Combined ab initio Quantum Mechanical 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, 718−730. (3) 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, 700−733. (4) Liu, M.; Wang, Y.; Chen, Y.; Field, M. J.; Gao, J. QM/MM through the 1990s: the First Twenty Years of Method Development and Applications. Isr. J. Chem. 2014, 54, 1250−1263. (5) Zheng, M.; Waller, M. P. Adaptive Quantum Mechanics/ Molecular Mechanics Methods. WIREs Comput. Mol. Sci. 2016, 6, 369−385. (6) Pezeshki, S.; Lin, H. Recent Developments in QM/MM Methods Towards Open-boundary Multi-scale Simulations. Mol. Simul. 2015, 41, 168−189. (7) Jiang, T.; Boereboom, J. M.; Michel, C.; Fleurat-Lessard, P.; Bulo, R. E. Proton Transfer in Aqueous Solution: Exploring the Boundaries of Adaptive QM/MM. In Quantum Modeling of Complex Molecular Systems; Rivail, J.-L., Ruiz-Lopez, M., Assfeld, X., Eds.; Springer International Publishing: Switzerland, 2015; Vol. 21, pp 51−91. (8) Pezeshki, S.; Lin, H. Recent Progress in Adaptive-partitioning QM/MM Methods for Born-Oppenheimer Molecular Dynamics. In Quantum Modeling of Complex Molecular Systems; Rivail, J.-L., RuizLopez, M., Assfeld, X., Eds.; Springer International Publishing: Switzerland, 2015; Vol. 21, pp 93−113. (9) Hofer, T. S.; Pribil, A. B.; Randolf, B. R.; Rode, B. M. Ab initio Quantum Mechanical Charge Field Molecular Dynamics  a Nonparametrized First-principle Approach to Liquids and Solutions. Adv. Quantum Chem. 2010, 59, 213−246. (10) Hofer, T. S.; Rode, B. M.; Pribil, A. B.; Randolf, B. R. Simulations of Liquids and Solutions based on Quantum Mechanical Forces. Adv. Inorg. Chem. 2010, 62, 143−175. 2351


MM Simulations - Journal of Chemical

Recommend Documents