Citation
Tools and methodology for the adaptive multiscale simulation of proton transport in complex environments

Material Information

Title:
Tools and methodology for the adaptive multiscale simulation of proton transport in complex environments
Creator:
Duster, Adam William ( author )
Language:
English
Physical Description:
1 electronic file (104 pages) : ;

Subjects

Subjects / Keywords:
Protons -- Physiological transport ( lcsh )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Review:
This thesis reports three developments in the combined quantum-mechanics/molecular-mechanics (QM/MM) methodology and software. First, an interface between the program QMMM and the MM package NAMD was implemented, substantially increasing the computational efficiency of MM calculations. Second, the permuted adaptive partitioning QM/MM algorithm implemented in the QMMM software package was parallelized using the OpenMP protocol. The algorithm can scale approximately linearly when the number of buffer groups vastly exceeds the number of threads requested. These two improvements will increase the scope and scale of future calculations conducted with QMMM. Finally, the feasibility of using a restrained proton indicator in biased (MD) simulations was investigated. A harmonic restraint potential was imposed on an existing proton indicator, the position of which approximates the location of an excess proton in water. When applied to a carbon nanotube model system, the restrained proton accelerated water recruitment into the nanotube when it was restrained above the nanotube opening, in agreement with previous simulations in the literature, where a different algorithm for proton representation was used.
Thesis:
Thesis (M.S.)--University of Colorado Denver.
Bibliography:
Includes bibliographical references.
System Details:
System requirements: Adobe Reader.
Statement of Responsibility:
by Adam William duster.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
on10003 ( NOTIS )
1000312229 ( OCLC )
on1000312229

Downloads

This item has the following downloads:


Full Text
TOOLS AND METHODOLOGY FOR THE ADAPTIVE MULTISCALE SIMULATION
OF PROTON TRANSPORT IN COMPLEX ENVIRONMENTS
by
ADAM WILLIAM DUSTER B.S., Metropolitan State University of Denver, 2014
A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Chemistry Program
2017


This thesis for the Master of Science degree by Adam William Duster has been approved for the Chemistry Program by
Hai Lin, Chair Haobin Wang Michael Crowley
Date: May 13, 2017


Duster, Adam William (M.S., Chemistry)
Tools and Methodology for the Adaptive Multiscale Simulation of Proton Transport in
Complex Environments
Thesis directed by Associate Professor Hai Lin
ABSTRACT
This thesis reports three developments in the combined quantum-mechanics/molecular-mechanics (QM/MM) methodology and software. First, an interface between the program QMMM and the MM package NAMD was implemented, substantially increasing the computational efficiency of MM calculations. Second, the permuted adaptive partitioning QM/MM algorithm implemented in the QMMM software package was parallelized using the OpenMP protocol. The algorithm can scale approximately linearly under appropriate conditions. These two improvements will increase the scope and scale of future calculations conducted with QMMM. Finally, the feasibility of using a restrained proton indicator in biased (MD) simulations was investigated. A harmonic restraint potential was imposed on an existing proton indicator, the position of which approximates the location of an excess proton in water. When applied to a carbon nanotube model system, the restrained proton accelerated water recruitment into the nanotube when it was restrained above the nanotube opening, in agreement with previous simulations in the literature, where a different algorithm for proton representation was used.
The form and content of this abstract are approved. I recommend its publication
m
Approved: Hai Lin


To my parents and sister for their unconditional love and support. To Nduka Odizor, for his mentorship in forming my values, character, and morals. To David and Lauren Suarez, Mary Rose Siason, Joshua Skotak, and Zachary Lowery for your friendship and support over many years. To Michael Jacobs and Andrew Bonham for making my Masters studies a reality. To Hai Lin and all members of the Lin Lab, for providing the opportunity to perform research and for making the past 3 years the happiest and most fulfilling of my life.
IV


ACKNOWLEDGEMENTS
I would like to thank my advisor, Dr. Hai Lin, for his guidance and help in preparing this thesis. I would like to thank Christina Garza for her help in extensively testing QMMM and Nara Chon for her helpful discussions. I would also like to thank Professor Walter Thiel, Dr. Xin Wu, and Dr. Pavlo Dral for their work on the OMx family of semiempirical methods, and their willingness to share their work and teach me about them.
This work is supported by National Science Foundation (CHE-1564349), Camille & Henry Dreyfus Foundation (TH-14-028), and NVIDIA Corporation. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) under grant CHE-140070, which is supported by National Science Foundation grant number ACI-1053575, and the National Energy Research Scientific Computing Centre (NERSC) under grant m2495, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This research is also supported by the University of Colorado Denver Chemistry Departmental Student Research Fellowship.


TABLE OF CONTENTS
CHAPTER
I. INTRODUCTION.............................................................1
1.1 QM/MM Multiscale Modeling...........................................1
1.1.1 What is Multiscale Molecular Modelling?........................1
1.1.2 What is QM/MM?..................................................1
1.2 Conventional QM/MM Methodology.......................................3
1.2.1 The MM Energy...................................................3
1.2.2 The QM Energy...................................................5
1.2.3 The OM3n Method.................................................9
1.2.4 Partitioning a QM/MM System....................................10
1.2.5 The QM/MM Hamiltonian..........................................11
1.2.4 QM/MM Partitions through Covalent Bonds........................12
1.3 Adaptive Partitioning QM/MM Methodology.............................13
1.3.1 Partitioning of an AP QM/MM System.............................14
1.3.2 The PAP Potential..............................................15
1.3.3 The PAP Gradients and the mPAP Scheme..........................17
1.4 Proton Transfer in Complex Environments.............................18
1.4.1 Grotthuss Shuttle Mechanism of Proton Transfer.................18
1.4.3 The Proton Indicator...........................................21
II. IMPLEMENTATION OF NAMD-QMMM INTERFACE...................................23
2.1 QMMM An Interfacing Program.......................................23
2.2 Motivation for Adding the NAMD Interface............................23
ii


2.3 Structure of the NAMD Interface
24
2.4 Topology Changes and the PSF....................................26
2.4.1 The Bond List..............................................27
2.4.2 The Angle List.............................................27
2.4.3 The Dihedrals List.........................................28
2.4.4 The Impropers and Cross Terms..............................28
2.5 Benchmarking....................................................29
III. IMPLEMENTATION OF PARALLEL ADAPTIVE PARTITIONING......................31
3.1 Theoretical Background..........................................31
3.2 Implementation of Parallel Permuted Adaptive Partitioning................31
3.3 Scaling and Benchmarks..........................................33
IV. IMPLEMENTATION OF INDICATOR RESTRAINT.................................38
4.1 Restraint Potential and Decomposed Gradients....................38
4.2 Behavior of the Indicator for the Zundel Cation.................40
4.3 Behavior of the Indicator for Collision Simulations.............43
4.4 MD Simulations of a Water Sphere................................45
V. APPLICATION OF INDICATOR RESTRAINT TO CARBON NANOTUBE..........................49
5.1 Introduction....................................................49
5.2 Simulations of Indicator Restrained at Nanotube Opening.........51
5.3 Simulations of Indicator Restrained at Different Distances from Nanotube.55
5.4 QM/MM Umbrella Sampling of PMF for Proton Transfer in Nanotube...........57
VI. CONCLUSION............................................................61
iii
REFERENCES
63


APPENDIX A: BENCHMARKING INFORMATION..............................................78
A.l Test System Specifications........................................78
A.2 Tinker/NAMD Benchmark Calculation.................................79
A. 3 PPAP Benchmarking of a Water Box for Various Systems.............79
Desktop Workstation 8..............................................79
XSEDE Stampede.....................................................80
Desktop Workstation 5..............................................80
APPENDEX B: DECOMPOSED GRADIENT INCLUDING HYDROGEN................................82
B. 1: Variables and Functions used in the Gradient Decomposition Equations.82
B.2: Decomposition of the Potential onto XAj...............................82
B.3 Decomposition of the Potential onto XD.................................83
B.4: Decomposition of the Potential onto XHm:.........................84
B. 5: Gradient Graphs.................................................85
APPENDIX C: Additional Details for MD Simulations.................................87
C. l NVE Collision Simulations........................................87
C.2 MD simulations of Water Sphere....................................89
C.3 MM System used for Nanotube Simulations...........................89
C.4 MD Simulations of Restrained Indicator at Pore Opening.................90
C.5 MD Simulations of Restraints at Different Locations....................92
C.6 Umbrella Sampling of Proton Transfer in CNT............................93
IV


LIST OF ABBREVIATIONS
AP Adaptive Partitioning
CPS Capped Primary Subsystem
CSH C-Shell
DBA Density Based Adaptive Partitioning QM/MM
DFT Density Functional Theory
ES Entire System
EVB Empirical Valence Bond (Theory)
HDD Hard Disk Drive
HF Hartree-Fock
HPC High Performance Comput(-er/-ing)
MD Molecular Dynamics
MM Molecular Mechanic(-s/-al)
MNDO Modified Neglect of Differential Overlap
MP2 2nd Order Moller-Plesset Perturbation Theory
mPAP Modified Permuted Adaptive Partitioning
MS-EVB Multistate Empirical Valence Bond (Theory)
NA N-adaptive
NAMD NAMD software package
NDDO Neglect of Differential Diatomic Overlap
OMx Orthogonalizati on-corrected Methods
PAP Permuted Adaptive Partitioning
PDB Protein Data Bank File
PMF Potential of Mean Force
post-HF Post-Hartree-Fock
PPAP Parallel Permuted Adaptive Partitioning
PS Primary Subsystem
PSF Protein Structure File
QM Quantum Mechanic(-s/al)
QMMM QMMM Software Package
QM/MM Combined Quantum-Mechanics/Molecular-Mechanics
RDF Radial Distribution Function
SCF Self-consistent-field
ss Secondary Subsystem
SSD Solid State Drive
SSE Strong Scaling Efficiency
vdW van der Waals
WHAM Weighted Histogram Analysis Method
v


CHAPTER I
INTRODUCTION
1.1 QM/MM Multiscale Modeling
1.1.1 What is Multiscale Molecular Modelling?
Multiscale models include different resolutions of length or time scales within the same model to ensure a model can represent a system of interest with sufficient accuracy while remaining computationally feasible.1 This has become an indispensable addition to the chemists toolbox for describing and studying chemical phenomena in complex environments, as evidenced by the proliferation of use of multiscale modelling methodology and by the awarding of the 2013 Nobel Prize in Chemistry to Martin Karplus, Michael Levitt, and Arieh Warshel for the development of multiscale models for complex chemical systems.2 Multiscale modelling has been applied to many chemical systems, including biological,3 materials science,4 and solid-state systems.5
Multiscale molecular modelling typically combines one or more of the following methods: quantum mechanics (QM), classical molecular mechanics (MM), coarse graining (CG) models, and finite element methods. These methods cover a wide range of scales from the sub-nanometer to beyond millimeters, and are employed in various types of calculations including single-point calculations, geometry optimizations, molecular dynamics (MD), and Monte Carlo (MC) simulations. This work is devoted to one specific multiscale modeling technique, the combined quantum-mechanics/molecular-mechanics (QM/MM) methodology.
1.1.2 What is QM/MM?
Quantum mechanics is essential to accurately describe many chemical phenomena such as electronic excitation, the making and breaking of bonds, extensive polarization, and


charge transfer. While semiempirical QM methods have been used for systems with >1000 heavy atoms,6 Density Functional Theory (DFT) and post-Hartree-Fock (post-HF) methods are typically computationally feasible for up to hundreds of atoms7 with dynamics simulation times in the hundreds of picoseconds.8 The computational expense of QM can be especially problematic if many atoms surrounding the system of interest must be considered for an accurate study.
Fortunately for many chemical processes, changes of interest in the electronic-structure are limited to a small number of atoms. Thus, sometimes it is not necessary to describe the whole system with QM to obtain an acceptable level of accuracy.1 The environment around the system of interest can be described by computationally efficient MM empirical force fields. MM system sizes can range to tens of millions of atoms9 or a millisecond timescale.10 Additionally, the total simulation time for MM studies with hundreds of thousands of atoms now routinely reaches one microsecond. By embedding a QM subsystem within an MM subsystem, the accuracy of QM can be combined with the computational efficiency of MM, allowing for simulations with hundreds of QM atoms surrounded by hundreds of thousands of MM atoms at the 100 ps timescale.11
Warshel and Levitt published the first paper on QM/MM in 197412, which described new methodology to explore the electrostatic stabilization of a carbonium ion intermediate in the cleavage of hexasaccharides in the enzyme lysozyme. Field et al. extended QM/MM past single-point calculations to include MD, Monte Carlo, and geometry optimization computations in 1990.13 The increase in computing power combined with this newly available methodology marked an explosion in QM/MM method development throughout the 90s.14-21
2


1.2 Conventional QM/MM Methodology
1.2.1 The MM Energy
The MM energy of a prototypical 1st generation force field can be calculated by an analytical function with a corresponding set of parameters called a force field. The energy can be expressed as the combination of the bonded and non-bonded interactions:
^MM ^bonded ~b f^non-bonded (1-1)
The bonded interactions can be expressed as the combination of the bond, angle, dihedral, improper, and cross term interactions:18
^bonded ^bond ~b Vang\e + V"cji]lecjrai H bVjmpr0per + V^r0ss (1.2)
^bond= Y, K^b~b0)2 bonds (1.3)
J'angle = ^ ~ 9 o) angles (1.4)
^dihedral ^ ^improper ^ ^>(0 0o) impropers (1.6)
Across / "b COs(ri(f) + TWp 5nm)) (1.7)
crosses
l^bond is the interaction energy between all pairs of bonded atoms (e.g. between the atoms labeled 1 and 2 in Fig 1.2, called 1-2), Fangie is the interaction energy between all connected triplets of atoms (e.g. atoms 1,2, and 3, called 1-3), and Vdihedrai is the interaction
3


energy from all sets of 4 atoms connected in a line (e.g. atoms 1,2,3 and 4, called 1-4). The symbols b, 6, and cp represent the bond length, angle, and dihedral respectively for a given geometry. Kb, Kg, and K(p are force constant parameters. The equilibrium bond distance and angle parameters are represented by b0 and 90, respectively. The dihedral expression has two additional terms, n and S, which represent the torsional multiplicity and phase, respectively, ^improper is the energy for out of plane bending, used to enforce planarity of atoms (for example, to keep atoms 6, 7, 8, and the hydrogen connected to 7 on the same plane). is the force constant parameter, 0O is the equilibrium out-of-plane angle, and

Fig 1.1: Atom indices for 1 -(1 //-imidazol-4-vl)propanc:
Non-bonded interactions are typically not computed between atoms included in 1-2 and 1-3 interactions, and may be excluded or scaled down for included in 1-4 interactions. A
4


minimal set of non-bonded interactions includes columbic interactions via point charges and van der Waals (vdW) interactions which implicitly model exchange-dispersion interactions with the empirical Lennard-Jones potential:
Vn
non-bonded ^coulomb rvan der Waals
+ Vv.
V,
coulomb
= x
non-bonded
1 RiRj 4ne0 rtj
(1.8)
(1.9)
V,
- I
van der Waals / '-ij
non-bonded
Krij,
12

l]
'll.
(1.10)
Here, qt is a parameter for the charge assigned to atom t, and e0 the vacuum permittivity of free space. In both equations, rtj is the distance between non-bonded atoms i and j, For the van der Waals equation (1. 10), atj, and e,y are parameters for the collision diameter between non-bonded atoms i and j, and the well depth of the interaction curve, respectively. These prototypical functions are similar to those from commonly used force fields,25 such as the CHARMM26"27, OPLS-AA,28 Amber,29 and GROMOS30 force fields.
1.2.2 The QM Energy
The QM energy for potentials in this work can be found by obtaining approximate solutions for the time-independent, non-relativistic Schrodinger equation, which can be expressed as:31"32
H|Y) = £|Y) (1.11)
where H is the Hamiltonian operator, V is the wave function, and E is the energy of the system. If the wave function is known, the energy may be found through the following expression:
5


(1.12)
The Hamiltonian operator consists of the nuclear and electronic kinetic and potential energy operators:
H Tn + Te + Vne + Vee + Vnn
(1.13)
where Tn and Te are the nuclear and electronic kinetic energy operators, respectively, and Vne, Vee, ar|d Vnn are the potential energy operators for the interaction energy between the nuclei and electrons, between the electrons, and between the nuclei, respectively.
To reduce the dimensionality of the problem, the Born-Oppenheimer approximation is invoked, which neglects the coupling between the nuclear and electronic motion by considering the electrons to be moving in a field of fixed nuclei. Therefore, the wave function only explicitly depends on the electronic coordinates, while parametrically depends on the nuclear coordinates. Thus, Tncan be set to 0, and Vnn is a constant. The total energy (in atomic units) becomes:
M M
£.. = Eelec + X X |R A_Br !
A=1B>A
Helec l^elec) = Ee\ec l^elec)
Helec Te + Vne + Vee
(1.14)
(1.15)
(1.16)
where M is the number of nuclei, ZA is the charge of nucleus A, R4 is the coordinate vector for nucleus A, V^gc is the electronic wave function which is explicitly dependent upon only
6


the electronic coordinates, and Heiec is the Hamiltonian operator acting on the electronic
wave function. Equation (1.1.5) is called the electronic Schrodinger equation.
To satisfy the antisymmetric requirement for the electronic wave function, the electronic wave function can be approximated by the Slater determinant of the spin orbitals of the individual electrons, which are assumed throughout this work to be orthonormal:
where Xk(n) is the spin orbital Xk occupied by electron n, N is the total number of electrons, and Sij is the Dirac delta function, which delineates the orthonormality between spin orbitals
approximated by one Slater determinant.
The electronic Hamiltonian can be decomposed into the core Hamiltonian, which depends only upon the coordinates of one electron, and the two-electron operator. The core Hamiltonian describes the kinetic energy of an electron and the potential energy of an electron in a field of nuclei. The two-electron operator acting on the electronic wave function describes the coulombic and exchange interactions between electrons. These operators may be related to the electron indices i and j through the following equations:
r/i(l) *2(1) " Xn(. 1)'
1 Zi(2) *2(2) Xn(. 2)
; (Xi\Xj) = 5ij (1.17)
-Xi (n) Xz(n) XnW-
Xt and Xj For the model systems studied within this work, the wave function is adequately
N
N N
(1.18)
i=l i=1 j>i
(1.19)
7


1
(1.20)
where h; is the core Hamiltonian for electron t, gtj is the two-electron operator for electrons i and j, and r, is the coordinate vector for electron i.
The electronic Schrodinger equation can be solved by evaluating the following equation:
Ee lec = <%lec|Helecl%lec) (1-21)
After many intermediate steps, an expression for the electronic energy can be obtained:
N N N
EeIec = 'YJ(Xi\hi\Xi) +^YJYj(.{Xj\]i\Xj) {Xj\Ki\Xj)) (1-22)
i=i i=ij=i
h\X)(2)) = IXj (2) > (1.23)
KiU,(2)) = Ui (1) | giz U/ (1)) Ui (2) > (1.24)
where Jj and Kj are the coulomb and exchange operators for electron i, respectively, which represent the average local potential at Xi from interactions with other electrons (note that Kj is a non-local operator). This treatment is known as the mean field approximation, and it neglects coupling between the motions of individual electrons of opposite spins.
As the spin orbitals of other electrons are required for Jj and Kj, the expression for Eeiec is non-linear. Approximate solutions can be obtained by iteratively improving an initial trail wave function until the energy change between two iterations falls below a specified threshold, a process called the self-consistent-field (SCF) method. According to the variational principle, the lowest energy found by varying the spin orbitals is an upper bound
8


of the exact ground state energy. The process of finding the QM energy using these equations is called the Hartree-Fock (HF) method.33"34 Post-Hartree-Fock (post-HF) methods like 2nd Order Moller-Plesset Perturbation Theory (MP2),35 are based upon the HF method and improve upon the HF energy at higher computational costs.
1.2.3 The OM3n Method
The OM3n method36 used extensively throughout this work, is a semiempirical QM method. In semiempirical methods, many of the integrals are approximated through parameterized expressions to reduce computational expense. Parameters are adjusted to reproduce higher-level QM calculations and experimental data. Thus, the accuracy of a given semiempirical method depends on the parameters used and the approximations made in the method.
The OM3n model is a reparametrized version of the OM337 method that accurately describes proton transfer in water. This method accurately reproduces the MP2 energy surface for proton transfer in the Zundel ion, yielding a barrier for proton transfer in agreement with path integral calculations that include nuclear quantum effects of the proton.38'39 The OM3n method is a member of the family of orthogonalization-corrected methods (OMx, x = 1,2, 3),40"41 which are based on the modified neglect of differential overlap (MNDO) model.42"43
In the OMx methods, core electrons are treated implicitly and valence electrons are treated explicitly. A minimum, orthogonalized basis set of Slater type atomic orbitals is used, called a Lowdin basis.44"45 Because the use of an orthogonalized basis can introduce significant errors in the one-electron integrals, additional empirical functions called orthogonalization corrections are added to the one-electron integrals in the OMx methods. As
9


the two-electron integrals are similar when evaluated in an orthogonal basis versus a nonorthogonal basis, no corrections are added to these terms.
Similar to other MNDO-type semiempirical models, many integrals are approximated in the OMx methods for computational efficiency.46 One-center electron integrals are calculated as the sum of a parameter which represents the attraction of the electron to the nucleus it is based on (taken from spectroscopic data) and the potential due to the other nuclei in the system. Two-center one-electron integrals are approximated based on an averaged pair of resonance parameters weighted by the corresponding overlap integral. The values of two-center two-electron integrals are based on semiempirical multipole-multipole interactions. The OMx methods invoke a modified version of the neglect of diatomic differential overlap (NDDO) approximation.44"45 Thus, one-electron integrals with three centers, and all three- and four-center two-electron integrals are neglected. Core-core electronic attractions and repulsions are evaluated through the corresponding integral, and summed with a parameterized function. A thorough description of the OMx formalism has been published by Dral et al.41
1.2.4 Partitioning a QM/MM System
In all QM/MM calculations, the system must be divided into a QM and MM partition (Fig 1.1). In conventional QM/MM methods, this partition is fixed and may not change throughout a simulation. The QM region, also called the primary subsystem (PS), is embedded in the region described with MM, the secondary subsystem (SS). The entire system (ES) consists of all atoms in the simulation. The choice of partitioning is one of the most important decisions in construction a QM/MM model, and great care must be taken to
10


ensure that all important atoms are designated as QM. The nature of the interaction between the PS and SS is one of the defining characteristics of a QM/MM method.
Fig 1.2: In a QM/MM simulation the system is partitioned into two interacting parts, the primary subsystem which is calculated using QM and the secondary subsystem which is calculated using MM.
1.2.5 The QM/MM Hamiltonian
The general form of the QM/MM potential comprises the QM potential of the PS, the MM potential of the SS, and the interaction energy between the PS and SS.12-16
^QM/MM;ES = ^QM;PS + ^MM;SS + ^QM/MM;PS|SS (1.25)
In this work, for simplicity, the mechanical embedding scheme is used, where the QM calculation for the PS is carried out in the gas phase. The subtractive definition of the QM/MM energy is employed, which has the following form:
^QM/MM;ES = ^MM;ES ^MM;PS + ^QM;PS (1.26)
In our calculations, the ES energy is computed at the MM level. Then, the MM energy of the PS is subtracted from the ES. Finally, a QM calculation of the PS is carried out in vacuum and added to obtain the final QM/MM energy. Electrostatic interactions between the PS and
11


SS are described at the MM level with the vdW and Coulombic nonbonded expressions.
More complicated formulations exist which describe mutual polarization18,47 and partial charge transfer48"50 between the PS and SS. However, they were not employed in this work due to their associated computational expenses.
1.2.4 QM/MM Partitions through Covalent Bonds
It is sometimes necessary to place the PS/SS boundary through a covalent bond in order incorporate parts of larger molecules, such as the sidechains of proteins, in the QM region. Many different methods for treating this situation can be found in the literature.18,51 However, only the link atom method with the scaled-bond-distance treatment52'53 is discussed, as it is the method implemented in our code, QMMM.
In this method, an SS boundary atom (designated M), the atom covalently bonded to a PS atom (designated Q), is replaced by a cap atom (designated L) in the PS. The cap atom is usually Hydrogen. The location of the link atom is placed on the line between atoms M and Q. The Q L bond length is then scaled by an adjustable parameter, CH. One recommendation for CH is to use the ratio of two equilibrium bond lengths:52
Ch = ^ (1.27)
kq-m
where Rq_h and Rq_m are the equilibrium bond distances taken from the MM parameters for the Q-H and Q-M bonds, respectively.
When PS is capped with a link atom, it is called the capped primary subsystem (CPS). However, when the term PS is used in this work, it encompasses primary subsystems with and without cap atoms. The addition of the link atom has serious implications for the
12


QM/MM system, as the topology of the CPS is different from that of the ES. Therefore, the topology of the CPS must be regenerated and the appropriate parameters must be available.
1.3 Adaptive Partitioning QM/MM Methodology
Adaptive partitioning (AP) methodology was developed to address two of the problems that arise when using QM/MM on diffusive systems: (1) molecules of interest may travel to the QM/MM boundary where their environment is not properly described with QM, and (2) frequently a large PS must be used to alleviate boundary effects, which increases the computational expense. AP methods address these problems by featuring the on-the-fly reclassifications of atoms as QM or MM so that molecules of interest are always in the QM subsystem. As such, smaller QM subsystems can be used with AP methods.
There are many different types of AP methodology in the literature, which include the hot-spot scheme,54 (HS) scheme, the our own n-layered integrated molecular orbital and molecular mechanics-exchange of solvent (ONIOM-XS) scheme,55"56 the permuted adaptive-partitioning (PAP) scheme57 and modified PAP (mPAP) scheme,58'60 the sorted-AP scheme (SAP),57'60 the difference-based adaptive solvation (DAS) scheme,61"63 the buffered-force (BF) scheme,64'67 the number-adaptive (NA) scheme,68"69 the size-consistent multipartitioning (SCMP) scheme,70 the density-based adaptive (DBA) scheme,71 and the time-adaptive switching interaction potential scheme.72 This work concentrates on the PAP and mPAP schemes, due to their demonstrated numerical stability.57
13


1.3.1 Partitioning of an AP QM/MM System
In many AP schemes, the system is divided into three parts (Fig 1.2): the active zone (QM part), environmental zone (MM part), and buffer zone (dual QM/MM characteristics). In this work, the partitions are distance based.
Fig 1.2 : An example of partitioning the system into three zones for adaptive partitioning. The zones are defined as a set of concentric spheres, defined by the radii Rmin and Rmax. The blue zone is the enviromnental zone and is treated with MM. The yellow zone is the buffer zone and lias dual QM/MM characteristics. The white zone in the middle is the active zone and is treated with QM only
First, a center for the active zone is chosen. This can be a stationary point, such as a dummy atom, or a moving point, such as a specific ion or binding site center (defined by the center of mass of several residues). Next, atoms are divided into groups which can be whole molecules or molecular fragments73 Groups are added to a zone based on the distance, ru from the active zone center to either the center of mass of the group or a designated reference point (such as the O atom in the H2O molecule) according to equation (1.3.1). All atoms within a sphere of radius /?min can be assigned to the active zone. All atoms within a spherical shell between /?min and Rmax are designated as buffer groups. This shell typically has a width of 1-2 A. All atoms greater than Rmax from the active zone center are assigned to
14


the environmental zone. If needed, a group can be designated as always QM or MM regardless of its location. The list of atoms belonging to each zone is updated on-the-fly throughout the course of the simulation, and atoms may change groups if bond-breaking and bond-forming reactions occur.
(Active Zone if rt < /?min
Buffer if^min < n < Rmax (1.28)
Environmental Zone if rt > Rmax
1.3.2 The PAP Potential
The PAP potential has the form of a many-body expansion and is expressed with the following equation:57
N-l N
V=VA+ £ P|(ViA-VA) + 2 Z ^ VM'
i=l i=l j=i+l \
N \ N-2 N-l N f
< > i >1 Vzz z w / ( -1 ( i 1 *1 1 r i 1 *1
1 -M J / i-l j-i+1 k-j+1 y
v-V
vi,j,k
vA+
N
£ (vA-vA)
r=i,j,k
N-1,N
- z
(p,q)=(i,j),(i,k),(j,k)
VM
A
/
(1.29)
Equation 1.3.2 can be formatted more concisely as:73
N N N N-l N N
v=va pt)+^ ]^[(i Pj)+z z p,piv n c1 ^ + (i3)
i=l i=l j*i i=l j=i+i kji
In these equations, VA represents the QM/MM energy with the active zone treated at the QM level, V[A represents the QM/MM energy with the tth buffer group and the active zone treated
15


at the QM level, V-j represents the QM/MM energy with the ith and jth buffer groups and the active zone treated at the QM level. This equation is expanded to the limit of V£2,...,n where all JV buffer groups and the active zone are treated at the QM level. Buffer groups are classified as QM in some partitioning configurations and MM in others which yield their dual QM and MM character. However, the active zone is always treated by QM, and the environmental zone is always treated by MM.
Pj is the smoothing function for the ith buffer group, which varies smoothly and monotonicaly in value between 0 and 1 depending on the position of the buffer group. As a group enters the buffer zone at the limit of the environmental zone, the smoothing function is equal to 0. As a group enters the buffer zone at the limit of the active zone, the smoothing function is equal to 1. In our implementation of PAP, a fifth-order spline is used for the smoothing function:
Pfaf) = -6a? + 15a? 10a? + 1 (1.31)
where at is a metric for the position of the buffer group calculated by:
at = Jl for Pmin < n < Pmax (1.32)
''max ''min
In total, 2N configurations must be evaluated to obtain the full PAP potential. However, this is rarely necessary as high-order many-body contributions are usually insignificant. Moreover, the product of the smoothing functions tends towards 0, further decreasing the contribution of higher order terms. Truncation of the potential at to the 3rd order has been shown to capture over 95% of the intermolecular interaction energy with other potentials, so the expansion can typically be truncated to increase computational
16


efficiency.74"76 Because many commonly used QM methods formally scale as 0(N4), where JV is the number of basis functions, it can be faster to compute multiple smaller QM calculations instead of one large one. In these cases, the use of AP methodology will lead to an increase in computational efficiency.
1.3.3 The PAP Gradients and the mPAP Scheme
The PAP scheme conserves both the total energy and momentum,57 and may be used in dynamics simulations under the microcanonical (NVE) and canonical ensembles (NVT). Unfortunately, forces due to the derivative of the smoothing function can be large enough to cause simulation artifacts such as erroneous solvation structures 73,77 For example, the full PAP potential for a system with one group in the buffer can be expressed as:
V = VA + PjVf P1VA (1.33)
The forces can then be found by computing the gradient of equation (1.33) with respect to the coordinates r:
dVA
dr
dVA
dr
Pi +
VA)
(1.34)
The following terms from equation (1.3.7) are deemed the physical forces of the system and called the adaptive forces:
dVA dVf dVA dr dr dr 1
(1.35)
The following term from equation (1.3.7) represents the unphysical forces of the system58,78 and is called the transition, extra, or drift force:
17


(1.36)
VA)
These forces arise due to the mismatch of potentials between the QM and MM systems, and cease to exist when the QM and MM potentials are equal. Currently, it is at least practically impossible to match QM and MM potentials with high levels of dimensionality, prohibiting potential matching for most systems. Therefore, a different approach is necessary to remove the transition forces. In the modified PAP (mPAP) scheme, the forces are simply deleted, which is equivalent to adding external forces that exactly cancel out the transition forces. In this treatment, the QM and MM systems are considered to be open systems in dynamic equilibrium with each other, while the entire system is a closed system. Because extra forces are added to the system, the entire system is now a non-Hamiltonian system. Therefore, the mPAP method cannot be used in microcanonical simulations, but can be coupled to thermostats and used under the canonical ensemble. The mPAP method has been shown to satisfactorily reproduce solvation structures observed in conventional QM/MM simulations with QM zones so large that boundary effects on the solvation become insignificant.58"59
1.4 Proton Transfer in Complex Environments
1.4.1 Grotthuss Shuttle Mechanism of Proton Transfer
Proton transfer is a ubiquitous process in many biological and chemical systems, but specific details of proton transfer, such as the stability of intermediary structures in the mechanism, are still a subject of debate.79'80 Proton diffusion in water is approximately seven times faster than for a Na+ cation, indicating that H3CC does not diffuse through solution as a discreet molecule. Mechanisms for proton diffusion are different from that of other
18


ordinary cations, such as Na+ or K+. and involve the hydrogen bonding81 network for water. In these mechanisms, the proton can be considered a charge defect which is shuttled through changes in the covalent and hydrogen bonding topology of water.
The most commonly accepted mechanism is the Grotthuss shuttling mechanism, named after von Grotthuss who published a mechanism in 1806 for the decomposition of water by galvanic electricity. Although significantly different from the modern mechanism, this early mechanism is the first published description detailing protons hopping from O to O atom through a water wire.82"83 Modem mechanisms describe proton transfer through water as a dynamic reorganization of the hydrogen-bonding environment of water surrounding the excess proton. The proton is typically found in two states, the Zundel84 (H5O2 ) and Eigen85 (H9O4 ) cations (Fig 1.3), which interconvert during transfer events.86 The most commonly accepted transferring mechanism in bulk water can be characterized as Eigen-Zundel-Eigen, where the proton spends most of the time in the Eigen resting state, and transfers through a Zundel transition state to another Eigen state.86'88 It has been suggested that these jumps sometimes occur in a quasiconcerted fashion over two to three waters, with two jumps more likely to occur than three or one.89 The transfer event is coupled to the dynamics of the water clusters surrounding the proton,90 and the reorganization of both the first and second
Eigen
Grotthuss Mechanism
Zundel
V j
p *, 1 ^%
^ xMt c *#- SJ
Fig 1.3: Eigen and Zundel structures of the solvated proton and the Grothuss Mechanism for proton transfer. The hydrogen bonding network of the water changes to shuttle an excess charge defect through solution.
19


solvation shells into specific arrangements,38,87 called presolvation,91"93 is crucial. Other mechanisms can be characterized as Zundel-Zundel, where molecules spend most of their time in Zundel states and transfer directly from Zundel to Zundel state.80,94 This mechanism has been observed for proton transfer in confined spaces, such as carbon nanotubes.95
The free energy barrier of proton transfer was calculated to be approximately 0.56 kcal mol"1 in Car-Parrinello MD simulations.38,96 Path integral molecular dynamics, which include quantum fluctuations of the nuclei, predict the barrier to be 0.15 kcal mol'1.38'39 These barriers are close to or less than kBT at 300 K ( approximately 0.6 kcal mol'1). Therefore, the hydrogen bonds between the donor and acceptor O atoms, and the transferring H at the Zundel transition state are classified as a low-barrier hydrogen bond,97"98 where nuclear tunneling effects are insignificant. This provides justification for approximating the nuclei in water as a classical particle in some dynamics simulations. However, when lower probability events, such as water autoprotolysis, are important to the chemical model, it may be necessary to incorporate nuclear quantum effects.99
Throughout the years, proton transfer has been described with theories and models at various levels. Full-quantum molecular dynamics are the most computationally expensive methods, but have provided the energetic data cited above and many mechanistic insights.38, 88,98-ioo Empiricai valence bond (EVB) theory by Warshel,101 and its extension, the multistate empirical valence bond (MS-EVB) theory have also been applied towards studying proton transfer in water and other environments.36,94,96,102-114 In addition, many reactive force-field methods115"124 and semi empirical QM methods,125"133 including the OM3n model134 used extensively in this work, have been developed.
20


1.4.3 The Proton Indicator
A major challenge in proton modelling is that the identity of the transferring proton is not known ahead of time and changes constantly throughout an MD simulation. This difficulty is compounded in biased simulations, where it can be unclear which atoms to apply restraints to. To overcome this problem, the proton indicator used in this work, first suggested by Hofer et al.,135 was developed by Wu et al.134 to track the approximate location of the excess proton over time. It was then extended by Pezeshki and Lin59 to update the active zone center centered at the indicator in adaptive QM/MM simulations. The indicator is used in this work as a collective variable in proton transfer reactions.
The indicator was constructed to fulfill certain criteria.6,59 The position of the indicator should be superimposed on the location of the shared proton in the limit of the Zundel cation, and on the donor O atom in the limit of the Eigen cation. The location of the indicator is a function of the coordinates of the donor oxygen D, all hydrogen atoms Hm bonded to D, and all possible acceptors Ay within a predefined enlisting radius rLISX from D. The most likely acceptor is denoted as A. When the proton is transferred, D and A switch.
The location of the indicator is determined by a linear combination of the coordinates of D and all possible acceptors Ay:
Here, gx is a normalization factor, g(pmj) is the weight function associated with the ratio of the projected donor-acceptor vector, pmy, and XD and XA are the Cartesian coordinates of the donor oxygen and the /-th acceptor oxygen respectively. M is the total
1
(1.37)
21


number of hydrogen atoms bonded to D, and / is the total number of acceptors within a radius of size rLISX. The ratio pmj is a metric for how close a given atom Hm is to donor D versus an acceptor Ay and is calculated according to the following equation:
Pmj ~
ID Hm ID Ay
(1.38)
where rDHm = XHmXD and rDA = XA XD. The linear combination of coordinates is
normalized by gh which is calculated according to equation (1.39) using the weight function described in equation (1.40):
J M
= l + g(Pmj) (1-39)
j=l m= 1
tO if 1 < x
g(x) = | 6x5 + 15x4 10x3 + 1 if 0 < x < 1 (.1 ifx < 0
The weight function g(pmj) depends on the reduced variable x:
x = x(pmy) = 1
Pm j Pmj 0-5 Pmj
(1.40)
(1.41)
where p?riy = (^dh/^list) XqH is a parameter set to a distance slightly larger than the equilibrium OH bond distance in a hydronium ion (1.0 A), and rLISX is the threshold distance for enlisting possible acceptors.
22


CHAPTER II
IMPLEMENTATION OF NAMD-QMMM INTERFACE
2.1 QMMM An Interfacing Program
QMMM136 is an interfacing program for conducting QM/MM single-point calculations, geometry minimizations, and molecular dynamics simulations. As a driver, QMMM calls an MM package (NAMD137 or Tinker138) and/or a QM package (Gaussian 03139 or 09,140 GAMESS,141 ORCA,142 or MNDO143) to calculate the energies, gradients, and/or hessians. QMMM parses the output from the selected MM or QM packages, synthesizes energies and energy derivatives, and propagates the trajectory (Fig 2.1). The source code files and associated scripts within the QMMM package now contain >70,000 lines, and are written in FORTRAN 95, FORTRAN 2003, C-shell (CSH), Perl, and Python. The program also requires a compiler that supports OpenMP version 3.0.
2.2 Motivation for Adding the NAMD Interface
Before the work described in this thesis was completed, Tinker was the only MM package which interfaced with QMMM. While Tinker is a flexible, open source, and free
QM Interface MM Interface
Figure 2.1: Structure of QMMM code. QMMM calls a QM package from the programs on the left and an MM package from the programs on the right. The interface to NAMD was added in this work.
23


MM package, the serial and parallel implementations of the program do not rival the speed of other MM packages. In QM/MM calculations for a very large ES and small PS, when employing a semiempirical QM potential and using Tinker, the MM calculation can be comparable to or slower than the QM calculation.
NAMD was specifically designed to run large MM systems quickly in high performance computing (HPC) environments. At the time of this publication, NAMD is available in the software stack for NERSC and XSEDE HPC resources. Moreover, it is freely available and open source. Tools for creating NAMD input files are also integrated within the widely used visualization software, VMD.144 VMD automates many of the steps for the creation of input files necessary for running a calculation within QMMM and cuts down on system-building time. Finally, the TCL scripting language can be used to extend the functionality of NAMD. Therefore, it is possible to add functionality (such as for specialized MM calculations needed for more advanced embedding schemes) without modifying the source code of NAMD.
2.3 Structure of the NAMD Interface
To use NAMD as the MM package for MM single-point gradient calculaitons, the user first creates a protein data bank (PDB) file of the entire system (ES) using the worldwide PDB format, version 3.30.145 Because QMMM uses the element symbol from the column described in the standards for assigning the element to an atom for QM calculations, a helper script, f ind_elsym. py, included within the QMMM distribution, assists the user with acquiescing to the PDB file standards by finding the elemental symbol for each atom listed in the PDB and appending it to the proper column. The user then needs to create an X-PLOR formatted CHARMM PSF file, which contains the topological information and atom types
24


for the system. Finally, the user supplies the force field to use for the MM calculation (Fig 2.2).
QMMM prepares NAMD calculations by writing 4 files: the coordinates binary file of the system, the PSF text file for the system, the input text file for NAMD, and the PDB text file for the system if the system is the primary system (PS). The ES PDB file is used to supply the ES coordinates to QMMM at the first step, and to fulfill the requirements that a PDB exists for an ES calculation in NAMD. Therefore, the ES PDB file is never rewritten through the course of a QM/MM calculation. However, the PS PDB, PS PSF, coordinate binary file, and the NAMD input file are rewritten at every step.
After the files are written, QMMM calls a CSH shuttle script, namdshuttle, to invoke the NAMD executable. The shuttle script is configurable by the user, and can
Figure 2.2: Schematic showing files used for communication in the QMMM interface with NAMD. Blue represents the files that QMMM creates, red represents the files that NAMD creates, and green represents the files that the user must supply to ran a QM/MM calculation.
25


accommodate the use of different environments (e.g. using ibrun to call the MPI/Charm++ parallelized version of NAMD). The user is also able to choose the number of MPI/OpenMP processes that NAMD will use a feature that is important for using parallel permuted adaptive partitioning and discussed in Chapter III. If the NAMD calculation is successful, NAMD will create several binary files and a text file with the output. QMMM only reads the binary force file for the gradient, and the text output file for the energy.
2.4 Topology Changes and the PSF
One of the primary challenges with all QM/MM methods is dealing with the topology of the MM system. If the QM partition cuts through covalent bonds, there will be new bonds, angles, and dihedrals that were not in the ES PSF. Parameters must be supplied for these new interactions if they are not already described in the forcefield.
Additionally, if bonds are made or broken within the QM region, the topology needs to be updated to avoid erroneous MM calculations. In principle, it is true that when using a subtractive QM/MM energy definition, the topology of the CPS should not matter as the MM energy of the CPS is subtracted from the MM energy of the ES. However, due to the structure of the NAMD program, this is not always the case. NAMD decomposes the MM system into discrete spatial units called patches for parallelization purposes. Patches comprise groups of heavy atoms as well as hydrogen atoms connected to them. If the topology of the system changes due to a chemical reaction, it is possible that atoms can move out of the spatial region of their patch. If this occurs, the MM energy of the CPS and/or ES will not be calculated as per the potential of the force field and the QM/MM calculation is no longer valid. Therefore, the topology of the system must be generated on the fly as the simulation progresses.
26


Updating the topology at every step is undesirable as it takes a significant amount of
time to write the ES PSF for large systems, a situation which is compounded when using optical hard drives (HDD) or when writing to a network file system (e.g. in HPC environments). Therefore, the topology for the ES is only updated when a bond is designated as formed/broken. Currently, the only automatic procedure for changing the topology implemented within QMMM deals with proton transfer and is discussed in Chapter IV. It is possible that there will be no universally applicable rule for handling topology changes, and new rules will need to be added on a case by case basis. When the topology is generated, lists of bonds, angles, dihedrals, impropers, and cross terms are created as described in the following sections.
2.4.1 The Bond List
The initial bond list is created by reading the list of bonds in the INBONDS section from the user-supplied PSF file. The bond list for the entire system is updated when a topology change occurs. Additionally, a separate bond list is created for the PS if the QM/MM partition cuts through covalent bonds.
2.4.2 The Angle List
The angle list is generated from the bond list. First the number of angles, iVangies, is calculated using number of bonds connected to the zth atom, bonds£, from the following formula:
Natoms
(2.1)
i=i
27


The angles list is equal to the set of all angles [/', i, k] defined by equation (2.4.2), where bondedi is the set of all atoms bonded to atom i:
l > j;Vj E bondedi, Vk e bondedi.Vi E atoms (2.2)
2.4.3 The Dihedrals List
The number of dihedrals can be calculated from the following formula:
Natoms bondsi
(bondsi 1) (bondsj l) (2.3)
i=i j=i
The dihedrals list is equal to the set of all dihedrals [k, /, /, l] defined by equation
(2.4.4):
j ^ i, k V- j, l ^ i, l ^ k\
Vi G atoms, Vy G bondedit\/k E bonded^, Vi G bondedj (2.4)
2.4.4 The Impropers and Cross Terms
No cross terms or impropers are automatically generated when the topology changes, as these are impossible to predict without specific knowledge of the molecule. These terms are read in from the user-supplied PSF and stored in memory. If the CPS contains all atoms of an improper term, with cap atoms corresponding to the atoms they substitute, that term is added to the CPS topology. Cross terms are added to the CPS topology in the same fashion.
If a section of the protein backbone is present in the CPS, cross terms which include the cap atom are not evaluated, and the CPS potential will not exactly cancel out the ES potential for subtractive schemes. This is an uncommon situation in QM/MM models.
28


2.5 Benchmarking
The following benchmark tests use QMMM version 2.2.7, a CUDA enabled version of NAMD 2.11, and Tinker 7.1.2. Each calculation consists of an MM MD calculation of a water box using the selected MM package with the TIP3P forcefield. Periodic boundaries are treated using the minimum image convention with the electrostatic interactions cut off at 14 A. The number of atoms in each benchmark varies. Full benchmarking details can be found in Appendix A.l.
The benchmarking data is shown in Table 2.1, with the NAMD speedup calculated as the average time per step of the NAMD calculation divided by the average time per step of the Tinker calculation (Equation 2.4.4),
NAMD Speedup
Average
Average
time
step
time
step
of NAMD simulation
of Tinker simulation
(2.5)
The benchmarking data indicates a linear speedup as the number of atoms in the system increases (Fig 2.3). The system sizes were chosen to approximate the number of atoms of various biomolecular systems that this methodology is currently applied to within the Lin Group. For small systems (<5,000 atoms), Tinker may still be a feasible solution for QM/MM dynamics. However, when conducting QM/MM dynamics simulations of larger systems, NAMD is the only feasible choice.
Table 2.1: Benchmarking Results for NAMD and Tinker Number of Atoms NAMD time/step (s) Tinker time/step (s) NAMD Speedup
4308 0.2 0.9 5
52233 0.6 35.2 64
96003 0.9 111.8 123
|l76454 1.6 330.9 221
29


Fig 2.3: Speedup when using NAMD versus the number of atoms. The speedup is calculated using equation (2.5)
30


CHAPTER III
IMPLEMENTATION OF PARALLEL ADAPTIVE PARTITIONING
3.1 Theoretical Background
As shown in section 1.3.2, the PAP potential can be expressed as thus:
N N N N-1 N N
= + + 2 £ Wf I~] U ft) + (3-1)
i=l i=l j*i i=l j=i+l kji
Where Pt is the smoothing function of the ith group based on the distance from the active zone center, VA is the QM/MM energy with only the active zone treated at the QM level, and VA n is the QM/MM energy with the active zone and the ith through nth groups treated at the QM/MM level. In the full PAP expansion, there are 2N terms where N is the number of buffer groups. The partitioning configurations are independent QM/MM calculations and may be calculated simultaneously.
3.2 Implementation of Parallel Permuted Adaptive Partitioning
The parallel PAP algorithm (PPAP) was implemented in our code, QMMM.
Currently, parallel energy and gradient calculations for MD simulations can be conducted using Tinker138 or NAMD137 for the MM package, and Gaussian 03,139 Gaussian 09,140 or MNDO143 for the QM package. QMMM was parallelized using the FORTRAN bindings from the OpenMP 3.0 protocol,146"147 a parallel language for writing parallel code on shared-memory multiprocessor architectures.
Fig 3.1 shows the flowchart for the implemented PPAP algorithm. After reading the initial coordinates and user-defined groups, QMMM obtains the entire system MM gradient for use in each of the individual QM/MM calculations. Next the program calculates at and
31


generates the list of groups in the active and buffer zones from this distance. Then, P£ is calculated for each of the buffer groups, and the list of partitioning configurations is generated based on the truncation order. The parallel region is then entered, and the main process spawns a number of threads based on an environmental variable set by the user. Each thread initializes a private data structure and calculates partitioning configurations at the QM/MM level. After that, the program enters a critical region, where the threads update the total energy and gradient in shared memory one at a time. Finally, the thread team is then destroyed and the program resumes normal function.
Fig 3.1: Scheme of parallel permuted adaptive partitioning implementation. Red boxes indicate serial instructions and green regions indicate parallel instructions
Work is assigned to the thread team using the dynamic scheduling option with a
chunk size of 1. This means that each thread is initially given a partitioning configuration to
calculate. When a thread finishes computing its assigned partitioning configuration, it is
32


given another until all partitioning configurations are calculated. Dynamic scheduling was chosen to allow for more flexibility in assigning work to the thread team, because the time to calculate a partitioning configuration is not known at runtime and may drastically vary between configurations. Due to the representation of numbers in finite precision in the computer, the computed gradients for a step may be slightly different if run with different numbers of threads or on different CPU architectures.
As the previous program design relied on FORTRAN common blocks for many energy and gradient functions, substantial sections of the entire program had to be rewritten to prevent data races in the parallel section. Additionally, the program was moved to dynamically allocating memory for most arrays at runtime, from using statically allocated arrays based on parameters in the code, which had to occasionally be modified. This change necessitated modifying almost every single subroutine within the program. As most of the static arrays for both the ES and CPS previously had one dimension determined by parameter set to be larger than the maximum number of atoms in the system, the memory profile of the program was reduced. This likely affected the PPAP algorithm the most, as these large arrays would have had to be continually reallocated for each member of the thread team.
3.3 Scaling and Benchmarks
Because the potential calculations in the PPAP calculation may be calculated simultaneously, the theoretical minimum time per step is proportional to the time of the longest QM calculation, since MM calculations usually take much less time. Additionally, this PPAP algorithm reaches its fullest potential when there are many buffer group molecules, as there are more partitioning configurations to assign to the thread team. Though
33


uncommon, it is possible that there are fewer QM calculations than available threads at a particular step. In this case, the extra threads do nothing.
The PPAP algorithm is benchmarked based on an MD simulation of 1435 H2O molecules and 1 H30+ molecule with the active zone center set to the location of the proton indicator. Rmin is set to 5 A and Rmax is set to 6 A yielding an active zone with a radius of 5 A encircled by a 1 A thick buffer zone shell. There were about 160 QM calculations per step during the simulation, with an average of 13.5 groups in the active zone and 17.5 groups in the buffer. The AP potential was truncated at the 2nd order. Due to the large number of QM calculations per step, this benchmark system is ideal for examining the parallel efficiency as the number of threads increase. Additional information about the benchmark systems can be found in Appendix A.2. Timing examples of ongoing experiments are also provided there, but the parallel efficiency is not rigorously benchmarked in them.
Under ideal calculation conditions, the PPAP algorithm scaled approximately linearly on a local workstation (Fig 3.2). These ideal conditions entail that there are enough calculations to spread across a thread team and that the calculations are run on a sufficiently fast disk drive. The strong scaling efficiency (SSE) can be used to measure how efficiently each thread in a thread team works compared to the serial calculation. The higher SSE values correspond to a more efficient thread team. The SSE is calculated using equation 3.2, where t1 is average wall time per step for the calculation with one thread, JV is the number of threads used in a calculation, and tN is the average wall time per step for the calculation with JV threads:
ti
Strong Scaling Efficiency = ---------* 100% (3.2)
N tN
34


Fig 3.2: Scaling of PPAP with increasing number of cores on local system described in Appendix A.3. Note that this computer only has 10 physical cores but 20 virtual cores. Left: Average wall time per step versus increasing number of cores. Right: Strong scaling efficiency (SSE) versus increasing number of cores.
The SSE can also be thought of as a metric for the parallel overhead for the calculation. The SSE in benchmark calculations ranges from 100% to 85% as the calculation was parallelized across all the physical cores on a local workstation (Fig 3.2). The time it takes to calculate one QM/MM potential is vastly greater than the cost of creating the thread team and initializing the private data structures. Because QMMM interfaces with the QM and MM packages by reading and writing files to disk, file EO congestion may contribute some overhead. This is evidenced by the difference in SSE between two calculations ran on the same computer but used file systems with 2 different speeds. In each system examined, the SSE is typically higher for the version ran on the faster filesystem than the slower file system (Fig 3.3)
The amount of file EO that occurs can even cause the simulation to crash if essential files become inaccessible by QMMM, or the MM or QM package due to file system traffic. This occurred frequently when using the LUSTRE network file system on the XSEDE Stampede HPC resource (Fig 3.3), and infrequently when using local file systems. The likelihood of a calculation crash increases as the filesystem speed decreases, i.e., calculations crash less on a SSD versus a HDD. This obstacle could be overcome by distributing the
35


Fig 3.3: Scaling of PPAP with increasing number of threads on the XSEDE Stampede HPC resource described in Appendix A.3. Note that this computer only lias 16 physical cores but 32 virtual cores. Tests on an internal node disk (blue) and network file system (red) were conducted. Tests with 32 threads crashed on the network file disk due to I/O congestion. Left: Average wall time per step versus increasing number of threads. Right: Strong scaling efficiency (SSE) versus increasing number of threads.
calculations and file operations across multiple nodes using the MPI protocol, by using multiple file systems on one node, or by developing an interface between packages that does not rely on files written to disk.
Interestingly, these benchmarks of the PAPP algorithm benefited from using virtual cores when they were run on an SSD. This may be attributed to the fact that many of the benchmarking QM calculations were on the order of 200 ms. Therefore, some threads were performing file system operations while others were performing QM calculations. For longer QM calculations, this may not be the case as QM calculations will be running on the all the virtual cores simultaneously. Typically, when running compute intensive processes, using more virtual cores than physical cores results in no speedup or slowdown.
Maximizing parallel performance also requires specific knowledge of the QM or MM package. If the SSE of QMMM at a certain thread count is less than the SSE of a package, then less threads of QMMM should be used and more threads of the package should be used. Timing data for MNDO using QMMM and Gaussian 09,148 indicate that QMMM should be run with as many threads as possible before disk EO speed becomes in issue, and that the
36


QM packages should use the remaining threads. For systems with few partitions, the QM or MM packages can be assigned more threads to avoid idle processors. Additionally, if the thread team has more threads than there are QM calculations, the SSE of the algorithm will suffer (Fig 3.4).
1.5 A Buffer 0.5 A Buffer
Fig 3.4: Strong scaling efficiency of two QM calculations with different buffer sizes. The buffer sizes were 1.5 A and 0.5 A which corresponded to an average of 115.0 and 9.3 QM calculations per step respectively. As there were not enough partitioning configurations to spread to the thread team, the SSE of the 0.5 A buffer size calculations quickly decreases.
37


CHAPTER IV
IMPLEMENTATION OF INDICATOR RESTRAINT
4.1 Restraint Potential and Decomposed Gradients
A harmonic restraint was chosen for the proton indicator for its intended usage with umbrella sampling and the Weighted Histogram Analysis Method (WHAM).149 The restraint potential can be expressed as:
Vhar is the harmonic restraint potential imposed on the indicator I, k is the force constant of the restraint potential, AIq is the origin of potential in Cartesian coordinates, and Xj is the vector for the current position in Cartesian coordinates of the proton indicator.
The position of the indicator, Xl5 is a function of the donor, D, and the possible acceptor atoms, Ay. The gradient is taken with respect to both of these categories of atoms. The calculation of the normalization factor for the weight function, gh and weight function g(Pmj) are treated as constants when taking the derivatives at any given time step of the simulation. When the gradient of these values was taken with respect to XD, XA and XHm, the gradients were too large, causing simulation artifacts. (Appendix B.5).
The gradient of the potential with respect to the /th acceptor Ay or donor D can be expressed by the following equations:
(4.1)
*(Xi X,0)

(4.2)
38


(4.3)
dVhar / \ 5Xj
= k(x, -x,0)
axr
aXr
dfr g{pmj)
dXA/ #i
_ J_
3Xd gi
(4.4)
(4.5)
The symbols XA and XD represent the location of a possible acceptor Ay or the donor D respectively. The direction of the gradients on all XA or XD are along to the vector (Fig 4.1):
(X,-XlQ) (4.6)
The form of g(pmj) (Equation 1.4.4) is designed such that atoms receive none of the gradient when they are further than rLISX away from the donor, smoothly receive gradient as they become an acceptor, and have the same gradient in the limit of the donor acceptor swap.
O Hydrogen ^ Indicator A Acceptor ^ Donor X 0r'gin
Fig 4.1: 3d depiction of gradients on the acceptor (red triangle), donor (yellow triangle), and indicator (green triangle) for the Zundel (Left) and Eigen (right) cations. The potential origin is the black X. The gradient unit vectors are represented by the arrows. The gradients all have the same direction but may have different magnitudes.
39


4.2 Behavior of the Indicator for the Zundel Cation
Single-point calculations were performed on the H5O2 (Zundel) structure to examine the gradients on the indicator for a simple system. The geometries were constructed such that D, A, and the transferring proton H are aligned along the x axis. Three different sets of singlepoint calculations were performed. In scan A (Fig 4.2), the shared proton was translated between the two oxygen atoms. In scan B, a H2O molecule was translated from rLISX towards an H30+ molecule to create the Zundel cation. In scan C, the entire system was translated with respect to the restraint potential origin, XIq.
Fig 4.2: Depictions of two scans. In scan A, a proton is translated across the x axis between two water molecules. In scan B, the blue water molecule is translated across the x axis from rLISX towards Oi. Two origin locations are chosen for each scan, one at Oi and one at H. Red spheres are oxygen and white are hydrogen
Scan A (Fig 4.3) shows the decomposed gradients with the respect to the translation of the proton. During the scan, the indicator changed position with the proton, but the potential origin stayed static. Two origins of potential were investigated, one placed on the location of the shared hydrogen (H), and one placed on an oxygen (Oi). The displacement Sx between the indicator and origin is calculated by:
8x Xorjgjn ^indicator (4-7)
The graphs show that when the donor and acceptor swap (at Sx = 0.0 A for the origin at H, and Sx = -1.2 A for the origin at Oi), the gradients are equivalent. Additionally, the
40


Fig. 4.3: For scan A: distance of indicator from origin (top) and decomposed gradients for the Zundel cation (bottom) with the indicator centered at the transferring proton (left) and the indicator centered at the donor oxygen (right). The indicator is in red, donor green, and acceptor blue. Note the y-axis is different for the gradient graphs. The force constant is 5 kcal mol"1 A'2.
Fig. 4.4: For scan B: distance of indicator from origin (top) and decomposed gradients for the Zundel cation (bottom) with the indicator centered at the transferring proton (left) and the indicator centered at the donor oxygen (right). The indicator is in red, donor green, and acceptor blue. Note the y-axis is different for the gradient graphs. rLISX is 3.5 A and the force constant is 5 kcal mol'1 A-2
41


Origin at Oi gradient graph (Fig 4.3) shows that this transition occurs smoothly when the origin potential is located away from the indicator.
Scan B (Fig 4.4) depicts the gradients as H2O is translated away from the complex from 1.5 A to past rLISX = 3.5 A. The complex starts as the donor at 1.5 A and becomes an acceptor at 2 A. Its gradient tends towards zero as it reaches rLISX. As the H2O translates, the g(Pmj) value of the acceptor constantly changes. The ratio of the gradients decomposed onto
the donor and acceptor are proportional to the magnitude of the g(pmj) f the acceptor, with a larger g(pmj) leading to more gradient decomposed onto the acceptor. The gradient of the donor will always be larger in magnitude than the gradient of the acceptor except for at a donor acceptor swap, when they are equal.
When the potential origin is located at Oi and the complex is translated past rLISX, the gradients on all atoms remain at 0 and the gradient is continuous. This is because the D H bond length is equal to the r^H parameter. If this is not the case, a discontinuity in the position of the indicator, and thus the gradients, occurs (Fig 4.5). While this happens to a minor extent throughout a simulation, a bond length would have to be substantially different from its equilibrium bond length to create a drastic discontinuity.
Scan C (Fig 4.6) shows the gradients decomposed onto a Zundel complex that is
Fig 4.5: A 0.4 A jump in the indicators position due to the OiH bond length of 1.2 A which is larger than the rLISX parameter of 1.0 A.
42


Sx [A]
^Acceptor Donor Indicator
Fig 4.6: For scan C: total gradient on indicator I and decomposed gradients on D and A with respect to the change in the origin of the restraint potential along the x axis Sx. Red spheres indicate oxygen, white spheres hydrogen, and the pink sphere is the indicator. The force constant is 5 kcal mol"1 A'2
translated with respect to a static restraint potential origin. In this scan, the g(pmj) value of the acceptor is constant. Thus, changing the position of the complex changes the magnitude of the donor and acceptor gradients, but the ratio of the gradients remains constant.
4.3 Behavior of the Indicator for Collision Simulations
NVE MD calculations were conducted on a system where one H2O molecule was placed 5.40 A from H?03+ cluster and given an initial velocity of 0.01 A fs"1 or 0.001 A fs"1 (Fig 4.7). The velocities of all other atoms were set to 0 A fs"1. All atoms were described using the OM3n potential. The simulations were run for 2 ps with 0.2 fs time steps. Four different restraint force constants were examined. The full set of simulation results can be found in Appendix C.l
The restraints keep the indicator closer to the origin of the potential after the collision than in the free simulation (Fig 4.8). When the restraint is active, the distance from the potential of origin is lower than in the absence of the restraint. Also, the total energy of the system increased when the restraint was active (Fig 4.9). Therefore, these restraints should
43


ro
-t'
w
b
6
4
Time (ps)
Fig 4.8: Distance (A) of indicator from origin of potential for Zundel (top) and Eigen(bottom) cations. Red, blue, green, and cyan correspond to force constants of 0 (no restraint), 5, 10, and 20 kcal mol"1 A'2 respectively. The initial velocity of the EEO was 0.001 A fs"1.
7
Time (ps)
k = o k = 5 k = 10 k = 20
Fig 4.9: Total energy of system for simulations of the Zundel (A) and Eigen (B) cations with different force constants. The blue, green, red, and cyan lines correspond to force constants of 0 (no restraint), 5, 10, and 20 kcal mol'1 A"2 respectively. The initial velocity of the EbO was 0.001 A fs'1.
44


only be used in simulations coupled to a thermostat. Additionally, in the set of 0.01 A fs"1 simulations for the Zundel cation, the decomposed gradients were strong enough to cause erroneous structures, leading to SCF convergence failures (Appendix C.l). However, these structures were not observed in the Eigen set. This is likely because there are more atoms to distribute the gradient to in the Eigen cluster, leading to reduced gradients on the individual atoms. Thus, low force constants should be used when the indicator can move far from the origin and there are few possible acceptors, for example, when studying proton diffusion in confined spaces.
4.4 MD Simulations of a Water Sphere
MD simulations of a sphere of H2O molecules surrounding 1 H3CE cation were conducted with and without the use of the indicator restraints to compare the difference in potential of mean force (PMF) for proton transfer and radial distribution functions (RDF) of water. Four simulations were conducted with the force constant for the indicator restraint set at 5 kcal mol A"2. For comparison, four simulations were conducted in the absence of the indicator restraint, resulting in a total of 8 simulations. The simulations used a droplet model similar to Wu et al.134 with the inner 12 A sphere described at the OM3n134 level and the outer 10 A shell described at the MM level with the SPC/Fw potential.150 Restraints with a force constant of 20 kcal mol"1 A'2 were applied to every MM oxygen atom. The full simulation details can be found in Appendix C.2.
The radial distribution functions (RDF) were calculated using saved snapshots from the simulations, excluding those from the first 0.5 ps. The radial distribution function was centered on the donor D. If the donor atom was located greater than 7 A away from the center of mass of the droplet (within 5 A of the QM/MM border), the snapshot was excluded
45


from the analysis to minimize the influence of MM waters on the RDF. In the simulations with no indicator restraint, the indicator quickly migrated far from the center (Fig 4.10), and only 7455 out of 20000 possible snapshots were usable. All snapshots were usable when the restraints were used. The D-0 RDF (Fig. 4.10) indicates close agreement between the simulations conducted with and without the restraint. The noise in the free line comes from the reduced amount of sampling due to the exclusion of snapshots. This agreement suggests the presence of the indicator restraint does not significantly alter the solvation shell surrounding the proton.
12
o
c
ro
-I'
(f)
b
D) <
A 1.6 1 1 B I -
J ^Restrained 1.2 1 ^Restrained
9 0.8 O) 0.4 1 Free
0.0 1 1
0
4 8 12 012345
Time (ps) r [A]
Fig 4.10: (A) Average distance of indicator from restraint potential origin (center of cluster) for restrained (blue) and unrestrained (red) simulations. (B) Area normalized radial distribution function go-o(r) for restrained (blue) and unrestrained (red) simulations, where r is the distance from the donor, D.
The potential of mean force (PMF) for the proton transfer between the donor oxygen and the most likely acceptor oxygen was calculated through Boltzmann inversion for two different reaction coordinates (Fig 4.11). The same set of snapshots was used for the PMF as the RDF. The most likely acceptor A was the acceptor with the highest pmj value at a given time step. The p reaction coordinate was the pmj value for the most likely acceptor in a snapshot calculated in equation (1.4.2). The Sr reaction coordinate was calculated as the
46


absolute value of the difference between the donor-proton distance DH, and acceptor-proton
distance AH (Equation 4.8).
Sr = |rDH -rAH|
(4.8)
Fig 4.11 Potential of mean force (PMF) for proton transfer for reaction coordinates (A) Sr and (B) p. The blue line corresponds to the PMF with the indicator restraint active, while the red line corresponds to a freely moving indicator.
The barrier heights for the PMF (Table 4.1) are measured at the reaction coordinates Sr = 0 A and p = 0.50, and correspond to a Zundel state. The barrier height for the restrained and free simulations are similar, with the difference in barrier height for p at 0.03 kcal mol"1 and for Sr at 0.11 kcal mol"1. The barrier heights similar to the previously reported barrier heights for the OM3n potential in the cluster model of 0.19 and 0.06 kcal mol"1 for p and Sr respectively by Wu et al.134 as well as for those obtained from path integral calculations of 0.16 kcal mol'1 for Sr,38'39 The minima at Sr = 0.25 A agrees with previous calculations38'39-%. 134 an(j represents the Eigen resting state. Thus, the use of the indicator restraints does not seem to significantly affect the free energy profile for proton transfer in water.
47


Table 4.1: Barrier heights for PMF calculations of proton transfer in a water cluster. The barrier height is the energy that corresponds to the reaction coordinate of Sr = 0 or p = 0.50
Reaction Coordinate Restrained/Free Barrier Height (kcal mot1)
Sr Restrained 0.15
Sr Free 0.26
P Restrained 0.22
P Free 0.25
48


CHAPTER V
APPLICATION OF INDICATOR RESTRAINT TO CARBON NANOTUBE
5.1 Introduction
A model system of a carbon nanotube embedded between two parallel graphene sheets similar to Peng et al.36 was constructed for three reasons: (1) to examine the behavior of the indicator in a greatly simplified model of a proton transport protein, (2) to explore and compare hydration dynamics of the tube with Peng et al,36 and (3) to examine the viability of obtaining a PMF for a proton transfer process using the indicator restraints. The carbon nanotube represents a simplistic model of a hydrophobic cavity surrounding a proton transfer relay,151 such as in the E. Coli CLC Chloride Ion Transport Protien.152 Thus, examining the hydration dynamics of the carbon nanotube may provide insights into the mechanism of water wire formation in hydrophobic cavities in biomolecular systems.
The nanotube system, based on Peng et al.,36 contains a single walled, armchair-type (6,6) carbon nanotube with a length of 28.2 A and diameter of 8.1 A. The nanotube bridges two graphene sheets with x andy dimensions of 42.99 A and 41.12 A, respectively (Fig 5.1). The nanotube dimensions allow for the formation of a 1-dimensional water wire along the tube axis, and have been used in several previous experiments to study proton transfer through a water wire.36,95l5K 153 A slab of water with 7 pairs of K /Cl ions was added to each side of the nanotube. The SPC F/w150 water model was used in conjunction with the CHARMM22154 parameters for the carbon atoms (atom type CA) and ions. Similarly to Peng et al.,36 the s parameter of the carbon nanotube atoms was changed to -0.0448 kcal mol"1 to scale the depth of the interpolated potential well for the vdW interaction between the oxygen and carbon atoms to 80% of its original value. However in this system, the 8 parameter of the
49


graphene sheets was also adjusted while in Peng et al.,36 the standard CHARMM22 parameter was used. This increased the hydrophobicity of the tube, which would otherwise be able to spontaneously hydrate. Visual inspection of the MM trajectories confirmed that only transient water wires of 2 4 waters could form. This system was then equilibrated as per the procedure described in Appendix C.3.
r
X
Fig 5.1: Carbon nanotube system used throughout the chapter. The nanotube axis is aligned with the z axis.
The initial geometries for the following QM/MM simulations obtained from the end of the equilibration runs. For the following QM/MM simulations, the MM waters and K+/CT ions were restrained with a force constant of 1 nN (approximately 14.39 kcal mol"1 A"2). No K /Cl ions were included in the QM region. The positions of all carbon atoms were fixed. All simulations were run under the NVT ensemble using the Nose-Hoover thermostat at 300 K, with no treatment of periodic boundaries, and with 0.25 fs time steps.
Peng et al.36 reported that an excess protonic charge restrained at the opening of the nanotube could recruit waters from the bulk into the tube. They demonstrated this behavior using the SPC/Fw water model and the multistate empirical valence bond version 3 (MS-EVB3) model110 with the multiscale reactive molecular dynamics (MS-RMD) method.103'104-
50


108, no, 155 jhgy use(j harmonic restraints on the Center of Excess Charge108 (CEC) of the system, which tracks the approximate location of the excess proton and differs from the indicator restraints described in Chapter IV. The CEC is calculated by the following equation:
where rlcoc is the center of charge vector of the H3CC ion in the ith MS-EVB state, and cf is the amplitude of the ith MS-EVB state.
Peng et al. described the water recruitment process as follows. Two waters are initially in the nanotube. When the proton nears the nanotube, it brings its solvation shell with it, filling the tube with two more waters for a total of four waters. The excess proton becomes trapped in an energy minima at the pore opening. Waters transfer through the proton into the nanotube as the proton enters the nanotube, and the tube is fully hydrated by the time the proton moves 2-3 A into the tube. This work aims to examine whether this process occurs in QM/MM simulations with the OM3n potential and with the previously described indicator restraints.
In this set of simulations, the restraint potential origin was placed 2.5 A above the center of the nanotube to look for evidence of water recruitment into the nanotube (Fig 5.2). The closest water to the restraint potential origin was protonated. To examine the effects of water in the tube, the number of waters initially in the tube was varied from zero to four. Four starting geometries with different atomic coordinates were created for each of these
N
t = l
5.2 Simulations of Indicator Restrained at Nanotube Opening
51


cases, for a total of 20 simulations. For each starting geometry, two simulations were conducted: one with the force constant set to 5 kcal mol"1 A"1, and one with the force constant set to 0 kcal mol'1 A'1. All waters within a 14 A radius of the middle of the pore opening were included in the QM subsystem. The simulations were run for up to 10.5 ps. The full simulation details can be found in Appendix C.4.
For simulations with the indicator restraint, waters begin to be recruited into the nanotube within the first few hundred femtoseconds (Fig 5.3). This begins with migration of an excess proton towards bulk. As the excess charge moves towards the bulk solution (in the positive z direction), the indicator moves with it and away from the restraint potential origin. Forces on donor and acceptor atoms, which were decomposed from the restraint force on the indicator, move the donor and acceptors towards the nanotube (in the negative z direction) in order to push the indicator towards the restraint origin. This moves waters into the tube. The proton then migrates back towards the bulk and the process repeats itself, eventually creating a water wire in the nanotube.
52


1 O 0 to Vjt

4 O / 0 5 n 5 w 3 4 9 Q 6 v. U
.. .
Wf
Fig 5.3: Water recruitment into the pore when restraints are applied to the indicator. The images depict the indicator (pink sphere), donor (blue sphere), most likely acceptor (yellow sphere), restraint potential origin (green sphere), water (ball and stick, red is oxygen and white is hydrogen), nanotube (cyan lines), and restraint forces (black arrows). Some waters are omitted for clarity. Restraint forces are shown in panels two and three.
(1-0 fs): The initial position of the indictor (pink) begins close to the restraint potential origin (green). Zero waters are initially inside the nanotube.
(2-2 fs): A Zundel cation is formed and the indicator migrates towards the acceptor (yellow). Forces (black arrows) are applied to the donor (blue) and acceptor in the direction of the nanotube.
(3-38 fs): The donor and acceptor have switched places. New waters have been recruited as possible acceptors. Forces are applied to a new list of possible acceptors.
(4-174 fs): The indicator continues to be drawn towards the bulk. The original donor has moved to the tube entrance.
(5 252 fs): The original donor has moved inside the tube. The acceptor from (4) is moving towards the tube opening.
(6 480 fs): Three waters are inside the nanotube. At the end of this simulation, 12 waters have been recruited into the tube.
53


For simulations with no restraint, the indicator quickly moved to the outside of the QM region (Fig 5.4), indicating a strong driving force for the indicator to move away from the tube opening towards bulk. Some simulations were terminated due to QM SCF convergence difficulties from a diffuse water wire which formed in the tube (Fig 5.4). This occurred due to water shooting into the nanotube, or initial waters in the tube drifting towards the other side of the nanotube. These simulations were truncated to a length less than 10.5 ps.
o 4 8 12 o
Time (ps)
4 8 12
Time (ps)
Fig 5.4: (A) Average distance of indicator from restraint potential origin versus simulation time for simulations that were restrained (blue) and unrestrained (red). (B) Picture showing diffuse water wire causing SCF convergence failure. (C, D) Average water occupancy of nanotube over simulation time for simulations with indicator restraints (C) and without restraints (D) for different numbers of initial waters in the tube.
The water occupancy (WO) of the nanotube, defined as the number of oxygen atoms within the tube at a given time step, was much higher for simulations with the indicator restraint than for simulations with no restraint (Fig 5.4). This indicates that the proton indicator restraints play a large role in recruiting waters into the nanotube. When the restraint
54


was applied, the water occupancy of the the nanotube increased within the first two picoseconds. The water occupancy plateaued at an average of 11 waters, which is close to wet minima for the tube of 13 reported by Peng et al.36. In the absence of indicator restraints, water recruitment into the nanotube was less frequent. When zero waters were initially in the tube, no water recruitment was observed. When one water was initially in the tube, on average, the waters exited the nanotube. However, when 2-4 waters were initially in the nanotube, water recruitment into the nanotube was observed, albiet at a much slower pace than when restraints were applied. As the average indicator position was approximately 10 A from the restraint potential origin and 12.5 A from the entrance of the tube after 4 ps of simulation time, it is unclear whether the presence of the excess proton around the nanotube entrance aids water recruitment into the tube.
5.3 Simulations of Indicator Restrained at Different Distances from Nanotube
The effect of different indicator restraint potential origins with respect to the entrance of the nanotube was examined to determine the effect of the restraint potential origin with respect to the water occupancy. The simulation setup was similar to that used in Section 5.2, with a QM hemisphere of water located above one nanotube entrance. However, the restraint potential origin was moved 3.5, 5.0, and 6.5 A above the center of the nanotube entrance (z direction). Two other potential origins with the origin shifted 2.5 A to the side of the nanotube opening (x direction) with heights of 3.5 and 5.0 A above the nanotube entrance (z direction) were also created for a total of five sets of simulations (Fig 5.5). Four different starting geometries were created for each set, for a total of 20 simulations over all sets. No further distances could be explored without the risk of significant boundary effects from the
55


QM/MM border. There were no waters in the tube at the start of the simulations. Complete simulation details can be found in Appendix C.5.
Fig 5.5: Graph of average water occupancy at different restraint locations. The color of the line corresponds with the sphere above the nanotube to indicate the location of the restraint potential origin. The legend indicates the distance above the nanotube the indicator was restrained, and the presence of an x indicates that the restraint was shifted 2.5 A in the x direction from the center of the tube.
For all the restraint potential origin locations, the average water occupancy over time increased (Fig 5.5). This indicates that the water molecules were recruited into the tube even when the indicator is restrained at a location of up to 6.5 A above the tube. As the average water occupancy without the restraint was 0 waters (Fig 5.4), the increase in water can be attributed solely to the presence of the restraint. The average water occupancy after 8 ps for these restraint potential origins (<8.75 waters) was less than the water occupancy when the restraint was 2.5 A above the tube (9.5-11 waters, Fig 5.4). No direct dependence of water occupancy on vertical distance above the tube can be drawn from this data. When the restraint potential origin was shifted away from the center of the tube, there was less water recruitment than when the restraint potential was above the tube. However, as the radius for enlisting possible acceptors was 3.5 A in the simulations, acceptors above the tube opening could be subjected to restraint forces and drawn into the nanotube.
56


5.4 QM/MM Umbrella Sampling of PMF for Proton Transfer in Nanotube
To provide proof-of-concept for using the indicator restraints for umbrella sampling, the nanotube was filled with water and a PMF for proton transfer was obtained using the z coordinate of the indicator as the reaction coordinate. Each umbrella window used different QM/MM partitions (Fig 5.6). Waters on both sides of the nanotube were described as QM. Force constants ranging from 2.5 7.5 kcal mol"1 A'2 in the z direction were used. The full set of simulation details can be found in Appendix C.6.
Fig 5.6: (A) Location of umbrella window origins along the reaction coordinate (z axis). Horizontal red lines indicate umbrella window locations and the black dot is the tube center. (B) Lowest umbrella window origin (z = 14) and (C) highest umbrella window origin (z = 30). The QM waters are depicted. The circle indicates the umbrella window origin.
The barrier height for proton transfer in the nanotube from z = 30.1 to z = 14.1 from the PMF (Fig 5.7) is 119.9 kcal mol"1, a value that is certainly too large. The PMF is relatively flat from z = 14 A (middle of nanotube) to r = 21 A. During exploratory calculations, no waters were recruited into the tube during this region. However, as the indicator moved past this region, waters were vigorously recruited into the nanotube, shuttling waters from the QM hemisphere in the positive z direction towards the negative z direction. This effect was so pronounced, that when z > 25 A, close to 20 waters were
57


16 20 24 28
z Coordinate [A]
Fig 5.7: (A) Histogram of samples along the umbrella reaction coordinate. Different colors indicate different umbrella window origins. (B) PMF for proton transfer from the middle of the nanotube towards the bulk. (C) Waters recruited through tube at different umbrella origins in preliminary umbrella sampling calculations. Positive water recruitment indicates net water movement from the positive z direction into the tube, negative water recruitment is in the opposite direction. (D) Restraints applied to the umbrella window centered at z = 14 A. Restrained and unrestrained QM waters are shown with blue and red respectively
58


recruited into the tube, leading to an unphysical density in the upper half of the tube. Therefore, the lower half of the QM subsystem was restrained to prevent this recruitment (Fig 5.7), increasing the PMF by approximately 40 kcal mol"1 in this region (Appendix C.6).
The mechanical embedding scheme may be another contributing factor for the steep PMF curve. As the indicator migrates into bulk, the PMF curve should flatten, not steepen. There is thus likely a strong driving force for the indicator to continue moving to the outside of the QM region. This is evidenced by the migration of the indicator to the outside of the QM region in unrestrained simulations for the water cluster (Fig 4.9) and nanotube systems (Fig 5.4). This driving force exacerbates the water recruitment mechanism observed in these simulations, but it is unknown what the cause of this force is. Further analysis and calculations using more accurate embedding schemes and different potentials will be necessary to determine the causes of this behavior.
Despite the problems in the model, these simulations do shed light on the basis of the restrained shooting mechanism. These results show that for the speedy recruitment of water into a dry nanotube, the excess proton must move away from the nanotube. This causes forces to be applied to the surrounding waters in the direction of the nanotube, shuttling the water into the tube. However, with this system, it appears that the probability of finding the excess charge around the nanotube opening is extremely low. Therefore, the recruitment of water into the nanotube with the excess charge defect may be more of a simulation artifact than feasible mechanism.
This mechanism becomes more feasible if there is an energy minima for the excess proton close to the nanotube opening, as is the case with the MS-EVB model,36,153 because a proton needs to be by the pore opening in order to recruit waters into it. Nevertheless, the
59


potential is asymmetric with respect to the minima at the pore opening, with a steeper barrier for motion into the tube and a shallower barrier towards bulk. Therefore, the excess charge will, on average, drift towards bulk, causing water to be recruited into the tube. For future studies, a DFT or ab initio potential could be used in order to confirm whether such a minima exists.
Nevertheless, these results indicate that this indicator restraint may become a viable option for obtaining the PMF for chemical processes involving proton transfer. As the distribution of forces is straightforward, it is easy to describe the bias applied to the simulation. However, further investigation is necessary to determine the exact causes behind the PMF generated from the restraints.
60


CHAPTER VI
CONCLUSION
In this work, three primary topics were considered. First, an interface between NAMD and QMMM was created to increase the speed of MM calculations. Benchmark calculations for MM MD simulations show that the speedup ranges from a factor of 5 for small calculations to >200 for larger calculations. In order to use the interface for QM/MM simulations where bonds were formed and broken, algorithms for the on-the-fly generation of the topology were implemented.
Next, the permuted adaptive partitioning algorithm present in QMMM was parallelized using the OMP protocol. The partitioning configurations may now be calculated simultaneously when using Tinker and NAMD for the MM package, and Gaussian 03 or 09 and MNDO for the QM package. The parallel speedup is approximately linear when there are enough QM/MM calculations to distribute to the thread team, and when the calculation is conducted on a sufficiently fast file system to minimize the overhead of file handling. The QMMM package was also updated to dynamically allocate memory for most arrays, making the program more user friendly and reducing the memory profile of the application.
Finally, a harmonic restraint potential was applied to the proton indicator. The gradients of the restraint potential indicator were decomposed onto the atoms which determine the location of the indicator. The gradients on the donor and acceptor atoms have the same direction of the gradient on the indicator. The ratio of the decomposed gradients depends on the ratio of the pmj values for the acceptors. The indicator restraints do not conserve energy. However, the PMF and radial distribution function for water was not significantly affected by the use of the restraints.


The indicator restraints were applied to a carbon nanotube system in order to provide proof of concept for using the restraints to determine the PMF for a process with a high energy barrier. These simulations confirmed results by Peng et al.,36 which indicated that excess protons restrained near the opening of a hydrophobic space cause water to enter that space. The recruitment of water by restraint of the excess charge has now been confirmed with two different simulation techniques: through our work with QM/MM and a proton indicator, and by Peng et al.36 with MS-EVB and restraints on the center of excess charge. The reasons for water recruitment are not yet fully understood, but it is apparent that the proton indicator felt a strong driving force to leave the nanotube opening. We suspect that the motion of the indicator coupled with the restraint forces recruited water into the tube, and the behavior could be a simulation artifact instead of a physical process. More research is needed to understand these observations.
62


REFERENCES
(1) Warshel, A., Multiscale modeling of biological functions: From enzymes to molecular machines (Nobel Lecture). Angew. Chem. Int. Ed. Engl. 2014, 53 (38), 10020-
10031.
(2) The Nobel Prize in Chemistry 2013.
https://www.nobelprize.org/nobel prizes/chemistry/laureates/2013/ (accessed Jan 2017).
(3) Senn, H. M.; Thiel, W., QM/MM methods for biological systems. In Atomistic Approaches in Modern Biology: From Quantum Chemistry to Molecular Simulations, Reiher, M., Ed. Springer Berlin Heidelberg: Berlin, Heidelberg, 2007; pp 173-290.
(4) Elliott, J. A., Novel approaches to multiscale modelling in materials science. Int Mater Rev 2011, 56 (4), 207-225.
(5) Rudd, R. E.; Broughton, J. Q., Concurrent coupling of length scales in solid state systems. Phys. Status SolidiB 2000, 277 (1), 251-291.
(6) Wu, X.; Koslowski, A.; Thiel, W., Semiempirical quantum chemical calculations accelerated on a hybrid multicore cpu-gpu computing platform. J. Chem. Theory. Comput. 2012,5(7), 2272-2281.
(7) Siegbahn, P. E.; Himo, F., Recent developments of the quantum chemical cluster approach for modeling enzyme reactions. J. Biol. Inorg. Chem. 2009,14 (5), 643-651.
(8) Ong, M. T.; Verners, O.; Draeger, E. W.; van Duin, A. C.; Lordi, V.; Pask, J. E., Lithium ion solvation and diffusion in bulk organic electrolytes from first-principles and classical reactive molecular dynamics. J. Phys. Chem. B. 2015,119 (4), 1535-1545.
(9) Zhao, G.; Perilla, J. R.; Yufenyuy, E. L.; Meng, X.; Chen, B.; Ning, J.; Ahn, J.; Gronenbom, A. M.; Schulten, K.; Aiken, C.; Zhang, P., Mature HIV-1 capsid structure by cryo-electron microscopy and all-atom molecular dynamics. Nature 2013, 497 (7451), 643-6.
(10) Shaw, D. E.; Maragakis, P.; Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Eastwood, M. P.; Bank, J. A.; Jumper, J. M.; Salmon, J. K.; Shan, Y.; Wriggers, W., Atomic-level characterization of the structural dynamics of proteins. Science 2010, 330 (6002), 341-346.
63


(11) Martins-Costa, M. T.; Ruiz-Lopez, M. F., Reaching multi-nanosecond timescales in combined QM/MM molecular dynamics simulations through parallel horsetail sampling../. Comput. Chem. 2017, 38 (10), 659-668.
(12) 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 (2), 227-249.
(13) 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.
(14) Gao, J. L., Hybrid quantum and molecular mechanical simulations: An alternative avenue to solvent effects in organic chemistry. Acc. Chem. Res. 1996, 29 (6), 298-
305.
(15) Gao, J.; Furlani, T. R., Simulating solvent effects in organic chemistry: Combining quantum and molecular mechanics. IEEE Comput. Sci. Eng. 1995, 2 (3), 24-33.
(16) 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 (2), 185-199.
(17) Mulholland, A. J., Chapter 14 The QM/MM approach to enzymatic reactions. In Theoretical and computational chemistry, Leif, A. E., Ed. Elsevier: 2001; Vol. Volume 9, pp 597-653.
(18) Senn, H. M.; Thiel, W., QM/MM methods for biomolecular systems. Angew. Chem. Int. Ed. Engl. 2009, 48 (7), 1198-1229.
(19) Riccardi, D.; Schaefer, P.; Yang, Y.; Yu, H.; Ghosh, N.; Prat-Resina, X.; Konig, P.; Li, G.; Xu, D.; Guo, H.; Elstner, M.; Cui, Q., Development of effective quantum mechanical/molecular mechanical (QM/MM) methods for complex biological processes. J. Rhys. Chem. B 2006,110 (13), 6458-6469.
(20) Lyne, P. D.; Hodoscek, M.; Karplus, M., A hybrid QM-MM potential employing Hartree-Fock or density functional methods in the quantum region. J. Phys. Chem. A 1999,103 (18), 3462-3471.
64


(21) Monard, G.; Prat-Resina, X.; Gonzalez-Lafont, A.; Lluch, J. M., Determination of enzymatic reaction pathways using QM/MM methods. Int. J. Quantum Chem. 2003, 93 (3), 229-244.
(22) Mackerell, A. D.; Feig, M.; Brooks, C. L., Extending the treatment of backbone energetics in protein force fields: Limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J. Comput. Chem. 2004, 25 (11), 1400-1415.
(23) MacKerell, A. D.; Feig, M.; Brooks, C. L., Improved treatment of the protein backbone in empirical force fields. J. Am. Chem. Soc. 2004,126 (3), 698-699.
(24) Maple, J. R.; Dinur, U.; Hagler, A. T., Derivation of force fields for molecular mechanics and dynamics from ab initio energy surfaces. Proc. Natl. Acad. Sci. U.S.A. 1988, 85 (15), 5350-5354.
(25) Mackerell, A. D., Jr., Empirical force fields for biological macromolecules: Overview and issues. J. Comput. Chem. 2004, 25 (13), 1584-1604.
(26) Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D., Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone c(), \\i and side-chain y\ and %2 dihedral angles. J. Chem. Theory Comput. 2012, 8 (9), 3257-3273.
(27) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M., CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 1983, 4 (2), 187-217.
(28) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J., Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996,118 (45), 11225-11236.
(29) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J.; Kollman, P., A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations. J. Comput. Chem. 2003, 24 (16), 1999-2012.
(30) Schmid, N.; Eichenberger, A. P.; Choutko, A.; Riniker, S.; Winger, M.; Mark, A. E.; van Gunsteren, W. F., Definition and testing of the GROMOS force-field versions 54a7 and 54b7. Eur. Biophys. J. 2011, 40 (7), 843-856.
65


(31) Jensen, F., Introduction to computational chemistry. John Wiley & Sons: 2016.
(32) Szabo, A.; Ostlund, N. S., Modern quantum chemistry: Introduction to advanced electronic structure theory. Courier Corporation: 2012.
(33) Fock, V., Naherungsmethode zur losung des quantenmechanischen mehrkorperproblems. Z. Phys. 1930, 61 (1), 126-148.
(34) Hartree, D. R., The wave mechanics of an atom with a non-coulomb central field. Part I. Theory and methods. Math. Proc. Cambridge 1928, 24 (01), 89-110.
(35) Moller, C.; Plesset, M. S., Note on an approximation treatment for many-electron systems. Phys. Rev. 1934, 46 (7), 618-622.
(36) 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 (29), 9212-9218.
(37) Scholten, M., Semiempirische verfahren mit orthogonalisierungskorrekturen: Die om3 methode. Ph. D. Dissertation, Max-Planck-Institut fur Kohlenforschung 2003.
(38) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M., The nature of the hydrated excess proton in water. Nature 1999, 397 (6720), 601-604.
(39) Marx, D.; Tuckerman, M. E.; Parrinello, M., Solvated excess protons in water: Quantum effects on the hydration structure. J. Phys. Condens. Matter 2000,12 (8a), A153-
A159.
(40) 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 (3), 1097-1120.
(41) Dral, P. O.; Wu, X.; Sporkel, L.; Koslowski, A.; Weber, W.; Steiger, R.; Scholten, M.; Thiel, W., Semiempirical quantum-chemical orthogonalization-corrected methods: Theory, implementation, and parameters. J. Chem. Theory Comput. 2016, 12 (3),
1082-1096.
(42) Dewar, M. J. S.; Thiel, W., Ground states of molecules. 38. The MNDO method. Approximations and parameters. J. Am. Chem. Soc. 1977, 99 (15), 4899-4907.
66


(43) Dewar, M. J. S.; Thiel, W., Ground states of molecules. 39. MNDO results for molecules containing hydrogen, carbon, nitrogen, and oxygen. J. Am. Chem. Soc. 1977, 99 (15), 4907-4917.
(44) Lowdin, P.-O., On the nonorthogonality problem. In A civ. Qunatum chem., Per-Olov, L., Ed. Academic Press: 1970; Vol. Volume 5, pp 185-199.
(45) Lowdin, P.-O., On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals. J. Chem. Phys. 1950, 18 (3),
365-375.
(46) Thiel, W., Semiempirical quantum-chemical methods. Wiley. Interdiscip. Rev. Comput. Mol. Sci. 2014, 4 (2), 145-157.
(47) Lin, H.; Truhlar, D. G., Redistributed charge and dipole schemes for combined quantum mechanical and molecular mechanical calculations. J. Phys. Chem. A 2005,109 (17), 3991-4004.
(48) Pezeshki, S.; Lin, H., Molecular dynamics simulations of ion solvation by flexible-boundary QM/MM: On-the-fly partial charge transfer between QM and MM subsystems. J. Comput. Chem. 2014, 35 (24), 1778-1788.
(49) Zhang, Y.; Lin, H., Flexible-boundary quantum-mechanical/molecular-mechanical calculations: Partial charge transfer between the quantum-mechanical and molecular-mechanical subsystems. J. Chem. Theory Comput. 2008, 4 (3), 414-425.
(50) Zhang, Y.; Lin, H., Flexible-boundary QM/MM calculations: Ii. Partial charge transfer across the QM/MM boundary that passes through a covalent bond. Theor. Chem. Acc. 2010,126(5-6), 315-322.
(51) Reuter, N.; Dejaegere, A.; Maigret, B.; Karplus, M., Frontier bonds in QM/MM methods: A comparison of different approaches. J. Phys. Chem. A 2000,104 (8), 1720-1735.
(52) Dapprich, S.; Komiromi, I.; Byun, K. S.; Morokuma, K.; Frisch, M. J., A new oniom implementation in Gaussian98. Part I. The calculation of energies, gradients, vibrational frequencies and electric field derivatives. Theochem 1999, 461-462, 1-21.
67


(53) 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(9), 1170-1179.
(54) 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 (1-3), 313-323.
(55) 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.
(56) Kerdcharoen, T.; Morokuma, K., Combined quantum mechanics and molecular mechanics simulation of Ca2+/ammonia solution based on the ONIOM-XS method: Octahedral coordination and implication to biology. J. Chem. Phys. 2003,118 (19), 8856-8862.
(57) 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.
(58) Pezeshki, S.; Davis, C.; Heyden, A.; Lin, H., Adaptive-partitioning QM/MM dynamics simulations: 3. Solvent molecules entering and leaving protein binding sites. J. Chem. Theory Comput. 2014,10 (11), 4765-4776.
(59) 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.
(60) Pezeshki, S.; Lin, H., Recent developments in QM/MM methods towards openboundary multi-scale simulations. Mol. Simul. 2015, 41 (1-3), 168-189.
(61) Bulo, R. E.; Ensing, B.; Sikkema, J.; Visscher, L., Toward a practical method for adaptive QM/MM simulations. J. Chem. Theory Comput. 2009, 5 (9), 2212-2221.
(62) 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.
68


(63) 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 (39), 12401-12414.
(64) 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.
(65) Mones, L.; Jones, A.; Gotz, A. W.; Laino, T.; Walker, R. C.; Leimkuhler, B.; Csanyi, G.; Bernstein, N., The adaptive buffered force QM/MM method in the CP2K and Amber software packages. J. Comput. Chem. 2015, 36 (9), 633-648.
(66) Peguiron, A.; Colombi Ciacchi, L.; De Vita, A.; Kermode, J. R.; Moras, G., Accuracy of buffered-force QM/MM simulations of silica. J. Chem. Phys. 2015,142 (6), 064116.
(67) Varnai, C.; Bernstein, N.; Mones, L.; Csanyi, G., Tests of an adaptive QM/MM calculation on free energy profiles of chemical reactions in solution. J. Phys. Chem. B. 2013, 7/7(40), 12202-12211.
(68) 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.
(69) Takenaka, N.; Kitamura, Y.; Koyano, Y.; Nagaoka, M., An improvement in quantum mechanical description of solute-solvent interactions in condensed systems via the number-adaptive multiscale quantum mechanical/molecular mechanical-molecular dynamics method: Application to zwitterionic glycine in aqueous solution. J Chem Phys 2012,137 (2), 024501.
(70) 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.
(71) Waller, M. P.; Kumbhar, S.; Yang, J., A density-based adaptive quantum mechanical/molecular mechanical method. ChemPhysChem 2014,15 (15), 3218-25.
(72) Bockmann, 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 (6), 2429-2439.
69


(73) 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 (11), 3625-3634.
(74) Deev, V.; Collins, M. A., Approximate ah initio energies by systematic molecular fragmentation. J. Chem. Phys. 2005,122 (15), 154102.
(75) Elrod, M. J.; Saykally, R. J., Many-body effects in intermolecular forces. Chem. Rev. 1994, 94 (7), 1975-1997.
(76) Hirata S.; Valiev, M.; Dupuis, M.; Xantheas, S. S.; Sugiki, S.; Sekino, EL, Fast electron correlation methods for molecular clusters in the ground and excited states. Mol. Phys. 2005,103 (15-16), 2255-2265.
(77) 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.
(78) 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.
(79) Agmon, N., Infrared spectroscopy: The acid test for water structure. Nat. Chem. 2016, 8 (3), 206-207.
(80) Thamer, M.; De Marco, L.; Ramasesha, K.; Mandal, A.; Tokmakoff, A., Ultrafast 2D IR spectroscopy of the excess proton in liquid water. Science 2015, 350 (6256), 78-82.
(81) Steiner, T., The hydrogen bond in the solid state. Angew. Chem. Int. Ed. Engl. 2002, 41 (1), 49-76.
(82) de Grotthuss, C. J., Memoir on the decomposition of water and of the bodies that it holds in solution by means of galvanic electricity. Biochim. Biophys. Acta. 2006,1757 (8),
871-875.
(83) Von Grotthuss, C. J. T., Sur la decomposition de l'eau et des corps qu'elle tient en dissolution a l'aide de l'electricite galvanique. Ann. Chim. 1806, 58, 54.
70


(84) Schuster, P.; Zundel, G.; Sandorfy, C., The hydrogen bond: Recent developments in theory and experiments: Structure and spectroscopy. II. North-Holland: 1976.
(85) Eigen, M.; de Maeyer, L., Self-dissociation and protonic charge transport in water and ice. Proc. Roy. Soc. London Ser. A 1958, 247 (1251), 505-533.
(86) Agmon, N., The Grotthuss mechanism. Chem. Phys. Lett. 1995, 244 (5-6), 456-
462.
(87) Markovitch, O.; Agmon, N., Structure and energetics of the hydronium hydration shells. J. Phys. Chem. A 2007, 111 (12), 2253-2256.
(88) Tuckerman, M.; Laasonen, K.; Sprik, M.; Parrinello, M., Ab initio molecular dynamics simulation of the solvation and transport of EbO+ and OH" ions in water. J. Phys. Chem. 1995, 99 (16), 5749-5752.
(89) Hassanali, A.; Giberti, F.; Cuny, J.; Kuhne, T. D.; Parrinello, M., Proton transfer through the water gossamer. Proc. Natl. Acad. Sci. U.S.A. 2013,110 (34), 13723-13728.
(90) Lapid, H.; Agmon, N.; Petersen, M. K.; Voth, G. A., A bond-order analysis of the mechanism for hydrated proton mobility in liquid water. J. Chem. Phys. 2005,122 (1), 14506.
(91) Marx, D., Proton transfer 200 years after von Grotthuss: Insights from ab initio simulations. ChemPhysChem 2006, 7 (9), 1848-1870.
(92) Chandra, A.; Tuckerman, M. E.; Marx, D., Connecting solvation shell structure to proton transport kinetics in hydrogen-bonded networks via population correlation functions. Phys. Rev. Lett. 2007, 99 (14), 145901.
(93) Marx, D.; Chandra, A.; Tuckerman, M. E., Aqueous basic solutions: Hydroxide solvation, structural diffusion, and comparison to the hydrated proton. Chem. Rev. 2010,110 (4), 2174-2216.
(94) Vuilleumier, R.; Borgis, D., Transport and spectroscopy of the hydrated proton: A molecular dynamics study. J. Chem. Phys. 1999, 111 (9), 4251-4266.
71


(95) Cao, Z.; Peng, Y.; Yan, T.; Li, S.; Li, A.; Voth, G. A., Mechanism of fast proton transport along one-dimensional water chains confined in carbon nanotubes. J. Am. Chem. Soc. 2010,132 (33), 11395-11397.
(96) Brancato, G.; Tuckerman, M. E., A polarizable multistate empirical valence bond model for proton transport in aqueous solution. J. Chem. Phys. 2005,122 (22), 224507.
(97) Cleland, W. W.; Kreevoy, M. M., Low-barrier hydrogen bonds and enzymic catalysis. Science 1994, 264 (5167), 1887-90.
(98) Tuckerman, M. E.; Marx, D.; Klein, M. L.; Parrinello, M., On the quantum nature of the shared proton in hydrogen bonds. Science 1997, 275 (5301), 817-820.
(99) Ceriotti, M.; Cuny, J.; Parrinello, M.; Manolopoulos, D. E., Nuclear quantum effects and hydrogen bond fluctuations in water. Proc. Natl. Acad. Sci. U.S.A. 2013,110 (39), 15591-15596.
(100) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M., Ab initio simulations of water and water ions. J. Phys. Condens. Matter 1994, 6 (23A), A93-A100.
(101) Warshel, A.; Weiss, R. M., An empirical valence bond approach for comparing reactions in solutions and in enzymes. J. Am. Chem. Soc. 1980,102 (20), 6218-6226.
(102) Sagnella, D. E.; Tuckerman, M. E., An empirical valence bond model for proton transfer in water. J. Chem. Phys. 1998,108 (5), 2073-2083.
(103) Schmitt, U. W.; Voth, G. A., Multistate empirical valence bond model for proton transport in water. J. Phys. Chem. B 1998,102 (29), 5547-5551.
(104) Schmitt, U. W.; Voth, G. A., The computer simulation of proton transport in water. J. Chem. Phys. 1999, 111 (20), 9361-9381.
(105) Vuilleumier, R.; Borgis, D., An extended empirical valence bond model for describing proton transfer in H+(H20)n clusters and liquid water. Chem. Phys. Lett. 1998, 284 (1-2), 71-77.
72


(106) Lefohn, A. E.; Ovchinnikov, M.; Voth, G. A., A multistate empirical valence bond approach to a polarizable and flexible water model. J. Phys. Chem. B 2001,105 (28), 6628-6637.
(107) Komyshev, A. A.; Kuznetsov, A. M.; Spohr, E.; Ulstrup, J., Kinetics of proton transport in water. J. Phys. Chem. B. 2003,107 (15), 3351-3366.
(108) Swanson, J. M.; Maupin, C. M.; Chen, EL; 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 (17), 4300-4314.
(109) Markovitch, O.; Chen, H.; Izvekov, S.; Paesani, F.; Voth, G. A.; Agmon, N., Special pair dance and partner selection: Elementary steps in proton transport in liquid water. J. Phys. Chem. B 2008,112 (31), 9456-9466.
(110) Wu, Y.; Chen, H.; 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 (2), 467-482.
(111) Park, K.; Lin, W.; Paesani, F., Fast and slow proton transfer in ice: The role of the quasi-liquid layer and hydrogen-bond network. J. Phys. Chem. B. 2014,118 (28), 8081-8089.
(112) Cao, Z.; Kumar, R.; Peng, Y.; Voth, G. A., Proton transport under external applied voltage. J. Phys. Chem. B 2014,118 (28), 8090-8098.
(113) Lee, S.; Swanson, J. M.; Voth, G. A., Multiscale simulations reveal key aspects of the proton transport mechanism in the CLC-ecl antiporter. Biophys. J. 2016,110 (6),
1334-1345.
(114) 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 (45), E6955-E6964.
(115) Billeter, S. R.; van Gunsteren, W. F., Protonizable water model for quantum dynamical simulations. J. Phys. Chem. A 1998,102 (24), 4669-4678.
(116) David, C. W., A variable charge central force model for water and its ionic dissociation products. J. Chem. Phys. 1996,104 (18), 7255-7260.
73


(117) Fogarty, J. C.; Aktulga, H. M.; Grama, A. Y.; van Duin, A. C.; Pandit, S. A., A reactive molecular dynamics simulation of the silica-water interface. J. Chem. Phys. 2010, 132 (17), 174704.
(118) Halley, J. W.; Rustad, J. R.; Rahman, A., A polarizable, dissociating, molecular-dynamics model for liquid water. J. Chem. Phys. 1993, 98 (5), 4110-4119.
(119) Kale, S.; Herzfeld, J.; Dai, S.; Blank, M., Lewis-inspired representation of dissociable water in clusters and grotthuss chains. J. Biol. Phys. 2012, 38 (1), 49-59.
(120) Lussetti, E.; Pastore, G.; Smargiassi, E., A fully polarizable and dissociable potential for water. Chem. Phys. Lett. 2003, 381 (3-4), 287-291.
(121) Mahadevan, T. S.; Garofalini, S. H., Dissociative water potential for molecular dynamics simulations. J. Phys. Chem. B 2007, 111 (30), 8919-8927.
(122) 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 (27), 11965-11976.
(123) Stillinger, F. H.; David, C. W., Polarization model for water and its ionic dissociation products. J. Chem. Phys. 1978, 69 (4), 1473-1484.
(124) Wolf, M. G.; Groenhof, G., Explicit proton transfer in classical molecular dynamics simulations. J. Comput. Chem. 2014, 35 (8), 657-671.
(125) Arillo-Flores, O. I.; Ruiz-Lopez, M. F.; Bernal-Uruchurtu, M. I., Can semi-empirical models describe HC1 dissociation in water? Theoretical Chemistry Accounts 2007, 118 (2), 425-435.
(126) Choi, T. H.; Jordan, K. D., Application of the SCC-DFTB method to H+(H20)6, H+(H20)2i, and H+(H20)22. J. Phys. Chem. B. 2010,114 (20), 6932-6936.
(127) 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.
74


(128) Goyal, P.; Qian, H. J.; Me, S.; Lu, X.; 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.
(129) Korth, M., Third-generation hydrogen-bonding corrections for semiempirical qm methods and force fields. J. Chem. Theory Comput. 2010, 6 (12), 3808-3816.
(130) Lin, Y.; Wynveen, A.; Halley, J. W.; Curtiss, L. A.; Redfern, P. C., Self consistent tight binding model for dissociable water. J. Chem. Phys. 2012,136 (17), 174507.
(131) Maupin, C. M.; Aradi, B.; Voth, G. A., The self-consistent charge density functional tight binding method applied to liquid water and the hydrated excess proton: Benchmark simulations. J. Phys. Chem. B. 2010,114 (20), 6922-6931.
(132) Riccardi, D.; Konig, P.; Prat-Resina, X.; Yu, H.; 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 (50), 16302-16311.
(133) Wang, S.; MacKay, L.; Lamoureux, G., Development of semiempirical models for proton transfer reactions in water. J. Chem. Theory. Comput. 2014,10 (8), 2881-2890.
(134) 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 (6), 2672-2686.
(135) Hofer, T. S.; Hitzenberger, M.; Randolf, B. R., Combining a dissociative water model with a hybrid QM/MM approacha simulation strategy for the study of proton transfer reactions in solution. J. Chem. Theory Comput. 2012, 8 (10), 3586-3595.
(136) Lin, H.; Zhang, Y.; Pezeshki, S.; Duster, A.; Wang, C. H.; Truhlar, D. G. Qmmm Version 2.2.7 CO, University of Minnesota, Minneapolis, 2017.
(137) 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(16), 1781-1802.
75


(138) Ponder, J. W.; Richards, F. M., An efficient Newton-like method for molecular mechanics energy minimization of large molecules. J. Comput. Chem. 1987, 8 (7), 1016-1024.
(139) Trucks, G.; Schlegel, H.; Scuseria, G.; Robb, M.; Cheeseman, J.; Montgomery Jr, J.; Vreven, T.; Kudin, K.; Burant, J.; Millam, J. Gaussian-03, RevisionE. 01, Gaussian Inc.: Wallingford, CT, 2004.
(140) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, FL; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery Jr., J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Gaussian, Inc.: Wallingford, CT, USA, 2009.
(141) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A., General atomic and molecular electronic structure system. J. Comput. Chem. 1993,14(11), 1347-1363.
(142) Neese, F., The ORCA program system. Wiley. Interdiscip. Rev. Comput. Mol. Sci. 2012, 2 (1), 73-78.
(143) Thiel, W. MND099, Max-Planck-Institut fur Kohlenforschung: Mulheim an der Ruhr, Germany, 2012.
(144) Humphrey, W.; Dalke, A.; Schulten, K., Vmd: Visual molecular dynamics. J. Mol. Graph. 1996,14 (1), 33-8, 27-8.
(145) Berman, H.; Henrick, K.; Nakamura, H., Announcing the worldwide protein data bank. Nat. Struct. Mol. Biol. 2003,10 (12), 980-980.
(146) Dagum, L.; Menon, R., Openmp: An industry standard api for shared-memory programming. IEEE Comput. Sci. Eng. 1998, 5 (1), 46-55.
76


(147) OpenMP, A., OpenMP application program interface, v. 3.0. OpenMP Architecture Review Board 2008.
(148) Chang, C. H.; Long, H.; Sides, S.; Vaidhynathan, D.; Jones, W., Parallel application performance on two generations of Intel Xeon HPC platforms. Laboratory, N. R. E., Ed. Golden, CO, 2015.
(149) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. EL; Kollman, P. A., The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. J. Comput. Chem. 1992,13 (8), 1011-1021.
(150) Wu, Y.; Tepper, H. L.; Voth, G. A., Flexible simple point-charge water model with improved liquid-state properties. J. Chem. Phys. 2006,124 (2), 024503.
(151) Zhu, F. Q.; Schulten, K., Water and proton conduction through carbon nanotubes as models for biological channels. Biophys. J. 2003, 85 (1), 236-244.
(152) Han, W.; Cheng, R. C.; Maduke, M. C.; Tajkhorshid, E., Water access points and hydration pathways in CLC H+/CL transporters. Proc. Natl. Acad. Sci. U.S.A. 2014, 111 (5), 1819-1824.
(153) Dellago, C.; Hummer, G., Kinetics and mechanism of proton transport across membrane nanopores. Phys. Rev. Lett. 2006, 97 (24), 245901.
(154) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; 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. E.; 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 (18), 3586-3616.
(155) Knight, C.; Lindberg, G. E.; Voth, G. A., Multiscale reactive molecular dynamics. J. Chem. Phys. 2012,137 (22), 22A525.
(156) Ponder, J. W. Tinker molecular modeling package Jay Ponder Lab. https://dasher.wustl.edu/tinker/ (accessed Feb. 3 2017).
77


APPENDIX A
BENCHMARKING INFORMATION
A.l Test System Specifications
Table A.1: System Specifications for Lenovo Laptop
Test system ID: laptop
CPU Intel Core i7-6700HQ 2.60 Ghz
Cores Physical/Virtual 4/8
GPU NVIDIA GeForce GTX 960M
RAM 16 GB
Disk 1 SanDisk Ultra II (X41200RL) SSD
Disk 2 WD My Passport Ultra USB 3.0 HDD
OS Ubuntu 16.10
Computer Type Laptop
Table A.2: System Specifications for compchemw5
Test System ID compchemw5
CPU Intel Xeon CPU E5-2687W @ 3.10 GHz
Cores Physical/Virtual 16/32
RAM 64 GB
GPU NVIDIA GTX 1080
Disk 1 Crucial BX200 SSD
Disk 2 Seagate Constellation ES.2 7200 RPM HDD
OS Centos 7
Computer Type Desktop
Table A.3: System Specifications for compchemw8
Test System ID compchemw8
CPU Intel Xeon CPU E5-2660 v2 a 2.20 GHz
Cores Physical/Virtual 10/20
RAM 64 GB
GPU NVIDIA GTX 980
GPU NVIDIA Tesla K40
Disk 1 Intel SSDSC2KI010X6 SSD
Disk 2 WDC WD2000FYYZ-01UL1B2 7200 RPM HDD
OS Centos 7
Computer Type Desktop


Table A.4: System Specifications for XSEDE Stampede
Test System ID Stampede
CPU 2x Intel XeonE5 2680 2.7 GHz
Cores Physical/Virtual 16/32
RAM 32 GB
GPU 61 Core XeonPhi SE10P 1.1GHz
GPU NVIDIA Tesla K40
Disk 1 250GB HDD (7.5 K RPM SATA)
Disk 2 LUSTRE Network Filey stem
OS Centos 7
Computer Type HPC
A.2 Tinker/NAMD Benchmark Calculation
Test system 1 described in Table A.l, was used for the following benchmark. Although the provided Tinker executables were claimed to be OpenMP capable,156 and the proper keyword was used in the .key file, we found that the testgrad binary rarely used more than 1 core in our calculations. Therefore, NAMD was ran with CUDA using 1 thread., and with CUDA utilizing 4 threads. The runs with NAMD using 1 thread were reported in the main body of the thesis. The additional timing data is highlighted in yellow in Table A.5:
Table A.5: Benchmarking Results for NAMD and Tinker.
Number of Atoms 4308 52233 96003 176454
NAMD time/step (s) 0.18 0.55 0.91 1.57
NAMD 4 cores time/step (s) 0.16 0.57 1.02 1.64
Tinker time/step (s) 0.9 35.2 111.8 330.9
1 Core: Speedup with NAMD 5 64 123 211
4 Cores: Speedup with NAMD 6 62 112 195
A.3 PPAP Benchmarking of a Water Box for Various Systems Desktop Workstation 8
These AP MD calculations were carried using the test system described in Table A.3 on a water box of 1435 H2O and 1 H3CU cation. The active zone has a 5.0 A radius and is surrounded by a 1.0 A buffer shell. On average, there were 13.4 active zone groups, 17.3
79


buffer zone groups, and 160.7 QM calculations per step for 100 steps. The AP potential was truncated at the 2nd order. The timing data can be found in Table A.6.
Table A.6: Tests with Desktop Workstation compchemw8 with an HDD and SSD. WT/S stands for wall time per step and SSE stands for strong scaling efficiency.
Threads SSD: WT/S (s) SSI): SSE HDD: WT/S (s) HDD: SSE
1 62.8 100.00% 65.8 100.00%
2 29.1 107.90% 28.8 114.24%
4 15.2 103.29% 16.0 102.81%
8 8.4 93.45% 9.5 86.58%
10 7.1 88.45% 8.5 77.41%
|20 5.4 58.15% 6.0 54.83%
XSEDE Stampede
The following tests used the normal compute nodes on the Stampede 1 system, described in Table A.3. The test model and parameters were the same as in the model in Table A.5, except N AMD 2.10 was used instead of2.11. There as an average of 13.5 active zone groups, 17.5 buffer zone groups, and 163.5 QM calculations per step. The timing data can be found in Table A.7.
Table A.7: Tests using XSEDE Stampede with an internal SSD and LUSTRE network filesystem. The calculation on the network fde system failed when 32 threads were used. WT/S stands for wall time per step and SSE stands for strong scaling efficiency
Threads Internal: WT/S (s) Internal: SSE Nehvork: WT/S (s) Nehvork: SSE
1 74.8 100.00% 100.3 100.00%
2 36.9 101.36% 60.1 83.44%
4 22.6 82.74% 38.8 64.63%
8 9.8 95.41% 21.1 59.42%
10 8.9 84.04% 19.5 51.44%
16 8.9 52.53% 21.7 28.89%
32 11.2 20.87% N/A N/A
Desktop Workstation 5
The following tests used the test system described in Table A.3 on a water box of 1435 H2O and 1 H3CE cation. Two active/buffer zone radii were investigated in order to
80


show the decrease in SSE due to a lack of available calculations to assign to the thread team. With an active zone radius of 4.0 A surrounded by a 0.5 A buffer shell, there was an average of 5.9 active zone groups, 3.4 buffer zone groups, and 9.3 QM calculations per step for 100 steps. With an active zone radius of 4.5 A surrounded by a 1.0 A buffer shell, there was an average of 7.5 active zone groups, 15.6 buffer zone groups, and 115.0 QM calculations per step for 100 steps. The timing data can be found in Table A.8.
Table A.8: Strong scaling efficiency (SSE) data for two active/buffer zone radii configurations: one with an active zone radius of 4.0 A surrounded by a 0.5 A buffer shell had 9.3 QM/MM calculations per step, the other with an active zone radius of 4.5 A surrounded by a 1.0 A buffer shell had 115.0 QM/MM calculations per step.
Threads 9.3 Calculations/Step SSE 115.0 Calculations /Step SSE
1 100.00% 100.00%
2 107.50% 99.49%
4 82.69% 93.81%
8 53.75% 84.91%
16 29.86% 70.36%
81


APPENDIX B
DECOMPOSED GRADIENT INCLUDING HYDROGEN
B.l: Variables and Functions used in the Gradient Decomposition Equations
Table B.l: Variables and Functions used in the Gradient Decomposition Equations Symbol Description
Restraint liannonic potential energy imposed on indicator I
k Force constant of the restraint potential
Xi Current position of indicator
Xi0 Targeted position of the indicator
V, Position of /-th possible acceptor oxygen A,
XD Position of donor oxygen D
XHm Position of w-th hydrogen bonded to donor D
rDHm Donor-hydrogen vector of w-th hydrogen
rDAy Donor-acceptor vector of y'-th possible acceptor
* DH Equilibrium bond distance parameter
rLIST Possible acceptor enlisting cutoff parameter
Pmj Ratio of vector projection for w-th H and /-th A
Pm} Cutoff parameter for pmj
Weight function for pmj
fli Nonnalization factor for the weight function
B.2: Decomposition of the Potential onto XA
Gradient of Fhar with respect to XA:
dVhar
~dXl~
fc|x,-xl0
axj
dxTj
Gradient of X( with respect to XA:
dXl
dxTj

M
'I
m= 1
/ da(Pmj)
] M
II
j=l m= 1
M
I
m= 1
d 9(Pmj)
dXAj
\
/
(5.2.1)
(5.2.2)
82


Gradient of g(pmj) with respect to XA :
dxAj
r0
(30x4 + 60x3 30x2)
dx
dxT
.0
if 1 < x t/ 0 < x < 1 if x < 0
(5.2.3)
Gradient of x with respect to XA:
dx 1 d pmj
dxZj = 0.5 p^j dXA.
(5.2.4)
Gradient of pmj with respect to XA
d Pmj
dXl~
rDH;
rDA,
2
(5.2.5)
B.3 Decomposition of the Potential onto XD
Gradient of V^ar with respect to XD:
d fh fly
= /c|Xj X
dX
D
dX!
IoldX^
(5.3.1)
Gradient of X( with respect to XD:
/
dX
dX
I -2
= 9\
J M
s.i i+ZZiXa
7 = 1 m= 1
J M
dg(Pmj)\\ \
1 dxD y
J M
V
XD + X I X I
/=! m=l
d^CPmy') ^ dXD /
7 = 1 m=l /
(5.3.2)
83


Gradient of g(pmj) with respect to XD:
dg{^Pmi)
dXD
= <(30x4 + 60x3 30x2)
dx
dX^
if 1 < x if 0 < x < 1 if x < 0
(5.3.3)
Gradient of x with respect to XD:
dx 1 d pmj
dX^ = 0.5 -pmj dXD
Gradient of pmj with respect to XD:
dPmj
dX
D
rDHm Ay
+ 2 (rDH]
(5. 3.4)
(5.3.5)
B.4: Decomposition of the Potential onto XHm:
Gradient of V^ar with respect to XHm:
dfftar / \ dXj
h ^(x, x,0) -
dX
dX

(5.4.1)
Gradient of X\ with respect to XHm:
] M
II
j=l m= 1
(5.4.2)
84


Gradient of g(pmj) f with respect to XHm:

dX
H
rO
= i(30x4 + 60x3 30x2)
^0
dx
dxiT
nM
if 1 < x if 0 < x < 1 if x < 0
Gradient of x with respect to XH :
dx
-1 dp.
my
dXH 0.5-p,dXH
nm r 771J nm
Gradients of pmj with respect to XHm:
dp
mj
dXi
rDA;
rDA;
(5.4.3)
(5.4.4)
(5.4.5)
B.5: Gradient Graphs
When the gradient was decomposed over all atoms, the gradient of the donor and acceptor atom was in the opposite direction of the gradient on the proton These gradients could be quite large, but cancelled each other out when summed to yield the total gradient on the indicator. The results of these exploratory calculations are shown in Fig B.l.
85


1,5 2.0 2.5 3.0 3.5
ro-o
o
£
o
j*'
c
0
03
0
70
60
50
40
30
20
10
0
-10
1.5
Donor and Acceptor Gradients
c A V V Acceptor ^B Donor

2.0
2.5
3.0
3.5
Fig B.l: Decomposed gradients for Zundel cation. (A) Indicator distance from origin and (B) gradient on indicator. Gradients on (B) donor and acceptor, and (D) indicator. The maximum gradient on the proton is around -130 kcal mol"1, and the maximum gradient on the acceptor and donor is around 55 and 65 kcal mol" 1 respectively.
86


APPENDIX C
Additional Details for MD Simulations
C.l NVE Collision Simulations
Only one half of the graphs were reported in the main text. The full sets of graphs are found below:
0.0 0.5 1.0 1.5 2.0
Time (ps)
Time (ps)
Time (ps)
0.0 0.5 1.0 1.5 2.0
Time (ps)
k = 0 k = 5 k = 10 k = 20
Fig C.l: Total energy of system for all simulations of the Zundel (A, C) and Eigen (B, D) cations with different force constants and initial EEO velocities of (A, B) 0.001 A fs_1 and (C, D) 0.01 A fs1. The blue, green, red, and cyan lines correspond to force constants of 0 (no restraint), 5, 10, and 20 kcal mol"1 A-2 respectively.
87


Distance [A] Distance [A]
Fig C.2: RMSD (A) of indicator distance from origin of potential for Zundel (rows 1,3) and Eigen (rows 2, 4) cations. Red, blue, green, and cyan correspond to force constants of 0 (no restraint), 5, 10, and 20 kcal mol"1 A-2 respectively. The initial velocity of the EEO was 0.001 A fs_1(rows 1,2) and 0.01 A fs'1 (rows 3,4).
88


C.l MD simulations of Water Sphere
Four simulations of a water sphere containing 1482 H2O molecules with one FbO+ at the center of the sphere were conducted. The force constant for indicator restraints were set at 5 kcal mol A"2. Four additional simulations were conducted in the absence of indicator restraints for comparison. The simulations used a droplet model with the inner 12 A sphere described at the OM3n level and the outer 10 A shell described at the MM level with the SPC/Fw potential. For the MM molecules, restraints of 20 kcal mol"1 applied to every oxygen atom. Other general simulation details can be found in Table C.l.
Table C.l General MD parameters for water sphere simulations
Atoms 4450
QM only atoms 796
Step length 0.25 fs
Equilibration time 0.5 ps
Production time 10.0 ps
Temperature 300 K
ES Charge Cutoff 14.0 A
ES Taper 13.0 A
Snapshot Save Frequency 2 fs
Density 1.0 g ml "1
Periodic Boundaries None
Average wall time per step ~10.2s
Simulation System compchemw5
MM Code NAMD2.10
QM Code MNDO
C.3 MM System used for Nanotube Simulations
The nanotube system was constructed and equilibrated via the following procedure. The epsilon parameter for the nanotube carbon, C A was changed to 80% of its original value to -0.056 to increase hydrophobicity of the tube. This was an error which would be later corrected, as the depth of the potential well between the oxygen water atom and carbon atom should have been scaled to 80% of its original value, not the epsilon parameter itself.
89


The simulation was equilibrated for 200 ps with 1 fs time steps at the NPTlevel. The x and >' dimensions were fixed, and unit cell and position changes due to pressure were only allowed to change in the z direction. Four snapshots were taken at 200, 190, 180, and 170 ps from the end of the equilibration trajectory and used as starting geometries for the next phase of the equilibration. Each of the four geometries were then equilibrated for an additional 5 ns at NPT with the x and >' dimensions fixed. Finally, the simulations were equilibrated for another 5 ns at NVT.
The well depth error was realized, and the epsilon parameter for CA was then changed to -0.0448. The final snapshot of each of the previous simulations was then equilibrated for an additional 5 ns at NPT with the x and >' dimensions fixed. The simulations were then equilibrated for another 5 ns at NVT.
C.4 MD Simulations of Restrained Indicator at Pore Opening
From each MM trajectory, a snapshot with zero, one, two, three, or four waters in the tube was obtained starting from the end of the trajectory, yielding a total of 20 starting geometries. All waters within a 14 A radius of the middle of the pore opening were included in the QM subsystem. The number QM atoms for these simulations ranges from 400 to 475, with an average of 437 QM atoms. Simulation lengths are listed in Table C.2. Also, Fig 5.3 in the main text is taken from the simulation named res Ow 2.
90


Table C.2: List of Simulations of Restrained Indicator at Pore Opening
Name Restraint Initial Waters Total Time (ps)
res Ow 1 Yes 0 2.66675
res Ow 2 Yes 0 10.50000
res Ow 3 Yes 0 10.50000
res Ow 4 Yes 0 3.5615
res lw 1 Yes 1 2.39225
res lw 2 Yes 1 4.2910
res lw 3 Yes 1 6.1550
res lw 4 Yes 1 10.50000
res 2w 1 Yes 2 10.50000
res 2w 2 Yes 2 4.34500
res 2w 3 Yes 2 10.50000
res 2w 4 Yes 2 6.56350
res 3w 1 Yes 3 10.50000
res 3w 2 Yes 3 10.50000
res 3w 3 Yes 3 10.50000
res 3w 4 Yes 3 8.22300
res 4w 1 Yes 4 10.50000
res 4w 2 Yes 4 10.50000
res 4w 3 Yes 4 10.50000
res 4w 4 Yes 4 1.96275
nor Ow 1 No 0 10.50000
nor Ow 2 No 0 10.50000
nor Ow 3 No 0 0.00000
nor Ow 4 No 0 8.94650
nor lw 1 No 1 10.50000
nor lw 2 No 1 10.50000
nor lw 3 No 1 10.50000
nor lw 4 No 1 10.50000
nor 2w 1 No 2 10.50000
nor 2w 2 No 2 10.50000
nor 2w 3 No 2 10.50000
nor 2w 4 No 2 10.50000
nor 3w 1 No 3 10.50000
nor 3w 2 No 3 7.79575
nor 3w 3 No 3 8.54700
nor 3w 4 No 3 10.50000
nor 4w 1 No 4 10.50000
nor 4w 2 No 4 10.50000
nor 4w 3 No 4 10.50000
nor 4w 4 No 4 3.64650
91


Full Text
xml version 1.0 encoding ISO-8859-1
DISS_submission publishing_option 0 embargo_code third_party_search Y
DISS_authorship
DISS_author type primary
DISS_name
DISS_surname Duster
DISS_fname Adam
DISS_middle William
DISS_suffix
DISS_affiliation University of Colorado at Denver
DISS_contact current
DISS_contact_effdt 04/23/2017
DISS_phone_fax P
DISS_cntry_cd 1
DISS_area_code 832
DISS_phone_num 248-6134
DISS_phone_ext
DISS_address
DISS_addrline 475 White Ash Dr
DISS_city Golden
DISS_st CO
DISS_pcode 80403
DISS_country US
DISS_email adam.duster@ucdenver.edu
future
04/23/2017
1
832
248-6134
475 White Ash Dr
Golden
CO
80403
US
adam.duster@ucdenver.edu
DISS_citizenship US
DISS_description page_count 104 masters external_id http:dissertations.umi.comucdenver:10847 apply_for_copyright no
DISS_title Tools and methodology for the adaptive multiscale simulation of proton transport in complex environments
DISS_dates
DISS_comp_date 2017
DISS_accept_date 01/01/2017
DISS_degree M.S.
DISS_institution
DISS_inst_code 0765
DISS_inst_name University of Colorado Denver
DISS_inst_contact Chemistry
DISS_processing_code N
DISS_advisor
Lin
Hai
DISS_cmte_member
Lin
Hai
Crowley
Michael
Wang
Haobin
DISS_categorization
DISS_category
DISS_cat_code 0485
DISS_cat_desc Chemistry
0786
Biophysics
0487
Biochemistry
DISS_keyword Adaptive Partitioning, Multiscale, PAP, Permuted Adaptive Partitioning, Proton transfer, QM/MM
DISS_language en
DISS_content
DISS_abstract
DISS_para This thesis reports three developments in the combined quantum-mechanics/molecular-mechanics (QM/MM) methodology and software. First, an interface between the program QMMM and the MM package NAMD was implemented, substantially increasing the computational efficiency of MM calculations. Second, the permuted adaptive partitioning QM/MM algorithm implemented in the QMMM software package was parallelized using the OpenMP protocol. The algorithm can scale approximately linearly when the number of buffer groups vastly exceeds the number of threads requested. These two improvements will increase the scope and scale of future calculations conducted with QMMM. Finally, the feasibility of using a restrained proton indicator in biased (MD) simulations was investigated. A harmonic restraint potential was imposed on an existing proton indicator, the position of which approximates the location of an excess proton in water. When applied to a carbon nanotube model system, the restrained proton accelerated water recruitment into the nanotube when it was restrained above the nanotube opening, in agreement with previous simulations in the literature, where a different algorithm for proton representation was used.
DISS_supp_abstract
DISS_binary PDF Duster_ucdenver_0765N_10847.pdf
DISS_restriction
DISS_repository
DISS_version 2011-11-08 15:37:33
DISS_agreement_decision_date 2017-04-23 00:59:34
DISS_acceptance 1
DISS_delayed_release
DISS_access_option



PAGE 1

TOOLS AND METHODOLOG Y FOR THE ADAPTIVE M ULTISCALE SIMULATION OF PROTON TRANSPORT IN COMPLEX ENVIRONME NTS by ADAM WILLIAM DUSTER B.S., Metropolitan State University of Denver 201 4 A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Chemistry Program 2017

PAGE 2

ii This thesis for the Mas ter of Science degree by Adam William Duster has been approved for the Chemistry Program by Ha i Lin, Chair Haobin Wang Michael Crowley Date: May 13, 2017

PAGE 3

iii Duster Adam William (M.S. Chemistry) Tools and Methodology for the Adaptive Multiscale Simulation of Proton Transport in Complex Environments Thesis directed by Associate Professor Hai Lin ABSTRACT This thesis reports three developments in the combined quantum mechanics/molecular mechanics (QM/MM) methodology and software. First, an interface between the program QMMM and the MM package NAMD was implemented, substantially increasing the compu tational efficiency of MM calculations. Second, t he permuted adaptive partitioning QM/MM algorithm implemented in the QMMM software package was parallelized using the O penMP protocol. The algorithm can scale approximately linearly under appropriate conditi ons These two improvements will increase the scope and scale of future calculations conducted with QMMM. Finally, the feasibility of using a restrained proton indica tor in biased (MD) simulations was investigated. A harmonic restraint potential was impose d on an existing proton indicato r, the position of which approximates the location of an excess proton in water. When applied to a carbon nanotube model system, the restrained proton accelerated water recruitment into the nanotube when it was restrained ab ove the nanotube opening in agreement with previous simulations in the literature, where a different algorithm for proton representation was used. The form and content of this abstract are approved. I recommend its publication Approved: Hai Lin

PAGE 4

iv To my pa rents and sister for their unconditional love and support. To Nduka Odizor, for his mentorship in forming my values, character, and morals. To David and Lauren Suarez Mary Rose Siason Joshua Skotak, and Zachary Lowery for your friendship and support over many years To Michael Jacobs and Andrew Bonham for making my Hai Lin and all members of the Lin Lab, for providing the opportunity to perform research and for making the past 3 years the happiest and most fulfilling of my l ife.

PAGE 5

ACKNOWLEDGEMENTS I would like to thank my advisor, Dr. Hai Lin for his guidance and help in preparing this thesis. I would like to thank Christina Garza for her help in extensively testing QMMM and Nara Chon for her helpful discussions. I would als o like to thank Professor Walter Thiel, Dr. Xin Wu, and Dr. Pavlo Dral for their work on the OM x family of semiempirical methods, and their willingness to share their work and teach me about them. This work is supported by National Science Foundation (CHE 1564349), Camille & Henry Dreyfus Foundation (TH 14 028), and NVIDIA Corporation This work used the Extreme Science and Engineering Discovery Environment (XSEDE) under grant CHE 140070, which is supported by National Science Foundation grant number ACI 10 53575, and the National Energy Research Scientific Computing Centre (NERSC) under grant m2495, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE AC02 05CH11231 This research is also supported by the University of Colorado Denver Chemistry Departmental Student Research Fellowship

PAGE 6

ii TABLE OF CONTENT S CHAPTER I. INTRODUCTION ................................ ................................ ................................ ........ 1 1.1 QM/MM Multiscale Modeling ................................ ................................ ............... 1 1.1.1 What is Multiscale Molecular Modelling? ................................ .................... 1 1.1.2 What is QM/MM? ................................ ................................ .......................... 1 1.2 Conventional QM/MM Methodology ................................ ................................ ..... 3 1.2.1 The MM Energy ................................ ................................ ............................. 3 1.2.2 The QM Energy ................................ ................................ ............................. 5 1.2.3 The OM3n Method ................................ ................................ ........................ 9 1.2.4 Partitioning a QM/MM System ................................ ................................ ... 10 1.2.5 The QM/M M Hamiltonian ................................ ................................ ........... 11 1.2.4 QM/MM Partitions through Covalent Bonds ................................ .............. 12 1.3 Adaptive Partitioning QM/MM Methodology ................................ ...................... 13 1.3.1 Partitioning of an AP QM/MM System ................................ ....................... 14 1.3.2 The PAP Potential ................................ ................................ ........................ 15 1.3.3 The PAP Gradients and the mPAP Scheme ................................ ................. 17 1.4 Proton Transfer in Complex Environments ................................ .......................... 18 1.4.1 Grotthuss Shuttle Mechanism of Proton Trans fer ................................ ....... 18 1.4.3 The Proton Indicator ................................ ................................ .................... 21 II. IMPLEMENTATION OF NA MD QMMM INTERFACE ................................ ........ 23 2.1 QMMM An Interfacing Program ................................ ................................ ....... 23 2.2 Motivation for Adding the NAMD Interface ................................ ........................ 23

PAGE 7

iii 2.3 Structure of the NAMD Interfac e ................................ ................................ ......... 24 2.4 Topology Changes and the PSF ................................ ................................ ............ 26 2.4.1 The Bond List ................................ ................................ .............................. 27 2.4 .2 The Angle List ................................ ................................ ............................. 27 2.4.3 The Dihedrals List ................................ ................................ ........................ 28 2.4.4 The Impropers and Cross Terms ................................ ................................ .. 28 2.5 Benchmarking ................................ ................................ ................................ ....... 29 III. IMPLEMENTATION OF PA RALLEL ADAPTIVE PART ITIONING ................... 31 3.1 Theoretical Background ................................ ................................ ........................ 31 3.2 Implementation of Parallel Permuted Adaptive Partitioning ................................ 31 3.3 Scaling and Benchmarks ................................ ................................ ....................... 33 IV. IMPLEMENTATION OF IN DICATOR RESTRAINT ................................ ............. 38 4.1 Restraint Potential and Decomposed Gradients ................................ .................... 38 4.2 Behavior of the Indicator for the Zundel Cation ................................ ................... 40 4.3 Behavior of the Indicator for Collision Simulations ................................ ............. 43 4.4 MD Simulations of a Water Sp here ................................ ................................ ...... 45 V. APPLICATION OF INDIC ATOR RESTRAINT TO CA RBON NANOTUBE ....... 49 5.1 Introduction ................................ ................................ ................................ ........... 49 5.2 Simulations of Indicator Restrained at Nanotube Opening ................................ .. 51 5.3 Simulations of Indicator Restrained at Different Distances from Nanotube ........ 5 5 5.4 QM/MM Umbrella Sampling of PMF for Proton Transfer in Nanotube .............. 57 VI. CONCLUSION ................................ ................................ ................................ ........... 61 REFERENCES ................................ ................................ ................................ ....................... 63

PAGE 8

iv APPENDIX A: BENCHMARKING INFORMATION ................................ ......................... 78 A.1 Test System Specifications ................................ ................................ .................. 78 A.2 Tinker/NAMD Benchmark Calculation ................................ ............................... 79 A.3 PPAP Benchmarking of a Water Box for Various Systems ................................ 79 Desktop Workstation 8 ................................ ................................ ......................... 79 XSEDE Stampede ................................ ................................ ................................ 80 Desktop Workstation 5 ................................ ................................ ......................... 80 APPENDEX B: DECOMPOSED GRADIENT INCLUDING HYDROGEN ....................... 82 B.1: Variables and Functions used in the Gradient Decomposition Equations .......... 82 B.2: Decomposition of the Potential onto ................................ ........................... 82 B.3 Decomposition of the Potential onto ................................ ............................. 83 B.4: Decomposition of the Potential onto : ................................ ....................... 84 B.5: Gradient Graphs ................................ ................................ ................................ .. 85 APPENDIX C: Additional Details for MD Simulations ................................ ........................ 87 C.1 NVE Collision Simulations ................................ ................................ ................. 87 C.2 MD simulations of Water Sphere ................................ ................................ ......... 89 C.3 MM System used for Nanotube Simulations ................................ ....................... 89 C.4 MD Simulations of Restrained Indicator at Pore Opening ................................ .. 90 C.5 MD Simulations of Restraints at Different Locations ................................ .......... 92 C.6 Umbrella Sampling of Proton Transfer in CNT ................................ ................... 93

PAGE 9

v LIST OF ABBREVIATION S AP A daptive P artitioning CPS Capped Primary Subsystem CSH C Shell DBA Density Based Adaptive Partitioning QM/MM DFT Density Functional Theory ES Enti re System EVB Empirical Valence Bond (Theory) HDD Hard Disk Drive HF Hartree Fock HPC High Performance Comput( er/ ing) MD M olecular D ynamics MM M olecular M echanic ( s/ al) MNDO Modified Negle ct of Differential Overlap MP2 2 nd Order Mller Plesset Perturbation T heory mPAP Modified Permuted Adaptive Partitioning MS EVB Multistate Empirical Valence Bond (Theory) NA N adaptive NAMD NAMD software package NDDO Neglect of Differential Diatomic Overlap OMx Orthogonalization corrected M ethods PAP Permuted Adaptive Partitioning PDB Protein Data Bank File PMF Potential of Mean Force post HF Post Hartree Fock PPAP Parallel Permuted Adaptive Partitioning PS Primary Subsystem PSF Protein Struct ure File QM Q uantum M echanic( s/al) QMMM QMMM Software Package QM/MM Combined Quantum M echanics/ M olecular M echanics RDF Radial Distribution Function SCF Self consistent f ield SS Secondary Subsystem SSD Solid State Drive SSE Strong Scaling Efficienc y vdW van der Wa a ls WHAM Weighted Histogram Analysis Method

PAGE 10

CHAPTER I INTRODUCTION 1.1 QM/MM Multiscale Modeling 1.1 .1 What is Multiscale Molecular Modelling? M ultiscale mode ls inclu de different resolutions of length or time scales within the same mo del to ensure a model can represent a system of interest with sufficient accuracy while remaining computationally feasible. 1 This has become an indispensable addition to the environments as evidenced by the proliferation of use of multiscale modelling methodology and by the awarding of the 2013 Nobel Prize in Chemistry to Martin K arplus, Michael Levitt, 2 Multiscale modelling has been applied to many chemical systems including biological, 3 materials science 4 and solid state systems. 5 Multiscale molecular modelling typically combines one or more of the fo llowing methods: quantum mechanics (QM), classical molecular mechanics (MM), coarse graining (CG) models and finite element meth ods. These methods cover a wide range of scales from the sub nanometer to beyond millimeters and are employed in various types of calculations including single point calculations, geometry optimizations, molecular dynamics (MD) and Monte Carlo (MC) simulations This work is devoted to one specific multiscale modeling technique, the combined quantum mechanics/molecular mechanics (QM/MM) methodology 1 1. 2 What is QM/MM? Quantum mechanics is essential to accurately describe many chemical phenomena such as electronic excitation, the making and breaking of bonds extensive polarization and

PAGE 11

2 charge transfer While semiempirical QM met hods have been used for systems with >1000 heavy atoms, 6 Density Functional Theory (DFT) and post Hartree Fock (post HF) methods are typically computationally feasible for up to hundreds of atoms 7 with d ynamics simulation times in the hundreds of picoseconds. 8 The computational expense of QM can be especially problematic if many atoms surrounding the system of interest must be considered for an accurate study. Fortunately for many chemical processes, changes of interest in the electronic structure are limited to a small number of atoms. Thus sometimes it is not necessary to describe the whole system with QM to obtain an acceptable level of accuracy 1 T he environment around the system of interest can be de scribed by computationally efficient MM empirical force fields MM system sizes can range to tens of millions of atoms 9 or a millisecond timescale 10 Additionally, the tota l simulation time for MM studies with h undreds of thousands of atoms now routinely reaches one microsecond By embedding a QM subsystem within a n MM subsystem, the accuracy of QM can be combined with the computational efficiency of MM, allowing for simulat ions with hundreds of QM atoms surrounded by hundreds of thousands of MM atoms at the 100 ps timescale 11 Warshel and Levitt published the first paper on QM/MM in 1974 12 which described new methodology to explore the electrostatic stabilization of a carbonium ion i ntermediate in the cleavage of hexasaccharides in the enzyme lysozyme. Field et al. extended QM/MM past single point calculations to include MD, M onte Carlo and geometry optimization computations in 19 90 13 The increase in computing power combined with this newly available methodology marked an explosi on in QM/MM method development throughout the 14 21

PAGE 12

3 1.2 Conventional QM/MM Methodology 1.2. 1 The MM Energy The MM energy of a prototypical 1 st generation force field can be calculated by an analytical fun c tion with a corresponding set of parameters called a force field. The energy c an be expressed as the combination of the bonded and non bonded interactions : T he bonded interactions can be expressed as the combination of the bond, angle, dihedral improper and cross term interactions : 18 is the interaction energy between all pair s of bonded atoms (e.g. between the atoms labeled 1 and 2 in Fig 1.2 called 1 2 ) is the interaction energy between all connected triplet s of atoms (e.g. atoms 1,2, and 3 called 1 3 ), and is the interaction

PAGE 13

4 energy from all sets of 4 atoms connected in a line (e.g. atoms 1,2,3 a nd 4 called 1 4 ). The symbols and represent the bond length, angle, and dihedral respectively for a given geometry and are force constant parameters The equilibrium bond distance and angle parameters are represented by and respectively. The dihedral expression has two additional terms, and which represent the torsional multiplicity and phase respectively is the energ y for out of plane bending, used to enforce planarity of atoms (f or example, to keep atoms 6, 7, 8, and the hydrogen connected to 7 on the same plane) is the force constant parameter, is the equilibrium out of plane angle, and is the measured out of plane angle. In Fig. 1.2 would represent the angl e between the hydrogen and the plane formed by atoms 6, 7, and 8. represents the interaction between two dihedrals, and where and are the multiplicity of and respectively, and where is the phase parameter for each cr oss term The cross term is commonly used to enforce specific angles between adjacent dihedrals in protein backbone s Additional bonded interactions can be added to the force field to increase the accuracy of the model. For example, atoms in the protein ba ckbone may be mapped to empirical energy maps, 22 23 and functions can be added which can account for the anharmon icity of bonds and angles. 24 Non bonded interactions are typically not computed between atoms included in 1 2 and 1 3 interactions, and may be excluded or scaled down for inclu ded in 1 4 interactions. A Fig 1. 1 : Atom indices for 1 ( 1 H imidazol 4 yl)propane:

PAGE 14

5 minimal set of non bonded interactions includes columbic interactions via point ch arges and van der Waals (vdW) interactions which implicitly model exchange dispersion interactions with the empirical Lennard Jones potential : Here, is a parameter for the charge assigned to atom and the vacuum permittivity of free space. In both equations, is the distance between non bonded atoms and For the van der Wa a ls equation (1. 10 ) and are paramete rs for the collision diameter between non bonded atoms and and the well depth of the interaction curve, respectively. These prototypical functions are similar to those from commonly used force fields, 25 such as the CHARMM 26 27 OPLS AA, 28 A mber 29 and GROMOS 30 force fields. 1.2. 2 The QM Energy The QM energy for potentials in this work can be found by obtaining app roximate solutions for the time independent non relativistic Schrdinger equ ation, which can be expressed as: 31 32 w here is the Hamiltonian operator, is the wave function, and is the energy of the system. If the wave function is known, the energy may be found through the following expression:

PAGE 15

6 The Hamiltonian opera tor consists of the nuclear and electronic kinetic and potential energy operators: where and are the nuclear and electronic kinetic energy operators respectively, and and are the potential energy operators for the interaction energy between the nuclei and electrons, between the electrons, and between the nuclei respectively. To reduce the dimensionality of the problem, t he Born Oppenheimer approximation is invoked, w hich neglects the coupling between the nuclear and electronic motion by considering the electrons to be moving in a field of fixed nuclei. Therefore, the wave function only explicitly depend s on the electronic coordinates, while parametrically depend s on t he nuclear coordinates. Thus can be set to 0 and is a constant. The total energy (in atomic units) becomes: where is the number of nuclei, is the charge of nucleus is the coordinate vector for nucle us is the electronic wave function which is explicitly dependent upon only

PAGE 16

7 the electronic coordinates, and is the Hamiltonian operator acting on the electronic wave function. Equation (1. 1.5 ) is called the electronic Schrdinger equation. To satisfy the anti symmetric requirement for the electronic wave function, the electronic wave function can be approximate d by the Slater determinant of the spin orbitals of the individual electrons, which are assumed throughout this work to be orthonormal: where is the spin orbital occupied by electron is the total number of electrons, and is the Dirac delta function which delineates the orthonormality between spin orbitals and For the model systems studied within this work, the wave function is adequately approximated by one Slater deter minant. The electronic Hamiltonian can be decomposed into the core Hamiltonian, which depends only upon the coordinates of one electron, and the two electron operator. The core Hamiltonian describes th e kinetic energy of an electron and the potential energ y of an electron in a field of nuclei. The two electron operator acting on the electronic wave function describes the coulombic and exchange interactions between electrons. These operators may be related to the electron indices and through the following equations:

PAGE 17

8 where is the core Hamiltonian for electron is the two electron operator for electrons and and is the coordinate vector for electron The electronic Schrdinger equation can be solved by e valuating the following equation: Aft er many intermediate steps, an expression for the electronic energy can be obtained: where and are the coulomb and exchange ope rators for electron respectively which represent the average local potential at from interactions with other electrons (note that is a non local operator) Th is treatment is known as the mean field approximation and it neglects coupli ng between the motions of individual electrons of opposite spins As the spin orbitals of other electrons are required for and the expression for is non linear. Approximate solutions can be obtained by iteratively improving an i nitial trail wave function until the energy change between two iterations falls below a specified threshold a process called the self consistent field (SCF) method. According to the variational principle, the lowest energy found by varying the spin orbita ls is an upper bound

PAGE 18

9 of the exact ground state energy. The process of finding the QM energy using these equations is called the Hartree Fock (HF) met hod 33 34 P ost H artree Fock (post HF) methods like 2 nd Order Mll er Plesset Perturbation Theory (MP2) 35 are based upon t he HF method and improve upon the HF energy at higher computational cost s 1.2. 3 The OM3n Method The OM3n method 36 used e xtensively throughout this work, is a semiempirical QM method. In s emiempirical methods many of the integrals are approximated through parameterized expressions to reduce computational expense Parameters are adjusted to reproduce higher level QM calculat ions and experimental data. Thus, the accuracy of a give n semiempirical method depends on the parameters used and the approximations made in the method. The OM3n model is a reparametrized version of the OM3 37 method that accurately describes proton transfer in water. This method accurately reproduces the MP2 energy surface for proton transfer in the Zundel ion, yielding a barrier for proton transfer in agreement with path integral calculations that in clude nuclear quantum effects of the proton 38 39 The OM3 n method is a member of the family of orthogonalization corrected methods (OM x x = 1, 2, 3 ), 40 41 which are based on the modified neglect of differential overlap (MNDO) model. 42 43 In the OM x methods core electrons are treated implicitly and valence electrons are treated explicitly A minimum, orthogonalized basis set of Slat er type atomic orbitals is used, called a Lwdin basis 44 45 Because the use of an orthogonalized basis can introduce significant errors in the one e lectron integrals, additional empirical functions called orthogonalization corrections are added to the one electron integrals in the OM x methods. As

PAGE 19

10 the two electron integrals are similar when evaluated in an orthogonal basis v er s us a nonorthogonal basis, no corrections are added to these terms. Similar to other MNDO type semiempirical models many integrals are approximated in the OM x methods for computational efficiency. 46 One center electron integrals are calculated as the sum of a parameter which represents the attraction of the electron to the nucleus it is based on (taken from spectroscopic data) and the potential due to the other nuclei in the system. Two ce nter one electron integrals are approximated based on an values of two center two electron integrals are based on semiempirical multipole multip o le interactions. Th e OM x methods invoke a modified version of the neglect of diatomic differential overlap (NDDO) approximation. 44 45 Thus, one electron integrals with three centers, and all three and four center two electron integrals are neglected. Core core electronic attractions and repulsions are evaluated through the corresponding integral, and summed with a parameterized function. A thorough description of the OM x formalism has been published by Dral et al. 41 1.2.4 Partitioning a QM/MM System In all QM/MM calculations, the system must be divi ded into a QM and MM partition ( Fig 1 .1 ). In conventional QM/MM methods, this partition is fixed and may not change throughout a simulation. The QM region, also called the primary subsystem (PS), is embedded in the region described with MM, the secondary s ubsystem (SS). The entire system (ES) consists of all atoms in the simulation. The choice of partitioning is one of the most important decisions in construction a QM/MM model, and great care must be taken to

PAGE 20

11 ensure that all important atoms are designated a s QM. The nature of the interaction between the PS and SS is one of the defining characteristics of a QM/MM method. 1.2. 5 The QM/MM Hamiltonian The general form of the QM/MM potential comprises the QM potential of the PS, the MM potential of the SS, and the interaction energy between the PS and SS. 12, 16 In this work, for simplicity, the mechanical embedding scheme is used where the QM calculation for the PS is carried out in the gas phase. The subtractive definition of the QM/MM energy is emplo yed, which has the following form: In our calculations the ES energy is computed at the MM level Then, the MM energy of the PS is subtracted from the ES. Finally, a QM calculation of the PS is carri ed out in vacuum and added to obtain the final QM/MM energy. Electrostatic interactions between the PS and Fig 1. 2 : In a QM/MM simulation the system is partitione d into two interacting parts the primary subsystem which is calculated using QM and the secondary subsystem which is calculated using MM

PAGE 21

12 SS are described at the MM level with the vdW and Coulombic nonbonded expressions. More complicated formulations exist which describe mutual polariza tion 18, 47 and partial charge transfer 48 50 between the PS and SS However, they were not employed in this work due to their associated computational expenses 1.2.4 QM/MM Partitions through Covalent Bonds It is sometimes necessary to place the PS/SS boundary through a covalent bond in order incorporate parts of larger molecules, such as the sidechains of proteins, in the QM region. Many different methods for treating this s ituation can be found in the literature 18, 51 However, only the link atom method with the scaled bond distance treatment 52 53 is discussed as it is the method implemented in our code, QMMM In this method, a n SS boundary atom (designated M) the atom covalently bonded to a PS atom (designated Q) is replaced by a cap atom (designated L) in the PS The cap atom is usually Hydrogen. The location of the link atom is plac ed o n the line between atom s M and Q. The Q L bond length is then scaled by an adjustable parameter, One recommendation for is to use the ratio of two equilibrium bond lengths: 52 where and are the equilibrium bond distances taken from the MM parameters for the Q H and Q M bonds, respectively When PS is capped with a link atom, it is called the capped primary subsystem (CPS). However, when the term PS is used in this work, it encompasses primary subsystems with and without cap atoms. The ad dition of the link atom has serious implications for th e

PAGE 22

13 QM/ MM system, as the topology of the CPS is different from that of the ES. Therefore, the topology of the CPS must be regenerated and the appropriate parameters must be available. 1.3 Adaptive Parti tioning QM/MM Methodology Adaptive partitioning (AP) methodology was developed to address two of the problems that arise when using QM/MM on diffusive systems : (1) molecules of interest may travel to the QM/MM boundary where their environment is not proper ly described with QM and (2) frequently a large PS must be used to alleviate boundary effects, which increases the computational expense AP methods address these problems by featur ing the on the fly reclassifications of atoms as QM or MM so that molecule s of interest are always in the QM subsystem As such, s maller QM subsyst ems can be used with AP methods There are many different types of AP methodology in the literature, which includ e the hot spot scheme, 54 layered integrated molecular orbital and molecular mechanics (ONIOM XS) scheme, 55 56 the permuted adaptive partiti oning (PAP) scheme 57 a nd modified PAP (mPAP) scheme 58 60 the sorted AP scheme (SAP), 57 60 the difference based adaptive solvation (DAS) scheme, 61 63 the buffered force (BF) scheme, 64 67 the number adaptive (NA) scheme, 68 69 the size consistent multi partitioning (SCMP) scheme, 7 0 the density based adaptive (DBA) scheme, 71 and the time adaptive switching interaction potential scheme. 72 Th is work concentrates o n the PAP and mPAP scheme s, due to their demonstrated numerical stability 57

PAGE 23

14 1.3.1 Partitioning of an AP QM/MM System In many AP schemes, the system is divided into three parts ( Fig 1.2 ): the active zone (QM part), environmental zone (MM part), and buffer zone (dual QM/MM characteristics). I n this work, the partitions are distance based. First, a center for the active zone is chosen. This can be a stationary point, such as a dummy atom, or a moving point, such as a specific ion or binding site center ( defined by the center of mass of several residues ) Next, atoms are divided into groups which can be whole molecules or molecular fragment s. 73 Groups are added to a zone based on the distance from the active zone center to either the center of mass of the group or a designated reference point (such as the O atom in the H 2 O molecule) according to equation (1. 3.1 ) All atoms within a sphere of radius can be assigned to the active zone. All atoms within a spherical shell between and are designated as buffer groups. This shell typically has a width of 1 2 All atoms greater than from the active zone center are assigned to Fig 1.2 : An example of partitioning the system into three zon es for adaptive partitioning. The zones are defined as a set of concentric spheres, defined by the radii and The blue zone is the environmental zone and is treated with MM. The yellow zone is the buffer zone and has dual QM/MM characteristics. The white zone in the middle is the active zone and is treated with QM only

PAGE 24

15 the environmental zone. If need ed, a group can be designated as always QM or MM regardless of its location. The list of atoms belonging to each zone is updated on the fly throughout the course of the simulation and atoms may change groups if bond breaking and bond forming reactions occ ur 1.3.2 The PAP Potential The PAP potential has the form of a many body expansion and is expressed with the following equa tion: 57 E quation 1.3.2 can be formatted more concise ly as: 73 In these equations, represents the QM/MM energy with the active zone treated at the QM level, represents the QM/MM energy with the th buffer group and the active zone treated

PAGE 25

16 at the QM level, represents the QM/MM energy with the th and th buffer groups and the active zone treated at the QM level. This equation is expanded to the limit of where all buffer groups and the active zon e are treated at the QM level. Buffer groups are classified as QM in some partitioning configurations and MM in others which yield their dual QM and MM character. However, the active zone is always treated by QM and the environmental zone is always treate d by MM. is the smoothing function for the th buffer group, which varies smoothly and monoton ical y in value between 0 and 1 depending on the position of the buffer group As a group enters the buffer zone at the limit of the environmental zone, the smoothi ng function is equal to 0. As a group enters the buffer zone at the limit of the active zone, the smoothing function is equal to 1. In our implementation of PAP, a fifth order spline is used for the smoothing function: where is a metric for the position of the buffer group calculated by: In total, configurations must be evaluated to obtain the fu ll PAP potential. However, this is rarely necessary as high order many body contributions are usually insignificant. Moreover, the product of the smoothing functions tends towards 0, further decreasing the contribution of higher order terms. T runcation of the potential at to the 3 rd order has been shown to capture over 95% of the intermolecular interaction energy with other potentials so the expansion can typically be truncated to increase computational

PAGE 26

17 efficiency 7 4 76 Because many commonly used QM methods formally scale as where is the number of basis functions, it can be faster to compute multiple smaller QM calculations instead of one large one. In these cases, the use of AP methodology will lead to an increase in computational efficiency. 1.3.3 The PAP Gradient s and the mPAP Scheme The PAP scheme conserves both the total energy and momentum, 57 and may be used in dynamics simulations under the microcanonical ( N VE ) and canonical ensembles ( N VT ). Unfortunately, forces due to the derivative of the smoothing function can be large enough to cause simulation artifacts such as erroneous solvation structures 73, 77 For example, the full PAP potential for a system with one group in the buffer can be expressed as: The forces can then be found by computing the gradient of equation (1 .33 ) with respect to the coordinates : The following terms fr om equation (1. 3.7 ) are deemed the physical forces of the system and : The following term from equation (1. 3.7 ) represents the unphysical forces of the system 58, 78 and is called the ri

PAGE 27

18 These forces arise due to the mismatch of potentials between the QM and MM systems, and cease to exist whe n the QM and MM potentials are equal. Current ly, it i s at least practically impossible to match QM and MM potentials with high levels of dimensionality prohibiting potential matching for most systems Therefore, a different approach is necessary to remove the transition forces. In the modified PAP (mPAP) scheme, the forces are simply deleted which is equivalent to adding external forces that exactly cancel out the transition forces In this treatment, the QM and MM systems are considered to be open system s in dynamic equilibrium with each other while the entire system is a closed system Because extra forces are added to the system, the entire system is no w a non Hamiltonian system. Therefore, the mPAP method cannot be used in microcanonical simulations, but can be coupled to thermostats and used under the canonical ensemble. The mPAP method has been shown to satisfactorily reproduce solvation structures observed in conventional QM/MM simulations with QM zones so large that boundary effects on the solvatio n become insignificant 58 59 1.4 Proton Transfer in Complex Environments 1.4.1 Grotthuss Shuttle Mechanism of Proton Transfer Proton transfer is a ubiquitous process in many biological and chemical systems but spe cific details of proton transfer, such as the stability of intermediary structures in the mechanism, are still a subject of debate 79 80 P roton diffusion in water is approximately seven times faster than for a Na + cation indicating that H 3 O + does not diffuse through solution as a discreet molecule M echanisms for proton diffusion are different from that of other

PAGE 28

19 + or K + and involve the hydrogen bonding 81 network for water. In these mechanisms, the proton can be considered a charge defect which is shuttled through changes in the covalent and hydrogen bonding topology of water. The most commonly accepted mechanism is t he Grotthuss shuttling mechanism, named after von Grotthus Although significantly different from the modern mechanism this early mechanism is the first published description detailing protons hopping from O to O atom thr ough a water wire. 82 83 Modern mechanism s describe proton transfer through water as a dynamic reorganization of t he hydrogen bonding environment of water surrounding the excess proton The proton is typically found in two states, the Z u ndel 84 (H 5 O 2 + ) and Eigen 85 (H 9 O 4 + ) cations ( Fig 1. 3 ) which interconvert during transfer events 86 The most commonly accepted transferring mechanism in bulk water c an be characterized as Eigen Zun del Eigen, where the proton spends most of the time in the Eigen resting s tate, and transfers through a Zu ndel transition state to another Eigen state. 86 88 It has been suggested that t hese jumps sometimes occur in a quasiconcerted fashion over two to three waters, with two jumps more likely to occur than three or one. 89 The transfer event is coupled to the dynamics of the water clusters surrounding the proton, 90 and the reorganization of both the first and second Fig 1.3 : Eigen and Zu ndel structures of the solvated proton and the Grothuss Mechanism for proton transfer. The hydrogen bonding network of the water changes to shuttle an excess charge defect through solution.

PAGE 29

20 solvation shells into specific arrangements 38, 87 called presolvation 91 93 is crucial Other mechani sms can be characterized as Zundel Zundel, where molecules spend most of their time in Zundel states and transfer directly from Zundel to Zundel state. 80, 94 This mechanism has been observed for proton transfer in confined spaces, such as carbon nanotubes. 95 The free energy barrier of proton transfer was calculated to be approximately 0.56 kca l mol 1 in Car Par r inello MD simulations 38, 96 P ath integral molecular dynamics which include quantum fluctuations of the nuclei predict the barrier to be 0. 15 kcal mol 1 38 39 T hese barriers are close to or less than at 300 K ( approximately 0.6 kcal mol 1 ) Therefore, the hydrogen bond s between the donor and acceptor O atoms and the transferring H at the Zundel transition state are classified as a low barrier hydrogen bond, 97 98 where nuclear tunneling effects are insignificant This provides justification for approximating the nuclei in water as a classical particle in some dynamics simulations. However, w hen lower probability events, such as water autoprotolysis, are important to the chemical model, it may be necessary to incorporate nuclear quantum effects. 99 T hroughout the years proton transfer has been described with theories and models at various levels Full quantum molecular dynamics are the most computationally expensive methods, but have provided the energetic data cited above and many mech anistic insights. 38, 88, 98 100 Empirical valence bond (EVB) theory by Warshel, 101 and its extension the multistate empirical valence bond ( MS EVB) theory have also been applied towards studying proton transfer in water and other envir onments 36, 94, 96, 102 114 In addition many r eactive force field methods 115 124 and semiempirical QM methods, 125 133 including t he OM3n model 134 used extensively in this work, have been dev eloped.

PAGE 30

21 1.4. 3 The Proton Indicator A major challenge in proton modelling is that the identity of the transferring proton is not known ahead of time and changes constantly throughout a n MD simulation. This difficulty is compounded in biased simulation s wh ere it can be unclear which atoms to apply restraints to To overcome this problem, the proton indicator used in this work first suggested by Hofer et al. 135 was developed by Wu et al. 134 to track the approximate location of the excess proton over time It was then extended by Pezeshki and Lin 59 to update the active zone center centered at the indicator in adaptive QM/MM simulations. The indicator is used in this work as a collective variable in proton transfer reactions The indicator wa s constructed to fulfill certain criteria. 6, 59 The position of the indicator should be superimpose d on the location of the shared proton in the limit of the Zundel cation, and on the donor O atom in the limit of t he Eigen cation. The location of the indicator is a function of the coordinates of the donor oxygen all hydrogen atoms bonded to and all possible acceptors within a predefined enlisting radius from D. The most likely acceptor is denoted as When the proton is transferred and switch The location of the indicator is determined by a linear comb ination of the coordinates of and all possible acceptors : Here, is a normalization factor, is the weight function associated with the ratio of the projected donor acceptor vector, and and are the Cartesian coordinates of the donor oxygen and the j th acceptor oxygen respectively. is the total

PAGE 31

22 number of hydrogen atoms bonded to and is the total number of acce ptors within a radius of size The ratio is a metric for how close a given atom is to donor versus an acceptor and is calculated according to the following equation: w here and The linear combination of coordinates is normalized by which is calculated according to equation (1.3 9 ) using the weight function described in equation (1. 40 ): The weight function depends on the reduced variable x : w here is a parameter set to a distance slightly larger than the equilibrium OH bond distance in a hydronium ion (1.0 ) and is the threshold distance for enlisting possible acceptors

PAGE 32

23 CHAPTER II IMPLE M ENTATION OF N AMD QMMM I NTERFACE 2.1 QMMM An Interfacing Program QMMM 136 is an interfacing program for conducting QM/MM sin gle point calculations geometry minimizations, and molecular dynamics simulations. As a driver, QMMM calls an MM package (NAMD 137 or Tinker 138 ) and/or a QM package ( Gaussian 03 139 or 09, 140 GAMESS, 141 ORCA, 142 or MNDO 143 ) to calculate the energ ies gradient s and/or hessian s Q MMM parses the output from the selected MM or QM packages, synthesizes e ne rgies and energy derivatives, and propagates the trajectory ( Fig 2.1 ). The source code files and associated scripts within the QMMM package now contain >70,000 lines, and are written in FORTRAN 95, FORTRAN 2003, C shell (CSH), Perl, and Python. The program also requires a compiler that supports OpenMP version 3.0. 2.2 Motivation for Adding the NAMD Interface Before th e work described in this thesis was completed, Tinker was the only MM package which interfaced with QMMM While Tinker is a flexible, open source, and free Figure 2.1 : Structure of QMMM code. QMMM calls a QM package from the programs on the left and an MM package from the programs on the right. The interface to NAMD was added in this work.

PAGE 33

24 MM package, the serial and parallel implementations of the program do not rival the speed of other MM packages. In QM/MM calculations for a very large ES and small PS, when employing a semiempirical QM potential and using Tinker the MM calculation can be comparable to or slower than the QM calculation. NAMD was specifically designed to run large MM systems quickly in high performance computing (HPC) environments A t th e time of this publication, NAMD is available in the software stack for NERSC and XSEDE HPC resources Moreover, it is freely available and open source Tools for creating NAMD input files are also integrated within the widely used visualization software, VMD 144 VMD automates many of the steps for the creation of input files necessary for running a calculation within QMMM and cuts down on syst em building time Finally, t he TCL scripting language can be used to extend the functionality of NAMD. Therefore, it is possible to add functionality (such as for specialized MM calculations needed for more advanced embedding schemes) without modifying the source code of NAMD. 2.3 Structure of the NAMD Interface To use NAMD as the MM package for MM single point gradient calculaitons the user first create s a protein data bank (PDB) file of the entire system (ES) using the worldwide PDB f ormat version 3.30 145 Because QMMM uses the element symbol from the column described in the standards for assigning the element to an atom for QM calculations, a helper script find_elsym .py included wi th in the QMMM distribution assists the user with acquiescing to the PDB file standards by finding the elemental symbol for each atom listed in the PDB and appending it to the proper column. The user then needs to create an X PLOR formatted CHARMM PSF file which contains the topological information and atom types

PAGE 34

25 for the system. Finally, the user suppl ies the force field to use for the MM calculation ( Fig 2.2 ) QMMM prepares NAMD calculations by writing 4 files: the coordinates binary file of the system, t he PSF text file for the system, the input text file for NAMD, and the PDB text file for the system if the system is the primary system (PS) The ES PDB file is used to supply the ES coordinates to QMMM at the first step and to fulfill the requirements th at a PDB exists for an ES calculation in NAMD. Therefore, the ES PDB file is never rewritten through the course of a QM/MM calculation. However, the PS PDB, PS PSF, coordinate binary file, and the NAMD input file are re written at every step. After the fil es are written, QMMM calls a CSH shuttle script, namdshuttle to invoke the NAMD executable The shuttle script is configurable by the user, and can Figure 2.2: Schematic showing files used for communication in the QMMM interface with NAMD. Blue represents the files that QMMM creat es, red represents the files that NAMD creates, and green represents the files that the user must supply to run a QM/MM calculation.

PAGE 35

26 accommodate the use of different environments (e.g. using ibrun to call the MPI/Charm++ parallelized vers ion of NAMD). The user is also able to choose the number of MPI/OpenMP processes that NAMD will use a feature that is important for using parallel permuted adaptive partitioning and discussed in Chapter III If the NAMD calculation is successful, NAMD w ill create several binary files and a text file with the output. QMMM only reads the binary force file for the gradient, and the text output file for the energy. 2.4 Topology Changes and the PSF One of the primary challenges with all QM/MM methods is deali ng with the topology of the MM system. If the QM partiti on cuts through covalent bonds, there will be new bonds, angles, and dihedrals that were not in the ES PSF. Parameters must be supplied for these new interactions if they are not already described in the forcefield Additionally, if bonds are made or broken within the QM region, the topology needs t o be updated to avoid erroneous MM calculation s In principle, i t is true that when using a subtractive QM/MM energy definition, the topology of the CPS sho uld not matter as the MM energy of the CPS is subtracted from the MM energy of the ES. However, due to the structure of the NAMD program, this is not always the case. NAMD decomposes the MM system into discrete spatial units called patches for parallelizat ion purposes Patches comprise groups of heavy atoms as well as hydrogen atoms connected to them. If the topology of the system changes due to a chemical reaction, it is possible that atoms can move out of the spatial region of their patch. If this occurs the MM energy of the CPS and/or ES will not be calculated as per the potential of the force field and the QM/MM calculation is no longer valid. Therefore, the topology of the system must be generated on the fly as the simulation progresses.

PAGE 36

27 Updating the topology at every step is undesirable as it takes a significant amount of time to write the ES PSF for large systems, a situation which is compounded when using optical hard drives (HDD) or when writing to a network file system (e.g. in HPC environments). Therefore, the topology for the ES is only updated when a bond is designated as formed/broken. Currently the only automatic procedure for changing the topology implemented within QMMM deals with proton transfer and is discussed in Chapter IV. It is possib le that there will be no universally applicable rule for handling topology changes and new rules will need to be added on a case by case basis. When the topology is generated, lists of bonds, angles, dihedrals, impropers, and cross terms are created as de scribed in the following sections 2.4.1 The Bond List The initial bond list is created by reading the list of bonds in the from the user supplied PSF file. The bond list for the entire system is updated when a topology change occurs Add itionally, a separate bond list is created for the PS if the QM/MM partition cuts through covalent bonds. 2.4.2 The Angle List The angle list is generated from the bond list. Firs t the number of angles is calculated using number of bonds conn ected to the i th atom, from the following formula :

PAGE 37

28 The angles list is equal to the set of all angles [ j, i, k ] defined by equation (2.4.2) where is the set of all atoms bonded to atom i : 2.4.3 The Dihedrals List The number of dihedrals can be calculated from the following formula: The dihedral s list is equal to the set of all dihedrals [ k i j l ] defined by equation (2.4.4): 2.4.4 Th e Impropers and Cross Terms No cross terms or impropers are automatically generated when the topology changes as these are impossible to predict without specific knowledge of the molecule These terms are read in from the user supplied PSF and stored in m emory. If the CPS contains all atoms of an improper term with cap atoms corresponding to the atoms they substitute that term is added to the CPS topology. Cross terms are added to the C PS topology in the same fashion. If a section of the protein backbone is present in the CPS, cross terms which include the cap atom are not evaluated and the CPS potential will not exactly cancel out the ES potential for subtractive schemes. This is an uncommon situation in QM/MM models.

PAGE 38

29 2.5 Benchmarking The following benc hmark tests use QMMM version 2.2.7, a CUDA enabled version of NAMD 2.11, and Tinker 7.1.2. Each calculation consists of a n MM MD calculation of a water box using the selected MM package with t he TIP3P forcefield. Periodic boundaries are treated using the m inimum image convention with the electrostatic interactions cut off at 14 The number of atoms in each benchmark var ies Full benchmarking details can be found in Appendix A .1 The benchmarking data is shown in Table 2.1 with the NAMD speedup calculated as the average time per step of the NAMD calculation divided by the average time per step of the Tinker calculation ( Equation 2.4.4 ) The benchmarking data indicates a linear speedup as the number of atoms i n the system increases ( Fig 2.3 ). The system sizes were chosen to approximate the number of atoms of various biomolecular system s that this methodology is currently applied to within the Lin Group. For small systems (<5,000 atoms), Tinker may still be a fe asible solution for QM/MM dynamics. However, when conducting QM/MM dynamics simulations of larger systems, NAMD is the only feasible choice Table 2.1 : Benchmar king Results for NAMD and Tinker Number of Atoms NAMD time/step (s) Tinker time/step (s) NAMD Speedup 4308 0.2 0.9 5 52233 0.6 35.2 64 96003 0.9 111.8 123 176454 1.6 330.9 221

PAGE 39

30 Fig 2.3: Speedup when using NAMD versus the number of atoms. The speedup is calculated using equation (2.5)

PAGE 40

31 CHAPTER III I MPLEMENTATION OF P ARALLEL A DAPTIVE P ARTITIONING 3.1 Theoretical Background As shown in section 1.3 .2, the PAP potential can be expressed as thus: Where is the smoothing function of the th group based on the distance from the act ive zone center, is the QM/MM energy with only the active zone treated at the QM level, and is the QM/MM energy with the active zone and the th through th groups treated at the QM/MM level. In the full PAP expansion, there are terms where is the number of buffer groups. The partitioning configuration s are independent QM/MM calculation s and may be calculated simultaneously 3.2 Implementation of Parallel Permuted Adaptive Partitioning T he parallel PAP algorithm (PPAP) wa s implemented in our code, QMMM Currently, parallel energy and gradient calculations for MD simulations can be conducted using Tinker 138 or NAMD 137 for the MM package, and Gaussian 03, 139 Gaussian 09, 140 or MNDO 143 for the QM package. QMMM was parallelized using the FORTRAN bindings from the OpenMP 3.0 protocol 146 147 a parallel language for writing par allel code on shared memory multiprocessor architectures Fig 3. 1 shows the flowchart for the implemented PPAP algorithm. After reading the initial coordinates and user defined groups, QMMM obtains the entire system MM gradient for use in each of the indi vidual QM/MM calculations. Next the program calculates and

PAGE 41

32 generates the list of groups in the active and buffer zones from this distance. Then, is calculated for each of the buffer groups and the list of partitioning configurations is gen erated based on the truncation order The parallel region is then entered, and the main process spawns a number of threads based on an environmental variable set by the user Each thread initializes a private data structure and calculates partitioning conf igurations at the QM/MM level. After that, t he program enters a critical region where the threads update the total energy and gradient in shared memory one at a time. Finally, t he thread t eam is then destroyed and the program resumes normal function. Wo rk is assigned to the thread team using t chunk size of 1. This means that each thread is initially given a partitioning configuration to calculate. When a thread finishes computing its assigned partitioning configurati on, it is Fig 3.1: Sch eme of parallel permuted adaptive partitioning implementation. Red boxes indicate serial instructions and green regions indicate parallel instructions

PAGE 42

33 given another until all partitioning configurations are calculated. D ynamic scheduling was chosen to allow for more flexibility in assigning work to the thread team because the time to calculate a partitioning configuration is not known at runti me and may drastically vary between configurations Due to the representation of numbers in finite precision in the computer, the computed gradients for a step may be slightly different if run with different numbers of threads or on different CPU architect ures A s the previous program design relied on FORTRAN common blocks for many energy and gradient functions, substantial sections of the entire program had to be rewritten to prevent data races in the parallel section. Additionally, the program was moved t o dynamically allocating memory for most arrays at runtime, from using statically allocated arrays based on parameters in the code which had to occasionally be modified. This change necessitated modifying almost every singl e subroutine within the program. As most of the static arrays for both the ES and CPS previously had one dimension determined by parameter set to be larger than the maximum number of atoms in the system, the memory profile of the program was reduced. This likely affected the PPAP algorit hm the most, as these large arrays would have had to be continually reallocated for each member of the thread team. 3.3 Scaling and Benchmarks Because the potential calculations in the PPAP calculation may be calculated simultaneously, t he theoretical mini mum time per step is proportional to the time of the longest QM calculation since MM calculations usually take much less time Additionally, t h is PPAP algorithm reaches its fullest potential when there are many buffer group molecules, as there are more pa rtition ing configurations to assign to the thread team. Though

PAGE 43

34 uncommon, i t is possible that there are fewer QM calculations than available threads at a particular step In this case, the extra threads do nothing. The PPAP algorithm is benchmarked based on a n MD simulation of 1435 H 2 O molecules and 1 H 3 O + molecule with the active zone center set to the location of the proton indicator. is set to 5 and is set to 6 yielding an active zone with a radius of 5 encircled by a 1 thick buffer zone shell. There were about 160 QM calculations per step during the simulation, with an average of 13. 5 groups in the active zone and 17. 5 groups in the buffer The AP potential was truncated at the 2 nd order. Due to the large number of QM calculations per step, this benchmark system is ideal for examining the parallel efficiency as the number of threads increase. Additional informa tion about the benchmark systems can be foun d in Appendix A .2 T iming examples of ongoing experiments are also provided there but the parallel efficiency is not rigorously benchmarked in them. Under ideal calculation conditions, the PPAP algorithm scale d approximately linearly on a local workstation ( Fig 3. 2 ). calculations to spread across a thread team and that the calculations are run on a sufficiently fast disk drive. The strong scaling efficiency (S SE) can be used to measure how efficiently each thread in a thread team works compared to the serial calculation T he higher SSE values correspond to a more efficient thread team The SSE is calculated using equation 3.2, where is average wall time per step for the calculation with one thread, is the number of threads used in a calculation, and is the average wall time per step for the calculation with threads :

PAGE 44

35 The SSE can also be thought of as a metric for the parallel overhead for the calculation. The SSE in benchmark calculations ranges from 100% to 85% as the calculation was parallelized across all the physical cores on a loc al workstation ( Fig 3.2 ) T he time it takes to calculate one QM/MM potential is vastly greater than the cost of creating the thread team and initializing the private data structures Because QMMM interfaces with the QM and MM packages by reading and writin g files to disk, file I/O congestion may contribute some overhead This is evidenced by the difference in SSE between two calculations ran on the same computer but us ed file systems with 2 different speeds In each system examined, the SSE is typically hig her for the version ran on the faster filesystem than the slower file system ( Fig 3.3 ) The amount of file I/O that occurs can even cause the simulation to crash if essential files become inaccessible by QMMM or the MM or QM package due to file system tra ffic This occurred frequently when using the LUSTRE network file system on the XSEDE Stampede HPC resource ( Fig 3.3 ) and infrequently when using local file systems. The likelihood of a calculation crash increases as the filesystem speed decreases i.e. calculations crash less on a SSD versus a HDD This obstacle could be overcome by distributing the Fig 3. 2 : Scaling of PPAP with increasing number of cores on local system described in Appendix A.3 N ot e that this computer only has 10 physical cores but 20 virtual cores. Left: Average wall time per step versus increasing number of cores. Right: Strong scaling efficiency (SSE) versus increasing number of cores.

PAGE 45

36 calculations and file operations across multiple nodes using the MPI protocol, by using multiple file systems on one node, or by developing an interface bet ween packages that does not rely on files written to disk. Interestingly, these benchmarks of the PAPP algorithm benefit ed from using virtual cores when they were r u n on an SSD. This may be attributed to the fact that many of the benchmarking QM calculatio ns were on the order of 200 ms. Therefore, some threads were performing file system operations while others were performing QM calculations. For longer QM calculations, this may not be the case as QM calculations will be running on the all the virtual core s simultaneously. Typically, when running compute intensive processes, using more virtual cores than physical cores results in no speedup or slowdown M axim izing parallel performance also requires specific knowledge of the QM or MM package. If the SSE of Q MMM at a certain thread count is less than the SSE of a package, then less threads of QMMM should be used and more threads of the package should be used. T iming data for MNDO using QMMM and Gaussian 09, 148 indicate that QMMM should be r u n with as many threads as possible before disk I/O speed becomes in issue and that the Fig 3. 3 : Scaling of PPAP with increasing nu mber of threads on the XSEDE Stampede HPC resource described in Appendix A.3 N ote that this computer only has 16 physical cores but 32 virtual cores. Tests on an internal node disk (blue) and network file system (red) were conducted. Tests with 32 threads crashed on the network file disk due to I/O congestion. Left: Average wall time per step versus increasing number of threads. Right: Strong scaling efficiency (SSE) versus increasing number of threads.

PAGE 46

37 QM packages should use the remaining thr eads. For systems with few partitions, the QM or MM packages can be assigned more threads to avoid idle processors. Additionally, if the thread team has more threads than there are QM calculations, the SSE of the algorithm will suffer ( Fig 3. 4 ) Fig 3.4 : Strong scaling efficiency of two QM calcula tions with different buffer sizes. The buffer sizes were 1.5 and 0.5 which corresponded to an average of 115.0 and 9.3 QM calculations per step respectively. As there were not enough partitioning configurations to spread to the thread team, the SSE of the 0.5 buffer size calculation s quickly decreases

PAGE 47

38 CHAPTER I V IMPLEMENTATION OF I NDICATOR R ESTRAINT 4 1 Restraint Potential and Decomposed Gradients A harmonic restraint was chosen for the proton indicator for its intended usage with umbrella sampling and the Weighted Histogram Analysis Method (WHAM). 149 The restraint potential can be expressed as: i s the harmonic restraint potential imposed on the indicator I, is the force constant of the restraint potential, is the origin of potential in Cartesian coordinates, and is th e vector for the current position in Cartesian coordinates of the proton indicator. T he position of the indicator, is a function of the donor and the possible acceptor atoms, T he gradient is taken with respect to both of these categorie s of atoms T he calculation of the normalization factor for the weight function, and weight function are treated as constants when taking the derivatives at any given time step of the simulation. When the gradient of these values was taken with respect to and the gradients were too large, causing simulation artifacts. ( Appendix B.5 ) The gradient of the potential with respect to the j th acceptor or donor can be expressed by the following equa tions:

PAGE 48

39 The symbols and represent the location of a possible acceptor or the donor respectively. The direction of the gradients on all or are along to the vector ( Fig 4.1 ) : The for m of ( Equation 1.4.4 ) is designed such that atoms receive none of the gradient when they are further than away from the donor, smoothly receive gradient as they become an acceptor, and have the same grad ient in the limit of the dono r acceptor swap Fig 4. 1 : 3d depiction of gradients on the acceptor (red triangle), donor (yellow triangle), and indicator (green triangle) for the Zundel (Left) and Eigen (right) cations. The potential origin is the b gradient unit vectors are represented by the arrows. The gradients all have the same direction but may have different magnitudes.

PAGE 49

40 4 2 Behavior of the Indicator for the Zundel Cation S ingle point calculations were performed on the H 5 O 2 + (Zundel) structure to examine the gradients on the indicator for a simple system The geometries were constructed such that and the transferring proton are aligned along the x axis. Three different sets of single point calculations were performed. In s can A ( Fig 4. 2 ) the shared proton was translated between the two oxygen atoms. In s can B, a H 2 O molecule was translated from towards an H 3 O + molecule to create the Zundel cation. In scan C the entire system was translated with respect to the restraint potential origin Scan A ( Fig 4. 3 ) shows the decomposed gradients with the respect to the translation of the proton. During the scan, the indicator changed position with the proton, but the potential origin stayed static. Two origins of potential were investigated, one placed on the location of the shared hydrogen (H) and one placed on an oxygen (O 1 ) The d isplacement between the indicator and origin is calculated by: The graphs show that when the donor and acceptor swap ( at = 0.0 for the origin at H, and = 1.2 for the origin at O 1 ), the gradients are equivalen t. Additionally, the Fig 4.2: Depictions of two scans. In s can A a proton is translated across the x axis between two water molecul es. In scan B the blue water mo lecule is translated across the x axis from towards O 1 Two origin locations are chosen for each scan, one at O 1 and one at H. Red spheres are oxygen and white are hydrogen

PAGE 50

41 Fig. 4.3 : For s can A : d istance of indica tor from origin (top) and decomposed gradients for the Zundel cation (bottom) with the indicator centered at the transferring proton (left) and the indicator centered at the donor oxygen (right). The indicator is in red, donor green, and acceptor blue. Not e the y axis is different for the gradient graphs. The force constant is 5 kcal mol 1 2 Fig. 4.4 : For s can B : d istance of indicator from origin (top) and decomposed gradients for the Zundel cation (bottom) with the indicator centered at the transferrin g proton (left) and the indicator centered at the donor oxygen (right). The indicator is in red, donor green, and acceptor blue. Note the y axis is different for the gradient graphs. is 3.5 and the force constant is 5 kcal mol 1 2

PAGE 51

42 1 Fig 4.3 ) shows that this transition occurs smoothly when the origin potential is located away from the indicator. Scan B ( Fig 4. 4 ) depicts the gradients as H 2 O is translated away from the complex fro m 1.5 to past = 3.5 The complex starts as the donor at 1.5 and becomes an acceptor at 2 Its gradient tends towards zero as it reaches As the H 2 O translates, the value of the acceptor constantly changes. The ratio of the gradient s decomposed onto the donor and acceptor are proportional to the magnitude of the of the acceptor, with a larger leading to more gradient decomposed onto the acceptor. The gradient of the donor will always be larger in magnitu de th an the gradient of the acceptor except for at a donor acceptor swap when they are equal. When the pot ential origin is located at O 1 and the complex is translated past the gradients on all atoms remain at 0 and the gradient is continuous This is because the bond length is equal to the parameter. If this is not the case, a discontinuity in the position of the indicator and thus the gradients, occurs ( Fig 4. 5 ). While this happens to a minor extent throughout a simulation, a bond length would have to be substantially different from its equilibrium bond length to create a drastic discontinuity Scan C ( Fig 4. 6 ) shows the gradients decomposed on to a Zundel complex that is Fig 4. 5 : A 0.4 1 H bond length of 1.2 which is larger than the parameter of 1.0

PAGE 52

43 translated with respect to a static restraint potenti al origin In this scan, the value of the acceptor is constant. Thus, changing the position of the complex changes the magnitude of the donor and acceptor gradients but the ratio of the gradients remains constant. 4 3 Behavior of the Indicator for Collision Simulations NVE MD calculations were conducted on a system where one H 2 O molecule was placed 5.40 from H 7 O 3 + cluster and given an initial velocity of 0.01 fs 1 or 0.001 fs 1 ( Fig 4. 7 ). The velocities of all other atoms were set to 0 fs 1 All atoms were descr ibed using the OM3n potential. The simulations were r u n for 2 ps with 0.2 fs time steps F our different restraint force constants were examined. The full set of simulation results can be found in Appendix C. 1 T he restraints keep the indicator closer to th e origin of the potential after the collision than in the free simulation ( Fig 4. 8 ) When the restraint is active, the distance from the potential of origin is lower than in the absence of the restraint. Also, the total energy of the system increased when the restraint was active ( Fig 4. 9 ). Therefore these restraints should Fig 4. 6 : For s can C: total gradient on indicator I and decomposed gradients on and with respect to the change in the origin of the restraint potential along the x axis Red spheres indicate oxygen, white spheres hydrogen, and the pink sphere is the indicator. The force constant is 5 kcal mol 1 2

PAGE 53

44 Fig 4. 9 : Total energy of system for simulations of the Zundel (A) and Eigen (B) cations with different force constants. The blue, green, red, and cyan lines correspond to force constants of 0 (no restraint), 5, 10, and 20 kcal mol 1 2 respectively. The initial velocity of the H 2 O was 0.001 fs 1 Fig 4. 8 : Distance ( ) of indicator from origin of potential for Zundel (top) and Eigen(bottom) cations. Red, bl ue, green, and cyan correspond to force constants of 0 (no restraint), 5, 10, and 20 kcal mol 1 2 respectively. The initial velocity of the H 2 O was 0.001 fs 1

PAGE 54

45 only be used in simulations coupled to a thermostat. Additionally, in the set of 0.01 fs 1 simulations for the Zundel cation, the decomposed gradients were strong enough to cause erroneous structures leading to SCF converge nce failures ( Appendix C.1 ) However, these structures were not observed in the E igen set. This is likely because there are more atoms to distribute the gradient to in the E igen cluster, leading to reduced grad ients on the individual atoms. Thus, low force constants should be used when the indicator can move far from the origin and there are few possible acceptors for example, when studying proton diffusion in confined spaces 4 4 MD Simulations of a Water Sphe re MD simulations of a sphere of H 2 O molecules surrounding 1 H 3 O + cation were conducted with and without the use of the indicator restraints to compare the difference in potential of mean force ( PMF ) for proton transfer and radial distribution functions (R DF) of water Four simulations were conducted with the force constant for the indicator restraint set at 5 kcal mol 2 For comparison, four simulations were conducted in the abse nce of the indicator restraint, resulting in a total of 8 simulations. The s imulations used a droplet model similar to Wu et al. 134 with the inner 12 sphere described at the OM3n 134 level and the outer 10 shell described at the MM level with the SPC/Fw potential 150 R estraints wi th a force constant of 20 kcal mol 1 2 were applied to every MM oxygen atom. The full simulation details can be found in Appendix C. 2 The radial distribution functions (RDF) were calculated using saved snapshots from the simulations excluding those fro m the first 0.5 ps The radial distribution function was centered on the donor D. If the donor atom was located greater than 7 away from the center of mass of the droplet (within 5 of the QM/MM border ) the snapshot was excluded

PAGE 55

46 from the analysis to mi nimize the influence of MM waters on the RDF. In the simulations with no indicator restraint the indicator quickly migrated far from the center ( Fig 4. 10 ) and only 7455 out of 20000 possible snapshots were usable. All snapshots were usable when the restr aints were used. The D O RDF ( Fig. 4 10 ) indicates close agreement between the simulations conducted with and without the restraint. The noise in the line comes from the reduced amount of sampling due to the exclusion of snapshots. This agreement su ggests the presence of the indicator restraint does not significantly alter the solvation shell surrounding the proton The potential of mean force (PMF) for the proton transfer between the donor oxygen and the most likely acceptor oxygen was calculated t hrough Boltzmann inversion for two different reaction coordinates ( Fig 4. 1 1 ) The same set of snapshots was used for the PMF as the RDF. The most likely acceptor was the acceptor with the highest value at a given time step The reaction coordinate was the value for the most likely acceptor in a snapshot calculated in equation (1.4.2). The reaction coordinate was c alculated as the Fig 4 10 : ( A ) Average distance of indicator from restraint potential origin (center of clus ter) for restrained (blue) and unrestrained (red) simulations. ( B ) Area normalized radial distribution function g D O (r) for restrained (blue) and unrestrained (red) simulations, where r is the distance from the donor,

PAGE 56

47 abso lute value of the difference between the donor proton distance and acceptor proton distance ( Equation 4.8 ). The barrier heights for the PMF ( Table 4. 1 ) are measured at the reaction coordinates = 0 and = 0.50, and correspond to a Zundel state. T he barrier height for the restrained and free simulations are similar, with the difference in barrier height for at 0.03 kcal mol 1 and for at 0.11 kcal mol 1 The barrier heights similar to the previously reported barrier heights for the OM3n potential in the cluster model of 0.19 and 0.06 kcal mol 1 for and respectively by Wu et al. 134 as well as for those obtained from path integral calculations of 0.16 kcal mol 1 for 38 39 The minima at = 0.25 agrees wi th previous calculations 38 39, 96, 134 and represents the Eigen resting state. Thus, the use of the indicator restraints does not seem to significantly affect the free energy profile for proton transfer in water. Fig 4. 11 Potential of mean force (PMF) for proton transfer for reaction coord inates (A) and (B) The blue line corresponds to the PMF with the indicator restraint active, while the red line corresponds to a freely moving indicator.

PAGE 57

48 Table 4.1: Barrier heights for PMF calculations of proton transfer in a water cluster. The barrier height is the energy that corresponds to the reaction coordinate of = 0 or = 0.50 Reaction Coordinate Restrained/Free Barrier Height (kcal mol 1 ) Restrained 0.15 Free 0.26 Restrained 0.22 Free 0.25

PAGE 58

49 CHAPTER V A PP LICATION OF INDICATO R R ESTRAINT TO C ARBON N ANOTUBE 5.1 Introduction A model syst em of a carbon nanotube embedded between two parallel graphene sheets similar to Peng et al. 36 was constructed for three reasons: ( 1 ) to examine the behavior of the indicator in a greatly simplified model of a proton transport protein, ( 2 ) to explore and compare hydration dynamics of the tube with Peng et al, 36 and ( 3 ) to examine the viability of obtaining a PMF for a proton transfer process using the indicator restraints The carbon nanotube represents a simplistic model of a hydrophobic cavi ty surrounding a proton transfer relay, 151 such as in the E. Coli CLC Chloride Ion Transport Protien 152 Thus, examining the hydration dynamics of the carbon nanotube may provide insights into the mechanism of water wire formation in hydrophobic cavities in biom olecular systems. The nanotube system based on Peng et al., 36 contains a single walled, armchair type (6,6) carbon nano tube with a length of 28.2 and diameter of 8.1 The nanotube bridges two graphene sheets with x and y dimensions of 42.99 and 41.12 respectively ( Fig 5.1 ) The nanotube dimensions allow for the formation of a 1 dimensional water wire along the tub e axis, and have been used in several previous experiments to study proton transfer through a water wire. 36, 95, 151, 153 A slab of water with 7 pairs of K + /Cl ions was added to each side of the nanotube. The SPC F/w 150 water model was used in conjunction with the CHARMM22 154 parameters for the carbon atoms (atom type CA) and ions. Similarly to Peng et al., 36 t he parameter of the carbon nanotube atoms was changed to 0.0448 kcal mol 1 to scale the depth of the interpolated potential well for the vdW interaction between the oxygen and carbon atoms to 80% of its original value. However in this system, the paramet er of the

PAGE 59

50 graphene sheets was also adjusted while in Peng et al., 36 the standard CHARMM22 parameter was used. This increa sed the hydrophobicity of the tube, which would otherwise be able to spontaneously hydrate. Visual inspection of the MM trajectories confirmed that only transient water wires of 2 4 waters could form. This system was then equilibrated as per the procedu re described in Appendix C. 3 The initial geometries for the following QM/MM simulations obtained from the end of the equilibration runs. For the following QM/MM simulations, the MM waters and K + /Cl ions were restrained with a force constant of 1 nN (ap proximately 14.39 kcal mol 1 2 ). No K + /Cl ions were included in the QM region. The positions of all carbon atoms were fixed. All simulations were r un under the NVT ensemble using the Nos Hoover thermostat at 300 K with no t reatment of periodic boundar ies and with 0.25 fs time steps. Peng et al. 36 reported that an excess protonic charge restrained at the opening of the nanotube could recruit waters from the bulk into the tube. They demonstrated this behavior using the SPC/Fw water model and the multistate empirical valence bond version 3 (MS EVB3) model 110 wit h the multiscale reactive molecular dynamics (MS RMD) method. 103 104, Fig 5.1 : Carbon nanotub e system used throughout the chapter. The nanotube axis is aligned with the z axis.

PAGE 60

51 108, 110, 155 108 (CEC) of the system, which tracks the approximate location of the excess proton and differs from the indicator restraints described in Chapter IV. The CEC is calculated by the following equation: where is the center of charge vector of the H 3 O + ion in the th MS EVB state, and is the amplitude of the th MS EVB state. Peng et al. described the water recruitment process as follows T wo waters are initially in the nanotube. When the proton nears the nanotube, it brings its solvation shell with it filling the tube with two more waters for a total of four waters. The excess proton becomes trapped in an energy minima at the pore opening. Waters transfer through the proton into the nanotube as the proton enters the nanotube, and the tube is fully hydrated by the time the proton moves 2 3 into the tube. This work aims to examine whether this process occurs in QM/MM simulations with the OM 3n potential and with the previously described indicator restraints. 5.2 Simulations of Indicator Restrained at Nanotube Opening In this set of simulations, the restraint potential origin was placed 2.5 above the center of the nano tube to look for eviden ce of water recruitment into the nanotube ( Fig 5.2 ). T he closest water to the restraint potential origin was protonated. To examine the effects of water in the tube, the number of waters initially in the tube was varied from zero to four. Four starting geo metries with different atomic coordinates were created for each of these

PAGE 61

52 cases for a total of 20 simulations For each starting geometry, two simulations were conducted: one with the force constant set to 5 kcal mol 1 1 and one with the force constant set to 0 kcal mol 1 1 All waters within a 14 radius of the middle of the pore opening we re included in the QM subsystem. The simulations w ere ru n for up to 10.5 ps The full simulation details can be found in Appendix C. 4 For simulations with the in dicator restraint, waters begin to be recruited into the nanotube within the first few hundred femtoseconds ( Fig 5.3 ). This begins with migration of an excess proton towards bulk. As the excess charge moves towards the bulk solution (in the positive z dire ction) the indicator moves with it and away from the restraint potential origin. Forces on donor and acceptor atoms which were decomposed from the restraint force on the indicator, move the donor and acceptors towards the nanotube (in the negative z dire ction) in order to push the indicator towards the restraint origin. This moves waters into the tube. The proton then migrates back towards the bulk and the process repeats itself, eventually creating a water wire in the nanotube.

PAGE 62

53 Fig 5.3: Water recruitment into the pore when restraints are applied to the indicator The images depict the indicator (pink sphere), donor (blue sphere), most likely acce ptor (yellow sphere), restraint potential origin (green sphere), water (ball and stick, red is oxygen and white is hydrogen), nanotube (cyan lines), and restraint forces (black arrows). Some waters are omitted for clarity. Restraint forces are shown in pan els two and three. ( 1 0 fs): The initial position of the indictor (pink) begins close to the restraint potential origin (green) Zero waters are initially inside the nanotube. ( 2 2 fs): A Zundel cation is formed and the indicator migrates towards the a cceptor (yellow). Forces (black arrows) are applied to the donor (blue) and acceptor in the direction of the nanotube. ( 3 38 fs): The donor and acceptor have switched places. New waters have been recruited as possible acceptors. Forces are applied to a n ew list of possible acceptors. ( 4 174 fs): The indicator continues to be drawn towards the bulk. The original donor has moved to the tube entrance. ( 5 252 fs): The original donor has moved inside the tube. The acceptor from (4) is moving towards the tu be opening. ( 6 480 fs): Three waters are inside the nanotube. At the end of this simulation, 12 waters have been recruited into the tube.

PAGE 63

54 For simulations with n o restraint, the indicator quickly moved to the outside of the QM region ( Fig 5 .4 ) indicating a strong driving force for the indicator to move away from the tube opening towards bulk Some s imulations were terminated due to QM SCF convergence difficulties from a diffuse water wire which formed in the tube ( Fig 5. 4 ) This occurred due to water shooting into the nanotube, or initial waters in the tube drifting towards the other side of the nanotube. These simulations were truncated to a length less than 10.5 ps. The water occupancy (WO) of the nano tube, defined as the number of oxygen atoms within the tube at a given time step, was much higher for simulations with the indicator restraint than for simulations with no restraint ( Fig 5.4 ). This indicates that t he proton indicator restraints play a large role in recruiting waters into the nanotube. When the restraint Fig 5. 4 : ( A ) Average distance of indicator from restraint potential origin v er s us simulation time for simulations t hat were restrained (blue) and unrestrained (red) ( B ) Picture showing diffuse water wire causing SCF convergence failure. ( C D ) Average water occupancy of nanotube over simulation time for simulations with indicator restraints ( C ) and without restraints ( D ) for different numbers of initial waters in the tube.

PAGE 64

55 was applied, the water occupancy of the the nano tube increased within the first two picoseconds The water occupancy platea ued at an average of 11 w aters which is close to reported by Peng et al. 36 In the absence of indicator restraint s, water recruitment into the nanotube was less frequent. When zero waters were initially in the tube, no water recruitment was observed. When one water was initially in the tube, on average, the waters exited the nanotube. However, when 2 4 waters were initially in the nanotube, water recruitment into the nanotube was observed, albiet at a much slower pace than when restraints were applied. As the average indicator position was approximately 10 from the restraint potential origin and 12.5 from the en trance of the tube after 4 ps of simulation time, it is unclear whether the presence of the excess proton around the nanotube entrance aids water recruitment into the tube. 5.3 Simulations of Indicator Restrained at Different Distances from Nanotube The ef fect of different indicator restraint potential origins with respect to the entrance of the nanotube was examined to determine the effect of the restraint potential origin with respect to the water occupancy. The simulation setup was similar to that used i n S ection 5.2 with a QM hemisphere of water located above one nanotube entrance. However, the restraint potential origin was moved 3.5, 5.0, and 6.5 above the center of the nanotube entrance ( z direction) Two other potential origins with the origin shi fted 2.5 to the side of the nano tube opening ( x direction) with heights of 3.5 and 5.0 above the nanotube entrance ( z direction) were also created for a total of five sets of simulations ( Fig 5.5 ). Four different starting geometries were created for ea ch set, for a total of 20 simulations over all sets. No further distances could be explored without the risk of significant boundary effects from the

PAGE 65

56 QM/MM border. There were no waters in the tube at the start of the simulations Complete simulation detai ls can be found in Appendix C.5 For all the restraint potential origin locations, the average water occupancy over time increased ( Fig 5.5 ). This indicates that the water molecules were recruited into the tube even when the indicator is restrained at a l ocation of up to 6.5 above the tube. As the average water occupancy without the restraint was 0 waters ( Fig 5.4 ) the increase in water can be attributed solely to the presence of the restraint. The average water occ upancy after 8 ps for these restraint potential origins (< 8.75 waters) was less than the water occupancy when the restraint was 2.5 above the tube (9.5 11 waters, Fig 5.4 ). N o direct dependence of water occupancy on vertical distance above the tube can be drawn from this data. When the rest raint potential origin was shifted away from the center of the tube, there was less water recruitment than when the restraint potential was above the tube. However, as the radius for enlisting possible acceptors was 3.5 in the simulations acceptors abov e the tube opening could be subjected to restraint forces and drawn into the nanotube. Fig 5. 5 : Graph of average water occupancy at different restraint locations. The color of the line corresponds with the sphere above the nanotube to indicate the location of the restraint potential o rigin The legend that the restraint was shifted 2. 5 in the x direction from the center of the tube.

PAGE 66

57 5. 4 QM/MM Umbrella Sampling of PMF for Proton Transfer in Nanotube To provide proof of concept for using the indicator restraints for umbrella sampling, the nanotube w as filled with water and a PMF for proton transfer was obtained using the z coordinate of the indicator as the reaction coordinate. Each umbrella window used different QM/MM partitions ( Fig 5.6 ) W aters on both sides of the nano tube were described as Q M. F orce constants ranging from 2.5 7.5 kcal mol 1 2 in the z direction were used. The full set of simulation details can be found in Appendix C.6 The barrier height for proton transfer in the nanotube from z = 30.1 to z = 14.1 from the PMF ( Fig 5.7 ) is 119.9 kcal mol 1 a value that is certainly too large The PMF is relatively flat from z = 14 (middle of nanotube) to z = 21 D uring exploratory calculations no waters were recruited into the tube during this region. However, as the indicator moved pa s t this region, waters were vigorously recruited into the nanotube, shuttling waters from the QM hemisphere in the positive z direction towards the negative z direction. This effect was so pronounced, that when z > 25 close to 20 waters were Fig 5.6 : (A) Location of umbrella window or igins along the reaction coordinate ( z axis). Horizontal red lines indicate umbrella window locations and the black dot is the tube center. (B) Lowest umbrella window origin ( z = 14) and (C) highest umbrella window origin ( z = 30). The QM waters are depict ed. The circle indicates the umbrella window origin.

PAGE 67

58 Fig 5.7 : ( A ) Histogram of samples along the umbrella reaction coordinate. Different colors indicate different umbrella window origins. ( B ) PMF for proton transfer from the middle of the nanotube towards the bulk. ( C ) Waters recruited through tube at different umbrella origins in preliminary umbrella sampling calculations. Positive water recruitment indicates net water movement from the positive z direction into the tube, negative water recruitment is in the opposite direction. ( D ) Restraints applied to the umbrella window centered at z = 14 Restrained and unrestrained QM waters are shown with blue and red respectively

PAGE 68

59 recruit ed into the tube leading to an unphysical density in the upper half of the tube Therefore, the lower half of the QM subsystem was restrained to prevent th is recruitment ( Fig 5.7 ) increasing the PMF by approximately 40 kcal mol 1 in this region ( Appendix C.6 ) T he mechanical embedding scheme may be another contributing factor for the steep PMF curve As the indicator migrates into bulk, the PMF curve should flatten, not steepen. There is thus likely a strong driving force for the indicator to continue mov ing to the outside of the QM region. This is evidenced by the migration of the indicator to the outside of the QM region in unrestrained simulations for the water cluster ( Fig 4. 9 ) and nanotube systems ( Fig 5.4 ) This driving force exacerbates the water re cruitment mechanism observed in these simulations but it is unknown what the cause of this force is. Further analysis and calculations using more accurate embedding schemes and different potentials w ill be necessary to determine the causes of this behavio r Despite the problems in the model, the se simulations do shed light on the basis of the restrained shooting mechanism. These results show that for the speedy recruitment of water into a dry nanotube, the excess proton must move away from the nanotube. Th is causes forces to be applied to the surrounding waters i n the direction of the nanotube, shuttling the water into the tube. However, with this system, it appears that the probability of finding the excess charge around the nanotube opening is extremely l ow. Therefore the recruitment of water into the nanotube with the excess charge defect may be more of a simulation artifact than feasible mechanism. This mechanism becomes more feasible if there is an energy minima for the excess proton close to the nano tube opening as is the case with the MS EVB model 36, 153 because a proton needs to be by the pore opening in order to recruit waters into it. Nevertheless, the

PAGE 69

60 potential is asymmetric with respect to the minima a t the pore opening, with a steeper barrier for motion into the tube and a shallower barrier towards bulk. Therefore, the excess charge will, on average, drift towards bulk causing water to be recruited into the tube For future studies, a DFT or ab initio potential could be used in order to confirm whether such a minima exists. Nevertheless, these results indicate that this indicator restraint may become a viable option for obtaining the PMF for chemical processes involving proton transfer As the distribu tion of forces is straightforward, it is easy to describe the bias applied to the simulation. However, further investigation is necessary to determine the exact causes behind the PMF generated from the restraints.

PAGE 70

CHAPTER VI CONCLUSION In this work, thre e primary topics were considered. First, an interface between NAMD and QMMM was created to increase the speed of MM calculations. Benchmark calculations for MM MD simulations show that the speedup ranges from a factor of 5 for small calculations to > 200 fo r larger calculations. In order to use the interface for QM/MM simulations where bonds were formed and broken, algorithms for the on the fly generation of the topology were implemented. Next, the permuted adaptive partitioning algorithm present in QMMM was parallelized using the OMP protocol. The partitioning configurations may now be calculated simultaneously when using Tinker and NAMD for the MM package and Gaussian 03 or 09 and MNDO for the QM package. The parallel speedup is approximately linear when t here are enough QM/MM calculations to distribute to the thread team and when the calculation is conducted on a sufficiently fast file system to minimize the overhead of file handling The QMMM package was also updated to dynamically allocate memory for mo st arrays making the program more user friendly and reducing the memory profile of the application. Finally, a harmonic restraint potential was applied to the proton indicator. The gradients of the restraint potential indicator were decomposed onto the at oms which determine the location of the indicator. The gradients on the donor and acceptor atoms have the same direction of the gradient on the indicator. The ratio of the decomposed gradients depend s on the ratio of the values for the acceptor s The indicator restraints do not conserve energy. However, the PMF and radial distribution function for water was not significantly affected by the use of the restraints.

PAGE 71

62 The indicator restraints were applied to a carbon nanotube system in order to provide proof of concept for using the restraints to determine the PMF for a process with a high energy barrier These simulations confirm ed results by Peng et al. 36 which indicated that excess protons restrained near the opening of a hydrophobic space cause water to enter that space. The recruitment of water by restraint of the excess charge has now been confirmed with two differ ent simulation techniques: through our work with QM/MM and a proton indicator, and by Peng et al. 36 with MS EVB and restr aints on the center of excess charge. The reasons for water recruitment are not yet fully understood but it is apparent that the proton indicator felt a strong driving force to leave the nanotube opening. We suspect that the motion of the indicator couple d with the restraint forces recruited water into the tube, and the behavior could be a simulation artifact instead of a physical process. More research is needed to understand these observations.

PAGE 72

63 R EFERENCES (1) Warshel, A., Multiscale modeling of biological functions: From enzymes to molecular machines ( N obel L ecture). Angew. Chem. Int. Ed. Engl. 2014, 53 (38), 10020 10031. (2) The N obel P rize in C hemistry 2013. https://www.nobelprize.org/nobel_prizes/chemistry/laureates/2013/ (accessed Jan 2017). (3) Senn, H. M.; Thiel, W., Q M / MM methods for biological systems. In Atomistic A pproaches in M odern B iology: From Q uantum C hemistry to M olecular S imulations Rei her, M., Ed. Springer Berlin Heidelberg: Berlin, Heidelberg, 2007; pp 173 290. (4) Elliott, J. A., Novel approaches to multiscale modelling in materials science. Int Mater Rev 2011, 56 (4), 207 225. (5) Rudd, R. E.; Broughton, J. Q., Concurrent coupling of length scales in solid state systems. Phys. Status Solidi B 2000, 217 (1), 251 291. (6) Wu, X.; Koslowski, A.; Thiel, W., Semiempirical quantum chemical calculations accelerated on a hybrid multicore cpu gpu computing platform. J. Chem. Theory. Comput. 20 12, 8 (7), 2272 2281. (7) Siegbahn, P. E.; Himo, F., Recent developments of the quantum chemical cluster approach for modeling enzyme reactions. J. Biol. Inorg. Chem. 2009, 14 (5), 643 651. (8) Ong, M. T.; Verners, O.; Draeger, E. W.; van Duin, A. C.; Lord i, V.; Pask, J. E., Lithium ion solvation and diffusion in bulk organic electrolytes from first principles and classical reactive molecular dynamics. J. Phys. Chem. B. 2015, 119 (4), 1535 1545. (9) Zhao, G.; Perilla, J. R.; Yufenyuy, E. L.; Meng, X.; Chen, B.; Ning, J.; Ahn, J.; Gronenborn, A. M.; Schulten, K.; Aiken, C.; Zhang, P., Mature HIV 1 capsid structure by cryo electron microscopy and all atom molecular dynamics. Nature 2013, 497 (7451), 643 6. (10) Shaw, D. E.; Maragakis, P.; Lindorff Larsen, K.; Piana, S.; Dror, R. O.; Eastwood, M. P.; Bank, J. A.; Jumper, J. M.; Salmon, J. K.; Shan, Y.; Wriggers, W., Atomic level characterization of the structural dynamics of proteins. Science 2010, 330 (6002), 341 346.

PAGE 73

64 (11) Martins Costa, M. T.; Ruiz Lopez, M. F ., Reaching multi nanosecond timescales in combined QM / MM molecular dynamics simulations through parallel horsetail sampling. J. Comput. Chem. 2017, 38 (10), 659 668. (12) Warshel, A.; Levitt, M., Theoretical studies of enzymic reactions: Dielectric, elect rostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103 (2), 227 249. (13) Field, M. J.; Bash, P. A.; Karplus, M., A combined quantum mechanical and molecular mechanical potential for molecular dynamics s imulations. J. Comput. Chem. 1990, 11 (6), 700 733. (14) Gao, J. L., Hybrid quantum and molecular mechanical simulations: An alternative avenue to solvent effects in organic chemistry. Acc. Chem. Res. 1996, 29 (6), 298 305. (15) Gao, J.; Furlani, T. R., Si mulating solvent effects in organic chemistry: Combining quantum and molecular mechanics. IEEE Comput. Sci. Eng. 1995, 2 (3), 24 33. (16) Lin, H.; Truhlar, D. G., Q M / MM : What have we learned, where are we, and where do we go from here? Theor. Chem. Acc. 20 07, 117 (2), 185 199. (17) Mulholland, A. J., Chapter 14 T he QM / MM approach to enzymatic reactions. In Theoretical and computational chemistry Leif, A. E., Ed. Elsevier: 2001; Vol. Volume 9, pp 597 653. (18) Senn, H. M.; Thiel, W., Q M / MM methods for bio molecular systems. Angew. Chem. Int. Ed. Engl. 2009, 48 (7), 1198 1229. (19) Riccardi, D.; Schaefer, P.; Yang, Y.; Yu, H.; Ghosh, N.; Prat Resina, X.; Konig, P.; Li, G.; Xu, D.; Guo, H.; Elstner, M.; Cui, Q., Development of effective quantum mechanical/mol ecular mechanical ( QM/MM ) methods for complex biological processes. J. Phys. Chem. B 2006, 110 (13), 6458 6469. (20) Lyne, P. D.; Hodoscek, M.; Karplus, M., A hybrid QM MM potential employing H F ock or density functional methods in the quantum region J. Phys. Chem. A 1999, 103 (18), 3462 3471.

PAGE 74

65 (21) Monard, G.; Prat Resina, X.; Gonzlez Lafont, A.; Lluch, J. M., Determination of enzymatic reaction pathways using QM/MM methods. Int. J. Quantum Chem. 2003, 93 (3), 229 244. (22) Mackerell, A. D.; Feig, M .; Brooks, C. L., Extending the treatment of backbone energetics in protein force fields: Limitations of gas phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J. Comput. Chem. 2004, 25 (11), 1400 1415. (23) MacKerell, A. D.; Feig, M.; Brooks, C. L., Improved treatment of the protein backbone in empirical force fields. J. Am. Chem. Soc. 2004, 126 (3), 698 699. (24) Maple, J. R.; Dinur, U.; Hagler, A. T., Derivation of force fields for molecular mec hanics and dynamics from ab initio energy surfaces. Proc. Natl. Acad. Sci. U.S.A. 1988, 85 (15), 5350 5354. (25) Mackerell, A. D., Jr., Empirical force fields for biological macromolecules: Overview and issues. J. Comput. Chem. 2004, 25 (13), 1584 1604. (2 6) Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D., Optimization of the additive CHARMM all atom protein force field targeting improved sampling of the backbone J. Chem. Theory Comput. 2012, 8 (9), 3257 3273. (27) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M., C HARMM : A program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 1983, 4 (2), 187 217. (28) Jorgensen, W. L.; Maxwell, D. S.; Tirado Rives, J., Development and testing of the OPLS all atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118 (45), 11225 11236. (29) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J.; Kollman, P., A point charge force field for molecular mechanics simulations of proteins based on condensed phase quantum mechanical calculati ons. J. Comput. Chem. 2003, 24 (16), 1999 2012. (30) Schmid, N.; Eichenberger, A. P.; Choutko, A.; Riniker, S.; Winger, M.; Mark, A. E.; van Gunsteren, W. F., Definition and testing of the GROMOS force field versions 54a7 and 54b7. Eur. Biophys. J. 2011, 4 0 (7), 843 856.

PAGE 75

66 (31) Jensen, F., Introduction to computational chemistry John Wiley & Sons: 2016. (32) Szabo, A.; Ostlund, N. S., Modern quantum chemistry: Introduction to advanced electronic structure theory Courier Corporation: 2012. (33) Fock, V., Nh erungsmethode zur lsung des quantenmechanischen mehrkrperproblems. Z. Phys. 1930, 61 (1), 126 148. (34) Hartree, D. R., The wave mechanics of an atom with a non coulomb central field. Part I Theory and methods. Math. Proc. Cambridge 1928, 24 (01), 89 11 0. (35) Mller, C.; Plesset, M. S., Note on an approximation treatment for many electron systems. Phys. Rev. 1934, 46 (7), 618 622. (36) Peng, Y.; Swanson, J. M. J.; Kang, S. g.; Zhou, R.; Voth, G. A., Hydrated excess protons can create their own water wir es. J. Phys. Chem. B. 2015, 119 (29), 9212 9218. (37) Scholten, M., Semiempirische verfahren mit orthogonalisierungskorrekturen: Die om3 methode. Ph. D. Dissertation, Max Planck Institut fr Kohlenforschung 2003 (38) Marx, D.; Tuckerman, M. E.; Hutter, J. ; Parrinello, M., The nature of the hydrated excess proton in water. Nature 1999, 397 (6720), 601 604. (39) Marx, D.; Tuckerman, M. E.; Parrinello, M., Solvated excess protons in water: Quantum effects on the hydration structure. J. Phys. Condens. Matter 2 000, 12 (8a), A153 A159. (40) Dral, P. O.; Wu, X.; Sprkel, L.; Koslowski, A.; Thiel, W., Semiempirical quantum chemical orthogonalization corrected methods: Benchmarks for ground state properties. J. Chem. Theory Comput. 2016, 12 (3), 1097 1120. (41) Dral P. O.; Wu, X.; Sprkel, L.; Koslowski, A.; Weber, W.; Steiger, R.; Scholten, M.; Thiel, W., Semiempirical quantum chemical orthogonalization corrected methods: Theory, implementation, and parameters. J. Chem. Theory Comput. 2016, 12 (3), 1082 1096. (42) Dewar, M. J. S.; Thiel, W., Ground states of molecules. 38. The MNDO method. Approximations and parameters. J. Am. Chem. Soc. 1977, 99 (15), 4899 4907.

PAGE 76

67 (43) Dewar, M. J. S.; Thiel, W., Ground states of molecules. 39. M NDO results for molecules containing h ydrogen, carbon, nitrogen, and oxygen. J. Am. Chem. Soc. 1977, 99 (15), 4907 4917. (44) Lwdin, P. O., On the nonorthogonality problem. In Adv. Qunatum chem. Per Olov, L., Ed. Academic Press: 1970; Vol. Volume 5, pp 185 199. (45) Lwdin, P. O., On the non atomic wave functions in the theory of molecules and crystals. J. Chem. Phys. 1950, 18 (3), 365 375. (46) Thiel, W., Semiempirical quantum chemical methods. Wiley. Interdiscip. Rev. Comput. Mol. Sci. 2014, 4 (2), 145 157. (47) Lin, H.; Truhlar, D. G., Redistributed charge and dipole schemes for combined quantum mechanical and molecular mechanical calculations. J. Phys. Chem. A 2005, 109 (17), 3991 4004. (48) Pezeshki, S.; Lin, H., Molecular dynamics simulatio ns of ion solvation by flexible boundary QM/MM : On the fly partial charge transfer between QM and MM subsystems. J. Comput. Chem. 2014, 35 (24), 1778 1788. (49) Zhang, Y.; Lin, H., Flexible boundary quantum mechanical/molecular mechanical calculations: Par tial charge transfer between the quantum mechanical and molecular mechanical subsystems. J. Chem. Theory Comput. 2008, 4 (3), 414 425. (50) Zhang, Y.; Lin, H., Flexible boundary QM/MM calculations: Ii. Partial charge transfer across the QM/MM boundary that passes through a covalent bond. Theor. Chem. Acc. 2010, 126 (5 6), 315 322. (51) Reuter, N.; Dejaegere, A.; Maigret, B.; Karplus, M., Frontier bonds in QM/MM methods: A comparison of different approaches. J. Phys. Chem. A 2000, 104 (8), 1720 1735. (52) Da pprich, S.; Komiromi, I.; Byun, K. S.; Morokuma, K.; Frisch, M. J., A new oniom implementation in G aussian98. Part I The calculation of energies, gradients, vibrational frequencies and electric field derivatives. Theochem 1999, 461 462 1 21.

PAGE 77

68 (53) 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 (9), 1170 1179. (54) Kerdcharoen, T.; Liedl, K. R.; Rode, B. M., A QM/MM sim ulation method applied to the solution of L i + in liquid ammonia. Chem. Phys. 1996, 211 (1 3), 313 323. (55) Kerdcharoen, T.; Morokuma, K., O NIOM XS : An extension of the ONIOM method for molecular simulation in condensed phase. Chem. Phys. Lett. 2002, 355 ( 3 4), 257 262. (56) Kerdcharoen, T.; Morokuma, K., Combined quantum mechanics and molecular mechanics simulation of Ca 2+ /ammonia solution based on the ONIOM XS method: Octahedral coordination and implication to biology. J. Chem. Phys. 2003, 118 (19), 8856 8862. (57) 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. (58) Pezeshki, S.; Davis, C.; Heyden, A.; Lin, H., Adaptive partitioning QM/MM dynamics simulations: 3. Solvent molecules entering and leaving protein binding sites. J. Chem. Theory Comput. 2014, 10 (11), 4765 4776. (59) 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. (60) Pezeshki, S.; Lin, H., Recent developments in QM/MM methods towards open boundary multi scale simulations. Mol. Simul. 2015, 41 (1 3 ), 168 189. (61) Bulo, R. E.; Ensing, B.; Sikkema, J.; Visscher, L., Toward a practical method for adaptive QM/MM simulations. J. Chem. Theory Comput. 2009, 5 (9), 2212 2221. (62) Bulo, R. E.; Michel, C.; Fleurat Lessard, P.; Sautet, P., Multiscale modelin g of chemistry in water: Are we there yet? J. Chem. Theory Comput. 2013, 9 (12), 5567 5577.

PAGE 78

69 (63) 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 (39), 12401 12414. (64) 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 (65) Mones, L.; Jones, A.; Gtz, A. W.; Laino, T.; Walker, R. C.; Leimkuhler, B.; Csnyi, G.; Bernstein, N., The adaptive buffered force QM/MM method in the CP2K and A mber software packages. J. Comput. Chem. 2015, 36 (9), 633 648. (66) Peguiron, A.; Colo mbi Ciacchi, L.; De Vita, A.; Kermode, J. R.; Moras, G., Accuracy of buffered force QM/MM simulations of silica. J. Chem. Phys. 2015, 142 (6), 064116. (67) Vrnai, C.; Bernstein, N.; Mones, L.; Csnyi, G., Tests of an adaptive QM/MM calculation on free ene rgy profiles of chemical reactions in solution. J. Phys. Chem. B. 2013, 117 (40), 12202 12211. (68) 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. (69) Takenaka, N.; Kitamura, Y.; Koyano, Y.; Nagaoka, M., An improvement in quantum mechanical description of solute solvent interactions in condensed systems via the number adaptive multiscale quantum mechanical/molecular mec hanical molecular dynamics method: Application to zwitterionic glycine in aqueous solution. J Chem Phys 2012, 137 (2), 024501. consistent multipartitioning QM/MM : A stable and efficient adaptive QM/MM meth od. J. Chem. Theory Comput. 2014, 10 (10), 4242 4252. (71) Waller, M. P.; Kumbhar, S.; Yang, J., A density based adaptive quantum mechanical/molecular mechanical method. ChemPhysChem 2014, 15 (15), 3218 25. (72) Bckmann, M.; Doltsinis, N. L.; Marx, D., Ad aptive 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 (6), 2429 2439.

PAGE 79

70 (73) Pezeshki, S.; Lin, H., Adaptive partitioni ng 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 (11), 3625 3634. (74) Deev, V.; Collins, M. A., Approximate ab initio energies by systematic molecular fragmentation. J. Chem. Phys. 2005, 122 (15), 154102. (75) Elrod, M. J.; Saykally, R. J., Many body effects in intermolecular forces. Chem. Rev. 1994, 94 (7), 1975 1997. (76) Hirata S.; Valiev, M.; Dupuis, M.; Xantheas, S. S.; Sug iki, S.; Sekino, H., Fast electron correlation methods for molecular clusters in the ground and excited states. Mol. Phys. 2005, 103 (15 16), 2255 2265. (77) Jiang, T.; Boereboom, J. M.; Michel, C.; Fleurat Lessard, P.; Bulo, R. E., Proton transfer in aque ous 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. (78) Boereboom, J. M.; Potestio, R.; Donad io, D.; Bulo, R. E., Toward H amiltonian adaptive QM/MM : Accurate solvent structures using many body potentials. J. Chem. Theory Comput. 2016, 12 (8), 3441 3448. (79) Agmon, N., Infrared spectroscopy: The acid test for water structure. Nat. Chem. 2016, 8 (3 ), 206 207. (80) Thmer, M.; De Marco, L.; Ramasesha, K.; Mandal, A.; Tokmakoff, A., Ultrafast 2 D IR spectroscopy of the excess proton in liquid water. Science 2015, 350 (6256), 78 82. (81) Steiner, T., The hydrogen bond in the solid state. Angew. Chem. In t. Ed. Engl. 2002, 41 (1), 49 76. (82) de Grotthuss, C. J., Memoir on the decomposition of water and of the bodies that it holds in solution by means of galvanic electricity. Biochim. Biophys. Acta. 2006, 1757 (8), 871 875. (83) Von Grotthuss, C. J. T., Su r la dcomposition de l'eau et des corps qu'elle tient en dissolution l'aide de l'lectricit galvanique. Ann. Chim. 1806, 58 54.

PAGE 80

71 (84) Schuster, P.; Zundel, G.; Sandorfy, C., The hydrogen bond: Recent developments in theory and experiments: Structure an d spectroscopy. I I North Holland: 1976. (85) Eigen, M.; de Maeyer, L., Self dissociation and protonic charge transport in water and ice. Proc. Roy. Soc. London Ser. A 1958, 247 (1251), 505 533. (86) Agmon, N., The G rotthuss mechanism. Chem. Phys. Lett. 19 95, 244 (5 6), 456 462. (87) Markovitch, O.; Agmon, N., Structure and energetics of the hydronium hydration shells. J. Phys. Chem. A 2007, 111 (12), 2253 2256. (88) Tuckerman, M.; Laasonen, K.; Sprik, M.; Parrinello, M., Ab initio molecular dynamics simula tion of the solvation and transport of H 3 O + and OH ions in water. J. Phys. Chem. 1995, 99 (16), 5749 5752. (89) Hassanali, A.; Giberti, F.; Cuny, J.; Khne, T. D.; Parrinello, M., Proton transfer through the water gossamer. Proc. Natl. Acad. Sci. U.S.A. 2 013, 110 (34), 13723 13728. (90) Lapid, H.; Agmon, N.; Petersen, M. K.; Voth, G. A., A bond order analysis of the mechanism for hydrated proton mobility in liquid water. J. Chem. Phys. 2005, 122 (1), 14506. (91) Marx, D., Proton transfer 200 years after vo n G rotthuss: Insights from ab initio simulations. ChemPhysChem 2006, 7 (9), 1848 1870. (92) Chandra, A.; Tuckerman, M. E.; Marx, D., Connecting solvation shell structure to proton transport kinetics in hydrogen bonded networks via population correlation fu nctions. Phys. Rev. Lett. 2007, 99 (14), 145901. (93) Marx, D.; Chandra, A.; Tuckerman, M. E., Aqueous basic solutions: Hydroxide solvation, structural diffusion, and comparison to the hydrated proton. Chem. Rev. 2010, 110 (4), 2174 2216. (94) Vuilleumier, R.; Borgis, D., Transport and spectroscopy of the hydrated proton: A molecular dynamics study. J. Chem. Phys. 1999, 111 (9), 4251 4266.

PAGE 81

72 (95) Cao, Z.; Peng, Y.; Yan, T.; Li, S.; Li, A.; Voth, G. A., Mechanism of fast proton transport along one dimensional water chains confined in carbon nanotubes. J. Am. Chem. Soc. 2010, 132 (33), 11395 11397. (96) Brancato, G.; Tuckerman, M. E., A polarizable multistate empirical valence bond model for proton transport in aqueous solution. J. Chem. Phys. 2005, 122 (22), 22 4507. (97) Cleland, W. W.; Kreevoy, M. M., Low barrier hydrogen bonds and enzymic catalysis. Science 1994, 264 (5167), 1887 90. (98) Tuckerman, M. E.; Marx, D.; Klein, M. L.; Parrinello, M., On the quantum nature of the shared proton in hydrogen bonds. Sci ence 1997, 275 (5301), 817 820. (99) Ceriotti, M.; Cuny, J.; Parrinello, M.; Manolopoulos, D. E., Nuclear quantum effects and hydrogen bond fluctuations in water. Proc. Natl. Acad. Sci. U.S.A. 2013, 110 (39), 15591 15596. (100) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M., Ab initio simulations of water and water ions. J. Phys. Condens. Matter 1994, 6 (23A), A93 A100. (101) Warshel, A.; Weiss, R. M., An empirical valence bond approach for comparing reactions in solutions and in enzymes. J. Am. Chem. Soc. 1980, 102 (20), 6218 6226. (102) Sagnella, D. E.; Tuckerman, M. E., An empirical valence bond model for proton transfer in water. J. Chem. Phys. 1998, 108 (5), 2073 2083. (103) Schmitt, U. W.; Voth, G. A., Multistate empirical valence bond model for proton transport in water. J. Phys. Chem. B 1998, 102 (29), 5547 5551. (104) Schmitt, U. W.; Voth, G. A., The computer simulation of proton transport in water. J. Chem. Phys. 1999, 111 (20), 9361 9381. (105) Vuilleumier, R.; Borgis, D., An extended em pirical valence bond model for describing proton transfer in H + ( H 2 O ) n clusters and liquid water. Chem. Phys. Lett. 1998, 284 (1 2), 71 77.

PAGE 82

73 (106) Lefohn, A. E.; Ovchinnikov, M.; Voth, G. A., A multistate empirical valence bond approach to a polarizable and flexible water model. J. Phys. Chem. B 2001, 105 (28), 6628 6637. (107) Kornyshev, A. A.; Kuznetsov, A. M.; Spohr, E.; Ulstrup, J., Kinetics of proton transport in water. J. Phys. Chem. B. 2003, 107 (15), 3351 3366. (108) Swanson, J. M.; Maupin, C. M.; Che n, 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 (17), 4300 4314. (109) Markovitch, O.; Chen, H.; Izvekov, S.; Paesani, F .; Voth, G. A.; Agmon, N., Special pair dance and partner selection: Elementary steps in proton transport in liquid water. J. Phys. Chem. B 2008, 112 (31), 9456 9466. (110) Wu, Y.; Chen, H.; Wang, F.; Paesani, F.; Voth, G. A., An improved multistate empiri cal valence bond model for aqueous proton solvation and transport. J. Phys. Chem. B 2008, 112 (2), 467 482. (111) Park, K.; Lin, W.; Paesani, F., Fast and slow proton transfer in ice: The role of the quasi liquid layer and hydrogen bond network. J. Phys. C hem. B. 2014, 118 (28), 8081 8089. (112) Cao, Z.; Kumar, R.; Peng, Y.; Voth, G. A., Proton transport under external applied voltage. J. Phys. Chem. B 2014, 118 (28), 8090 8098. (113) Lee, S.; Swanson, J. M.; Voth, G. A., Multiscale simulations reveal key a spects of the proton transport mechanism in the CLC ec 1 antiporter. Biophys. J. 2016, 110 (6), 1334 1345. (114) Liang, R.; Swanson, J. M. J.; Madsen, J. J.; Hong, M.; DeGrado, W. F.; Voth, G. A., Acid activation mechanism of the influenza A M 2 proton chann el. Proc. Natl. Acad. Sci. U.S.A. 2016, 113 (45), E6955 E6964. (115) Billeter, S. R.; van Gunsteren, W. F., Protonizable water model for quantum dynamical simulations. J. Phys. Chem. A 1998, 102 (24), 4669 4678. (116) David, C. W., A variable charge centra l force model for water and its ionic dissociation products. J. Chem. Phys. 1996, 104 (18), 7255 7260.

PAGE 83

74 (117) Fogarty, J. C.; Aktulga, H. M.; Grama, A. Y.; van Duin, A. C.; Pandit, S. A., A reactive molecular dynamics simulation of the silica water interfac e. J. Chem. Phys. 2010, 132 (17), 174704. (118) Halley, J. W.; Rustad, J. R.; Rahman, A., A polarizable, dissociating, molecular dynamics model for liquid water. J. Chem. Phys. 1993, 98 (5), 4110 4119. (119) Kale, S.; Herzfeld, J.; Dai, S.; Blank, M., Lewi s inspired representation of dissociable water in clusters and grotthuss chains. J. Biol. Phys. 2012, 38 (1), 49 59. (120) Lussetti, E.; Pastore, G.; Smargiassi, E., A fully polarizable and dissociable potential for water. Chem. Phys. Lett. 2003, 381 (3 4) 287 291. (121) Mahadevan, T. S.; Garofalini, S. H., Dissociative water potential for molecular dynamics simulations. J. Phys. Chem. B 2007, 111 (30), 8919 8927. (122) Selvan, M. E.; Keffer, D. J.; Cui, S.; Paddison, S. J., A reactive molecular dynamics a lgorithm for proton transport in aqueous systems. J. Phys. Chem. C 2010, 114 (27), 11965 11976. (123) Stillinger, F. H.; David, C. W., Polarization model for water and its ionic dissociation products. J. Chem. Phys. 1978, 69 (4), 1473 1484. (124) Wolf, M. G.; Groenhof, G., Explicit proton transfer in classical molecular dynamics simulations. J. Comput. Chem. 2014, 35 (8), 657 671. (125) Arillo Flores, O. I.; Ruiz Lpez, M. F.; Bernal Uruchurtu, M. I., Can semi empirical models describe HC l dissociation in w ater? Theoretical Chemistry Accounts 2007, 118 (2), 425 435. (126) Choi, T. H.; Jordan, K. D., Application of the SCC DFTB method to H + ( H 2 O ) 6 H + ( H 2 O ) 21 and H + ( H 2 O ) 22 J. Phys. Chem. B. 2010, 114 (20), 6932 6936. (127) Goyal, P.; Elstner, M.; Cui, Q., App lication of the SCC DFTB method to neutral and protonated water clusters and bulk water. J. Phys. Chem. B. 2011, 115 (20), 6790 6805.

PAGE 84

75 (128) Goyal, P.; Qian, H. J.; Irle, S.; Lu, X.; 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. (129) Korth, M., Third generation hydrogen bonding corrections for semiempirical qm methods and force fields. J. Chem. Theory Comput. 2010, 6 (12), 3808 3816. (130) Lin, Y.; Wynveen, A.; Halley, J. W.; Curtiss, L. A.; Redfern, P. C., Self consistent tight binding model for dissociable water. J. Chem. Phys. 2012, 136 (17), 174507. (131) Maupin, C. M.; Aradi, B.; V oth, G. A., The self consistent charge density functional tight binding method applied to liquid water and the hydrated excess proton: Benchmark simulations. J. Phys. Chem. B. 2010, 114 (20), 6922 6931. (132) Riccardi, D.; Konig, P.; Prat Resina, X.; Yu, H .; 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 (50), 16302 16311. (133) Wang, S.; MacKay, L.; Lamoureux, G., Development of semiem pirical models for proton transfer reactions in water. J. Chem. Theory. Comput. 2014, 10 (8), 2881 2890. (134) 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 (6), 2672 2686. (135) 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. T heory Comput. 2012, 8 (10), 3586 3595. (136) Lin, H.; Zhang, Y.; Pezeshki, S.; Duster, A.; Wang, C. H.; Truhlar, D. G. Qmmm V ersion 2.2.7 CO University of Minnesota, Minneapolis, 2017. (137) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kal, L.; Schulten, K., Scalable molecular dynamics with NAMD J. Comput. Chem. 2005, 26 (16), 1781 1802.

PAGE 85

76 (138) Ponder, J. W.; Richards, F. M., An efficient N ewton like method for molecular mechanics energy minimiza tion of large molecules. J. Comput. Chem. 1987, 8 (7), 1016 1024. (139) Trucks, G.; Schlegel, H.; Scuseria, G.; Robb, M.; Cheeseman, J.; Montgomery Jr, J.; Vreven, T.; Kudin, K.; Burant, J.; Millam, J. Gaussian 03, R evision E 01 Gaussian Inc.: Wallingfor d, CT, 2004. (140) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery Jr., J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomper ts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, .; Foresman, J. B.; Ortiz, J V.; Cioslowski, J.; Fox, D. J. Gaussian 09 Gaussian, Inc.: Wallingford, CT, USA, 2009. (141) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Du puis, M.; Montgomery, J. A., General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14 (11), 1347 1363. (142) Neese, F., The ORCA program system. Wiley. Interdiscip. Rev. Comput. Mol. Sci. 2012, 2 (1), 73 78. (143) Thiel, W. M NDO 9 9 Max Planck Institut f r Kohlenforschung: M lheim an der Ruhr, Germany, 2012. (144) Humphrey, W.; Dalke, A.; Schulten, K., Vmd: Visual molecular dynamics. J. Mol. Graph. 1996, 14 (1), 33 8, 27 8. (145) Berman, H.; Henrick, K.; Nakamura, H., Announcing the worldwide protein data bank. Nat. Struct. Mol. Biol. 2003, 10 (12), 980 980. (146) Dagum, L.; Menon, R., Openmp: An industry standard api for shared memory programming. IEEE Comput. Sci. Eng. 1998, 5 (1), 46 55.

PAGE 86

77 (147) OpenMP, A., Open MP application pro gram interface, v. 3.0. OpenMP Architecture Review Board 2008 (148) Chang, C. H.; Long, H.; Sides, S.; Vaidhynathan, D.; Jones, W., Parallel application performance on two generations of I ntel X eon HPC platforms. Laboratory, N. R. E., Ed. Golden, CO, 2015 (149) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A., The weighted histogram analysis method for free energy calculations on biomolecules. I. The method. J. Comput. Chem. 1992, 13 (8), 1011 1021. (150) Wu, Y.; Tepper, H. L.; Vo th, G. A., Flexible simple point charge water model with improved liquid state properties. J. Chem. Phys. 2006, 124 (2), 024503. (151) Zhu, F. Q.; Schulten, K., Water and proton conduction through carbon nanotubes as models for biological channels. Biophys J. 2003, 85 (1), 236 244. (152) Han, W.; Cheng, R. C.; Maduke, M. C.; Tajkhorshid, E., Water access points and hydration pathways in CLC H + / C l transporters. Proc. Natl. Acad. Sci. U.S.A. 2014, 111 (5), 1819 1824. (153) Dellago, C.; Hummer, G., Kinetics and mechanism of proton transport across membrane nanopores. Phys. Rev. Lett. 2006, 97 (24), 245901. (154) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; 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. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wirkiewicz Kuczera, J.; Yin, D.; Karplus, M., All atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B. 1998, 102 (18), 3586 3616. (155) Knight, C.; Lindberg, G. E.; Voth, G. A., Multiscale reactive molecular dynamics. J. Chem. Phys. 2012, 137 (22), 22A525. (156) Ponder, J. W. Tinker molecular modeling package J ay P onder L ab. https://dasher.wustl.edu/tinker/ (accessed Feb. 3 2017).

PAGE 87

APPENDIX A BENCHMARKING INFORMATION A.1 Test System Specifications Table A.1 : Sys tem Specifications for Lenovo Laptop Test system ID: l aptop CPU Intel Core i7 6700HQ 2.60 Ghz Cores Physical/Virtual 4/8 GPU NVIDIA GeForce GTX 960M RAM 16 GB D isk 1 SanDisk Ultra II (X41200RL) SSD D isk 2 WD My Passport Ultra USB 3.0 HDD OS Ubuntu 16.10 Computer Type Laptop Table A.2 : System Specifications for compchemw5 Test System ID compchemw5 CPU Intel Xeon CPU E5 2687W @ 3.10 GHz Cores Physical/Virtual 16/32 RAM 64 GB GPU NVIDIA GTX 10 80 Disk 1 C rucial BX200 SSD Disk 2 Seagate Constel lation ES.2 7200 RPM HDD OS Centos 7 Computer Type Desktop Table A. 3 : System Specifications for compchemw8 Test System ID compchemw8 CPU Intel Xeon CPU E 5 2660 v2 @ 2.20 GHz Cores Physical/Virtual 10/20 RAM 64 GB GPU NVIDIA GTX 980 GPU NVIDIA Tes la K40 D isk 1 Intel SSDSC2KIO1OX6 SSD D isk 2 WDC WD2000FYYZ 01UL1B2 7200 RPM HDD OS Centos 7 Computer Type Desktop

PAGE 88

79 Table A. 4 : System Specifications for XSEDE Stampede Test System ID Stampede CPU 2x Intel Xeon E5 2680 2.7 GHz Cores Physical/Virt ual 16/32 RAM 32 GB GPU 61 Core Xeon Phi SE10P 1.1GHz GPU NVIDIA Tesla K40 Disk 1 250GB HDD (7.5 K RPM SATA) Disk 2 LUSTRE Network Fileystem OS Centos 7 Computer Type HPC A 2 Tinker/NAMD Benchmark Calculation Test system 1 described in Table A.1 was used for the following benchmark. 156 and the we found that the testgrad binary rarely used more than 1 core in our calculations Therefore, NAMD was ran with CUDA using 1 thread. and with CUDA util izing 4 threads. The runs with NAMD using 1 thread were reported in the main body of the thesis. The additional timing data is highlighted in yellow in Table A. 5 : Table A. 5 : Benchmarking Results for NAMD and Tinker Number of Atoms 4308 52233 96003 176454 NAMD time/step (s) 0.18 0. 55 0. 91 1. 57 NAMD 4 cores time/step (s) 0.16 0.57 1.0 2 1.64 Tinker time/step (s) 0.9 35.2 111.8 330.9 1 Core: Speedup with NAMD 5 64 123 211 4 Cores: Speedup with NAMD 6 62 112 195 A. 3 PPAP Benchmarking of a Water Box for Various Systems Desktop Workstation 8 These AP MD calculations were carried u sing the test system described in Table A. 3 on a water box of 1435 H 2 O and 1 H 3 O + cation The active zone has a 5.0 radius and is surrounded by a 1 .0 buffer shell. On average, t here were 13.4 active zone groups, 17.3

PAGE 89

80 buffer zone groups, and 160.7 QM calculations per step for 100 steps The AP potential was truncated at the 2 nd order. The timing data can be found in Table A. 6 Table A.6: Tests with Desktop Workstation compchemw8 with an HDD and SSD WT/S stands for wall time per step and SSE stands for strong scaling efficiency Threads SSD: WT/S (s) SSD: SSE HDD: WT/S (s) HDD: SSE 1 62.8 100.00% 65.8 100.00% 2 29.1 107.90% 28.8 114.24% 4 15.2 103.29% 16.0 102.81% 8 8.4 93. 45% 9.5 86.58% 10 7.1 88.45% 8.5 77.41% 20 5.4 58.15% 6.0 54.83% XSEDE Stampede The following described in Table A.3 The test model and parameters were the same as in the model in Table A .5 except NAMD 2.10 was used instead of 2.11 There as an average of 13. 5 active zone groups, 17. 5 buffer zone groups, and 16 3.5 QM calculations per step. The timing data can be found in Table A. 7 Table A. 7 : Tests using XSEDE Stampede with an internal SS D and LUSTRE network filesystem. The calculation on the network file system failed when 32 threads were used. WT/S stands for wall time per step and SSE stands for strong scaling efficiency Threads Internal: WT/S (s) Internal: SSE Network: WT/S (s) Netwo rk: SSE 1 74.8 100.00% 100.3 100.00% 2 36.9 101.36% 60.1 83.44% 4 22.6 82.74% 38.8 64.63% 8 9.8 95.41% 21.1 59.42% 10 8.9 84.04% 19.5 51.44% 16 8.9 52.53% 21.7 28.89% 32 11.2 20.87% N/A N/A Desktop Workstation 5 The following tests used the test s ystem described in Table A.3 on a water box of 1435 H 2 O and 1 H 3 O + cation Two active/buffer zone radii were investigated in order to

PAGE 90

81 show the decrease in SSE due to a lack of available calculations to assign to the thread team With an active zone radius of 4.0 surrounded by a 0.5 buffer shell, there was an average of 5.9 active zone groups, 3.4 buffer zone groups, and 9.3 QM calculations per step for 100 steps. With an active zone radius of 4.5 surrounded by a 1.0 buffer shell, there was an averag e of 7.5 active zone groups, 15.6 buffer zone groups, and 115.0 QM calculations per step for 100 steps. The timing data can be found in Table A.8 Table A.8: Strong scaling efficiency (SSE) data for two active/buffer zone radii configurations : one with an active zone radius of 4.0 surrounded by a 0.5 buffer shell had 9.3 QM/MM calculations per step, the other with an active zone radius of 4.5 surrounded by a 1.0 buffer shell had 115.0 QM/MM calculations per step. Threads 9.3 Calculations / Step SS E 115.0 Calculations / Step SSE 1 100.00% 100.00% 2 107.50% 99.49% 4 82.69% 93.81% 8 53.75% 84.91% 16 29.86% 70.36%

PAGE 91

82 APPENDIX B DECOMPOSED GRADIENT INCLUDING HYDROGEN B .1: Variables and Functions used in the Gradient Decomposition Equations Table B. 1: Variables and Functions used in the Gradient Decomposition Equations Symbol Description Restraint harmonic potential energy imposed on indicator I Force constant of the restraint potential Current position of indicator Targeted position of the indicator Position of j th possible acceptor oxy gen A j Position of donor oxygen D Position of m th hydrogen bonded to donor D Donor hydrogen vector of m th hydrogen Donor acceptor vector of j th possible acceptor Equilibrium bond dista nce parameter Possible acceptor enlisting cutoff parameter Ratio of vector projection for m th H and j th A Cutoff parameter for mj Weight function for mj Normalization factor for the weight function B .2 : Decomposition of the Potential onto Gradient of with respect to : Gradient of with respect to :

PAGE 92

83 Gradient of with respect to : Gradient of with respect to : Gradient of with respect to B .3 Decomposition of the Potential onto Grad ient of with respect to : Gradient of with respect to :

PAGE 93

84 Gradient of with respect to : Gradient of with respect to : Gradient of with respect to : B .4: Decomposition of the Potential onto : Gradient of with respect to : Gradient of with respect to :

PAGE 94

85 Gradient of of with respect to : Gradient of with respect to : Gradients of with respect to : B.5: Gradient Graph s When the gradient was decomposed over all atoms, the gradient of the donor and acceptor atom was in the opposite direction of the gradient on the proton These gradients could be quite large, but cancelled each other out when summed to yield the total gradi ent on the indicator. The results of these exploratory calculations are shown in Fig B.1

PAGE 95

86 Fig B.1 : Decomposed gradients for Zundel cation. ( A ) I ndicator distance from origin and ( B ) gradient on indicator. G radients on ( B ) donor and acceptor, and ( D ) indicator. The maximum gradient on the proton is around 130 kcal mol 1 and the maximum gradient on the acceptor and donor is around 55 and 65 kcal mol 1 respectively.

PAGE 96

87 APPENDIX C Additional Details for MD Simulations C.1 NVE Collision Simulations Only one half of the graphs were reported in the main text. The full sets of gr aphs are found below: Fig C.1: T otal energy of system for all simulations of the Zundel ( A C ) and Eigen ( B D ) cations with different force constants and initial H 2 O velocities of ( A B ) 0.001 fs 1 and ( C D ) 0.01 fs 1 The blue, green, red, and cyan lines correspond to force consta nts of 0 (no restraint), 5, 10, and 20 kcal mol 1 2 respectively.

PAGE 97

88 Fig C.2 : RMSD ( ) of indicator distance from origin of potential for Zundel ( rows 1, 3 ) and Eigen ( rows 2, 4 ) cations. Red, blue, green, and cyan correspond to force constants of 0 (no re straint), 5, 10, and 20 kcal mol 1 2 respectively. The initial velocity of the H 2 O was 0.001 fs 1 (rows 1,2) and 0.01 fs 1 (rows 3,4).

PAGE 98

89 C. 2 MD simulations of Water Sphere Four simulations of a water sphere containing 1482 H 2 O molecules with one H 3 O + at the center of the sphere were conducted T he force constant for indicator restraints were set at 5 kcal mol 2 Fo ur additional simulations were conducted in the absence of indicator restrai nts for comparison The simulations used a droplet model with the inner 12 sphere described at the OM3n level and the outer 10 shell described at the MM level with the SPC /Fw potential For the MM molecules, restraints of 20 kcal mol 1 applied to every oxygen atom. Other general simulation details can be found in Table C.1 Table C.1 General MD parameters for water sphere simulations Atoms 4450 QM only atoms 796 Ste p length 0.25 fs Equilibration time 0.5 ps Production time 10.0 ps Temperature 300 K ES Charge Cutoff 14.0 ES Taper 13.0 Snapshot Save Frequency 2 fs Density 1.0 g ml 1 Periodic Boundaries None Average wall time per step ~ 10.2 s Simulation Sys tem c ompchemw5 MM Code NAMD 2.10 QM Code MNDO C. 3 MM System used for Nanotube Simulations The nanotube system was constructed and equilibrated via the following procedure. The epsilon parameter for the nanotube carbon, CA was changed to 80% of its orig inal value to 0.056 to increase hydrophobicity of the tube. This was an error which would be later corrected as the depth of the potential well between the oxygen water atom and carbon atom should have been scaled to 80% of its original value, not the ep silon parameter itself.

PAGE 99

90 The simulation was equilibrated for 200 ps with 1 fs time steps at the NPT level. The x and y dimensions were fixed, and unit cell and position changes due to pressure were only allowed to change in the z direction. Four snapshots w ere taken at 200, 190, 180, and 170 ps from the end of the equilibration trajectory and used as starting geometries for the next phase of the equilibration. Each of the four geometries were then equilibrated for an additional 5 ns at NPT with the x and y d imensions fixed. Finally, the simulations were equilibrated for another 5 ns at NVT The well depth error was realized, and the epsilon parameter for CA was then changed to 0.0448. The final snapshot of each of the previous simulations was then equilibrat ed for an additional 5 ns at NPT with the x and y dimensions fixed. The simulations were then equilibrated for another 5 ns at NVT C.4 MD Simulations of Restrained Indicator at Pore Opening From each MM trajectory, a snapshot with zero, one, two, three, o r four waters in the tube was obtained starting from the end of the trajectory, yielding a total of 20 starting geometries. All waters within a 14 radius of the middle of the pore opening were included in the QM subsystem. The number QM atoms for these s imulations ranges from 400 to 475, wi th an average of 437 QM atoms. Simulation lengths are listed in Table C. 2 Also, Fig 5.3 in the main text is taken from the simulation named res_0w_2

PAGE 100

91 Table C.2 : List of Simulations of Restrained Indicator at Pore O pening Name Restraint Initial Waters Total Time (ps) res_0w_1 Yes 0 2.66675 res_0w_2 Yes 0 10.50000 res_0w_3 Yes 0 10.50000 res_0w_4 Yes 0 3.5615 res_1w_1 Yes 1 2.39225 res_1w_2 Yes 1 4.2910 res_1w_3 Yes 1 6.1550 res_1w_4 Yes 1 10.50000 res_2w_1 Yes 2 10.50000 res_2w_2 Yes 2 4.34500 res_2w_3 Yes 2 10.50000 res_2w_4 Yes 2 6.56350 res_3w_1 Yes 3 10.50000 res_3w_2 Yes 3 10.50000 res_3w_3 Yes 3 10.50000 res_3w_4 Yes 3 8.22300 res_4w_1 Yes 4 10.50000 res_4w_2 Yes 4 10.50000 res_4w_3 Yes 4 10. 50000 res_4w_4 Yes 4 1.96275 nor_0w_1 No 0 10.50000 nor_0w_2 No 0 10.50000 nor_0w_3 No 0 0.00000 nor_0w_4 No 0 8.94650 nor_1w_1 No 1 10.50000 nor_1w_2 No 1 10.50000 nor_1w_3 No 1 10.50000 nor_1w_4 No 1 10.50000 nor_2w_1 No 2 10.50000 nor_2w_2 No 2 10.50000 nor_2w_3 No 2 10.50000 nor_2w_4 No 2 10.50000 nor_3w_1 No 3 10.50000 nor_3w_2 No 3 7.79575 nor_3w_3 No 3 8.54700 nor_3w_4 No 3 10.50000 nor_4w_1 No 4 10.50000 nor_4w_2 No 4 10.50000 nor_4w_3 No 4 10.50000 nor_4w_4 No 4 3.64650

PAGE 101

92 C 5 MD Simulations of Restraints at Different Locations One snapshot was taken from the end of a n MM simulation and a water close to 3.5 A above the opening of the tube was protonated. A harmonic restraint with a force constant of 20 kcal mol 1 2 was applied to restrain the hydronium oxygen to 3.5 above the center of the tube. The system was then equilibrated at NVT for 50 ps. Snapshots were taken from the steps at 50 ps, 47.5 ps, 45 ps, and 42.5 ps and used as input for QM/MM simulations w h ere the indicato r was restrained 3.5 above the center of the tube opening. The last snapshot of th is simulation was used as the initial geometry for a simulation where the hydronium oxygen atom was restrained at 5 above the pore opening and a simulation where the hyd ronium oxygen atom was restrained 3.5 above the pore opening and 2 in the positive x direction from the center of the pore. The last step of the 5 simulation was used to generate an initial geometry for a simulation where the O was restrained 6.5 a bove the center of the pore opening and for a simulation where the O was restrained 5 above and 2 in the x direction away from the center of the pore. Simulation names and lengths are listed in Table C.3

PAGE 102

93 Table C.3: List of Simulations for Restrain ts at Different Locations Name Distance Above Tube ( ) Distance to Side of Tube ( ) Simulation Time ( ps ) 35_1 3.5 0.0 6.09300 35_2 3.5 0.0 10.50000 35_3 3.5 0.0 4.00825 35_4 3.5 0.0 10.50000 35x_1 3.5 2.5 7.55000 35x_2 3.5 2.5 10.50000 35x_3 3.5 2. 5 10.50000 35x_4 3.5 2.5 6.72950 50_1 5.0 0.0 3.76350 50_2 5.0 0.0 7.66200 50_3 5.0 0.0 10.50000 50_4 5.0 0.0 10.50000 50x_1 5.0 2.5 10.50000 50x_2 5.0 2.5 10.50000 50x_3 5.0 2.5 10.50000 50x_4 5.0 2.5 4.37750 65_1 6.5 0.0 10.50000 65_2 6.5 0.0 10.50000 65_3 6.5 0.0 7.15500 65_4 6.5 0.0 8.25275 C.6 Umbrella Sampling of Proton Transfer in CNT One snapshot was taken from the end of an MM simulation and 14 waters were manually added to the tube. The system was minimized for 50 steps. Then a wate r molecule inside the tube was protonated, and the model was minimized for 5 steps. The simulations were ran at the QM/MM level with a truncated PS for 10 ps, and snapshots from the simulation were used to generate the starting geometries for 18 umbrella w indows ( Table C. 4 )

PAGE 103

94 Table C. 4 : List of umbrella windows Window ID Origin ( ) Force Constant s ( kcal mol 1 1 ) Simulation Length ( ps ) x y z x y z 1_k5 0.000 0.000 14.122 0.0 0.0 5.0 10.69450 1_k2.5 0.000 0.000 14.122 0.0 0.0 2.5 10.50425 2_k5 0.000 0.000 15.131 0.0 0.0 5.0 10.72925 2_k2.5 0.000 0.000 15.131 0.0 0.0 2.5 10.56050 3 0.000 0.000 16.140 0.0 0.0 5.0 10.95000 4 0.000 0.000 17.148 0.0 0.0 5.0 10.95100 5 0.000 0.000 18.157 0.0 0.0 5.0 10.63175 6 0.000 0.000 19.166 0.0 0.0 5.0 10.62225 7 0.000 0.000 20.175 0.0 0.0 5.0 10.50175 8 0.000 0.000 21.814 0.0 0.0 5.0 10.29750 9 0.000 0.000 22.192 0.0 0.0 5.0 10.49675 10 0.000 0.000 23.201 0.0 0.0 5.0 10.45875 11 0.000 0.000 24.210 0.0 0.0 5.0 9.01500 12 20.995 20.560 25.219 2.0 2.0 7.5 9.05 325 13 20.995 20.560 26.227 2.0 2.0 7.5 12.55225 14 20.995 20.560 27.236 2.0 2.0 7.5 8.22800 15 20.995 20.560 28.245 2.0 2.0 7.5 7.61000 16 20.995 20.560 30.197 2.0 2.0 7.5 13.68400 17 20.995 20.560 25.219 2.0 2.0 7.5 4.76150 18 20.995 20.560 25.201 2.0 2.0 7.5 4.99700 As mentioned in the main body of the text, no QM waters were restrained in the first the first set of umbrella sampling simulations. This led to vigorous recruitment of water into the tube The PMF before restraints and the final PMF are shown side by side for comparison ( Fig C. 3 ). The PMF shown in the main text ( Fig 5.7 ) was truncated at the z coordinate of the furthest window away from the center of the nanotube ( z = 30.197). The full sampling and PMF are included for completeness ( F ig C. 4 ), but there are certainly sampling errors past the furthest umbrella window origin.

PAGE 104

95 Fig C. 4 : ( A ) Histogram of samples along the umbrella reaction coordinate. Different colors indicate different umbr ella window origins. ( B ) PMF for proton transfer from the middle of the nanotube towards the bulk. Note that the PMF past z = 30.197 may be inaccurate because this is the last umbrella window origin Fig C.3 : Final PMF obtained from umbrella sampling (blue) compared with PMF from Umbrella sampling before QM region restraints.