| 1 |
< |
\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND} |
| 1 |
> |
\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND} |
| 2 |
|
|
| 3 |
< |
\section{\label{sec:IntroIntegrate}Rigid Body Motion} |
| 3 |
> |
The following dissertation presents the primary aspects of the |
| 4 |
> |
research I have performed and been involved with over the last several |
| 5 |
> |
years. Rather than presenting the topics in a chronological fashion, |
| 6 |
> |
they are arranged to form a series where the later topics apply and |
| 7 |
> |
extend the findings of the former topics. This layout does lead to |
| 8 |
> |
occasional situations where knowledge gleaned from earlier chapters |
| 9 |
> |
(particularly chapter \ref{chap:electrostatics}) is not applied |
| 10 |
> |
outright in the later work; however, I feel that this organization is |
| 11 |
> |
more instructive and provides a more cohesive progression of research |
| 12 |
> |
efforts. |
| 13 |
> |
|
| 14 |
> |
This introductory chapter gives a general overview of molecular |
| 15 |
> |
simulations, with particular emphasis on considerations that need to |
| 16 |
> |
be made in order to apply the technique properly. This leads quite |
| 17 |
> |
naturally into chapter \ref{chap:electrostatics}, where we investigate |
| 18 |
> |
correction techniques for proper handling of long-ranged electrostatic |
| 19 |
> |
interactions in molecular simulations. In particular we develop and |
| 20 |
> |
evaluate some new simple pairwise methods. These techniques make an |
| 21 |
> |
appearance in the later chapters, as they are applied to specific |
| 22 |
> |
systems of interest, showing how it they can improve the quality of |
| 23 |
> |
various molecular simulations. |
| 24 |
> |
|
| 25 |
> |
In chapter \ref{chap:water}, we focus on simple water models, |
| 26 |
> |
specifically the single-point soft sticky dipole (SSD) family of water |
| 27 |
> |
models. These single-point models are more efficient than the common |
| 28 |
> |
multi-point partial charge models and, in many cases, better capture |
| 29 |
> |
the dynamic properties of water. We discuss improvements to these |
| 30 |
> |
models in regards to long-range electrostatic corrections and show |
| 31 |
> |
that these models can work well with the techniques discussed in |
| 32 |
> |
chapter \ref{chap:electrostatics}. By investigating and improving |
| 33 |
> |
simple water models, we are extending the limits of computational |
| 34 |
> |
efficiency for systems that depend heavily on water calculations. |
| 35 |
> |
|
| 36 |
> |
In chapter \ref{chap:ice}, we study a unique polymorph of ice that we |
| 37 |
> |
discovered while performing water simulations with the fast simple |
| 38 |
> |
water models discussed in the previous chapter. This form of ice, |
| 39 |
> |
which we call ``imaginary ice'' (Ice-$i$), has a low-density structure |
| 40 |
> |
which is different from any known polymorph characterized in either |
| 41 |
> |
experiment or other simulations. In this work, we perform a free |
| 42 |
> |
energy analysis and see that this structure is in fact the |
| 43 |
> |
thermodynamically preferred form of ice for both the single-point and |
| 44 |
> |
commonly used multi-point water models under the chosen simulation |
| 45 |
> |
conditions. We also consider electrostatic corrections, again |
| 46 |
> |
including the techniques discussed in chapter |
| 47 |
> |
\ref{chap:electrostatics}, to see how they influence the free energy |
| 48 |
> |
results. This work, to some degree, addresses the appropriateness of |
| 49 |
> |
using these simplistic water models outside of the conditions for |
| 50 |
> |
which they were developed. |
| 51 |
> |
|
| 52 |
> |
Finally, in chapter \ref{chap:conclusion}, we summarize the work |
| 53 |
> |
presented in the previous chapters and connect ideas together with |
| 54 |
> |
some final comments. The supporting information follows in the |
| 55 |
> |
appendix, and it gives a more detailed look at systems discussed in |
| 56 |
> |
chapter \ref{chap:electrostatics}. |
| 57 |
> |
|
| 58 |
> |
\section{On Molecular Simulation} |
| 59 |
> |
|
| 60 |
> |
In order to advance our understanding of natural chemical and physical |
| 61 |
> |
processes, researchers develop explanations for events observed |
| 62 |
> |
experimentally. These hypotheses, supported by a body of corroborating |
| 63 |
> |
observations, can develop into accepted theories for how these |
| 64 |
> |
processes occur. This validation process, as well as testing the |
| 65 |
> |
limits of these theories, is essential in developing a firm foundation |
| 66 |
> |
for our knowledge. Developing and validating theories involving |
| 67 |
> |
molecular scale systems is particularly difficult because their size |
| 68 |
> |
and often fast motions make them difficult to observe |
| 69 |
> |
experimentally. One useful tool for addressing these difficulties is |
| 70 |
> |
computer simulation. |
| 71 |
> |
|
| 72 |
> |
In computer simulations, we can develop models from either the theory |
| 73 |
> |
or our experimental knowledge and then test them in controlled |
| 74 |
> |
environments. Work done in this manner allows us to further refine |
| 75 |
> |
theories, more accurately represent what is happening in experimental |
| 76 |
> |
observations, and even make predictions about what one will see in |
| 77 |
> |
future experiments. Thus, computer simulations of molecular systems |
| 78 |
> |
act as a bridge between theory and experiment. |
| 79 |
> |
|
| 80 |
> |
Depending on the system of interest, there are a variety of different |
| 81 |
> |
computational techniques that can used to test and gather information |
| 82 |
> |
from the developed models. In the study of classical systems, the two |
| 83 |
> |
most commonly used techniques are Monte Carlo and molecular |
| 84 |
> |
dynamics. Both of these methods operate by calculating interactions |
| 85 |
> |
between particles of our model systems; however, the progression of |
| 86 |
> |
the simulation under the different techniques is vastly |
| 87 |
> |
different. Monte Carlo operates through random configuration changes |
| 88 |
> |
that that follow rules adhering to a specific statistical mechanics |
| 89 |
> |
ensemble, while molecular dynamics is chiefly concerned with solving |
| 90 |
> |
the classic equations of motion to move between configurations within |
| 91 |
> |
an ensemble. Thermodynamic properties can be calculated from both |
| 92 |
> |
techniques; but because of the random nature of Monte Carlo, only |
| 93 |
> |
molecular dynamics can be used to investigate dynamical |
| 94 |
> |
quantities. The research presented in the following chapters utilized |
| 95 |
> |
molecular dynamics near exclusively, so we will present a general |
| 96 |
> |
introduction to molecular dynamics. There are several resources |
| 97 |
> |
available for those desiring a more in-depth presentation of either of |
| 98 |
> |
these techniques.\cite{Allen87,Frenkel02,Leach01} |
| 99 |
> |
|
| 100 |
> |
\section{\label{sec:MolecularDynamics}Molecular Dynamics} |
| 101 |
> |
|
| 102 |
> |
As stated above, in molecular dynamics we focus on evolving |
| 103 |
> |
configurations of molecules over time. In order to use this as a tool |
| 104 |
> |
for understanding experiments and testing theories, we want the |
| 105 |
> |
configuration to evolve in a manner that mimics real molecular |
| 106 |
> |
systems. To do this, we start with clarifying what we know about a |
| 107 |
> |
given configuration of particles at time $t_1$, basically that each |
| 108 |
> |
particle in the configuration has a position ($\mathbf{q}$) and velocity |
| 109 |
> |
($\dot{\mathbf{q}}$). We now want to know what the configuration will be at |
| 110 |
> |
time $t_2$. To find out, we need the classical equations of |
| 111 |
> |
motion, and one useful formulation of them is the Lagrangian form. |
| 112 |
> |
|
| 113 |
> |
The Lagrangian ($L$) is a function of the positions and velocities that |
| 114 |
> |
takes the form, |
| 115 |
> |
\begin{equation} |
| 116 |
> |
L = K - V, |
| 117 |
> |
\label{eq:lagrangian} |
| 118 |
> |
\end{equation} |
| 119 |
> |
where $K$ is the kinetic energy and $V$ is the potential energy. We |
| 120 |
> |
can use Hamilton's principle, which states that the integral of the |
| 121 |
> |
Lagrangian over time has a stationary value for the correct path of |
| 122 |
> |
motion, to say that the variation of the integral of the Lagrangian |
| 123 |
> |
over time is zero,\cite{Tolman38} |
| 124 |
> |
\begin{equation} |
| 125 |
> |
\delta\int_{t_1}^{t_2}L(\mathbf{q},\dot{\mathbf{q}})dt = 0. |
| 126 |
> |
\end{equation} |
| 127 |
> |
The variation can be transferred to the variables that make up the |
| 128 |
> |
Lagrangian, |
| 129 |
> |
\begin{equation} |
| 130 |
> |
\int_{t_1}^{t_2}\sum_{i=1}^{3N}\left( |
| 131 |
> |
\frac{\partial L}{\partial \mathbf{q}_i}\delta \mathbf{q}_i |
| 132 |
> |
+ \frac{\partial L}{\partial \dot{\mathbf{q}}_i}\delta |
| 133 |
> |
\dot{\mathbf{q}}_i\right)dt = 0. |
| 134 |
> |
\end{equation} |
| 135 |
> |
Using the fact that $\dot{\mathbf{q}}$ is the derivative of |
| 136 |
> |
$\mathbf{q}$ with respect to time and integrating the second partial |
| 137 |
> |
derivative in the parenthesis by parts, this equation simplifies to |
| 138 |
> |
\begin{equation} |
| 139 |
> |
\int_{t_1}^{t_2}\sum_{i=1}^{3N}\left( |
| 140 |
> |
\frac{d}{dt}\frac{\partial L}{\partial \dot{\mathbf{q}}_i} |
| 141 |
> |
- \frac{\partial L}{\partial \mathbf{q}_i}\right) |
| 142 |
> |
\delta {\mathbf{q}}_i dt = 0, |
| 143 |
> |
\end{equation} |
| 144 |
> |
and since each variable is independent, we can separate the |
| 145 |
> |
contribution from each of the variables: |
| 146 |
> |
\begin{equation} |
| 147 |
> |
\frac{d}{dt}\frac{\partial L}{\partial \dot{\mathbf{q}}_i} |
| 148 |
> |
- \frac{\partial L}{\partial \mathbf{q}_i} = 0 |
| 149 |
> |
\quad\quad(i=1,2,\dots,N). |
| 150 |
> |
\label{eq:formulation} |
| 151 |
> |
\end{equation} |
| 152 |
> |
To obtain the classical equations of motion for the particles, we can |
| 153 |
> |
substitute equation (\ref{eq:lagrangian}) into the above equation with |
| 154 |
> |
$m\dot{\mathbf{r}}^2/2$ for the kinetic energy, giving |
| 155 |
> |
\begin{equation} |
| 156 |
> |
\frac{d}{dt}(m\dot{\mathbf{r}})+\frac{dV}{d\mathbf{r}}=0, |
| 157 |
> |
\end{equation} |
| 158 |
> |
or more recognizably, |
| 159 |
> |
\begin{equation} |
| 160 |
> |
\mathbf{f} = m\mathbf{a}, |
| 161 |
> |
\end{equation} |
| 162 |
> |
where $\mathbf{f} = -dV/d\mathbf{r}$ and $\mathbf{a} = |
| 163 |
> |
d^2\mathbf{r}/dt^2$. This Lagrangian formulation shown in equation |
| 164 |
> |
(\ref{eq:formulation}) is generalized, and it can be used to determine |
| 165 |
> |
equations of motion in forms outside of the typical Cartesian case |
| 166 |
> |
shown here. |
| 167 |
> |
|
| 168 |
> |
\subsection{\label{sec:Verlet}Verlet Integration} |
| 169 |
> |
|
| 170 |
> |
In order to perform molecular dynamics, we need an algorithm that |
| 171 |
> |
integrates the equations of motion described above. Ideal algorithms |
| 172 |
> |
are both simple in implementation and highly accurate. There have been |
| 173 |
> |
a large number of algorithms developed for this purpose; however, for |
| 174 |
> |
reasons discussed below, we are going to focus on the Verlet class of |
| 175 |
> |
integrators.\cite{Gear66,Beeman76,Berendsen86,Allen87,Verlet67,Swope82} |
| 176 |
> |
|
| 177 |
> |
In Verlet's original study of computer ``experiments'', he directly |
| 178 |
> |
integrated the Newtonian second order differential equation of motion, |
| 179 |
> |
\begin{equation} |
| 180 |
> |
m\frac{d^2\mathbf{r}_i}{dt^2} = \sum_{j\ne i}\mathbf{f}(r_{ij}), |
| 181 |
> |
\end{equation} |
| 182 |
> |
with the following algorithm: |
| 183 |
> |
\begin{equation} |
| 184 |
> |
\mathbf{r}_i(t+\delta t) = -\mathbf{r}_i(t-\delta t) + 2\mathbf{r}_i(t) |
| 185 |
> |
+ \sum_{j\ne i}\mathbf{f}(r_{ij}(t))\delta t^2, |
| 186 |
> |
\label{eq:verlet} |
| 187 |
> |
\end{equation} |
| 188 |
> |
where $\delta t$ is the time step of integration.\cite{Verlet67} It is |
| 189 |
> |
interesting to note that equation (\ref{eq:verlet}) does not include |
| 190 |
> |
velocities, and this makes sense since they are not present in the |
| 191 |
> |
differential equation. The velocities are necessary for the |
| 192 |
> |
calculation of the kinetic energy and can be calculated at each time |
| 193 |
> |
step with the equation: |
| 194 |
> |
\begin{equation} |
| 195 |
> |
\mathbf{v}_i(t) = \frac{\mathbf{r}_i(t+\delta t)-\mathbf{r}_i(t-\delta t)} |
| 196 |
> |
{2\delta t}. |
| 197 |
> |
\end{equation} |
| 198 |
> |
|
| 199 |
> |
Like the equation of motion it solves, the Verlet algorithm has the |
| 200 |
> |
beneficial property of being time-reversible, meaning that you can |
| 201 |
> |
integrate the configuration forward and then backward and end up at |
| 202 |
> |
the original configuration. Some other methods for integration, like |
| 203 |
> |
predictor-corrector methods, lack this property in that they require |
| 204 |
> |
higher order information that is discarded after integrating |
| 205 |
> |
steps. Another interesting property of this algorithm is that it is |
| 206 |
> |
symplectic, meaning that it preserves area in phase-space. Symplectic |
| 207 |
> |
algorithms keep the system evolving in the region of phase-space |
| 208 |
> |
dictated by the ensemble and enjoy a greater degree of energy |
| 209 |
> |
conservation.\cite{Frenkel02} |
| 210 |
> |
|
| 211 |
> |
While the error in the positions calculated using the Verlet algorithm |
| 212 |
> |
is small ($\mathcal{O}(\delta t^4)$), the error in the velocities is |
| 213 |
> |
substantially larger ($\mathcal{O}(\delta t^2)$).\cite{Allen87} Swope |
| 214 |
> |
{\it et al.} developed a corrected version of this algorithm, the |
| 215 |
> |
`velocity Verlet' algorithm, which improves the error of the velocity |
| 216 |
> |
calculation and thus the energy conservation.\cite{Swope82} This |
| 217 |
> |
algorithm involves a full step of the positions, |
| 218 |
> |
\begin{equation} |
| 219 |
> |
\mathbf{r}(t+\delta t) = \mathbf{r}(t) + \delta t\mathbf{v}(t) |
| 220 |
> |
+ \frac{1}{2}\delta t^2\mathbf{a}(t), |
| 221 |
> |
\end{equation} |
| 222 |
> |
and a half step of the velocities, |
| 223 |
> |
\begin{equation} |
| 224 |
> |
\mathbf{v}\left(t+\frac{1}{2}\delta t\right) = \mathbf{v}(t) |
| 225 |
> |
+ \frac{1}{2}\delta t\mathbf{a}(t). |
| 226 |
> |
\end{equation} |
| 227 |
> |
After forces are calculated at the new positions, the velocities can |
| 228 |
> |
be updated to a full step, |
| 229 |
> |
\begin{equation} |
| 230 |
> |
\mathbf{v}(t+\delta t) = \mathbf{v}\left(t+\frac{1}{2}\delta t\right) |
| 231 |
> |
+ \frac{1}{2}\delta t\mathbf{a}(t+\delta t). |
| 232 |
> |
\end{equation} |
| 233 |
> |
By integrating in this manner, the error in the velocities reduces to |
| 234 |
> |
$\mathcal{O}(\delta t^3)$. It should be noted that the error in the |
| 235 |
> |
positions increases to $\mathcal{O}(\delta t^3)$, but the resulting |
| 236 |
> |
improvement in the energies coupled with the maintained simplicity, |
| 237 |
> |
time-reversibility, and symplectic nature make it an improvement over |
| 238 |
> |
the original form. |
| 239 |
> |
|
| 240 |
> |
\subsection{\label{sec:IntroIntegrate}Rigid Body Motion} |
| 241 |
> |
|
| 242 |
> |
Rigid bodies are non-spherical particles or collections of particles |
| 243 |
> |
that have a constant internal potential and move |
| 244 |
> |
collectively.\cite{Goldstein01} To move these particles, the |
| 245 |
> |
translational and rotational motion can be integrated |
| 246 |
> |
independently. Discounting iterative constraint procedures like {\sc |
| 247 |
> |
shake} and {\sc rattle} for approximating rigid body dynamics, rigid |
| 248 |
> |
bodies are not included in most simulation packages because of the |
| 249 |
> |
algorithmic complexity involved in propagating the orientational |
| 250 |
> |
degrees of freedom.\cite{Ryckaert77,Andersen83,Krautler01} Integrators |
| 251 |
> |
which propagate orientational motion with an acceptable level of |
| 252 |
> |
energy conservation for molecular dynamics are relatively new |
| 253 |
> |
inventions. |
| 254 |
> |
|
| 255 |
> |
Moving a rigid body involves determination of both the force and |
| 256 |
> |
torque applied by the surroundings, which directly affect the |
| 257 |
> |
translational and rotational motion in turn. In order to accumulate |
| 258 |
> |
the total force on a rigid body, the external forces and torques must |
| 259 |
> |
first be calculated for all the internal particles. The total force on |
| 260 |
> |
the rigid body is simply the sum of these external forces. |
| 261 |
> |
Accumulation of the total torque on the rigid body is more complex |
| 262 |
> |
than the force because the torque is applied to the center of mass of |
| 263 |
> |
the rigid body. The space-fixed torque on rigid body $i$ is |
| 264 |
> |
\begin{equation} |
| 265 |
> |
\boldsymbol{\tau}_i= |
| 266 |
> |
\sum_{a}\biggl[(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia} |
| 267 |
> |
+ \boldsymbol{\tau}_{ia}\biggr], |
| 268 |
> |
\label{eq:torqueAccumulate} |
| 269 |
> |
\end{equation} |
| 270 |
> |
where $\boldsymbol{\tau}_i$ and $\mathbf{r}_i$ are the torque on and |
| 271 |
> |
position of the center of mass respectively, while $\mathbf{f}_{ia}$, |
| 272 |
> |
$\mathbf{r}_{ia}$, and $\boldsymbol{\tau}_{ia}$ are the force on, |
| 273 |
> |
position of, and torque on the component particles of the rigid body. |
| 274 |
> |
|
| 275 |
> |
The summation of the total torque is done in the body fixed axis. In |
| 276 |
> |
order to move between the space fixed and body fixed coordinate axes, |
| 277 |
> |
parameters describing the orientation must be maintained for each |
| 278 |
> |
rigid body. At a minimum, the rotation matrix ($\mathsf{A}$) can be |
| 279 |
> |
described by the three Euler angles ($\phi, \theta,$ and $\psi$), |
| 280 |
> |
where the elements of $\mathsf{A}$ are composed of trigonometric |
| 281 |
> |
operations involving $\phi, \theta,$ and $\psi$.\cite{Goldstein01} |
| 282 |
> |
Direct propagation of the Euler angles has a known $1/\sin\theta$ |
| 283 |
> |
divergence in the equations of motion for $\phi$ and $\psi$, leading |
| 284 |
> |
to numerical instabilities any time one of the directional atoms or |
| 285 |
> |
rigid bodies has an orientation near $\theta=0$ or |
| 286 |
> |
$\theta=\pi$.\cite{Allen87} One of the most practical ways to avoid |
| 287 |
> |
this ``gimbal point'' is to switch to another angular set defining the |
| 288 |
> |
orientation of the rigid body near this point.\cite{Barojas73} This |
| 289 |
> |
procedure results in additional book-keeping and increased algorithm |
| 290 |
> |
complexity. In the search for more elegant alternative methods, Evans |
| 291 |
> |
proposed the use of quaternions to describe and propagate |
| 292 |
> |
orientational motion.\cite{Evans77} |
| 293 |
> |
|
| 294 |
> |
The quaternion method for integration involves a four dimensional |
| 295 |
> |
representation of the orientation of a rigid |
| 296 |
> |
body.\cite{Evans77,Evans77b,Allen87} Thus, the elements of |
| 297 |
> |
$\mathsf{A}$ can be expressed as arithmetic operations involving the |
| 298 |
> |
four quaternions ($q_0, q_1, q_2,$ and $q_3$), |
| 299 |
> |
\begin{equation} |
| 300 |
> |
\mathsf{A} = \left( \begin{array}{l@{\quad}l@{\quad}l} |
| 301 |
> |
q_0^2+q_1^2-q_2^2-q_3^2 & 2(q_1q_2+q_0q_3) & 2(q_1q_3-q_0q_2) \\ |
| 302 |
> |
2(q_1q_2-q_0q_3) & q_0^2-q_1^2+q_2^2-q_3^2 & 2(q_2q_3+q_0q_1) \\ |
| 303 |
> |
2(q_1q_3+q_0q_2) & 2(q_2q_3-q_0q_1) & q_0^2-q_1^2-q_2^2+q_3^2 \\ |
| 304 |
> |
\end{array}\right). |
| 305 |
> |
\end{equation} |
| 306 |
> |
Integration of the equations of motion involves a series of arithmetic |
| 307 |
> |
operations involving the quaternions and angular momenta and leads to |
| 308 |
> |
performance enhancements over Euler angles, particularly for very |
| 309 |
> |
small systems.\cite{Evans77} This integration method works well for |
| 310 |
> |
propagating orientational motion in the canonical ensemble ($NVT$); |
| 311 |
> |
however, energy conservation concerns arise when using the simple |
| 312 |
> |
quaternion technique under the microcanonical ($NVE$) ensemble. An |
| 313 |
> |
earlier implementation of our simulation code utilized quaternions for |
| 314 |
> |
propagation of rotational motion; however, a detailed investigation |
| 315 |
> |
showed that they resulted in a steady drift in the total energy, |
| 316 |
> |
something that had also been observed by Kol {\it et al.} (see |
| 317 |
> |
section~\ref{sec:waterSimMethods} for information on this |
| 318 |
> |
investigation).\cite{Kol97} |
| 319 |
> |
|
| 320 |
> |
Because of these issues involving integration of the orientational |
| 321 |
> |
motion using both Euler angles and quaternions, we decided to focus on |
| 322 |
> |
a relatively new scheme that propagates the entire nine parameter |
| 323 |
> |
rotation matrix. This techniques is a velocity-Verlet version of the |
| 324 |
> |
symplectic splitting method proposed by Dullweber, Leimkuhler and |
| 325 |
> |
McLachlan ({\sc dlm}).\cite{Dullweber97} When there are no directional |
| 326 |
> |
atoms or rigid bodies present in the simulation, this {\sc dlm} |
| 327 |
> |
integrator reduces to the standard velocity-Verlet integrator, which |
| 328 |
> |
is known to effectively sample the constant energy $NVE$ |
| 329 |
> |
ensemble.\cite{Frenkel02} |
| 330 |
> |
|
| 331 |
> |
The key aspect of the integration method proposed by Dullweber |
| 332 |
> |
\emph{et al.} is that the entire $3 \times 3$ rotation matrix is |
| 333 |
> |
propagated from one time step to the next. In the past, this would not |
| 334 |
> |
have been as feasible, since the rotation matrix for a single body has |
| 335 |
> |
nine elements compared with the more memory-efficient methods (using |
| 336 |
> |
three Euler angles or four quaternions). Computer memory has become |
| 337 |
> |
much less costly in recent years, and this can be translated into |
| 338 |
> |
substantial benefits in energy conservation. |
| 339 |
> |
|
| 340 |
> |
The integration of the equations of motion is carried out in a |
| 341 |
> |
velocity-Verlet style two-part algorithm.\cite{Swope82} The first part |
| 342 |
> |
({\tt moveA}) consists of a half-step ($t + \delta t/2$) of the linear |
| 343 |
> |
velocity (${\bf v}$) and angular momenta ({\bf j}) and a full-step ($t |
| 344 |
> |
+ \delta t$) of the positions ({\bf r}) and rotation matrix, |
| 345 |
> |
\begin{equation*} |
| 346 |
> |
{\tt moveA} = \left\{\begin{array}{r@{\quad\leftarrow\quad}l} |
| 347 |
> |
{\bf v}\left(t + \delta t / 2\right) & {\bf v}(t) |
| 348 |
> |
+ \left( {\bf f}(t) / m \right)(\delta t/2), \\ |
| 349 |
> |
% |
| 350 |
> |
{\bf r}(t + \delta t) & {\bf r}(t) |
| 351 |
> |
+ {\bf v}\left(t + \delta t / 2 \right)\delta t, \\ |
| 352 |
> |
% |
| 353 |
> |
{\bf j}\left(t + \delta t / 2 \right) & {\bf j}(t) |
| 354 |
> |
+ \boldsymbol{\tau}^b(t)(\delta t/2), \\ |
| 355 |
> |
% |
| 356 |
> |
\mathsf{A}(t + \delta t) & \mathrm{rotate}\left( {\bf j} |
| 357 |
> |
(t + \delta t / 2)\delta t \cdot |
| 358 |
> |
\overleftrightarrow{\mathsf{I}}^{-1} \right), |
| 359 |
> |
\end{array}\right. |
| 360 |
> |
\end{equation*} |
| 361 |
> |
where $\overleftrightarrow{\mathsf{I}}^{-1}$ is the inverse of the |
| 362 |
> |
moment of inertia tensor. The $\mathrm{rotate}$ function is the |
| 363 |
> |
product of rotations about the three body-fixed axes, |
| 364 |
> |
\begin{equation} |
| 365 |
> |
\mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot |
| 366 |
> |
\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y / |
| 367 |
> |
2) \cdot \mathsf{G}_x(a_x /2), |
| 368 |
> |
\label{eq:dlmTrot} |
| 369 |
> |
\end{equation} |
| 370 |
> |
where each rotational propagator, $\mathsf{G}_\alpha(\theta)$, rotates |
| 371 |
> |
both the rotation matrix ($\mathsf{A}$) and the body-fixed angular |
| 372 |
> |
momentum (${\bf j}$) by an angle $\theta$ around body-fixed axis |
| 373 |
> |
$\alpha$, |
| 374 |
> |
\begin{equation} |
| 375 |
> |
\mathsf{G}_\alpha( \theta ) = \left\{ |
| 376 |
> |
\begin{array}{l@{\quad\leftarrow\quad}l} |
| 377 |
> |
\mathsf{A}(t) & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^\textrm{T},\\ |
| 378 |
> |
{\bf j}(t) & \mathsf{R}_\alpha(\theta) \cdot {\bf j}(0). |
| 379 |
> |
\end{array} |
| 380 |
> |
\right. |
| 381 |
> |
\end{equation} |
| 382 |
> |
$\mathsf{R}_\alpha$ is a quadratic approximation to the single-axis |
| 383 |
> |
rotation matrix. For example, in the small-angle limit, the rotation |
| 384 |
> |
matrix around the body-fixed x-axis can be approximated as |
| 385 |
> |
\begin{equation} |
| 386 |
> |
\mathsf{R}_x(\theta) \approx \left( |
| 387 |
> |
\begin{array}{ccc} |
| 388 |
> |
1 & 0 & 0 \\ |
| 389 |
> |
0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} & -\frac{\theta}{1+\theta^2 / 4} \\ |
| 390 |
> |
0 & \frac{\theta}{1+\theta^2 / 4} & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} |
| 391 |
> |
\end{array} |
| 392 |
> |
\right). |
| 393 |
> |
\end{equation} |
| 394 |
> |
The remaining rotations follow in a straightforward manner. As seen |
| 395 |
> |
from the form of equation~(\ref{eq:dlmTrot}), the {\sc dlm} method |
| 396 |
> |
uses a Trotter factorization of the orientational |
| 397 |
> |
propagator.\cite{Trotter59} This has three effects: |
| 398 |
> |
\begin{enumerate} |
| 399 |
> |
\item the integrator is area-preserving in phase space (i.e. it is |
| 400 |
> |
{\it symplectic}), |
| 401 |
> |
\item the integrator is time-{\it reversible}, and |
| 402 |
> |
\item the error for a single time step is of order |
| 403 |
> |
$\mathcal{O}\left(\delta t^3\right)$ for time steps of length $\delta t$. |
| 404 |
> |
\end{enumerate} |
| 405 |
> |
|
| 406 |
> |
After the initial half-step ({\tt moveA}), the forces and torques are |
| 407 |
> |
evaluated for all of the particles. Once completed, the velocities can |
| 408 |
> |
be advanced to complete the second-half of the two-part algorithm |
| 409 |
> |
({\tt moveB}), resulting an a completed full step of both the |
| 410 |
> |
positions and momenta, |
| 411 |
> |
\begin{equation*} |
| 412 |
> |
{\tt moveB} = \left\{\begin{array}{r@{\quad\leftarrow\quad}l} |
| 413 |
> |
{\bf v}\left(t + \delta t \right) & |
| 414 |
> |
{\bf v}\left(t + \delta t / 2 \right) |
| 415 |
> |
+ \left({\bf f}(t + \delta t) / m \right)(\delta t/2), \\ |
| 416 |
> |
% |
| 417 |
> |
{\bf j}\left(t + \delta t \right) & |
| 418 |
> |
{\bf j}\left(t + \delta t / 2 \right) |
| 419 |
> |
+ \boldsymbol{\tau}^b(t + \delta t)(\delta t/2). |
| 420 |
> |
\end{array}\right. |
| 421 |
> |
\end{equation*} |
| 422 |
> |
|
| 423 |
> |
The matrix rotations used in the {\sc dlm} method end up being more |
| 424 |
> |
costly computationally than the simpler arithmetic quaternion |
| 425 |
> |
propagation. With the same time step, a 1024-molecule water simulation |
| 426 |
> |
incurs approximately a 10\% increase in computation time using the |
| 427 |
> |
{\sc dlm} method in place of quaternions. This cost is more than |
| 428 |
> |
justified when comparing the energy conservation achieved by the two |
| 429 |
> |
methods. Figure \ref{fig:quatdlm} provides a comparative analysis of |
| 430 |
> |
the {\sc dlm} method versus the traditional quaternion scheme. |
| 431 |
> |
|
| 432 |
> |
\begin{figure} |
| 433 |
> |
\centering |
| 434 |
> |
\includegraphics[width=3.5in]{./figures/dlmVsQuat.pdf} |
| 435 |
> |
\caption[Energy conservation analysis of the {\sc dlm} and quaternion |
| 436 |
> |
integration methods]{Analysis of the energy conservation of the {\sc |
| 437 |
> |
dlm} and quaternion integration methods. $\delta \mathrm{E}_1$ is the |
| 438 |
> |
linear drift in energy over time and $\delta \mathrm{E}_0$ is the |
| 439 |
> |
standard deviation of energy fluctuations around this drift. All |
| 440 |
> |
simulations were of a 1024 SSD water system at 298 K starting from the |
| 441 |
> |
same initial configuration. Note that the {\sc dlm} method provides |
| 442 |
> |
more than an order-of-magnitude improvement in both the energy drift |
| 443 |
> |
and the size of the energy fluctuations when compared with the |
| 444 |
> |
quaternion method at any given time step. At time steps larger than 4 |
| 445 |
> |
fs, the quaternion scheme resulted in rapidly rising energies which |
| 446 |
> |
eventually lead to simulation failure. Using the {\sc dlm} method, |
| 447 |
> |
time steps up to 8 fs can be taken before this behavior is evident.} |
| 448 |
> |
\label{fig:quatdlm} |
| 449 |
> |
\end{figure} |
| 450 |
> |
|
| 451 |
> |
In figure \ref{fig:quatdlm}, $\delta \mbox{E}_1$ is a measure of the |
| 452 |
> |
linear energy drift in units of $\mbox{kcal mol}^{-1}$ per particle |
| 453 |
> |
over a nanosecond of simulation time, and $\delta \mbox{E}_0$ is the |
| 454 |
> |
standard deviation of the energy fluctuations in units of $\mbox{kcal |
| 455 |
> |
mol}^{-1}$ per particle. In the top plot, it is apparent that the |
| 456 |
> |
energy drift is reduced by a significant amount (2 to 3 orders of |
| 457 |
> |
magnitude improvement at all tested time steps) by choosing the {\sc |
| 458 |
> |
dlm} method over the simple non-symplectic quaternion integration |
| 459 |
> |
method. In addition to this improvement in energy drift, the |
| 460 |
> |
fluctuations in the total energy are also dampened by 1 to 2 orders of |
| 461 |
> |
magnitude by utilizing the {\sc dlm} method. |
| 462 |
> |
|
| 463 |
> |
\begin{figure} |
| 464 |
> |
\centering |
| 465 |
> |
\includegraphics[width=\linewidth]{./figures/compCost.pdf} |
| 466 |
> |
\caption[Energy drift as a function of required simulation run |
| 467 |
> |
time]{Energy drift as a function of required simulation run time. |
| 468 |
> |
$\delta \mathrm{E}_1$ is the linear drift in energy over time. |
| 469 |
> |
Simulations were performed on a single 2.5 GHz Pentium 4 |
| 470 |
> |
processor. Simulation time comparisons can be made by tracing |
| 471 |
> |
horizontally from one curve to the other. For example, a simulation |
| 472 |
> |
that takes 24 hours using the {\sc dlm} method will take roughly |
| 473 |
> |
210 hours using the simple quaternion method if the same degree of |
| 474 |
> |
energy conservation is desired.} |
| 475 |
> |
\label{fig:cpuCost} |
| 476 |
> |
\end{figure} |
| 477 |
> |
Although the {\sc dlm} method is more computationally expensive than |
| 478 |
> |
the traditional quaternion scheme for propagating of a time step, |
| 479 |
> |
consideration of the computational cost for a long simulation with a |
| 480 |
> |
particular level of energy conservation is in order. A plot of energy |
| 481 |
> |
drift versus computational cost was generated |
| 482 |
> |
(Fig.~\ref{fig:cpuCost}). This figure provides an estimate of the CPU |
| 483 |
> |
time required under the two integration schemes for 1 nanosecond of |
| 484 |
> |
simulation time for the model 1024-molecule system. By choosing a |
| 485 |
> |
desired energy drift value it is possible to determine the CPU time |
| 486 |
> |
required for both methods. If a $\delta \mbox{E}_1$ of |
| 487 |
> |
0.001~kcal~mol$^{-1}$ per particle is desired, a nanosecond of |
| 488 |
> |
simulation time will require ~19 hours of CPU time with the {\sc dlm} |
| 489 |
> |
integrator, while the quaternion scheme will require ~154 hours of CPU |
| 490 |
> |
time. This demonstrates the computational advantage of the {\sc dlm} |
| 491 |
> |
integration scheme. |
| 492 |
> |
|
| 493 |
> |
\section{Accumulating Interactions} |
| 494 |
> |
|
| 495 |
> |
In the force calculation between {\tt moveA} and {\tt moveB} mentioned |
| 496 |
> |
in section \ref{sec:IntroIntegrate}, we need to accumulate the |
| 497 |
> |
potential and forces (and torques if the particle is a rigid body or |
| 498 |
> |
multipole) on each particle from their surroundings. This can quickly |
| 499 |
> |
become a cumbersome task for large systems since the number of pair |
| 500 |
> |
interactions increases by $\mathcal{O}(N(N-1)/2)$ when accumulating |
| 501 |
> |
interactions between all particles in the system. (Note the |
| 502 |
> |
utilization of Newton's third law to reduce the interaction number |
| 503 |
> |
from $\mathcal{O}(N^2)$.) Using periodic boundary conditions (section |
| 504 |
> |
\ref{sec:periodicBoundaries}) further complicates matters by turning |
| 505 |
> |
the finite system into an infinitely repeating one. Fortunately, we |
| 506 |
> |
can reduce the scale of this problem by using spherical cutoff |
| 507 |
> |
methods. |
| 508 |
> |
|
| 509 |
> |
\begin{figure} |
| 510 |
> |
\centering |
| 511 |
> |
\includegraphics[width=3.5in]{./figures/sphericalCut.pdf} |
| 512 |
> |
\caption{When using a spherical cutoff, only particles within a chosen |
| 513 |
> |
cutoff radius distance, $R_\textrm{c}$, of the central particle are |
| 514 |
> |
included in the pairwise summation. This reduces a problem that scales |
| 515 |
> |
by $\sim\mathcal{O}(N^2)$ to one that scales by $\sim\mathcal{O}(N)$.} |
| 516 |
> |
\label{fig:sphereCut} |
| 517 |
> |
\end{figure} |
| 518 |
> |
With spherical cutoffs, rather than accumulating the full set of |
| 519 |
> |
interactions between all particles in the simulation, we only |
| 520 |
> |
explicitly consider interactions between particles separated by less |
| 521 |
> |
than a specified cutoff radius distance, $R_\textrm{c}$, (see figure |
| 522 |
> |
\ref{fig:sphereCut}). This reduces the scaling of the interaction to |
| 523 |
> |
$\mathcal{O}(N\cdot\textrm{c})$, where `c' is a value that depends on |
| 524 |
> |
the size of $R_\textrm{c}$ (c $\approx R_\textrm{c}^3$). Determination |
| 525 |
> |
of which particles are within the cutoff is also an issue, because |
| 526 |
> |
this process requires a full loop over all $N(N-1)/2$ pairs. To reduce |
| 527 |
> |
the this expense, we can use neighbor lists.\cite{Verlet67,Thompson83} |
| 528 |
> |
|
| 529 |
> |
When using neighbor lists, we build a second list of particles built |
| 530 |
> |
from a list radius $R_\textrm{l}$, which is larger than |
| 531 |
> |
$R_\textrm{c}$. Once any particle within $R_\textrm{l}$ moves half the |
| 532 |
> |
distance of $R_\textrm{l}-R_\textrm{c}$ (the ``skin'' thickness), we |
| 533 |
> |
rebuild the list with the full $N(N-1)/2$ loop.\cite{Verlet67} With an |
| 534 |
> |
appropriate skin thickness, these updates are only performed every |
| 535 |
> |
$\sim$20 time steps, significantly reducing the time spent on |
| 536 |
> |
pair-list bookkeeping operations.\cite{Allen87} If these neighbor |
| 537 |
> |
lists are utilized, it is important that these list updates occur |
| 538 |
> |
regularly. Incorrect application of this technique leads to |
| 539 |
> |
non-physical dynamics, such as the ``flying block of ice'' behavior |
| 540 |
> |
for which improper neighbor list handling was identified a one of the |
| 541 |
> |
possible causes.\cite{Harvey98,Sagui99} |
| 542 |
> |
|
| 543 |
> |
\subsection{Correcting Cutoff Discontinuities} |
| 544 |
> |
\begin{figure} |
| 545 |
> |
\centering |
| 546 |
> |
\includegraphics[width=3.5in]{./figures/ljCutoffSquare.pdf} |
| 547 |
> |
\caption{The common methods to smooth the potential discontinuity |
| 548 |
> |
introduced when using a cutoff include a shifted potential, a shifted |
| 549 |
> |
force, and a switching function. The shifted potential and shifted |
| 550 |
> |
force both lift the whole potential so that it zeroes at |
| 551 |
> |
$R_\textrm{c}$, thereby reducing the strength of the interaction. The |
| 552 |
> |
(cubic) switching function only alters the potential in the switching |
| 553 |
> |
region in order to smooth out the discontinuity.} |
| 554 |
> |
\label{fig:ljCutoff} |
| 555 |
> |
\end{figure} |
| 556 |
> |
As the distance between a pair of particles fluctuates around |
| 557 |
> |
$R_\textrm{c}$, there will be sudden discontinuous jumps in the |
| 558 |
> |
potential (and forces) due to their inclusion and exclusion from the |
| 559 |
> |
interaction loop. In order to prevent heating and poor energy |
| 560 |
> |
conservation in the simulations, this discontinuity needs to be |
| 561 |
> |
smoothed out. There are several ways to modify the potential function |
| 562 |
> |
to eliminate these discontinuities, and the easiest method is shifting |
| 563 |
> |
the potential. To shift the potential, we simply subtract out the |
| 564 |
> |
value of the function at $R_\textrm{c}$ from the whole potential. The |
| 565 |
> |
shifted form of the Lennard-Jones potential is |
| 566 |
> |
\begin{equation} |
| 567 |
> |
V_\textrm{SLJ} = \left\{\begin{array}{l@{\quad\quad}l} |
| 568 |
> |
V_\textrm{LJ}(r_{ij}) - V_\textrm{LJ}(R_\textrm{c}) & r_{ij} < R_\textrm{c}, \\ |
| 569 |
> |
0 & r_{ij} \geqslant R_\textrm{c}, |
| 570 |
> |
\end{array}\right. |
| 571 |
> |
\end{equation} |
| 572 |
> |
where |
| 573 |
> |
\begin{equation} |
| 574 |
> |
V_\textrm{LJ}(r_{ij}) = |
| 575 |
> |
4\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12} - |
| 576 |
> |
\left(\frac{\sigma}{r_{ij}}\right)^6\right]. |
| 577 |
> |
\end{equation} |
| 578 |
> |
In figure \ref{fig:ljCutoff}, the shifted form of the potential |
| 579 |
> |
reaches zero at the cutoff radius at the expense of the correct |
| 580 |
> |
magnitude inside $R_\textrm{c}$. This correction method also does |
| 581 |
> |
nothing to correct the discontinuity in the forces. We can account for |
| 582 |
> |
this force discontinuity by applying the shifting in the forces as |
| 583 |
> |
well as in the potential via |
| 584 |
> |
\begin{equation} |
| 585 |
> |
V_\textrm{SFLJ} = \left\{\begin{array}{l@{\quad\quad}l} |
| 586 |
> |
V_\textrm{LJ}({r_{ij}}) - V_\textrm{LJ}(R_\textrm{c}) - |
| 587 |
> |
\left(\frac{d V(r_{ij})}{d r_{ij}}\right)_{r_{ij}=R_\textrm{c}} |
| 588 |
> |
(r_{ij} - R_\textrm{c}) & r_{ij} < R_\textrm{c}, \\ |
| 589 |
> |
0 & r_{ij} \geqslant R_\textrm{c}. |
| 590 |
> |
\end{array}\right. |
| 591 |
> |
\end{equation} |
| 592 |
> |
The forces are continuous with this potential; however, the inclusion |
| 593 |
> |
of the derivative term distorts the potential even further than |
| 594 |
> |
shifting only the potential as shown in figure \ref{fig:ljCutoff}. |
| 595 |
> |
|
| 596 |
> |
The method for correcting these discontinuities which results in the |
| 597 |
> |
smallest perturbation in both the potential and the forces is the use |
| 598 |
> |
of a switching function. The cubic switching function, |
| 599 |
> |
\begin{equation} |
| 600 |
> |
S(r) = \left\{\begin{array}{l@{\quad\quad}l} |
| 601 |
> |
1 & r_{ij} \leqslant R_\textrm{sw}, \\ |
| 602 |
> |
\frac{(R_\textrm{c} + 2r_{ij} - 3R_\textrm{sw}) |
| 603 |
> |
(R_\textrm{c} - r_{ij})^2}{(R_\textrm{c} - R_\textrm{sw})^3} |
| 604 |
> |
& R_\textrm{sw} < r_{ij} \leqslant R_\textrm{c}, \\ |
| 605 |
> |
0 & r_{ij} > R_\textrm{c}, |
| 606 |
> |
\end{array}\right. |
| 607 |
> |
\end{equation} |
| 608 |
> |
is sufficient to smooth the potential (again, see figure |
| 609 |
> |
\ref{fig:ljCutoff}) and the forces by only perturbing the potential in |
| 610 |
> |
the switching region. If smooth second derivatives are required, a |
| 611 |
> |
higher order polynomial switching function (i.e. fifth order |
| 612 |
> |
polynomial) can be used.\cite{Andrea83,Leach01} It should be noted |
| 613 |
> |
that the higher the order of the polynomial, the more abrupt the |
| 614 |
> |
switching transition. |
| 615 |
> |
|
| 616 |
> |
\subsection{\label{sec:LJCorrections}Long-Range Interaction Corrections} |
| 617 |
> |
|
| 618 |
> |
While a good approximation, accumulating interactions only from nearby |
| 619 |
> |
particles can lead to some issues, because particles at distances |
| 620 |
> |
greater than $R_\textrm{c}$ still have a small effect. For instance, |
| 621 |
> |
while the strength of the Lennard-Jones interaction is quite weak at |
| 622 |
> |
$R_\textrm{c}$ in figure \ref{fig:ljCutoff}, we are discarding all of |
| 623 |
> |
the attractive interactions that extend out to extremely long |
| 624 |
> |
distances. Thus, the potential will be a little too high and the |
| 625 |
> |
pressure on the central particle from the surroundings will be a |
| 626 |
> |
little too low. For homogeneous Lennard-Jones systems, we can correct |
| 627 |
> |
for this effect by assuming a uniform density ($\rho$) and integrating |
| 628 |
> |
the missing part, |
| 629 |
> |
\begin{equation} |
| 630 |
> |
V_\textrm{full}(r_{ij}) \approx V_\textrm{c} |
| 631 |
> |
+ 2\pi N\rho\int^\infty_{R_\textrm{c}}r^2V_\textrm{LJ}(r)dr, |
| 632 |
> |
\end{equation} |
| 633 |
> |
where $V_\textrm{c}$ is the truncated Lennard-Jones |
| 634 |
> |
potential.\cite{Allen87} Like the potential, other properties can be |
| 635 |
> |
corrected by integration over the relevant function. Note that with |
| 636 |
> |
heterogeneous systems, this correction breaks down because the density |
| 637 |
> |
is no longer uniform. |
| 638 |
> |
|
| 639 |
> |
Correcting long-range electrostatic interactions is a topic of great |
| 640 |
> |
importance in the field of molecular simulations. There have been |
| 641 |
> |
several techniques developed to address this issue, like the Ewald |
| 642 |
> |
summation and the reaction field technique. An in-depth analysis of |
| 643 |
> |
this problem, as well as useful correction techniques, is presented in |
| 644 |
> |
chapter \ref{chap:electrostatics}. |
| 645 |
> |
|
| 646 |
> |
\subsection{\label{sec:periodicBoundaries}Periodic Boundary Conditions} |
| 647 |
> |
|
| 648 |
> |
In typical molecular dynamics simulations there are no restrictions |
| 649 |
> |
placed on the motion of particles outside of what the inter-particle |
| 650 |
> |
interactions dictate. This means that if a particle collides with |
| 651 |
> |
other particles, it is free to move away from the site of the |
| 652 |
> |
collision. If we consider the entire system as a collection of |
| 653 |
> |
particles, they are not confined by walls of the ``simulation box'' |
| 654 |
> |
and can freely move away from the other particles. With no boundary |
| 655 |
> |
considerations, particles moving outside of the simulation box |
| 656 |
> |
enter a vacuum. This is correct behavior for cluster simulations in a |
| 657 |
> |
vacuum; however, if we want to simulate bulk or spatially infinite |
| 658 |
> |
systems, we need to use periodic boundary conditions. |
| 659 |
> |
|
| 660 |
> |
\begin{figure} |
| 661 |
> |
\centering |
| 662 |
> |
\includegraphics[width=4.5in]{./figures/periodicImage.pdf} |
| 663 |
> |
\caption{With periodic boundary conditions imposed, when particles |
| 664 |
> |
move out of one side the simulation box, they wrap back in the |
| 665 |
> |
opposite side. In this manner, a finite system of particles behaves as |
| 666 |
> |
an infinite system.} |
| 667 |
> |
\label{fig:periodicImage} |
| 668 |
> |
\end{figure} |
| 669 |
> |
In periodic boundary conditions, as a particle moves outside one wall |
| 670 |
> |
of the simulation box, the coordinates are remapped such that the |
| 671 |
> |
particle enters the opposing side of the box. This process is easy to |
| 672 |
> |
visualize in two dimensions as shown in figure \ref{fig:periodicImage} |
| 673 |
> |
and can occur in three dimensions, though it is not as easy to |
| 674 |
> |
visualize. Remapping the actual coordinates of the particles can be |
| 675 |
> |
problematic in that the we are restricting the distance a particle can |
| 676 |
> |
move from it's point of origin to a diagonal of the simulation |
| 677 |
> |
box. Thus, even though we are not confining the system with hard |
| 678 |
> |
walls, we are confining the particle coordinates to a particular |
| 679 |
> |
region in space. To avoid this ``soft'' confinement, it is common |
| 680 |
> |
practice to allow the particle coordinates to expand in an |
| 681 |
> |
unrestricted fashion while calculating interactions using a wrapped |
| 682 |
> |
set of ``minimum image'' coordinates. These coordinates need not be |
| 683 |
> |
stored because they are easily calculated while determining particle |
| 684 |
> |
distances. |
| 685 |
> |
|
| 686 |
> |
\section{Calculating Properties} |
| 687 |
> |
|
| 688 |
> |
In order to use simulations to model experimental processes and |
| 689 |
> |
evaluate theories, we need to be able to extract properties from the |
| 690 |
> |
results. In experiments, we can measure thermodynamic properties such |
| 691 |
> |
as the pressure and free energy. In computer simulations, we can |
| 692 |
> |
calculate properties from the motion and configuration of particles in |
| 693 |
> |
the system and make connections between these properties and the |
| 694 |
> |
experimental thermodynamic properties through statistical mechanics. |
| 695 |
> |
|
| 696 |
> |
The work presented in the later chapters use the canonical ($NVT$), |
| 697 |
> |
isobaric-isothermal ($NPT$), and microcanonical ($NVE$) statistical |
| 698 |
> |
mechanical ensembles. The different ensembles lend themselves to more |
| 699 |
> |
effectively calculating specific properties. For instance, if we |
| 700 |
> |
concerned ourselves with the calculation of dynamic properties, which |
| 701 |
> |
are dependent upon the motion of the particles, it is better to choose |
| 702 |
> |
an ensemble that does not add artificial motions to keep properties |
| 703 |
> |
like the temperature or pressure constant. In this case, the $NVE$ |
| 704 |
> |
ensemble would be the most appropriate choice. In chapter |
| 705 |
> |
\ref{chap:ice}, we discuss calculating free energies, which are |
| 706 |
> |
non-mechanical thermodynamic properties, and these calculations also |
| 707 |
> |
depend on the chosen ensemble.\cite{Allen87} The Helmholtz free energy |
| 708 |
> |
($A$) depends on the $NVT$ partition function ($Q_{NVT}$), |
| 709 |
> |
\begin{equation} |
| 710 |
> |
A = -k_\textrm{B}T\ln Q_{NVT}, |
| 711 |
> |
\end{equation} |
| 712 |
> |
while the Gibbs free energy ($G$) depends on the $NPT$ partition |
| 713 |
> |
function ($Q_{NPT}$), |
| 714 |
> |
\begin{equation} |
| 715 |
> |
G = -k_\textrm{B}T\ln Q_{NPT}. |
| 716 |
> |
\end{equation} |
| 717 |
> |
It is also useful to note that the conserved quantities of the $NVT$ |
| 718 |
> |
and $NPT$ ensembles are related to the Helmholtz and Gibbs free |
| 719 |
> |
energies respectively.\cite{Melchionna93} |
| 720 |
> |
|
| 721 |
> |
Integrating the equations of motion is a simple method to obtain a |
| 722 |
> |
sequence of configurations that sample the chosen ensemble. For each |
| 723 |
> |
of these configurations, we can calculate an instantaneous value for a |
| 724 |
> |
chosen property like the density in the $NPT$ ensemble, where the |
| 725 |
> |
volume is allowed to fluctuate. The density for the simulation is |
| 726 |
> |
calculated from an average over the densities for the individual |
| 727 |
> |
configurations. Since the configurations from the integration process |
| 728 |
> |
are related to one another by the time evolution of the interactions, |
| 729 |
> |
this average is technically a time average. In calculating |
| 730 |
> |
thermodynamic properties, we would actually prefer an ensemble average |
| 731 |
> |
that is representative of all accessible states of the system. We can |
| 732 |
> |
calculate thermodynamic properties from the time average by taking |
| 733 |
> |
advantage of the ergodic hypothesis, which states that for a |
| 734 |
> |
sufficiently chaotic system, and over a long enough period of time, |
| 735 |
> |
the time and ensemble averages are the same. |
| 736 |
> |
|
| 737 |
> |
In addition to the average, the fluctuations of a particular property |
| 738 |
> |
can be determined via the standard deviation. Not only are |
| 739 |
> |
fluctuations useful for determining the spread of values around the |
| 740 |
> |
average and the error in the calculation of the value, they are also |
| 741 |
> |
useful for measuring various thermodynamic properties in computer |
| 742 |
> |
simulations. In section \ref{sec:t5peThermo}, we use fluctuations in |
| 743 |
> |
properties like the enthalpy and volume to calculate other |
| 744 |
> |
thermodynamic properties, such as the constant pressure heat capacity |
| 745 |
> |
and the isothermal compressibility. |
| 746 |
> |
|
| 747 |
> |
\section{OOPSE} |
| 748 |
> |
|
| 749 |
> |
In the following chapters, the above techniques are all utilized in |
| 750 |
> |
the study of molecular systems. There are a number of excellent |
| 751 |
> |
simulation packages available, both free and commercial, which |
| 752 |
> |
incorporate many of these |
| 753 |
> |
methods.\cite{Brooks83,MacKerell98,Pearlman95,Berendsen95,Lindahl01,Smith96,Ponder87} |
| 754 |
> |
Because of our interest in rigid body dynamics, point multipoles, and |
| 755 |
> |
systems where the orientational degrees of freedom cannot be handled |
| 756 |
> |
using the common constraint procedures (like {\sc shake}), the |
| 757 |
> |
majority of the following work was performed using {\sc oopse}, the |
| 758 |
> |
object-oriented parallel simulation engine.\cite{Meineke05} The {\sc |
| 759 |
> |
oopse} package started out as a collection of separate programs |
| 760 |
> |
written within our group, and has developed into one of the few |
| 761 |
> |
parallel molecular dynamics packages capable of accurately integrating |
| 762 |
> |
rigid bodies and point multipoles. This simulation package is |
| 763 |
> |
open-source software, available to anyone interested in performing |
| 764 |
> |
molecular dynamics simulations. More information about {\sc oopse} can |
| 765 |
> |
be found in reference \cite{Meineke05} or at the {\tt |
| 766 |
> |
http://oopse.org} website. |
| 767 |
> |
|
| 768 |
> |
|
