| 1 |
|
| 2 |
\section{\label{sec:mechanics}Mechanics} |
| 3 |
|
| 4 |
\subsection{\label{integrate}Integrating the Equations of Motion: the |
| 5 |
DLM method} |
| 6 |
|
| 7 |
The default method for integrating the equations of motion in {\sc |
| 8 |
oopse} is a velocity-Verlet version of the symplectic splitting method |
| 9 |
proposed by Dullweber, Leimkuhler and McLachlan |
| 10 |
(DLM).\cite{Dullweber1997} When there are no directional atoms or |
| 11 |
rigid bodies present in the simulation, this integrator becomes the |
| 12 |
standard velocity-Verlet integrator which is known to sample the |
| 13 |
microcanonical (NVE) ensemble.\cite{} |
| 14 |
|
| 15 |
Previous integration methods for orientational motion have problems |
| 16 |
that are avoided in the DLM method. Direct propagation of the Euler |
| 17 |
angles has a known $1/\sin\theta$ divergence in the equations of |
| 18 |
motion for $\phi$ and $\psi$,\cite{AllenTildesley} leading to |
| 19 |
numerical instabilities any time one of the directional atoms or rigid |
| 20 |
bodies has an orientation near $\theta=0$ or $\theta=\pi$. More |
| 21 |
modern quaternion-based integration methods have relatively poor |
| 22 |
energy conservation. While quaternions work well for orientational |
| 23 |
motion in other ensembles, the microcanonical ensemble has a |
| 24 |
constant energy requirement that is quite sensitive to errors in the |
| 25 |
equations of motion. An earlier implementation of {\sc oopse} |
| 26 |
utilized quaternions for propagation of rotational motion; however, a |
| 27 |
detailed investigation showed that they resulted in a steady drift in |
| 28 |
the total energy, something that has been observed by |
| 29 |
others.\cite{Laird97} |
| 30 |
|
| 31 |
The key difference in the integration method proposed by Dullweber |
| 32 |
\emph{et al.} is that the entire rotation matrix is propagated from |
| 33 |
one time step to the next. In the past, this would not have been |
| 34 |
feasible, since that the rotation matrix for a single body has nine |
| 35 |
elements compared to 3 or 4 elements for Euler angles and quaternions |
| 36 |
respectively. System memory has become less costly in recent times, |
| 37 |
and this can be translated into substantial benefits in energy |
| 38 |
conservation. |
| 39 |
|
| 40 |
The basic equations of motion being integrated are derived from the |
| 41 |
Hamiltonian for conservative systems containing rigid bodies, |
| 42 |
\begin{equation} |
| 43 |
H = \sum_{i} \frac{{\bf v}_i^T m_i {\bf v}_i}{2} + \sum_i |
| 44 |
\frac{{\bf j}_i^T \cdot {\bf j}_i}{2 \overleftrightarrow{\mathsf{I}}_i} + |
| 45 |
V\left(\left\{{\bf r}\right\}, \left\{\mathsf{A}\right\}\right) |
| 46 |
\end{equation} |
| 47 |
where ${\bf r}_i$ and ${\bf v}_i$ are the cartesian position vector |
| 48 |
and velocity of the center of mass of particle $i$, and ${\bf j}_i$ |
| 49 |
and $\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular |
| 50 |
momentum and moment of inertia tensor, respectively. $\mathsf{A}_i$ |
| 51 |
is the $3 \times 3$ rotation matrix describing the instantaneous |
| 52 |
orientation of the particle. $V$ is the potential energy function |
| 53 |
which may depend on both the positions $\left\{{\bf r}\right\}$ and |
| 54 |
orientations $\left\{\mathsf{A}\right\}$ of all particles. The |
| 55 |
equations of motion for the particle centers of mass are derived from |
| 56 |
Hamilton's equations and are quite simple, |
| 57 |
\begin{eqnarray} |
| 58 |
\dot{{\bf r}} & = & {\bf v} \\ |
| 59 |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} |
| 60 |
\end{eqnarray} |
| 61 |
where ${\bf f}$ is the instantaneous force on the center of mass |
| 62 |
of the particle, |
| 63 |
\begin{equation} |
| 64 |
{\bf f} = \frac{\partial}{\partial |
| 65 |
{\bf r}} V(\left\{{\bf r}(t)\right\}, \left\{\mathsf{A}(t)\right\}). |
| 66 |
\end{equation} |
| 67 |
|
| 68 |
|
| 69 |
The equations of motion for the orientational degrees of freedom are |
| 70 |
\begin{eqnarray} |
| 71 |
\dot{\mathsf{A}} & = & \mathsf{A} |
| 72 |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\ |
| 73 |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 74 |
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \frac{\partial |
| 75 |
V}{\partial \mathsf{A}} \right) |
| 76 |
\end{eqnarray} |
| 77 |
In these equations of motion, the $\mbox{skew}$ matrix of a vector |
| 78 |
${\bf v} = \left( v_1, v_2, v_3 \right)$ is defined: |
| 79 |
\begin{equation} |
| 80 |
\mbox{skew}\left( {\bf v} \right) := \left( |
| 81 |
\begin{array}{ccc} |
| 82 |
0 & v_3 & - v_2 \\ |
| 83 |
-v_3 & 0 & v_1 \\ |
| 84 |
v_2 & -v_1 & 0 |
| 85 |
\end{array} |
| 86 |
\right) |
| 87 |
\end{equation} |
| 88 |
The $\mbox{rot}$ notation refers to the mapping of the $3 \times 3$ |
| 89 |
rotation matrix to a vector of orientations by first computing the |
| 90 |
skew-symmetric part $\left(\mathsf{A} - \mathsf{A}^{T}\right)$ and |
| 91 |
then associating this with a length 3 vector by inverting the |
| 92 |
$\mbox{skew}$ function above: |
| 93 |
\begin{equation} |
| 94 |
\mbox{rot}\left(\mathsf{A}\right) := \mbox{ skew}^{-1}\left(\mathsf{A} |
| 95 |
- \mathsf{A}^{T} \right) |
| 96 |
\end{equation} |
| 97 |
|
| 98 |
|
| 99 |
|
| 100 |
In the integration method, the |
| 101 |
orientational propagation involves a sequence of matrix evaluations to |
| 102 |
update the rotation matrix.\cite{Dullweber1997} These matrix rotations |
| 103 |
end up being more costly computationally than the simpler arithmetic |
| 104 |
quaternion propagation. With the same time step, a 1000 SSD particle |
| 105 |
simulation shows an average 7\% increase in computation time using the |
| 106 |
symplectic step method in place of quaternions. This cost is more than |
| 107 |
justified when comparing the energy conservation of the two methods as |
| 108 |
illustrated in figure |
| 109 |
\ref{timestep}. |
| 110 |
|
| 111 |
\begin{figure} |
| 112 |
\epsfxsize=6in |
| 113 |
\epsfbox{timeStep.epsi} |
| 114 |
\caption{Energy conservation using quaternion based integration versus |
| 115 |
the symplectic step method proposed by Dullweber \emph{et al.} with |
| 116 |
increasing time step. For each time step, the dotted line is total |
| 117 |
energy using the symplectic step integrator, and the solid line comes |
| 118 |
from the quaternion integrator. The larger time step plots are shifted |
| 119 |
up from the true energy baseline for clarity.} |
| 120 |
\label{timestep} |
| 121 |
\end{figure} |
| 122 |
|
| 123 |
In figure \ref{timestep}, the resulting energy drift at various time |
| 124 |
steps for both the symplectic step and quaternion integration schemes |
| 125 |
is compared. All of the 1000 SSD particle simulations started with the |
| 126 |
same configuration, and the only difference was the method for |
| 127 |
handling rotational motion. At time steps of 0.1 and 0.5 fs, both |
| 128 |
methods for propagating particle rotation conserve energy fairly well, |
| 129 |
with the quaternion method showing a slight energy drift over time in |
| 130 |
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the |
| 131 |
energy conservation benefits of the symplectic step method are clearly |
| 132 |
demonstrated. Thus, while maintaining the same degree of energy |
| 133 |
conservation, one can take considerably longer time steps, leading to |
| 134 |
an overall reduction in computation time. |
| 135 |
|
| 136 |
Energy drift in these SSD particle simulations was unnoticeable for |
| 137 |
time steps up to three femtoseconds. A slight energy drift on the |
| 138 |
order of 0.012 kcal/mol per nanosecond was observed at a time step of |
| 139 |
four femtoseconds, and as expected, this drift increases dramatically |
| 140 |
with increasing time step. To insure accuracy in the constant energy |
| 141 |
simulations, time steps were set at 2 fs and kept at this value for |
| 142 |
constant pressure simulations as well. |
| 143 |
|
| 144 |
|
| 145 |
\subsection{\label{sec:extended}Extended Systems for other Ensembles} |
| 146 |
|
| 147 |
{\sc oopse} implements a number of extended system integrators for |
| 148 |
sampling from other ensembles relevant to chemical physics. The |
| 149 |
integrator can selected with the {\tt ensemble} keyword in the |
| 150 |
{\tt .bass} file: |
| 151 |
|
| 152 |
\begin{center} |
| 153 |
\begin{tabular}{lll} |
| 154 |
{\bf Integrator} & {\bf Ensemble} & {\bf {\tt .bass} line} \\ |
| 155 |
NVE & microcanonical & {\tt ensemble = ``NVE''; } \\ |
| 156 |
NVT & canonical & {\tt ensemble = ``NVT''; } \\ |
| 157 |
NPTi & isobaric-isothermal (with isotropic volume changes) & {\tt |
| 158 |
ensemble = ``NPTi'';} \\ |
| 159 |
NPTf & isobaric-isothermal (with changes to box shape) & {\tt |
| 160 |
ensemble = ``NPTf'';} \\ |
| 161 |
NPTxyz & approximate isobaric-isothermal & {\tt ensemble = |
| 162 |
``NPTxyz'';} \\ |
| 163 |
& (with separate barostats on each box dimension) & |
| 164 |
\end{tabular} |
| 165 |
\end{center} |
| 166 |
|
| 167 |
The relatively well-known Nos\'e-Hoover thermostat is implemented in |
| 168 |
{\sc oopse}'s NVT integrator. This method couples an extra degree of |
| 169 |
freedom (the thermostat) to the kinetic energy of the system, and has |
| 170 |
been shown to sample the canonical distribution in the system degrees |
| 171 |
of freedom while conserving a quantity that is, to within a constant, |
| 172 |
the Helmholtz free energy. |
| 173 |
|
| 174 |
NPT algorithms attempt to maintain constant pressure in the system by |
| 175 |
coupling the volume of the system to a barostat. {\sc oopse} contains |
| 176 |
three different constant pressure algorithms. The first two, NPTi and |
| 177 |
NPTf have been shown to conserve a quantity that is, to within a |
| 178 |
constant, the Gibbs free energy. The Melchionna modification to the |
| 179 |
Hoover barostat is implemented in both NPTi and NPTf. NPTi allows |
| 180 |
only isotropic changes in the simulation box, while box {\it shape} |
| 181 |
variations are allowed in NPTf. The NPTxyz integrator has {\it not} |
| 182 |
been shown to sample from the isobaric-isothermal ensemble. It is |
| 183 |
useful, however, in that it maintains orthogonality for the axes of |
| 184 |
the simulation box while attempting to equalize pressure along the |
| 185 |
three perpendicular directions in the box. |
| 186 |
|
| 187 |
Each of the extended system integrators requires additional keywords |
| 188 |
to set target values for the thermodynamic state variables that are |
| 189 |
being held constant. Keywords are also required to set the |
| 190 |
characteristic decay times for the dynamics of the extended |
| 191 |
variables. |
| 192 |
|
| 193 |
\begin{tabular}{llll} |
| 194 |
{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf |
| 195 |
default value} \\ |
| 196 |
$T_{\mathrm{target}}$ & {\tt targetTemperature = 300;} & K & none \\ |
| 197 |
$P_{\mathrm{target}}$ & {\tt targetPressure = 1;} & atm & none \\ |
| 198 |
$\tau_T$ & {\tt tauThermostat = 1e3;} & fs & none \\ |
| 199 |
$\tau_B$ & {\tt tauBarostat = 5e3;} & fs & none \\ |
| 200 |
& {\tt resetTime = 200;} & fs & none \\ |
| 201 |
& {\tt useInitialExtendedSystemState = ``true'';} & logical & |
| 202 |
false |
| 203 |
\end{tabular} |
| 204 |
|
| 205 |
Two additional keywords can be used to either clear the extended |
| 206 |
system variables periodically ({\tt resetTime}), or to maintain the |
| 207 |
state of the extended system variables between simulations ({\tt |
| 208 |
useInitialExtendedSystemState}). More details on these variables |
| 209 |
and their use in the integrators follows below. |
| 210 |
|
| 211 |
\subsubsection{\label{sec:noseHooverThermo}Nos\'{e}-Hoover Thermostatting} |
| 212 |
|
| 213 |
The Nos\'e-Hoover equations of motion are given by\cite{Hoover85} |
| 214 |
\begin{eqnarray} |
| 215 |
\dot{{\bf r}} & = & {\bf v} \\ |
| 216 |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} \\ |
| 217 |
\dot{\mathsf{A}} & = & \\ |
| 218 |
\dot{{\bf j}} & = & - \chi {\bf j} |
| 219 |
\label{eq:nosehoovereom} |
| 220 |
\end{eqnarray} |
| 221 |
|
| 222 |
$\chi$ is an ``extra'' variable included in the extended system, and |
| 223 |
it is propagated using the first order equation of motion |
| 224 |
\begin{equation} |
| 225 |
\dot{\chi} = \frac{1}{\tau_{T}^2} \left( \frac{T}{T_{\mathrm{target}}} - 1 \right). |
| 226 |
\label{eq:nosehooverext} |
| 227 |
\end{equation} |
| 228 |
|
| 229 |
The instantaneous temperature $T$ is proportional to the total kinetic |
| 230 |
energy (both translational and orientational) and is given by |
| 231 |
\begin{equation} |
| 232 |
T = \frac{2 K}{f k_B} |
| 233 |
\end{equation} |
| 234 |
Here, $f$ is the total number of degrees of freedom in the system, |
| 235 |
\begin{equation} |
| 236 |
f = 3 N + 3 N_{\mathrm{orient}} - N_{\mathrm{constraints}} |
| 237 |
\end{equation} |
| 238 |
and $K$ is the total kinetic energy, |
| 239 |
\begin{equation} |
| 240 |
K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i + |
| 241 |
\sum_{i=1}^{N_{\mathrm{orient}}} \sum_{\alpha=x,y,z} \frac{{\bf |
| 242 |
j}_{i\alpha}^T \cdot {\bf j}_{i\alpha}}{2 |
| 243 |
\overleftrightarrow{\mathsf{I}}_{i,\alpha \alpha}} |
| 244 |
\end{equation} |
| 245 |
|
| 246 |
In eq.(\ref{eq:nosehooverext}), $\tau_T$ is the time constant for |
| 247 |
relaxation of the temperature to the target value. To set values for |
| 248 |
$\tau_T$ or $T_{\mathrm{target}}$ in a simulation, one would use the |
| 249 |
{\tt tauThermostat} and {\tt targetTemperature} keywords in the {\tt |
| 250 |
.bass} file. The units for {\tt tauThermostat} are fs, and the units |
| 251 |
for the {\tt targetTemperature} are degrees K. The integration of |
| 252 |
the equations of motion is carried out in a velocity-Verlet style 2 |
| 253 |
part algorithm: |
| 254 |
|
| 255 |
{\tt moveA:} |
| 256 |
\begin{eqnarray} |
| 257 |
T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\ |
| 258 |
{\bf v}\left(t + \delta t / 2\right) & \leftarrow & {\bf |
| 259 |
v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) |
| 260 |
\chi(t)\right) \\ |
| 261 |
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t {\bf |
| 262 |
v}\left(t + \delta t / 2 \right) \\ |
| 263 |
{\bf j}\left(t + \delta t / 2 \right) & \leftarrow & {\bf |
| 264 |
j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 265 |
\chi(t) \right) \\ |
| 266 |
\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}({\bf j}(t + |
| 267 |
\delta t / 2) \overleftrightarrow{\mathsf{I}}^{b}, |
| 268 |
\delta t) \\ |
| 269 |
\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) + |
| 270 |
\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1 |
| 271 |
\right) |
| 272 |
\end{eqnarray} |
| 273 |
|
| 274 |
Here $\mathrm{rot}( {\bf j} / {\bf I} )$ is the same symplectic Trotter |
| 275 |
factorization of the three rotation operations that was discussed in |
| 276 |
the section on the DLM integrator. Note that this operation modifies |
| 277 |
both the rotation matrix $\mathsf{A}$ and the angular momentum ${\bf |
| 278 |
j}$. {\tt moveA} propagates velocities by a half time step, and |
| 279 |
positional degrees of freedom by a full time step. The new positions |
| 280 |
(and orientations) are then used to calculate a new set of forces and |
| 281 |
torques. |
| 282 |
|
| 283 |
{\tt doForces:} |
| 284 |
\begin{eqnarray} |
| 285 |
{\bf f}(t + \delta t) & \leftarrow & - \frac{\partial V}{\partial {\bf |
| 286 |
r}(t + \delta t)} \\ |
| 287 |
{\bf \tau}^{s}(t + \delta t) & \leftarrow & {\bf u}(t + \delta t) |
| 288 |
\times \frac{\partial V}{\partial {\bf u}} \\ |
| 289 |
{\bf \tau}^{b}(t + \delta t) & \leftarrow & \mathsf{A}(t + \delta t) |
| 290 |
\cdot {\bf \tau}^s(t + \delta t) |
| 291 |
\end{eqnarray} |
| 292 |
|
| 293 |
Here ${\bf u}$ is a unit vector pointing along the principal axis of |
| 294 |
the rigid body being propagated, ${\bf \tau}^s$ is the torque in the |
| 295 |
space-fixed (laboratory) frame, and ${\bf \tau}^b$ is the torque in |
| 296 |
the body-fixed frame. ${\bf u}$ is automatically calculated when the |
| 297 |
rotation matrix $\mathsf{A}$ is calculated in {\tt moveA}. |
| 298 |
|
| 299 |
Once the forces and torques have been obtained at the new time step, |
| 300 |
the velocities can be advanced to the same time value. |
| 301 |
|
| 302 |
{\tt moveB:} |
| 303 |
\begin{eqnarray} |
| 304 |
T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\}, |
| 305 |
\left\{{\bf j}(t + \delta t)\right\} \\ |
| 306 |
\chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t / |
| 307 |
2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta |
| 308 |
t)}{T_{\mathrm{target}}} - 1 \right) \\ |
| 309 |
{\bf v}\left(t + \delta t \right) & \leftarrow & {\bf |
| 310 |
v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( |
| 311 |
\frac{{\bf f}(t + \delta t)}{m} - {\bf v}(t + \delta t) |
| 312 |
\chi(t \delta t)\right) \\ |
| 313 |
{\bf j}\left(t + \delta t \right) & \leftarrow & {\bf |
| 314 |
j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf |
| 315 |
\tau}^b(t + \delta t) - {\bf j}(t + \delta t) |
| 316 |
\chi(t + \delta t) \right) |
| 317 |
\end{eqnarray} |
| 318 |
|
| 319 |
Since ${\bf v}(t + \delta t)$ and ${\bf j}(t + \delta t)$ are required |
| 320 |
to caclculate $T(t + \delta t)$ as well as $\chi(t + \delta t)$, the |
| 321 |
indirectly depend on their own values at time $t + \delta t$. {\tt |
| 322 |
moveB} is therefore done in an iterative fashion until $\chi(t + |
| 323 |
\delta t)$ becomes self-consistent. The relative tolerance for the |
| 324 |
self-consistency check defaults to a value of $\mbox{10}^{-6}$, but |
| 325 |
{\sc oopse} will terminate the iteration after 4 loops even if the |
| 326 |
consistency check has not been satisfied. |
| 327 |
|
| 328 |
The Nos\'e-Hoover algorithm is known to conserve a Hamiltonian for the |
| 329 |
extended system that is, to within a constant, identical to the |
| 330 |
Helmholtz free energy, |
| 331 |
\begin{equation} |
| 332 |
H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left( |
| 333 |
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t\prime) dt\prime |
| 334 |
\right) |
| 335 |
\end{equation} |
| 336 |
Poor choices of $\delta t$ or $\tau_T$ can result in non-conservation |
| 337 |
of $H_{\mathrm{NVT}}$, so the conserved quantity is maintained in the |
| 338 |
last column of the {\tt .stat} file to allow checks on the quality of |
| 339 |
the integration. |
| 340 |
|
| 341 |
Bond constraints are applied at the end of both the {\tt moveA} and |
| 342 |
{\tt moveB} portions of the algorithm. Details on the constraint |
| 343 |
algorithms are given in section \ref{sec:rattle}. |
| 344 |
|
| 345 |
\subsubsection{\label{sec:NPTi}Constant-pressure integration (isotropic box)} |
| 346 |
|
| 347 |
|
| 348 |
\subsection{\label{Sec:zcons}Z-Constraint Method} |
| 349 |
|
| 350 |
Based on fluctuatin-dissipation theorem,\bigskip\ force auto-correlation |
| 351 |
method was developed to investigate the dynamics of ions inside the ion |
| 352 |
channels.\cite{Roux91} Time-dependent friction coefficient can be calculated |
| 353 |
from the deviation of the instaneous force from its mean force. |
| 354 |
|
| 355 |
% |
| 356 |
|
| 357 |
\begin{equation} |
| 358 |
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T |
| 359 |
\end{equation} |
| 360 |
|
| 361 |
|
| 362 |
where% |
| 363 |
\begin{equation} |
| 364 |
\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle |
| 365 |
\end{equation} |
| 366 |
|
| 367 |
|
| 368 |
If the time-dependent friction decay rapidly, static friction coefficient can |
| 369 |
be approximated by% |
| 370 |
|
| 371 |
\begin{equation} |
| 372 |
\xi^{static}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt |
| 373 |
\end{equation} |
| 374 |
|
| 375 |
|
| 376 |
Hence, diffusion constant can be estimated by |
| 377 |
\begin{equation} |
| 378 |
D(z)=\frac{k_{B}T}{\xi^{static}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty |
| 379 |
}\langle\delta F(z,t)\delta F(z,0)\rangle dt}% |
| 380 |
\end{equation} |
| 381 |
|
| 382 |
|
| 383 |
\bigskip Z-Constraint method, which fixed the z coordinates of the molecules |
| 384 |
with respect to the center of the mass of the system, was proposed to obtain |
| 385 |
the forces required in force auto-correlation method.\cite{Marrink94} However, |
| 386 |
simply resetting the coordinate will move the center of the mass of the whole |
| 387 |
system. To avoid this problem, a new method was used at {\sc oopse}. Instead of |
| 388 |
resetting the coordinate, we reset the forces of z-constraint molecules as |
| 389 |
well as subtract the total constraint forces from the rest of the system after |
| 390 |
force calculation at each time step. |
| 391 |
\begin{verbatim} |
| 392 |
$F_{\alpha i}=0$ |
| 393 |
$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}$ |
| 394 |
$F_{\alpha i}=F_{\alpha i}-\frac{M_{_{\alpha i}}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}\sum\limits_{\beta}F_{\beta}$ |
| 395 |
$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}}}$ |
| 396 |
\end{verbatim} |
| 397 |
|
| 398 |
At the very beginning of the simulation, the molecules may not be at its |
| 399 |
constraint position. To move the z-constraint molecule to the specified |
| 400 |
position, a simple harmonic potential is used% |
| 401 |
|
| 402 |
\begin{equation} |
| 403 |
U(t)=\frac{1}{2}k_{Harmonic}(z(t)-z_{cons})^{2}% |
| 404 |
\end{equation} |
| 405 |
where $k_{Harmonic}$\bigskip\ is the harmonic force constant, $z(t)$ is |
| 406 |
current z coordinate of the center of mass of the z-constraint molecule, and |
| 407 |
$z_{cons}$ is the restraint position. Therefore, the harmonic force operated |
| 408 |
on the z-constraint molecule at time $t$ can be calculated by% |
| 409 |
\begin{equation} |
| 410 |
F_{z_{Harmonic}}(t)=-\frac{\partial U(t)}{\partial z(t)}=-k_{Harmonic}% |
| 411 |
(z(t)-z_{cons}) |
| 412 |
\end{equation} |
| 413 |
Worthy of mention, other kinds of potential functions can also be used to |
| 414 |
drive the z-constraint molecule. |
