1 |
chrisfen |
2987 |
\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND} |
2 |
|
|
|
3 |
chrisfen |
3001 |
In order to advance our understanding of natural chemical and physical |
4 |
|
|
processes, researchers develop explanations for events observed |
5 |
|
|
experimentally. These hypotheses, supported by a body of corroborating |
6 |
|
|
observations, can develop into accepted theories for how these |
7 |
|
|
processes occur. This validation process, as well as testing the |
8 |
|
|
limits of these theories, is essential in developing a firm foundation |
9 |
|
|
for our knowledge. Theories involving molecular scale systems cause |
10 |
|
|
particular difficulties in this process because their size and often |
11 |
|
|
fast motions make them difficult to observe experimentally. One useful |
12 |
|
|
tool for addressing these difficulties is computer simulation. |
13 |
|
|
|
14 |
|
|
In computer simulations, we can develop models from either the theory |
15 |
|
|
or experimental knowledge and then test them in a controlled |
16 |
|
|
environment. Work done in this manner allows us to further refine |
17 |
|
|
theories, more accurately represent what is happening in experimental |
18 |
|
|
observations, and even make predictions about what one will see in |
19 |
|
|
experiments. Thus, computer simulations of molecular systems act as a |
20 |
|
|
bridge between theory and experiment. |
21 |
|
|
|
22 |
|
|
Depending on the system of interest, there are a variety of different |
23 |
|
|
computational techniques that can used to test and gather information |
24 |
|
|
from the developed models. In the study of classical systems, the two |
25 |
|
|
most commonly used techniques are Monte Carlo and molecular |
26 |
|
|
dynamics. Both of these methods operate by calculating interactions |
27 |
|
|
between particles of our model systems; however, the progression of |
28 |
|
|
the simulation under the different techniques is vastly |
29 |
|
|
different. Monte Carlo operates through random configuration changes |
30 |
|
|
that that follow rules adhering to a specific statistical mechanics |
31 |
|
|
ensemble, while molecular dynamics is chiefly concerned with solving |
32 |
|
|
the classic equations of motion to move between configurations within |
33 |
|
|
an ensemble. Thermodynamic properties can be calculated from both |
34 |
|
|
techniques; but because of the random nature of Monte Carlo, only |
35 |
|
|
molecular dynamics can be used to investigate dynamical |
36 |
|
|
quantities. The research presented in the following chapters utilized |
37 |
|
|
molecular dynamics near exclusively, so we will present a general |
38 |
|
|
introduction to molecular dynamics and not Monte Carlo. There are |
39 |
|
|
several resources available for those desiring a more in-depth |
40 |
|
|
presentation of either of these |
41 |
|
|
techniques.\cite{Allen87,Frenkel02,Leach01} |
42 |
|
|
|
43 |
|
|
\section{\label{sec:MolecularDynamics}Molecular Dynamics} |
44 |
|
|
|
45 |
|
|
\section{\label{sec:MovingParticles}Propagating Particle Motion} |
46 |
|
|
|
47 |
|
|
\subsection{\label{sec:IntroIntegrate}Rigid Body Motion} |
48 |
|
|
|
49 |
|
|
Rigid bodies are non-spherical particles or collections of particles |
50 |
|
|
(e.g. $\mbox{C}_{60}$) that have a constant internal potential and |
51 |
|
|
move collectively.\cite{Goldstein01} Discounting iterative constraint |
52 |
|
|
procedures like {\sc shake} and {\sc rattle}, they are not included in |
53 |
|
|
most simulation packages because of the algorithmic complexity |
54 |
|
|
involved in propagating orientational degrees of |
55 |
|
|
freedom.\cite{Ryckaert77,Andersen83,Krautler01} Integrators which |
56 |
|
|
propagate orientational motion with an acceptable level of energy |
57 |
|
|
conservation for molecular dynamics are relatively new inventions. |
58 |
|
|
|
59 |
|
|
Moving a rigid body involves determination of both the force and |
60 |
|
|
torque applied by the surroundings, which directly affect the |
61 |
|
|
translational and rotational motion in turn. In order to accumulate |
62 |
|
|
the total force on a rigid body, the external forces and torques must |
63 |
|
|
first be calculated for all the internal particles. The total force on |
64 |
|
|
the rigid body is simply the sum of these external forces. |
65 |
|
|
Accumulation of the total torque on the rigid body is more complex |
66 |
|
|
than the force because the torque is applied to the center of mass of |
67 |
|
|
the rigid body. The space-fixed torque on rigid body $i$ is |
68 |
|
|
\begin{equation} |
69 |
|
|
\boldsymbol{\tau}_i= |
70 |
|
|
\sum_{a}\biggl[(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia} |
71 |
|
|
+ \boldsymbol{\tau}_{ia}\biggr], |
72 |
|
|
\label{eq:torqueAccumulate} |
73 |
|
|
\end{equation} |
74 |
|
|
where $\boldsymbol{\tau}_i$ and $\mathbf{r}_i$ are the torque on and |
75 |
|
|
position of the center of mass respectively, while $\mathbf{f}_{ia}$, |
76 |
|
|
$\mathbf{r}_{ia}$, and $\boldsymbol{\tau}_{ia}$ are the force on, |
77 |
|
|
position of, and torque on the component particles. |
78 |
|
|
|
79 |
|
|
The summation of the total torque is done in the body fixed axis. In |
80 |
|
|
order to move between the space fixed and body fixed coordinate axes, |
81 |
|
|
parameters describing the orientation must be maintained for each |
82 |
|
|
rigid body. At a minimum, the rotation matrix ($\mathsf{A}$) can be |
83 |
|
|
described by the three Euler angles ($\phi, \theta,$ and $\psi$), |
84 |
|
|
where the elements of $\mathsf{A}$ are composed of trigonometric |
85 |
|
|
operations involving $\phi, \theta,$ and $\psi$.\cite{Goldstein01} |
86 |
|
|
Direct propagation of the Euler angles has a known $1/\sin\theta$ |
87 |
|
|
divergence in the equations of motion for $\phi$ and $\psi$, leading |
88 |
|
|
to numerical instabilities any time one of the directional atoms or |
89 |
|
|
rigid bodies has an orientation near $\theta=0$ or |
90 |
|
|
$\theta=\pi$.\cite{Allen87} One of the most practical ways to avoid |
91 |
|
|
this ``gimbal point'' is to switch to another angular set defining the |
92 |
|
|
orientation of the rigid body near this point.\cite{Barojas73} This |
93 |
|
|
procedure results in additional book-keeping and increased algorithm |
94 |
|
|
complexity. In the search for more elegant alternative methods, Evans |
95 |
|
|
proposed the use of quaternions to describe and propagate |
96 |
|
|
orientational motion.\cite{Evans77} |
97 |
|
|
|
98 |
|
|
The quaternion method for integration involves a four dimensional |
99 |
|
|
representation of the orientation of a rigid |
100 |
|
|
body.\cite{Evans77,Evans77b,Allen87} Thus, the elements of |
101 |
|
|
$\mathsf{A}$ can be expressed as arithmetic operations involving the |
102 |
|
|
four quaternions ($q_w, q_x, q_y,$ and $q_z$), |
103 |
|
|
\begin{equation} |
104 |
|
|
\mathsf{A} = \left( \begin{array}{l@{\quad}l@{\quad}l} |
105 |
|
|
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) \\ |
106 |
|
|
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) \\ |
107 |
|
|
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 \\ |
108 |
|
|
\end{array}\right). |
109 |
|
|
\end{equation} |
110 |
|
|
Integration of the equations of motion involves a series of arithmetic |
111 |
|
|
operations involving the quaternions and angular momenta and leads to |
112 |
|
|
performance enhancements over Euler angles, particularly for very |
113 |
|
|
small systems.\cite{Evans77} This integration methods works well for |
114 |
|
|
propagating orientational motion in the canonical ensemble ($NVT$); |
115 |
|
|
however, energy conservation concerns arise when using the simple |
116 |
|
|
quaternion technique under the microcanonical ($NVE$) ensemble. An |
117 |
|
|
earlier implementation of {\sc oopse} utilized quaternions for |
118 |
|
|
propagation of rotational motion; however, a detailed investigation |
119 |
|
|
showed that they resulted in a steady drift in the total energy, |
120 |
|
|
something that has been observed by Kol {\it et al.} (also see |
121 |
|
|
section~\ref{sec:waterSimMethods}).\cite{Kol97} |
122 |
|
|
|
123 |
|
|
Because of the outlined issues involving integration of the |
124 |
|
|
orientational motion using both Euler angles and quaternions, we |
125 |
|
|
decided to focus on a relatively new scheme that propagates the entire |
126 |
|
|
nine parameter rotation matrix. This techniques is a velocity-Verlet |
127 |
|
|
version of the symplectic splitting method proposed by Dullweber, |
128 |
|
|
Leimkuhler and McLachlan ({\sc dlm}).\cite{Dullweber97} When there are |
129 |
|
|
no directional atoms or rigid bodies present in the simulation, this |
130 |
|
|
integrator becomes the standard velocity-Verlet integrator which is |
131 |
|
|
known to effectively sample the microcanonical ($NVE$) |
132 |
|
|
ensemble.\cite{Frenkel02} |
133 |
|
|
|
134 |
|
|
The key aspect of the integration method proposed by Dullweber |
135 |
|
|
\emph{et al.} is that the entire $3 \times 3$ rotation matrix is |
136 |
|
|
propagated from one time step to the next. In the past, this would not |
137 |
|
|
have been as feasible, since the rotation matrix for a single body has |
138 |
|
|
nine elements compared with the more memory-efficient methods (using |
139 |
|
|
three Euler angles or four quaternions). Computer memory has become |
140 |
|
|
much less costly in recent years, and this can be translated into |
141 |
|
|
substantial benefits in energy conservation. |
142 |
|
|
|
143 |
|
|
The integration of the equations of motion is carried out in a |
144 |
|
|
velocity-Verlet style two-part algorithm.\cite{Swope82} The first-part |
145 |
|
|
({\tt moveA}) consists of a half-step ($t + \delta t/2$) of the linear |
146 |
|
|
velocity (${\bf v}$) and angular momenta ({\bf j}) and a full-step ($t |
147 |
|
|
+ \delta t$) of the positions ({\bf r}) and rotation matrix, |
148 |
|
|
\begin{equation*} |
149 |
|
|
{\tt moveA} = \left\{\begin{array}{r@{\quad\leftarrow\quad}l} |
150 |
|
|
{\bf v}\left(t + \delta t / 2\right) & {\bf v}(t) |
151 |
|
|
+ \left( {\bf f}(t) / m \right)(\delta t/2), \\ |
152 |
|
|
% |
153 |
|
|
{\bf r}(t + \delta t) & {\bf r}(t) |
154 |
|
|
+ {\bf v}\left(t + \delta t / 2 \right)\delta t, \\ |
155 |
|
|
% |
156 |
|
|
{\bf j}\left(t + \delta t / 2 \right) & {\bf j}(t) |
157 |
|
|
+ \boldsymbol{\tau}^b(t)(\delta t/2), \\ |
158 |
|
|
% |
159 |
|
|
\mathsf{A}(t + \delta t) & \mathrm{rotate}\left( {\bf j} |
160 |
|
|
(t + \delta t / 2)\delta t \cdot |
161 |
|
|
\overleftrightarrow{\mathsf{I}}^{-1} \right), |
162 |
|
|
\end{array}\right. |
163 |
|
|
\end{equation*} |
164 |
|
|
where $\overleftrightarrow{\mathsf{I}}^{-1}$ is the inverse of the |
165 |
|
|
moment of inertia tensor. The $\mathrm{rotate}$ function is the |
166 |
|
|
product of rotations about the three body-fixed axes, |
167 |
|
|
\begin{equation} |
168 |
|
|
\mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot |
169 |
|
|
\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y / |
170 |
|
|
2) \cdot \mathsf{G}_x(a_x /2), |
171 |
|
|
\label{eq:dlmTrot} |
172 |
|
|
\end{equation} |
173 |
|
|
where each rotational propagator, $\mathsf{G}_\alpha(\theta)$, rotates |
174 |
|
|
both the rotation matrix ($\mathsf{A}$) and the body-fixed angular |
175 |
|
|
momentum (${\bf j}$) by an angle $\theta$ around body-fixed axis |
176 |
|
|
$\alpha$, |
177 |
|
|
\begin{equation} |
178 |
|
|
\mathsf{G}_\alpha( \theta ) = \left\{ |
179 |
|
|
\begin{array}{l@{\quad\leftarrow\quad}l} |
180 |
|
|
\mathsf{A}(t) & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^\textrm{T},\\ |
181 |
|
|
{\bf j}(t) & \mathsf{R}_\alpha(\theta) \cdot {\bf j}(0). |
182 |
|
|
\end{array} |
183 |
|
|
\right. |
184 |
|
|
\end{equation} |
185 |
|
|
$\mathsf{R}_\alpha$ is a quadratic approximation to the single-axis |
186 |
|
|
rotation matrix. For example, in the small-angle limit, the rotation |
187 |
|
|
matrix around the body-fixed x-axis can be approximated as |
188 |
|
|
\begin{equation} |
189 |
|
|
\mathsf{R}_x(\theta) \approx \left( |
190 |
|
|
\begin{array}{ccc} |
191 |
|
|
1 & 0 & 0 \\ |
192 |
|
|
0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} & -\frac{\theta}{1+\theta^2 / 4} \\ |
193 |
|
|
0 & \frac{\theta}{1+\theta^2 / 4} & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} |
194 |
|
|
\end{array} |
195 |
|
|
\right). |
196 |
|
|
\end{equation} |
197 |
|
|
The remaining rotations follow in a straightforward manner. As seen |
198 |
|
|
from the form of equation~(\ref{eq:dlmTrot}), the {\sc dlm} method |
199 |
|
|
uses a Trotter factorization of the orientational |
200 |
|
|
propagator.\cite{Trotter59} This has three effects: |
201 |
|
|
\begin{enumerate} |
202 |
|
|
\item the integrator is area-preserving in phase space (i.e. it is |
203 |
|
|
{\it symplectic}), |
204 |
|
|
\item the integrator is time-{\it reversible}, and |
205 |
|
|
\item the error for a single time step is of order |
206 |
|
|
$\mathcal{O}\left(\delta t^4\right)$ for time steps of length $\delta t$. |
207 |
|
|
\end{enumerate} |
208 |
|
|
|
209 |
|
|
After the initial half-step ({\tt moveA}), the forces and torques are |
210 |
|
|
evaluated for all of the particles. Once completed, the velocities can |
211 |
|
|
be advanced to complete the second-half of the two-part algorithm |
212 |
|
|
({\tt moveB}), resulting an a completed full step of both the |
213 |
|
|
positions and momenta, |
214 |
|
|
\begin{equation*} |
215 |
|
|
{\tt moveB} = \left\{\begin{array}{r@{\quad\leftarrow\quad}l} |
216 |
|
|
{\bf v}\left(t + \delta t \right) & |
217 |
|
|
{\bf v}\left(t + \delta t / 2 \right) |
218 |
|
|
+ \left({\bf f}(t + \delta t) / m \right)(\delta t/2), \\ |
219 |
|
|
% |
220 |
|
|
{\bf j}\left(t + \delta t \right) & |
221 |
|
|
{\bf j}\left(t + \delta t / 2 \right) |
222 |
|
|
+ \boldsymbol{\tau}^b(t + \delta t)(\delta t/2). |
223 |
|
|
\end{array}\right. |
224 |
|
|
\end{equation*} |
225 |
|
|
|
226 |
|
|
The matrix rotations used in the {\sc dlm} method end up being more |
227 |
|
|
costly computationally than the simpler arithmetic quaternion |
228 |
|
|
propagation. With the same time step, a 1024-molecule water simulation |
229 |
|
|
incurs approximately a 10\% increase in computation time using the |
230 |
|
|
{\sc dlm} method in place of quaternions. This cost is more than |
231 |
|
|
justified when comparing the energy conservation achieved by the two |
232 |
|
|
methods. Figure \ref{fig:quatdlm} provides a comparative analysis of |
233 |
|
|
the {\sc dlm} method versus the traditional quaternion scheme. |
234 |
|
|
|
235 |
|
|
\begin{figure} |
236 |
|
|
\centering |
237 |
|
|
\includegraphics[width=3.5in]{./figures/dlmVsQuat.pdf} |
238 |
|
|
\caption[Energy conservation analysis of the {\sc dlm} and quaternion |
239 |
|
|
integration methods]{Analysis of the energy conservation of the {\sc |
240 |
|
|
dlm} and quaternion integration methods. $\delta \mathrm{E}_1$ is the |
241 |
|
|
linear drift in energy over time and $\delta \mathrm{E}_0$ is the |
242 |
|
|
standard deviation of energy fluctuations around this drift. All |
243 |
|
|
simulations were of a 1024 SSD water system at 298 K starting from the |
244 |
|
|
same initial configuration. Note that the {\sc dlm} method provides |
245 |
|
|
more than an order-of-magnitude improvement in both the energy drift |
246 |
|
|
and the size of the energy fluctuations when compared with the |
247 |
|
|
quaternion method at any given time step. At time steps larger than 4 |
248 |
|
|
fs, the quaternion scheme resulted in rapidly rising energies which |
249 |
|
|
eventually lead to simulation failure. Using the {\sc dlm} method, |
250 |
|
|
time steps up to 8 fs can be taken before this behavior is evident.} |
251 |
|
|
\label{fig:quatdlm} |
252 |
|
|
\end{figure} |
253 |
|
|
|
254 |
|
|
In figure \ref{fig:quatdlm}, $\delta \mbox{E}_1$ is a measure of the |
255 |
|
|
linear energy drift in units of $\mbox{kcal mol}^{-1}$ per particle |
256 |
|
|
over a nanosecond of simulation time, and $\delta \mbox{E}_0$ is the |
257 |
|
|
standard deviation of the energy fluctuations in units of $\mbox{kcal |
258 |
|
|
mol}^{-1}$ per particle. In the top plot, it is apparent that the |
259 |
|
|
energy drift is reduced by a significant amount (2 to 3 orders of |
260 |
|
|
magnitude improvement at all tested time steps) by choosing the {\sc |
261 |
|
|
dlm} method over the simple non-symplectic quaternion integration |
262 |
|
|
method. In addition to this improvement in energy drift, the |
263 |
|
|
fluctuations in the total energy are also dampened by 1 to 2 orders of |
264 |
|
|
magnitude by utilizing the {\sc dlm} method. |
265 |
|
|
|
266 |
|
|
\begin{figure} |
267 |
|
|
\centering |
268 |
|
|
\includegraphics[width=\linewidth]{./figures/compCost.pdf} |
269 |
|
|
\caption[Energy drift as a function of required simulation run |
270 |
|
|
time]{Energy drift as a function of required simulation run time. |
271 |
|
|
$\delta \mathrm{E}_1$ is the linear drift in energy over time. |
272 |
|
|
Simulations were performed on a single 2.5 GHz Pentium 4 |
273 |
|
|
processor. Simulation time comparisons can be made by tracing |
274 |
|
|
horizontally from one curve to the other. For example, a simulation |
275 |
|
|
that takes 24 hours using the {\sc dlm} method will take roughly |
276 |
|
|
210 hours using the simple quaternion method if the same degree of |
277 |
|
|
energy conservation is desired.} |
278 |
|
|
\label{fig:cpuCost} |
279 |
|
|
\end{figure} |
280 |
|
|
Although the {\sc dlm} method is more computationally expensive than |
281 |
|
|
the traditional quaternion scheme for propagating of a time step, |
282 |
|
|
consideration of the computational cost for a long simulation with a |
283 |
|
|
particular level of energy conservation is in order. A plot of energy |
284 |
|
|
drift versus computational cost was generated |
285 |
|
|
(Fig.~\ref{fig:cpuCost}). This figure provides an estimate of the CPU |
286 |
|
|
time required under the two integration schemes for 1 nanosecond of |
287 |
|
|
simulation time for the model 1024-molecule system. By choosing a |
288 |
|
|
desired energy drift value it is possible to determine the CPU time |
289 |
|
|
required for both methods. If a $\delta \mbox{E}_1$ of |
290 |
|
|
0.001~kcal~mol$^{-1}$ per particle is desired, a nanosecond of |
291 |
|
|
simulation time will require ~19 hours of CPU time with the {\sc dlm} |
292 |
|
|
integrator, while the quaternion scheme will require ~154 hours of CPU |
293 |
|
|
time. This demonstrates the computational advantage of the integration |
294 |
|
|
scheme utilized in {\sc oopse}. |
295 |
|
|
|
296 |
|
|
\subsection{Periodic Boundary Conditions} |
297 |
|
|
|
298 |
|
|
\begin{figure} |
299 |
|
|
\centering |
300 |
|
|
\includegraphics[width=4.5in]{./figures/periodicImage.pdf} |
301 |
|
|
\caption{With periodic boundary conditions imposed, when particles |
302 |
|
|
move out of one side the simulation box, they wrap back in the |
303 |
|
|
opposite side. In this manner, a finite system of particles behaves as |
304 |
|
|
an infinite system.} |
305 |
|
|
\label{fig:periodicImage} |
306 |
|
|
\end{figure} |
307 |
|
|
|
308 |
|
|
\section{Accumulating Interactions} |
309 |
|
|
|
310 |
|
|
In the force calculation between {\tt moveA} and {\tt moveB} mentioned |
311 |
|
|
in section \ref{sec:IntroIntegrate}, we need to accumulate the |
312 |
|
|
potential and forces (and torques if the particle is a rigid body or |
313 |
|
|
multipole) on each particle from their surroundings. This can quickly |
314 |
|
|
become a cumbersome task for large systems since the number of pair |
315 |
|
|
interactions increases by $\mathcal{O}(N(N-1)/2)$ if you accumulate |
316 |
|
|
interactions between all particles in the system, note the utilization |
317 |
|
|
of Newton's third law to reduce the interaction number from |
318 |
|
|
$\mathcal{O}(N^2)$. The case of periodic boundary conditions further |
319 |
|
|
complicates matters by turning the finite system into an infinitely |
320 |
|
|
repeating one. Fortunately, we can reduce the scale of this problem by |
321 |
|
|
using spherical cutoff methods. |
322 |
|
|
|
323 |
|
|
\begin{figure} |
324 |
|
|
\centering |
325 |
|
|
\includegraphics[width=3.5in]{./figures/sphericalCut.pdf} |
326 |
|
|
\caption{When using a spherical cutoff, only particles within a chosen |
327 |
|
|
cutoff radius distance, $R_\textrm{c}$, of the central particle are |
328 |
|
|
included in the pairwise summation. This reduces a problem that scales |
329 |
|
|
by $\sim\mathcal{O}(N^2)$ to one that scales by $\sim\mathcal{O}(N)$.} |
330 |
|
|
\label{fig:sphereCut} |
331 |
|
|
\end{figure} |
332 |
|
|
With spherical cutoffs, rather than accumulating the full set of |
333 |
|
|
interactions between all particles in the simulation, we only |
334 |
|
|
explicitly consider interactions between local particles out to a |
335 |
|
|
specified cutoff radius distance, $R_\textrm{c}$, (see figure |
336 |
|
|
\ref{fig:sphereCut}). This reduces the scaling of the interaction to |
337 |
|
|
$\mathcal{O}(N\cdot\textrm{c})$, where `c' is a value that depends on |
338 |
|
|
the size of $R_\textrm{c}$ (c $\approx R_\textrm{c}^3$). Determination |
339 |
|
|
of which particles are within the cutoff is also an issue, because |
340 |
|
|
this process requires a full loop over all $N(N-1)/2$ pairs. To reduce |
341 |
|
|
the this expense, we can use neighbor lists.\cite{Verlet67,Thompson83} |
342 |
|
|
With neighbor lists, we have a second list of particles built from a |
343 |
|
|
list radius $R_\textrm{l}$, which is larger than $R_\textrm{c}$. Once |
344 |
|
|
any of the particles within $R_\textrm{l}$ move a distance of |
345 |
|
|
$R_\textrm{l}-R_\textrm{c}$ (the ``skin'' thickness), we rebuild the |
346 |
|
|
list with the full $N(N-1)/2$ loop.\cite{Verlet67} With an appropriate |
347 |
|
|
skin thickness, these updates are only performed every $\sim$20 |
348 |
|
|
time steps, significantly reducing the time spent on pair-list |
349 |
|
|
bookkeeping operations.\cite{Allen87} If these neighbor lists are |
350 |
|
|
utilized, it is important that these list updates occur |
351 |
|
|
regularly. Incorrect application of this technique leads to |
352 |
|
|
non-physical dynamics, such as the ``flying block of ice'' behavior |
353 |
|
|
for which improper neighbor list handling was identified a one of the |
354 |
|
|
possible causes.\cite{Harvey98,Sagui99} |
355 |
|
|
|
356 |
|
|
\subsection{Correcting Cutoff Discontinuities} |
357 |
|
|
\begin{figure} |
358 |
|
|
\centering |
359 |
|
|
\includegraphics[width=3.5in]{./figures/ljCutoffSquare.pdf} |
360 |
|
|
\caption{The common methods to smooth the potential discontinuity |
361 |
|
|
introduced when using a cutoff include a shifted potential, a shifted |
362 |
|
|
force, and a switching function. The shifted potential and shifted |
363 |
|
|
force both lift the whole potential so that it zeroes at |
364 |
|
|
$R_\textrm{c}$, thereby reducing the strength of the interaction. The |
365 |
|
|
(cubic) switching function only alters the potential in the switching |
366 |
|
|
region in order to smooth out the discontinuity.} |
367 |
|
|
\label{fig:ljCutoff} |
368 |
|
|
\end{figure} |
369 |
|
|
As particle move in and out of $R_\textrm{c}$, there will be sudden |
370 |
|
|
discontinuous jumps in the potential (and forces) due to their |
371 |
|
|
appearance and disappearance. In order to prevent heating and poor |
372 |
|
|
energy conservation in the simulations, this discontinuity needs to be |
373 |
|
|
smoothed out. There are several ways to modify the function so that it |
374 |
|
|
crosses $R_\textrm{c}$ in a continuous fashion, and the easiest |
375 |
|
|
methods is shifting the potential. To shift the potential, we simply |
376 |
|
|
subtract out the value we calculate at $R_\textrm{c}$ from the whole |
377 |
|
|
potential. For the shifted form of the Lennard-Jones potential is |
378 |
|
|
\begin{equation} |
379 |
|
|
V_\textrm{SLJ} = \left\{\begin{array}{l@{\quad\quad}l} |
380 |
|
|
V_\textrm{LJ}(r_{ij}) - V_\textrm{LJ}(R_\textrm{c}) & r_{ij} < R_\textrm{c}, \\ |
381 |
|
|
0 & r_{ij} \geqslant R_\textrm{c}, |
382 |
|
|
\end{array}\right. |
383 |
|
|
\end{equation} |
384 |
|
|
where |
385 |
|
|
\begin{equation} |
386 |
|
|
V_\textrm{LJ} = 4\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12} - |
387 |
|
|
\left(\frac{\sigma}{r_{ij}}\right)^6\right]. |
388 |
|
|
\end{equation} |
389 |
|
|
In figure \ref{fig:ljCutoff}, the shifted form of the potential |
390 |
|
|
reaches zero at the cutoff radius at the expense of the correct |
391 |
|
|
magnitude below $R_\textrm{c}$. This correction method also does |
392 |
|
|
nothing to correct the discontinuity in the forces. We can account for |
393 |
|
|
this force discontinuity by shifting in the by applying the shifting |
394 |
|
|
in the forces rather than just the potential via |
395 |
|
|
\begin{equation} |
396 |
|
|
V_\textrm{SFLJ} = \left\{\begin{array}{l@{\quad\quad}l} |
397 |
|
|
V_\textrm{LJ}({r_{ij}}) - V_\textrm{LJ}(R_\textrm{c}) - |
398 |
|
|
\left(\frac{d V(r_{ij})}{d r_{ij}}\right)_{r_{ij}=R_\textrm{c}} |
399 |
|
|
(r_{ij} - R_\textrm{c}) & r_{ij} < R_\textrm{c}, \\ |
400 |
|
|
0 & r_{ij} \geqslant R_\textrm{c}. |
401 |
|
|
\end{array}\right. |
402 |
|
|
\end{equation} |
403 |
|
|
The forces are continuous with this potential; however, the inclusion |
404 |
|
|
of the derivative term distorts the potential even further than the |
405 |
|
|
simple shifting as shown in figure \ref{fig:ljCutoff}. The method for |
406 |
|
|
correcting the discontinuity which results in the smallest |
407 |
|
|
perturbation in both the potential and the forces is the use of a |
408 |
|
|
switching function. The cubic switching function, |
409 |
|
|
\begin{equation} |
410 |
|
|
S(r) = \left\{\begin{array}{l@{\quad\quad}l} |
411 |
|
|
1 & r_{ij} \leqslant R_\textrm{sw}, \\ |
412 |
|
|
\frac{(R_\textrm{c} + 2r_{ij} - 3R_\textrm{sw}) |
413 |
|
|
(R_\textrm{c} - r_{ij})^2}{(R_\textrm{c} - R_\textrm{sw})^3} |
414 |
|
|
& R_\textrm{sw} < r_{ij} \leqslant R_\textrm{c}, \\ |
415 |
|
|
0 & r_{ij} > R_\textrm{c}, |
416 |
|
|
\end{array}\right. |
417 |
|
|
\end{equation} |
418 |
|
|
is sufficient to smooth the potential (again, see figure |
419 |
|
|
\ref{fig:ljCutoff}) and the forces by only perturbing the potential in |
420 |
|
|
the switching region. If smooth second derivatives are required, a |
421 |
|
|
higher order polynomial switching function (i.e. fifth order |
422 |
|
|
polynomial) can be used.\cite{Andrea83,Leach01} It should be noted |
423 |
|
|
that the higher the order of the polynomial, the more abrupt the |
424 |
|
|
switching transition. |
425 |
|
|
|
426 |
|
|
\subsection{Long-Range Interaction Corrections} |
427 |
|
|
|
428 |
|
|
While a good approximation, accumulating interaction only from the |
429 |
|
|
nearby particles can lead to some issues, because the further away |
430 |
|
|
surrounding particles do still have a small effect. For instance, |
431 |
|
|
while the strength of the Lennard-Jones interaction is quite weak at |
432 |
|
|
$R_\textrm{c}$ in figure \ref{fig:ljCutoff}, we are discarding all of |
433 |
|
|
the attractive interaction that extends out to extremely long |
434 |
|
|
distances. Thus, the potential is a little too high and the pressure |
435 |
|
|
on the central particle from the surroundings is a little too low. For |
436 |
|
|
homogeneous Lennard-Jones systems, we can correct for this neglect by |
437 |
|
|
assuming a uniform density and integrating the missing part, |
438 |
|
|
\begin{equation} |
439 |
|
|
V_\textrm{full}(r_{ij}) \approx V_\textrm{c} |
440 |
|
|
+ 2\pi N\rho\int^\infty_{R_\textrm{c}}r^2V_\textrm{LJ}(r)dr, |
441 |
|
|
\end{equation} |
442 |
|
|
where $V_\textrm{c}$ is the truncated Lennard-Jones |
443 |
|
|
potential.\cite{Allen87} Like the potential, other Lennard-Jones |
444 |
|
|
properties can be corrected by integration over the relevant |
445 |
|
|
function. Note that with heterogeneous systems, this correction begins |
446 |
|
|
to break down because the density is no longer uniform. |
447 |
|
|
|
448 |
|
|
Correcting long-range electrostatic interactions is a topic of great |
449 |
|
|
importance in the field of molecular simulations. There have been |
450 |
|
|
several techniques developed to address this issue, like the Ewald |
451 |
|
|
summation and the reaction field technique. An in-depth analysis of |
452 |
|
|
this problem, as well as useful corrective techniques, is presented in |
453 |
|
|
chapter \ref{chap:electrostatics}. |
454 |
|
|
|
455 |
|
|
\section{Measuring Properties} |
456 |
|
|
|
457 |
|
|
\section{Application of the Techniques} |
458 |
|
|
|
459 |
|
|
In the following chapters, the above techniques are all utilized in |
460 |
|
|
the study of molecular systems. There are a number of excellent |
461 |
|
|
simulation packages available, both free and commercial, which |
462 |
|
|
incorporate many of these |
463 |
|
|
methods.\cite{Brooks83,MacKerell98,Pearlman95,Berendsen95,Lindahl01,Smith96,Ponder87} |
464 |
|
|
Because of our interest in rigid body dynamics, point multipoles, and |
465 |
|
|
systems where the orientational degrees cannot be handled using the |
466 |
|
|
common constraint procedures (like {\sc shake}), the majority of the |
467 |
|
|
following work was performed using {\sc oopse}, the object-oriented |
468 |
|
|
parallel simulation engine.\cite{Meineke05} The {\sc oopse} package |
469 |
|
|
started out as a collection of separate programs written within our |
470 |
|
|
group, and has developed into one of the few parallel molecular |
471 |
|
|
dynamics packages capable of accurately integrating rigid bodies and |
472 |
|
|
point multipoles. |
473 |
|
|
|
474 |
|
|
In chapter \ref{chap:electrostatics}, we investigate correction |
475 |
|
|
techniques for proper handling of long-ranged electrostatic |
476 |
|
|
interactions. In particular we develop and evaluate some new pairwise |
477 |
|
|
methods which we have incorporated into {\sc oopse}. These techniques |
478 |
|
|
make an appearance in the later chapters, as they are applied to |
479 |
|
|
specific systems of interest, showing how it they can improve the |
480 |
|
|
quality of various molecular simulations. |
481 |
|
|
|
482 |
|
|
In chapter \ref{chap:water}, we focus on simple water models, |
483 |
|
|
specifically the single-point soft sticky dipole (SSD) family of water |
484 |
|
|
models. These single-point models are more efficient than the common |
485 |
|
|
multi-point partial charge models and, in many cases, better capture |
486 |
|
|
the dynamic properties of water. We discuss improvements to these |
487 |
|
|
models in regards to long-range electrostatic corrections and show |
488 |
|
|
that these models can work well with the techniques discussed in |
489 |
|
|
chapter \ref{chap:electrostatics}. By investigating and improving |
490 |
|
|
simple water models, we are extending the limits of computational |
491 |
|
|
efficiency for systems that heavily depend on water calculations. |
492 |
|
|
|
493 |
|
|
In chapter \ref{chap:ice}, we study a unique polymorph of ice that we |
494 |
|
|
discovered while performing water simulations with the fast simple |
495 |
|
|
water models discussed in the previous chapter. This form of ice, |
496 |
|
|
which we called ``imaginary ice'' (Ice-$i$), has a low-density |
497 |
|
|
structure which is different from any known polymorph from either |
498 |
|
|
experiment or other simulations. In this study, we perform a free |
499 |
|
|
energy analysis and see that this structure is in fact the |
500 |
|
|
thermodynamically preferred form of ice for both the single-point and |
501 |
|
|
commonly used multi-point water models under the chosen simulation |
502 |
|
|
conditions. We also consider electrostatic corrections, again |
503 |
|
|
including the techniques discussed in chapter |
504 |
|
|
\ref{chap:electrostatics}, to see how they influence the free energy |
505 |
|
|
results. This work, to some degree, addresses the appropriateness of |
506 |
|
|
using these simplistic water models outside of the conditions for |
507 |
|
|
which they were developed. |
508 |
|
|
|
509 |
|
|
In the final chapter we summarize the work presented previously. We |
510 |
|
|
also close with some final comments before the supporting information |
511 |
|
|
presented in the appendices. |