Restrained Proton Indicator in Combined Quantum-Mechanics

Article pubs.acs.org/JPCB

Restrained Proton Indicator in Combined Quantum-Mechanics/ Molecular-Mechanics Dynamics Simulations of Proton Transfer through a Carbon Nanotube Adam W. Duster and Hai Lin* Chemistry Department, CB 194, University of Colorado Denver, Denver, Colorado 80217, United States S Supporting Information *

ABSTRACT: Recently, a collective variable “proton indicator” was purposed for tracking an excess proton solvated in bulk water in molecular dynamics simulations. In this work, we demonstrate the feasibility of utilizing the position of this proton indicator as a reaction coordinate to model an excess proton migrating through a hydrophobic carbon nanotube in combined quantum-mechanics/molecular-mechanics simulations. Our results indicate that applying a harmonic restraint to the proton indicator in the bulk solvent near the nanotube pore entrance leads to the recruitment of water molecules into the pore. This is consistent with an earlier study that employed a multistate empirical valence bond potential and a different representation (center of excess charge) of the proton. We attribute this water recruitment to the delocalized nature of the solvated proton, which prefers to be in high-dielectric bulk solvent. While water recruitment into the pore is considered an artifact in the present simulations (because of the artificially imposed restraint on the proton), if the proton were naturally restrained, it could assist in building water wires prior to proton transfer through the pore. The potential of mean force for a proton translocation through the water-filled pore was computed by umbrella sampling, where the bias potentials were applied to the proton indicator. The free energy curve and barrier heights agree reasonably with those in the literature. The results suggest that the proton indicator can be used as a reaction coordinate in simulations of proton transport in confined environments.

1. INTRODUCTION A hydrated H+ can exist as a hydronium ion (H3O+) or a Zundel ion (H5O2+, or a H+ shared by two water molecules). The hydronium ion is often hydrogen-bonded by three water molecules to form an Eigen ion (H9O4+). A solvated H+ can diffuse quickly via a water wire formed by water molecules (or a proton wire formed by water molecules and titratable functional groups such as carboxyl groups) through reorganizing the network of covalent bonds and hydrogen bonds. This H+ “hopping” mechanism is known as the Grotthuss mechanism,1,2 which is crucial to many biological and chemical processes. The Grotthuss mechanism presents a challenge in molecular dynamics simulations of H+ solvation and transfer. Although many reactive molecular mechanical (MM) models3−12 and multistate empirical valence bond (MS-EVB) models13−21 have been developed and applied to simulate H+ solvation and transport, the involved breaking and forming covalent bonds are ideally described by quantum mechanics (QM). This is especially the case when various titratable groups are part of the H+ relay path, for which elaborate and time-consuming parametrizations of (specific) MM models may be required. However, QM calculations are more expensive, thus limiting the QM simulations to relatively small model systems (tens to © XXXX American Chemical Society

hundreds of atoms, depending on the chosen level of theory) and to relatively short simulation times (usually less than 100 ps). Combined QM/MM22−26 can reduce computational costs, but special treatments of the QM/MM boundary are needed to ensure that the diffusing excess H+ and its solvation shell remain in the QM region. This can be done in two ways: (1) by imposing necessary restraints on the QM/MM boundary to prevent QM atoms from moving into the MM region27,28 or (2) by using an adaptive QM/MM29−50 boundary to relocate the QM region on the fly following the hopping H+ wherever it goes.43 Both methods have their own merits and pitfalls. The first “restraint” approach is straightforward and easy to implement, but the restraints at the QM/MM boundary may restrict the time lengths and scopes of the simulations. On the other hand, the second “adaptive” approach is more involved but permits the use of a smaller, mobile QM subsystem, in principle allowing for simulations of much longer time scales. When using the adaptive-partitioning QM/MM algorithms for H+ transfer, it is necessary to track the delocalized H+ on the fly. The simplest solution is to represent the H+ by the donor D Received: July 6, 2017 Revised: August 18, 2017 Published: August 18, 2017 A

DOI: 10.1021/acs.jpcb.7b06657 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

enon observed by Peng et al.57 and to compare the potential of mean force via umbrella sampling for H+ transfer through the water-filled pore with previous calculations.

(the hydronium O atom), but the abrupt changes in the donor identity make it difficult to update the QM subsystem in a gradual and smooth manner. Therefore, an explicit H+ indicator was introduced, the position of which represents the approximate location of the excess H+. Its coordinates are a function of the coordinates of D, the H atoms bonded to D, and the possible acceptors (the O atoms of the water molecules within a threshold distance rLIST from D). The list of possible acceptors is updated at every step, and at a given time, the “most likely” acceptor is denoted A. In the limit of Eigen ion, the indicator superimposes with D, whereas in the limit of the Zundel ion, it superimposes with the shared H atom at the midpoint between the two water O atoms. When the structure varies between the Eigen and Zundel limits, the indicator’s position changes accordingly. When proton transfer occurs, A and D switch identities; the switch is triggered if the shared H+ becomes closer to A than D, as indicated by the value of a parameter ρ becoming larger than a threshold value (the definition of ρ is given in the next section). These changes should lead to minimal perturbation to the indicator’s position, which is necessary for smooth variations in the indicator’s position. The H+ indicator has been successfully applied in adaptive-partitioning QM/MM simulations of an excess H+ in bulk water, serving as the center of the mobile QM subsystem.43 There are other forms of collective variables for depicting a delocalized excess H+ in the literature. For example, Pomès and Roux51 tracked the location of excess H+ using the center of excess charge (CEC), the position of which is given by projecting the total dipole moment of the water wire on the axis along the water wire (taken to be the z-axis). Chakrabarti et al.52 adopted an average z coordinate over the O and H atoms in a single-file water wire, weighted by the number of excess H atom in close vicinity of each O atom. Integrating the above two schemes, König et al.53 proposed a new collective variable called modified CEC (mCEC), which was generalized to threedimensional space. The feasibility of using the position of mCEC as a reaction coordinate in H+ translocation has been nicely demonstrated in the simulations of a variety of model systems.53−56 An alternative definition of CEC was used by Day et al.16 in the MS-EVB model for H+ transport in aqueous solution, where the CEC is expressed as the sum of the “pivot” hydronium centers of involved valence states, each weighted by the EVB amplitude of that state. Peng et al.57 recently applied the MS-EVB/CEC to study H+ migration through a carbon nanotube. They concluded that water molecules move through a hydrated excess H+ charge defect that resides at the interface between bulk water and a hydrophobic pore in order to wet the path ahead for subsequent H+ charge migration. As the wetting of the hydrophobic pore helps to lower the barrier of the H+ translocation, the authors proposed that this could be an important mechanism for H+ transport in biological systems. Given the success of employing the H+ indicator as the reference for on-the-fly relocation of the QM subsystems in adaptive QM/MM simulations, it is natural to ask if the H+ indicator, just like the CEC and mCEC, can also be utilized as a reaction coordinate in biased simulations. This requires the application of restraint potentials on the H+ indicator, which has not been tested in dynamics simulations previously. In this work, we attempt to answer this question by imposing a harmonic restraint on the H+ indicator in QM/MM simulations of an excess H+ moving through a confined nanotube. Specifically, we want to verify the “water-shooting” phenom-

2. THEORY AND CALCULATIONS Details about the H+ indicator can be found in ref 43. Here we only give a brief introduction. As illustrated in Figure 1, the

Figure 1. Position of the H+ indicator (I, highlighted in green) is a function of the coordinates of the donor (D, highlighted in yellow) and all possible acceptors (Aj, j = 1, 2, ..., J, one of which is highlighted in blue) that are within an enlisting cutoff distance rLIST from D. The H atoms covalently bonded to D are denoted Hm (m = 1, 2, ..., M, one of which is highlighted in red). The topology (bonding connectivity) is updated on the fly in dynamics simulations.

position of the indicator (I) is a linear combination of the coordinates of the donor (D) and possible acceptors (Aj, j = 1, 2, ..., J) surrounding the donor: XI =

1 (X D + gI J

gI = 1 +

J

M

∑ ∑ g(ρmj )X A ) j

(1)

j=1 m=1 M

∑ ∑ g(ρmj ) (2)

j=1 m=1

Here, XI, XD, and XAj are the Cartesian coordinates of I, D, and Aj, respectively, M the total number of hydrogen atoms Hm bonded to D, and g(ρmj) a weight function that depends on ρmj, a ratio indicating how close Hm is to D versus to Aj, and gI is the normalization factor. The ratio ρmj is defined using projected donor−acceptor vectors:58 rDHm·rDA j ρmj = |rDA j|2 (3) where rDHm = XHm − XD and rDAj = XAj − XD. The weight function is chosen to be a fifth-order polynomial that depends on a reduced variable x: ⎧0 if 1 ≤ x ⎪ 5 4 3 g (x) = ⎨−6x + 15x − 10x + 1 if 0 ≤ x < 1 ⎪ ⎩1 if x < 0 x = x(ρmj ) = 1 −

(4)

ρmj − ρmj0 0.5 − ρmj0

(5)

By design, g(ρmj) should be 1 at ρmj = 0.5, which corresponds to equal sharing of a H+ in the Zundel structure, and fall quickly to B


Article

The Journal of Physical Chemistry B 0 as ρmj decreases approaching ρ0mj = r0DH/rLIST, where r0DH = 1.00 Å is a parameter taken to be slightly larger than the equilibrium OH bond length in the hydronium ion. The smooth and fast decrease of g(ρmj) to 0 ensures that XI is dominated by the XD and XA (A denotes the most likely acceptor with the largest ρmj) and that changes in XI are minimal when D and A switch. Any forces that act on I are decomposed proportionally to D and Aj. A single-walled carbon nanotube model is employed in the present study. The nanotube was placed along the z-axis (Figure 2), which is 28.2 Å long with a diameter of 8 Å. On

out under the NpT ensemble (p = 1 bar and T = 300 K). The x and y dimensions were fixed, and volume changes due to pressure only occurred along the z dimension. Next, the model system was equilibrated by 5 ns under the NVT ensemble at 300 K. A Nosé−Hoover thermostat62 was employed for temperature control and a Langevin-piston63 barostat for pressure control. The C atoms were held fixed. As expected, water molecules refrained from entering deep into the pore, with typically 0−4 water molecules seen in the pore near each entrance; a water molecule is classified as in the pore if its O atom is within the pore. Occasionally (0.5 Å) when the donor and acceptor switch identities. This happens when a third water molecule is in rather close distance from the acceptor before the switching (i.e., the donor after the switching), with the arrangement of these three water molecules in a geometry such that concerted H+ hopping becomes possible. Figures 6a−e illustrate such long-distance

Figure 5. Potential of mean force for proton transfer through the nanotube from the umbrella sampling using the dumbbell model. Because of the symmetry, sampling was carried out for only half of the proton translocation path. Figure 6. Long-distance movement of the H+ indicator during converted H+ transfer. The donor and acceptor are highlighted in yellow and blue, respectively, and the indicator is represented by a green circle.

pore of nanotube. Oscillations centered on the umbrella windows in the curve are clearly visible. These “bumps” are most likely due to the restrains imposed on the pore water O atoms, as the proton is most likely spending its time in the Zundel state formed by two adjacent restrained water molecules. Nevertheless, overall the potential of mean force shows a gradually increasing free energy profile for translocation of one H+ from the bulk solvent into the pore. The barrier height of 15 kcal/mol is compared favorably with those in earlier calculations.51,53,57,71 For example, electrostatic calculations by Burykin and Warshel yielded a barrier of 15 kcal/mol.71 Pomès and Roux51 determined a barrier of 8 kcal/ mol for proton translocation along a water wire described by a polarizable and dissociable MM water model PM6.3,72 The authors attributed the barrier to the reorientation of the water molecules in the chain (propagation of the bonding defect), as the translocation of the H+ (propagation of the ionic defect) is essentially activationless. Using a simplified dumbbell model, König et al.53 carried out dynamics simulations employing the self-consistent charge density functional tight binding (SCC-

migration if concerted H+ hopping occurs. As the donor changes identities from the left water to the central water (Figures 6b,c), the indicator moves from near the shared H between the left and central water molecules to close to the O atoms of the central water because both the left and right water molecules are contributing roughly the same weights in eq 1. Similarly, the indicator will move across a long distance again when the donor is switched from the central to the right water (Figures 6c,d). The long-distance movements of the indicator appear to be in line with the concerted H+ migration process.

4. SUMMARY In this contribution, we investigated the restrained H+ indicator in combined QM/MM molecular dynamics simulations of H+ transport through a hydrophobic C nanotube, where the indicator represented the hydrated H+. First, our simulations E



■

reveal that a restrained H+ near the nanotube entrance can recruit water molecules into the hydrophobic pore, confirming the recent observations by Peng et al.,57 who employed totally different computational methods. We attribute this to the delocalization nature of the solvated H+, which prefers to reside in high dielectric medium (the bulk solvent). Because the restraint imposed on the excess H+ prevents the H+ from moving away from the interface between the low- and highdielectric media, the water molecules over which the H+ is delocalized are led into the pore to form a more favorable solvation environment for the H+. If the H+ is restrained naturally, the water-recruitment mechanism may be important in preparing the water wire for H+ transfer through the pore. Second, potential of mean force for H+ translocation through the water-filled pore was computed by umbrella sampling, where the bias potentials were applied to the H+ indicator. The estimated free energy curve and barrier heights agree reasonably well with those in the literature. These results suggest that the restrained H+ indicator can be used in QM/ MM simulations of H+ transport in confined environments.

■

REFERENCES

(1) Agmon, N. The Grotthuss mechanism. Chem. Phys. Lett. 1995, 244, 456−462. (2) Marx, D. Proton transfer 200 years after von Grotthuss: Insights from ab initio simulations. ChemPhysChem 2006, 7, 1848−1870. (3) Stillinger, F. H.; David, C. W. Polarization model for water and its ionic dissociation products. J. Chem. Phys. 1978, 69, 1473−1484. (4) David, C. W. A variable charge central force model for water and its ionic dissociation products. J. Chem. Phys. 1996, 104, 7255−7260. (5) Halley, J. W.; Rustad, J. R.; Rahman, A. A polarizable, dissociating molecular dynamics model for liquid water. J. Chem. Phys. 1993, 98, 4110−4119. (6) Billeter, S. R.; van Gunsteren, W. F. Protonizable water model for quantum dynamical simulations. J. Phys. Chem. A 1998, 102, 4669− 4678. (7) Lussetti, E.; Pastore, G.; Smargiassi, E. A fully polarizable and dissociable potential for water. Chem. Phys. Lett. 2003, 381, 287−291. (8) Mahadevan, T. S.; Garofalini, S. H. Dissociative water potential for molecular dynamics simulations. J. Phys. Chem. B 2007, 111, 8919− 8927. (9) Fogarty, J. C.; Aktulga, H. M.; Grama, A. Y.; van Duin, A. C. T.; Pandit, S. A. A reactive molecular dynamics simulation of the silicawater interface. J. Chem. Phys. 2010, 132, 174704. (10) Selvan, M. E.; Keffer, D. J.; Cui, S.; Paddison, S. J. A reactive molecular dynamics algorithm for proton transport in aqueous systems. J. Phys. Chem. C 2010, 114, 11965−11976. (11) Kale, S.; Herzfeld, J.; Dai, S.; Blank, M. Lewis-inspired representation of dissociable water in clusters and Grotthuss chains. J. Biol. Phys. 2012, 38, 49−59. (12) Wolf, M. G.; Groenhof, G. Explicit proton transfer in classical molecular dynamics simulations. J. Comput. Chem. 2014, 35, 657−671. (13) Sagnella, D. E.; Tuckerman, M. E. An empirical valence bond model for proton transfer in water. J. Chem. Phys. 1998, 108, 2073− 2083. (14) Schmitt, U. W.; Voth, G. A. Multistate empirical valence bond model for proton transport in water. J. Phys. Chem. B 1998, 102, 5547−5551. (15) Vuilleumier, R.; Borgis, D. An extended empirical valence bond model for describing proton transfer in H+(H2O)n clusters and liquid water. Chem. Phys. Lett. 1998, 284, 71−77. (16) Day, T. J. F.; Soudackov, A. V.; Č uma, M.; Schmitt, U. W.; Voth, G. A. A second generation multistate empirical valence bond model for proton transport in aqueous systems. J. Chem. Phys. 2002, 117, 5839− 5849. (17) Kornyshev, A. A.; Kuznetsov, A. M.; Spohr, E.; Ulstrup, J. Kinetics of proton transport in water. J. Phys. Chem. B 2003, 107, 3351−3366. (18) Brancato, G.; Tuckerman, M. E. A polarizable multistate empirical valence bond model for proton transport in aqueous solution. J. Chem. Phys. 2005, 122, 224507. (19) Swanson, J. M. J.; Maupin, C. M.; Chen, H.; Petersen, M. K.; Xu, J.; Wu, Y.; Voth, G. A. Proton solvation and transport in aqueous and biomolecular systems: Insights from computer simulations. J. Phys. Chem. B 2007, 111, 4300−4314. (20) Wu, Y. J.; Chen, H. N.; Wang, F.; Paesani, F.; Voth, G. A. An improved multistate empirical valence bond model for aqueous proton solvation and transport. J. Phys. Chem. B 2008, 112, 467−482. (21) Park, K.; Lin, W.; Paesani, F. A refined MS-EVB model for proton transport in aqueous environments. J. Phys. Chem. B 2012, 116, 343−352. (22) 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. (23) 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. (24) Gao, J. Methods and applications of combined quantum mechanical and molecular mechanical potentials. Rev. Comput. Chem. 1996, 7, 119−185.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b06657. Figure S1: the number of water molecules in the pore as a function of simulation time; Figure S2: compares how far the H+ is from the pore entrance in simulations with and without the restraint imposed on the H+ indicator (PDF) Video showing the recruitment of water molecules into the nanotube by a restrained H+ in the bulk near the pore entrance (MPG)

■

Article

AUTHOR INFORMATION

Corresponding Author

*(H.L.) E-mail [email protected]. ORCID

Hai Lin: 0000-0002-3525-9122 Author Contributions

H.L. designed the project. A.D. carried out the code implementation and calculations. A.D. and H.L. analyzed the results and wrote the manuscript. Both authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is supported by the NSF (CHE-1564349), Camille & Henry Dreyfus Foundation (TH-14-028), and NVIDIA Corporation. This work used XSEDE under Grant CHE140070, supported by NSF Grant ACI-1053575, and NERSC under Grant m2495. We thank Danielle Miller and Christina Garza for helpful discussions.

■

ABBREVIATIONS CEC, center of excess charge; mCEC, modified center of excess charge; MM, molecular mechanics; MS-EVB, multistate empirical valence bond; QM, quantum-mechanics; QM/MM, quantum-mechanics/molecular-mechanics; SCC-DFTB, selfconsistent charge density functional tight binding. F


Article


using many-body potentials. J. Chem. Theory Comput. 2016, 12, 3441− 3448. (46) Duster, A.; Garza, C.; Lin, H. Adaptive partitioning QM/MM dynamics simulations for substrate uptake, product release, and solvent exchange. In Computational Approaches for Studying Enzyme Mechanism; Voth, G. A., Ed.; Elsevier: Cambridge, MA, 2016; pp 342−358. (47) Zheng, M.; Waller, M. P. Adaptive quantum mechanics/ molecular mechanics methods. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2016, 6, 369−385. (48) Duster, A. W.; Wang, C.-H.; Garza, C. M.; Miller, D. E.; Lin, H. Adaptive quantum/molecular mechanics: what have we learned, where are we, and where do we go from here. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2017, 7, e1310. (49) Field, M. J. An algorithm for adaptive QC/MM simulations. J. Chem. Theory Comput. 2017, 13, 2342−2351. (50) Zheng, M.; Kuriappan, J. A.; Waller, M. P. Toward more efficient density-based adaptive QM/MM methods. Int. J. Quantum Chem. 2017, 117, e25336. (51) Pomès, R.; Roux, B. Free energy profiles for H+ conduction along hydrogen-bonded chains of water molecules. Biophys. J. 1998, 75, 33−40. (52) Chakrabarti, N.; Tajkhorshid, E.; Roux, B.; Pomès, R. Molecular basis of proton blockage in aquaporins. Structure 2004, 12, 65−74. (53) Konig, P. H.; Ghosh, N.; Hoffmann, M.; Elstner, M.; Tajkhorshid, E.; Frauenheim, T.; Cui, Q. Toward theoretical analyis of long-range proton transfer kinetics in biomolecular pumps. J. Phys. Chem. A 2006, 110, 548−563. (54) Riccardi, D.; Konig, P.; Prat-Resina, X.; Yu, H. B.; Elstner, M.; Frauenheim, T.; Cui, Q. ″Proton holes” in long-range proton transfer reactions in solution and enzymes: A theoretical analysis. J. Am. Chem. Soc. 2006, 128, 16302−16311. (55) Liang, R.; Swanson, J. M. J.; Madsen, J. J.; Hong, M.; DeGrado, W. F.; Voth, G. A. Acid activation mechanism of the influenza A M2 proton channel. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, E6955− E6964. (56) Lee, S.; Swanson, J. M. J.; Voth, G. A. Multiscale simulations reveal key aspects of the proton transport mechanism in the ClC-ec1 antiporter. Biophys. J. 2016, 110, 1334−1345. (57) Peng, Y.; Swanson, J. M. J.; Kang, S.-g.; Zhou, R.; Voth, G. A. Hydrated excess protons can create their own water wires. J. Phys. Chem. B 2015, 119, 9212−9218. (58) Hofer, T. S.; Hitzenberger, M.; Randolf, B. R. Combining a dissociative water model with a hybrid QM/MM approach-a simulation strategy for the study of proton transfer reactions in solution. J. Chem. Theory Comput. 2012, 8, 3586−3595. (59) Wu, Y.; Tepper, H. L.; Voth, G. A. Flexible simple point-charge water model with improved liquid-state properties. J. Chem. Phys. 2006, 124, 024503−12. (60) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (61) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (62) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (63) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. Constant pressure molecular dynamics simulation: The Langevin piston method. J. Chem. Phys. 1995, 103, 4613−4621. (64) Roux, B. The calculation of the potential of mean force using computer simulations. Comput. Phys. Commun. 1995, 91, 275−282. (65) Lin, H.; Zhang, Y.; Pezeshki, S.; Duster, A.; Truhlar, D. G. QMMM, Version 2.2.8.CO; University of Minnesota: Minneapolis, 2017. (66) Thiel, W. MNDO2005, Version 7.0; Max-Planck-Institut für Kohlenforschung: Müheim an der Ruhr, Germany, 2005.

(25) Lin, H.; Truhlar, D. G. QM/MM: What have we learned, where are we, and where do we go from here? Theor. Chem. Acc. 2007, 117, 185−199. (26) Senn, H. M.; Thiel, W. QM/MM methods for biological systems. Top. Curr. Chem. 2007, 268, 173−290. (27) Wu, X.; Thiel, W.; Pezeshki, S.; Lin, H. Specific reaction path Hamiltonian for proton transfer in water: Reparameterized semiempirical models. J. Chem. Theory Comput. 2013, 9, 2672−2686. (28) 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, 3526−3535. (29) Kerdcharoen, T.; Liedl, K. R.; Rode, B. M. A QM/MM simulation method applied to the solution of Li+ in liquid ammonia. Chem. Phys. 1996, 211, 313−323. (30) Kerdcharoen, T.; Morokuma, K. ONIOM-XS: an extension of the ONIOM method for molecular simulation in condensed phase. Chem. Phys. Lett. 2002, 355, 257−262. (31) 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. (32) 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. (33) Nielsen, S. O.; Bulo, R. E.; Moore, P. B.; Ensing, B. Recent progress in adaptive multiscale molecular dynamics simulations of soft matter. Phys. Chem. Chem. Phys. 2010, 12, 12401−12414. (34) 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. (35) 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. (36) 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, 5567−5577. (37) Pezeshki, S.; Davis, C.; Heyden, A.; Lin, H. Adaptivepartitioning QM/MM dynamics simulations: 3. Solvent molecules entering and leaving protein binding sites. J. Chem. Theory Comput. 2014, 10, 4765−4776. (38) Pezeshki, S.; Lin, H. Recent developments in QM/MM methods towards open-boundary multi-scale simulations. Mol. Simul. 2015, 41, 168−189. (39) Waller, M. P.; Kumbhar, S.; Yang, J. A density-based adaptive quantum mechanical/molecular mechanical method. ChemPhysChem 2014, 15, 3218−3225. (40) 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. (41) Böckmann, M.; Doltsinis, N. L.; Marx, D. Adaptive switching of interaction potentials in the time domain: An extended lagrangian approach tailored to transmute force field to QM/MM simulations and back. J. Chem. Theory Comput. 2015, 11, 2429−2439. (42) 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: Cham, 2015; pp 51−91. (43) Pezeshki, S.; Lin, H. Adaptive-partitioning QM/MM for molecular dynamics simulations: 4. Proton hopping in bulk water. J. Chem. Theory Comput. 2015, 11, 2398−2411. (44) Pezeshki, S.; Lin, H. Recent developments in adaptive QM/ MM. In Quantum Modeling of Complex Molecular Systems; Rivail, J.-L., Ruiz-Lopez, M., Assfeld, X., Eds.; Springer: New York, 2015; pp 93− 113. (45) Boereboom, J. M.; Potestio, R.; Donadio, D.; Bulo, R. E. Toward Hamiltonian adaptive QM/MM: Accurate solvent structures G


Article

The Journal of Physical Chemistry B (67) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. The nature of the hydrated excess proton in water. Nature 1999, 397, 601− 604. (68) Marx, D.; Tuckerman, M. E.; Parrinello, M. Solvated excess protons in water: quantum effects on the hydration structure. J. Phys.: Condens. Matter 2000, 12, A153−A159. (69) Bakowies, D.; Thiel, W. Hybrid models for combined quantum mechanical and molecular mechanical approaches. J. Phys. Chem. 1996, 100, 10580−10594. (70) Maseras, F.; Morokuma, K. IMOMM: a new integrated ab initio + molecular mechanics geometry optimization scheme of equilibrium structures and transition states. J. Comput. Chem. 1995, 16, 1170− 1179. (71) Burykin, A.; Warshel, A. What really prevents proton transport through aquaporin? Charge self-energy versus proton wire proposals. Biophys. J. 2003, 85, 3696−3706. (72) Weber, T. A.; Stillinger, F. H. Reactive collisions of hydronium and hydroxide ions studied with the polarization model. J. Phys. Chem. 1982, 86, 1314−1318.

H


Restrained Proton Indicator in Combined Quantum-Mechanics

Recommend Documents