| 1 |
|
| 2 |
\section{\label{sec:mechanics}Mechanics} |
| 3 |
|
| 4 |
\subsection{\label{integrate}Integrating the Equations of Motion: the Symplectic Step Integrator} |
| 5 |
|
| 6 |
Integration of the equations of motion was carried out using the |
| 7 |
symplectic splitting method proposed by Dullweber \emph{et |
| 8 |
al.}.\cite{Dullweber1997} The reason for this integrator selection |
| 9 |
deals with poor energy conservation of rigid body systems using |
| 10 |
quaternions. While quaternions work well for orientational motion in |
| 11 |
alternate ensembles, the microcanonical ensemble has a constant energy |
| 12 |
requirement that is quite sensitive to errors in the equations of |
| 13 |
motion. The original implementation of this code utilized quaternions |
| 14 |
for rotational motion propagation; however, a detailed investigation |
| 15 |
showed that they resulted in a steady drift in the total energy, |
| 16 |
something that has been observed by others.\cite{Laird97} |
| 17 |
|
| 18 |
The key difference in the integration method proposed by Dullweber |
| 19 |
\emph{et al.} is that the entire rotation matrix is propagated from |
| 20 |
one time step to the next. In the past, this would not have been as |
| 21 |
feasible a option, being that the rotation matrix for a single body is |
| 22 |
nine elements long as opposed to 3 or 4 elements for Euler angles and |
| 23 |
quaternions respectively. System memory has become much less of an |
| 24 |
issue in recent times, and this has resulted in substantial benefits |
| 25 |
in energy conservation. There is still the issue of 5 or 6 additional |
| 26 |
elements for describing the orientation of each particle, which will |
| 27 |
increase dump files substantially. Simply translating the rotation |
| 28 |
matrix into its component Euler angles or quaternions for storage |
| 29 |
purposes relieves this burden. |
| 30 |
|
| 31 |
The symplectic splitting method allows for Verlet style integration of |
| 32 |
both linear and angular motion of rigid bodies. In the integration |
| 33 |
method, the orientational propagation involves a sequence of matrix |
| 34 |
evaluations to update the rotation matrix.\cite{Dullweber1997} These |
| 35 |
matrix rotations end up being more costly computationally than the |
| 36 |
simpler arithmetic quaternion propagation. With the same time step, a |
| 37 |
1000 SSD particle simulation shows an average 7\% increase in |
| 38 |
computation time using the symplectic step method in place of |
| 39 |
quaternions. This cost is more than justified when comparing the |
| 40 |
energy conservation of the two methods as illustrated in figure |
| 41 |
\ref{timestep}. |
| 42 |
|
| 43 |
\begin{figure} |
| 44 |
\epsfxsize=6in |
| 45 |
\epsfbox{timeStep.epsi} |
| 46 |
\caption{Energy conservation using quaternion based integration versus |
| 47 |
the symplectic step method proposed by Dullweber \emph{et al.} with |
| 48 |
increasing time step. For each time step, the dotted line is total |
| 49 |
energy using the symplectic step integrator, and the solid line comes |
| 50 |
from the quaternion integrator. The larger time step plots are shifted |
| 51 |
up from the true energy baseline for clarity.} |
| 52 |
\label{timestep} |
| 53 |
\end{figure} |
| 54 |
|
| 55 |
In figure \ref{timestep}, the resulting energy drift at various time |
| 56 |
steps for both the symplectic step and quaternion integration schemes |
| 57 |
is compared. All of the 1000 SSD particle simulations started with the |
| 58 |
same configuration, and the only difference was the method for |
| 59 |
handling rotational motion. At time steps of 0.1 and 0.5 fs, both |
| 60 |
methods for propagating particle rotation conserve energy fairly well, |
| 61 |
with the quaternion method showing a slight energy drift over time in |
| 62 |
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the |
| 63 |
energy conservation benefits of the symplectic step method are clearly |
| 64 |
demonstrated. Thus, while maintaining the same degree of energy |
| 65 |
conservation, one can take considerably longer time steps, leading to |
| 66 |
an overall reduction in computation time. |
| 67 |
|
| 68 |
Energy drift in these SSD particle simulations was unnoticeable for |
| 69 |
time steps up to three femtoseconds. A slight energy drift on the |
| 70 |
order of 0.012 kcal/mol per nanosecond was observed at a time step of |
| 71 |
four femtoseconds, and as expected, this drift increases dramatically |
| 72 |
with increasing time step. To insure accuracy in the constant energy |
| 73 |
simulations, time steps were set at 2 fs and kept at this value for |
| 74 |
constant pressure simulations as well. |
| 75 |
|
| 76 |
|
| 77 |
\subsection{\label{sec:extended}Extended Systems for other Ensembles} |
| 78 |
|
| 79 |
|
| 80 |
{\sc oopse} implements a |
| 81 |
|
| 82 |
|
| 83 |
\subsubsection{\label{sec:noseHooverThermo}Nose-Hoover Thermostatting} |
| 84 |
|
| 85 |
To mimic the effects of being in a constant temperature ({\sc nvt}) |
| 86 |
ensemble, {\sc oopse} uses the Nose-Hoover extended system |
| 87 |
approach.\cite{Hoover85} In this method, the equations of motion for |
| 88 |
the particle positions and velocities are |
| 89 |
\begin{eqnarray} |
| 90 |
\dot{{\bf r}} & = & {\bf v} \\ |
| 91 |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} |
| 92 |
\label{eq:nosehoovereom} |
| 93 |
\end{eqnarray} |
| 94 |
|
| 95 |
$\chi$ is an ``extra'' variable included in the extended system, and |
| 96 |
it is propagated using the first order equation of motion |
| 97 |
\begin{equation} |
| 98 |
\dot{\chi} = \frac{1}{\tau_{T}} \left( \frac{T}{T_{target}} - 1 \right) |
| 99 |
\label{eq:nosehooverext} |
| 100 |
\end{equation} |
| 101 |
where $T_{target}$ is the target temperature for the simulation, and |
| 102 |
$\tau_{T}$ is a time constant for the thermostat. |
| 103 |
|
| 104 |
To select the Nose-Hoover {\sc nvt} ensemble, the {\tt ensemble = NVT;} |
| 105 |
command would be used in the simulation's {\sc bass} file. There is |
| 106 |
some subtlety in choosing values for $\tau_{T}$, and it is usually set |
| 107 |
to values of a few ps. Within a {\sc bass} file, $\tau_{T}$ could be |
| 108 |
set to 1 ps using the {\tt tauThermostat = 1000; } command. |
| 109 |
|
| 110 |
|
| 111 |
\subsection{\label{Sec:zcons}Z-Constraint Method} |
| 112 |
|
| 113 |
Based on fluctuatin-dissipation theorem,\bigskip\ force auto-correlation |
| 114 |
method was developed to investigate the dynamics of ions inside the ion |
| 115 |
channels.\cite{Roux91} Time-dependent friction coefficient can be calculated |
| 116 |
from the deviation of the instaneous force from its mean force. |
| 117 |
|
| 118 |
% |
| 119 |
|
| 120 |
\begin{equation} |
| 121 |
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T |
| 122 |
\end{equation} |
| 123 |
|
| 124 |
|
| 125 |
where% |
| 126 |
\begin{equation} |
| 127 |
\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle |
| 128 |
\end{equation} |
| 129 |
|
| 130 |
|
| 131 |
If the time-dependent friction decay rapidly, static friction coefficient can |
| 132 |
be approximated by% |
| 133 |
|
| 134 |
\begin{equation} |
| 135 |
\xi^{static}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt |
| 136 |
\end{equation} |
| 137 |
|
| 138 |
|
| 139 |
Hence, diffusion constant can be estimated by |
| 140 |
\begin{equation} |
| 141 |
D(z)=\frac{k_{B}T}{\xi^{static}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty |
| 142 |
}\langle\delta F(z,t)\delta F(z,0)\rangle dt}% |
| 143 |
\end{equation} |
| 144 |
|
| 145 |
|
| 146 |
\bigskip Z-Constraint method, which fixed the z coordinates of the molecules |
| 147 |
with respect to the center of the mass of the system, was proposed to obtain |
| 148 |
the forces required in force auto-correlation method.\cite{Marrink94} However, |
| 149 |
simply resetting the coordinate will move the center of the mass of the whole |
| 150 |
system. To avoid this problem, a new method was used at {\sc oopse}. Instead of |
| 151 |
resetting the coordinate, we reset the forces of z-constraint molecules as |
| 152 |
well as subtract the total constraint forces from the rest of the system after |
| 153 |
force calculation at each time step. |
| 154 |
\begin{verbatim} |
| 155 |
$F_{\alpha i}=0$ |
| 156 |
$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}$ |
| 157 |
$F_{\alpha i}=F_{\alpha i}-\frac{M_{_{\alpha i}}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}\sum\limits_{\beta}F_{\beta}$ |
| 158 |
$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}$ |
| 159 |
\end{verbatim} |
| 160 |
|
| 161 |
At the very beginning of the simulation, the molecules may not be at its |
| 162 |
constraint position. To move the z-constraint molecule to the specified |
| 163 |
position, a simple harmonic potential is used% |
| 164 |
|
| 165 |
\begin{equation} |
| 166 |
U(t)=\frac{1}{2}k_{Harmonic}(z(t)-z_{cons})^{2}% |
| 167 |
\end{equation} |
| 168 |
where $k_{Harmonic}$\bigskip\ is the harmonic force constant, $z(t)$ is |
| 169 |
current z coordinate of the center of mass of the z-constraint molecule, and |
| 170 |
$z_{cons}$ is the restraint position. Therefore, the harmonic force operated |
| 171 |
on the z-constraint molecule at time $t$ can be calculated by% |
| 172 |
\begin{equation} |
| 173 |
F_{z_{Harmonic}}(t)=-\frac{\partial U(t)}{\partial z(t)}=-k_{Harmonic}% |
| 174 |
(z(t)-z_{cons}) |
| 175 |
\end{equation} |
| 176 |
Worthy of mention, other kinds of potential functions can also be used to |
| 177 |
drive the z-constraint molecule. |
