ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/tengDissertation/Methodology.tex
(Generate patch)

Comparing trunk/tengDissertation/Methodology.tex (file contents):
Revision 2804 by tim, Tue Jun 6 19:47:27 2006 UTC vs.
Revision 2914 by tim, Fri Jun 30 14:35:34 2006 UTC

# Line 2 | Line 2 | In order to mimic the experiments, which are usually p
2  
3   \section{\label{methodSection:rigidBodyIntegrators}Integrators for Rigid Body Motion in Molecular Dynamics}
4  
5 < In order to mimic the experiments, which are usually performed under
5 > In order to mimic experiments which are usually performed under
6   constant temperature and/or pressure, extended Hamiltonian system
7   methods have been developed to generate statistical ensembles, such
8 < as canonical ensemble and isobaric-isothermal ensemble \textit{etc}.
9 < In addition to the standard ensemble, specific ensembles have been
10 < developed to account for the anisotropy between the lateral and
11 < normal directions of membranes. The $NPAT$ ensemble, in which the
12 < normal pressure and the lateral surface area of the membrane are
13 < kept constant, and the $NP\gamma T$ ensemble, in which the normal
14 < pressure and the lateral surface tension are kept constant were
15 < proposed to address this issue.
8 > as the canonical and isobaric-isothermal ensembles. In addition to
9 > the standard ensemble, specific ensembles have been developed to
10 > account for the anisotropy between the lateral and normal directions
11 > of membranes. The $NPAT$ ensemble, in which the normal pressure and
12 > the lateral surface area of the membrane are kept constant, and the
13 > $NP\gamma T$ ensemble, in which the normal pressure and the lateral
14 > surface tension are kept constant were proposed to address the
15 > issues.
16  
17 < Integration schemes for rotational motion of the rigid molecules in
18 < microcanonical ensemble have been extensively studied in the last
19 < two decades. Matubayasi and Nakahara developed a time-reversible
20 < integrator for rigid bodies in quaternion representation. Although
21 < it is not symplectic, this integrator still demonstrates a better
22 < long-time energy conservation than traditional methods because of
23 < the time-reversible nature. Extending Trotter-Suzuki to general
24 < system with a flat phase space, Miller and his colleagues devised an
25 < novel symplectic, time-reversible and volume-preserving integrator
26 < in quaternion representation, which was shown to be superior to the
27 < time-reversible integrator of Matubayasi and Nakahara. However, all
28 < of the integrators in quaternion representation suffer from the
17 > Integration schemes for the rotational motion of the rigid molecules
18 > in the microcanonical ensemble have been extensively studied over
19 > the last two decades. Matubayasi developed a
20 > time-reversible integrator for rigid bodies in quaternion
21 > representation\cite{Matubayasi1999}. Although it is not symplectic, this integrator still
22 > demonstrates a better long-time energy conservation than Euler angle
23 > methods because of the time-reversible nature. Extending the
24 > Trotter-Suzuki factorization to general system with a flat phase
25 > space, Miller\cite{Miller2002} and his colleagues devised a novel
26 > symplectic, time-reversible and volume-preserving integrator in the
27 > quaternion representation, which was shown to be superior to the
28 > Matubayasi's time-reversible integrator. However, all of the
29 > integrators in the quaternion representation suffer from the
30   computational penalty of constructing a rotation matrix from
31   quaternions to evolve coordinates and velocities at every time step.
32 < An alternative integration scheme utilizing rotation matrix directly
33 < proposed by Dullweber, Leimkuhler and McLachlan (DLM) also preserved
34 < the same structural properties of the Hamiltonian flow. In this
35 < section, the integration scheme of DLM method will be reviewed and
36 < extended to other ensembles.
32 > An alternative integration scheme utilizing the rotation matrix
33 > directly proposed by Dullweber, Leimkuhler and McLachlan (DLM) also
34 > preserved the same structural properties of the Hamiltonian
35 > propagator\cite{Dullweber1997}. In this section, the integration
36 > scheme of DLM method will be reviewed and extended to other
37 > ensembles.
38  
39   \subsection{\label{methodSection:DLM}Integrating the Equations of Motion: the
40   DLM method}
# Line 44 | Line 46 | Monte Carlo applications, and
46   {\it symplectic}),
47   \item the integrator is time-{\it reversible}, making it suitable for Hybrid
48   Monte Carlo applications, and
49 < \item the error for a single time step is of order $\mathcal{O}\left(h^4\right)$
49 > \item the error for a single time step is of order $\mathcal{O}\left(h^3\right)$
50   for timesteps of length $h$.
51   \end{enumerate}
50
52   The integration of the equations of motion is carried out in a
53   velocity-Verlet style 2-part algorithm, where $h= \delta t$:
54  
# Line 62 | Line 63 | velocity-Verlet style 2-part algorithm, where $h= \del
63   {\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t)
64      + \frac{h}{2} {\bf \tau}^b(t), \\
65   %
66 < \mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
66 > \mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
67      (t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right).
68   \end{align*}
68
69   In this context, the $\mathrm{rotate}$ function is the reversible
70   product of the three body-fixed rotations,
71   \begin{equation}
# Line 74 | Line 74 | rotates both the rotation matrix ($\mathsf{A}$) and th
74   / 2) \cdot \mathsf{G}_x(a_x /2),
75   \end{equation}
76   where each rotational propagator, $\mathsf{G}_\alpha(\theta)$,
77 < rotates both the rotation matrix ($\mathsf{A}$) and the body-fixed
78 < angular momentum (${\bf j}$) by an angle $\theta$ around body-fixed
77 > rotates both the rotation matrix $\mathsf{Q}$ and the body-fixed
78 > angular momentum ${\bf j}$ by an angle $\theta$ around body-fixed
79   axis $\alpha$,
80   \begin{equation}
81   \mathsf{G}_\alpha( \theta ) = \left\{
82   \begin{array}{lcl}
83 < \mathsf{A}(t) & \leftarrow & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
83 > \mathsf{Q}(t) & \leftarrow & \mathsf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
84   {\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf
85   j}(0).
86   \end{array}
# Line 100 | Line 100 | All other rotations follow in a straightforward manner
100   \end{array}
101   \right).
102   \end{equation}
103 < All other rotations follow in a straightforward manner.
103 > All other rotations follow in a straightforward manner. After the
104 > first part of the propagation, the forces and body-fixed torques are
105 > calculated at the new positions and orientations
106  
105 After the first part of the propagation, the forces and body-fixed
106 torques are calculated at the new positions and orientations
107
107   {\tt doForces:}
108   \begin{align*}
109   {\bf f}(t + h) &\leftarrow
110      - \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)}, \\
111   %
112   {\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h)
113 <    \times \frac{\partial V}{\partial {\bf u}}, \\
113 >    \times (\frac{\partial V}{\partial {\bf u}})_{u(t+h)}, \\
114   %
115 < {\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h)
115 > {\bf \tau}^{b}(t + h) &\leftarrow \mathsf{Q}(t + h)
116      \cdot {\bf \tau}^s(t + h).
117   \end{align*}
118 <
119 < {\sc oopse} automatically updates ${\bf u}$ when the rotation matrix
121 < $\mathsf{A}$ is calculated in {\tt moveA}.  Once the forces and
118 > ${\bf u}$ is automatically updated when the rotation matrix
119 > $\mathsf{Q}$ is calculated in {\tt moveA}.  Once the forces and
120   torques have been obtained at the new time step, the velocities can
121   be advanced to the same time value.
122  
# Line 132 | Line 130 | be advanced to the same time value.
130   \right)
131      + \frac{h}{2} {\bf \tau}^b(t + h) .
132   \end{align*}
135
133   The matrix rotations used in the DLM method end up being more costly
134   computationally than the simpler arithmetic quaternion propagation.
135   With the same time step, a 1000-molecule water simulation shows an
136   average 7\% increase in computation time using the DLM method in
137   place of quaternions. This cost is more than justified when
138   comparing the energy conservation of the two methods as illustrated
139 < in Fig.~\ref{timestep}.
139 > in Fig.~\ref{methodFig:timestep} where the resulting energy drifts at
140 > various time steps for both the DLM and quaternion integration
141 > schemes are compared. All of the 1000 molecule water simulations
142 > started with the same configuration, and the only difference was the
143 > method for handling rotational motion. At time steps of 0.1 and 0.5
144 > fs, both methods for propagating molecule rotation conserve energy
145 > fairly well, with the quaternion method showing a slight energy
146 > drift over time in the 0.5 fs time step simulation. At time steps of
147 > 1 and 2 fs, the energy conservation benefits of the DLM method are
148 > clearly demonstrated. Thus, while maintaining the same degree of
149 > energy conservation, one can take considerably longer time steps,
150 > leading to an overall reduction in computation time.
151  
152   \begin{figure}
153   \centering
# Line 153 | Line 161 | from the true energy baseline for clarity.}
161   from the true energy baseline for clarity.}
162   \label{methodFig:timestep}
163   \end{figure}
156
157 In Fig.~\ref{methodFig:timestep}, the resulting energy drift at
158 various time steps for both the DLM and quaternion integration
159 schemes is compared. All of the 1000 molecule water simulations
160 started with the same configuration, and the only difference was the
161 method for handling rotational motion. At time steps of 0.1 and 0.5
162 fs, both methods for propagating molecule rotation conserve energy
163 fairly well, with the quaternion method showing a slight energy
164 drift over time in the 0.5 fs time step simulation. At time steps of
165 1 and 2 fs, the energy conservation benefits of the DLM method are
166 clearly demonstrated. Thus, while maintaining the same degree of
167 energy conservation, one can take considerably longer time steps,
168 leading to an overall reduction in computation time.
164  
165   \subsection{\label{methodSection:NVT}Nos\'{e}-Hoover Thermostatting}
166  
# Line 173 | Line 168 | The Nos\'e-Hoover equations of motion are given by\cit
168   \begin{eqnarray}
169   \dot{{\bf r}} & = & {\bf v}, \\
170   \dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} ,\\
171 < \dot{\mathsf{A}} & = & \mathsf{A} \cdot
171 > \dot{\mathsf{Q}} & = & \mathsf{Q} \cdot
172   \mbox{ skew}\left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right), \\
173   \dot{{\bf j}} & = & {\bf j} \times \left(
174   \overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j} \right) - \mbox{
175 < rot}\left(\mathsf{A}^{T} \cdot \frac{\partial V}{\partial
176 < \mathsf{A}} \right) - \chi {\bf j}. \label{eq:nosehoovereom}
175 > rot}\left(\mathsf{Q}^{T} \cdot \frac{\partial V}{\partial
176 > \mathsf{Q}} \right) - \chi {\bf j}. \label{eq:nosehoovereom}
177   \end{eqnarray}
183
178   $\chi$ is an ``extra'' variable included in the extended system, and
179   it is propagated using the first order equation of motion
180   \begin{equation}
181   \dot{\chi} = \frac{1}{\tau_{T}^2} \left(
182   \frac{T}{T_{\mathrm{target}}} - 1 \right). \label{eq:nosehooverext}
183   \end{equation}
184 <
185 < The instantaneous temperature $T$ is proportional to the total
186 < kinetic energy (both translational and orientational) and is given
193 < by
184 > where $\tau_T$ is the time constant for relaxation of the
185 > temperature to the target value, and the instantaneous temperature
186 > $T$ is given by
187   \begin{equation}
188 < T = \frac{2 K}{f k_B}
188 > T = \frac{2 K}{f k_B}.
189   \end{equation}
190   Here, $f$ is the total number of degrees of freedom in the system,
191   \begin{equation}
192   f = 3 N + 3 N_{\mathrm{orient}} - N_{\mathrm{constraints}},
193   \end{equation}
194 < and $K$ is the total kinetic energy,
195 < \begin{equation}
196 < K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i +
204 < \sum_{i=1}^{N_{\mathrm{orient}}}  \frac{1}{2} {\bf j}_i^T \cdot
205 < \overleftrightarrow{\mathsf{I}}_i^{-1} \cdot {\bf j}_i.
206 < \end{equation}
194 > where $N_{\mathrm{orient}}$ is the number of molecules with
195 > orientational degrees of freedom. The integration of the equations of motion
196 > is carried out in a velocity-Verlet style 2 part algorithm:
197  
208 In eq.(\ref{eq:nosehooverext}), $\tau_T$ is the time constant for
209 relaxation of the temperature to the target value.  To set values
210 for $\tau_T$ or $T_{\mathrm{target}}$ in a simulation, one would use
211 the {\tt tauThermostat} and {\tt targetTemperature} keywords in the
212 {\tt .bass} file.  The units for {\tt tauThermostat} are fs, and the
213 units for the {\tt targetTemperature} are degrees K.   The
214 integration of the equations of motion is carried out in a
215 velocity-Verlet style 2 part algorithm:
216
198   {\tt moveA:}
199   \begin{align*}
200   T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\
# Line 229 | Line 210 | T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\b
210      + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
211      \chi(t) \right) ,\\
212   %
213 < \mathsf{A}(t + h) &\leftarrow \mathrm{rotate}
214 <    \left(h * {\bf j}(t + h / 2)
213 > \mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}
214 >    \left(h {\bf j}(t + h / 2)
215      \overleftrightarrow{\mathsf{I}}^{-1} \right) ,\\
216   %
217   \chi\left(t + h / 2 \right) &\leftarrow \chi(t)
218      + \frac{h}{2 \tau_T^2} \left( \frac{T(t)}
219      {T_{\mathrm{target}}} - 1 \right) .
220   \end{align*}
240
221   Here $\mathrm{rotate}(h * {\bf j}
222 < \overleftrightarrow{\mathsf{I}}^{-1})$ is the same symplectic
223 < Trotter factorization of the three rotation operations that was
224 < discussed in the section on the DLM integrator.  Note that this
225 < operation modifies both the rotation matrix $\mathsf{A}$ and the
226 < angular momentum ${\bf j}$.  {\tt moveA} propagates velocities by a
227 < half time step, and positional degrees of freedom by a full time
228 < step.  The new positions (and orientations) are then used to
229 < calculate a new set of forces and torques in exactly the same way
230 < they are calculated in the {\tt doForces} portion of the DLM
231 < integrator.
232 <
253 < Once the forces and torques have been obtained at the new time step,
254 < the temperature, velocities, and the extended system variable can be
222 > \overleftrightarrow{\mathsf{I}}^{-1})$ is the same symplectic Strang
223 > factorization of the three rotation operations that was discussed in
224 > the section on the DLM integrator.  Note that this operation
225 > modifies both the rotation matrix $\mathsf{Q}$ and the angular
226 > momentum ${\bf j}$.  {\tt moveA} propagates velocities by a half
227 > time step, and positional degrees of freedom by a full time step.
228 > The new positions (and orientations) are then used to calculate a
229 > new set of forces and torques in exactly the same way they are
230 > calculated in the {\tt doForces} portion of the DLM integrator. Once
231 > the forces and torques have been obtained at the new time step, the
232 > temperature, velocities, and the extended system variable can be
233   advanced to the same time value.
234  
235   {\tt moveB:}
# Line 273 | Line 251 | T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
251      \left( {\bf \tau}^b(t + h) - {\bf j}(t + h)
252      \chi(t + h) \right) .
253   \end{align*}
276
254   Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to
255   caclculate $T(t + h)$ as well as $\chi(t + h)$, they indirectly
256   depend on their own values at time $t + h$.  {\tt moveB} is
257   therefore done in an iterative fashion until $\chi(t + h)$ becomes
258 < self-consistent.  The relative tolerance for the self-consistency
259 < check defaults to a value of $\mbox{10}^{-6}$, but {\sc oopse} will
260 < terminate the iteration after 4 loops even if the consistency check
284 < has not been satisfied.
285 <
286 < The Nos\'e-Hoover algorithm is known to conserve a Hamiltonian for
287 < the extended system that is, to within a constant, identical to the
288 < Helmholtz free energy,\cite{Melchionna1993}
258 > self-consistent. The Nos\'e-Hoover algorithm is known to conserve a
259 > Hamiltonian for the extended system that is, to within a constant,
260 > identical to the Helmholtz free energy,\cite{Melchionna1993}
261   \begin{equation}
262   H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left(
263   \frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime)
264   dt^\prime \right).
265   \end{equation}
266   Poor choices of $h$ or $\tau_T$ can result in non-conservation of
267 < $H_{\mathrm{NVT}}$, so the conserved quantity is maintained in the
268 < last column of the {\tt .stat} file to allow checks on the quality
297 < of the integration.
267 > $H_{\mathrm{NVT}}$, so the conserved quantity should be checked
268 > periodically to verify the quality of the integration.
269  
270   \subsection{\label{methodSection:NPTi}Constant-pressure integration with
271 < isotropic box deformations (NPTi)}
271 > isotropic box (NPTi)}
272  
273 < To carry out isobaric-isothermal ensemble calculations {\sc oopse}
274 < implements the Melchionna modifications to the
275 < Nos\'e-Hoover-Andersen equations of motion,\cite{Melchionna1993}
305 <
273 > We can used an isobaric-isothermal ensemble integrator which is
274 > implemented using the Melchionna modifications to the
275 > Nos\'e-Hoover-Andersen equations of motion\cite{Melchionna1993}
276   \begin{eqnarray}
277   \dot{{\bf r}} & = & {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right), \\
278   \dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\eta + \chi) {\bf v}, \\
279 < \dot{\mathsf{A}} & = & \mathsf{A} \cdot
279 > \dot{\mathsf{Q}} & = & \mathsf{Q} \cdot
280   \mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right),\\
281   \dot{{\bf j}} & = & {\bf j} \times \left(
282   \overleftrightarrow{I}^{-1} \cdot {\bf j} \right) - \mbox{
283 < rot}\left(\mathsf{A}^{T} \cdot \frac{\partial
284 < V}{\partial \mathsf{A}} \right) - \chi {\bf j}, \\
283 > rot}\left(\mathsf{Q}^{T} \cdot \frac{\partial
284 > V}{\partial \mathsf{Q}} \right) - \chi {\bf j}, \\
285   \dot{\chi} & = & \frac{1}{\tau_{T}^2} \left(
286   \frac{T}{T_{\mathrm{target}}} - 1 \right) ,\\
287   \dot{\eta} & = & \frac{1}{\tau_{B}^2 f k_B T_{\mathrm{target}}} V
# Line 319 | Line 289 | P_{\mathrm{target}} \right), \\
289   P_{\mathrm{target}} \right), \\
290   \dot{\mathcal{V}} & = & 3 \mathcal{V} \eta . \label{eq:melchionna1}
291   \end{eqnarray}
322
292   $\chi$ and $\eta$ are the ``extra'' degrees of freedom in the
293   extended system.  $\chi$ is a thermostat, and it has the same
294   function as it does in the Nos\'e-Hoover NVT integrator.  $\eta$ is
# Line 352 | Line 321 | P(t) = \frac{1}{3} \mathrm{Tr} \left(
321   the Pressure tensor,
322   \begin{equation}
323   P(t) = \frac{1}{3} \mathrm{Tr} \left(
324 < \overleftrightarrow{\mathsf{P}}(t). \right)
324 > \overleftrightarrow{\mathsf{P}}(t) \right) .
325   \end{equation}
326 <
327 < In eq.(\ref{eq:melchionna1}), $\tau_B$ is the time constant for
359 < relaxation of the pressure to the target value.  To set values for
360 < $\tau_B$ or $P_{\mathrm{target}}$ in a simulation, one would use the
361 < {\tt tauBarostat} and {\tt targetPressure} keywords in the {\tt
362 < .bass} file.  The units for {\tt tauBarostat} are fs, and the units
363 < for the {\tt targetPressure} are atmospheres.  Like in the NVT
326 > In Eq.~\ref{eq:melchionna1}, $\tau_B$ is the time constant for
327 > relaxation of the pressure to the target value. Like in the NVT
328   integrator, the integration of the equations of motion is carried
329   out in a velocity-Verlet style 2 part algorithm:
330  
# Line 378 | Line 342 | P(t) &\leftarrow \left\{{\bf r}(t)\right\}, \left\{{\b
342      + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
343      \chi(t) \right), \\
344   %
345 < \mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h *
345 > \mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}\left(h *
346      {\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1}
347      \right) ,\\
348   %
# Line 398 | Line 362 | P(t) &\leftarrow \left\{{\bf r}(t)\right\}, \left\{{\b
362   \mathsf{H}(t + h) &\leftarrow e^{-h \eta(t + h / 2)}
363      \mathsf{H}(t).
364   \end{align*}
401
365   Most of these equations are identical to their counterparts in the
366   NVT integrator, but the propagation of positions to time $t + h$
367 < depends on the positions at the same time.  {\sc oopse} carries out
368 < this step iteratively (with a limit of 5 passes through the
369 < iterative loop).  Also, the simulation box $\mathsf{H}$ is scaled
370 < uniformly for one full time step by an exponential factor that
371 < depends on the value of $\eta$ at time $t + h / 2$.  Reshaping the
409 < box uniformly also scales the volume of the box by
367 > depends on the positions at the same time. The simulation box
368 > $\mathsf{H}$ is scaled uniformly for one full time step by an
369 > exponential factor that depends on the value of $\eta$ at time $t +
370 > h / 2$.  Reshaping the box uniformly also scales the volume of the
371 > box by
372   \begin{equation}
373   \mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)}.
374   \mathcal{V}(t)
375   \end{equation}
414
376   The {\tt doForces} step for the NPTi integrator is exactly the same
377   as in both the DLM and NVT integrators.  Once the forces and torques
378   have been obtained at the new time step, the velocities can be
# Line 443 | Line 404 | P(t + h) &\leftarrow  \left\{{\bf r}(t + h)\right\},
404      \tau}^b(t + h) - {\bf j}(t + h)
405      \chi(t + h) \right) .
406   \end{align*}
446
407   Once again, since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required
408   to caclculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, and $\eta(t +
409   h)$, they indirectly depend on their own values at time $t + h$.
410   {\tt moveB} is therefore done in an iterative fashion until $\chi(t
411 < + h)$ and $\eta(t + h)$ become self-consistent.  The relative
452 < tolerance for the self-consistency check defaults to a value of
453 < $\mbox{10}^{-6}$, but {\sc oopse} will terminate the iteration after
454 < 4 loops even if the consistency check has not been satisfied.
411 > + h)$ and $\eta(t + h)$ become self-consistent.
412  
413   The Melchionna modification of the Nos\'e-Hoover-Andersen algorithm
414   is known to conserve a Hamiltonian for the extended system that is,
# Line 461 | Line 418 | Poor choices of $\delta t$, $\tau_T$, or $\tau_B$ can
418   \frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime)
419   dt^\prime \right) + P_{\mathrm{target}} \mathcal{V}(t).
420   \end{equation}
421 < Poor choices of $\delta t$, $\tau_T$, or $\tau_B$ can result in
422 < non-conservation of $H_{\mathrm{NPTi}}$, so the conserved quantity
423 < is maintained in the last column of the {\tt .stat} file to allow
467 < checks on the quality of the integration.  It is also known that
468 < this algorithm samples the equilibrium distribution for the enthalpy
469 < (including contributions for the thermostat and barostat),
421 > It is also known that this algorithm samples the equilibrium
422 > distribution for the enthalpy (including contributions for the
423 > thermostat and barostat),
424   \begin{equation}
425   H_{\mathrm{NPTi}} = V + K + \frac{f k_B T_{\mathrm{target}}}{2}
426   \left( \chi^2 \tau_T^2 + \eta^2 \tau_B^2 \right) +
427   P_{\mathrm{target}} \mathcal{V}(t).
428   \end{equation}
429  
476 Bond constraints are applied at the end of both the {\tt moveA} and
477 {\tt moveB} portions of the algorithm.  Details on the constraint
478 algorithms are given in section \ref{oopseSec:rattle}.
479
430   \subsection{\label{methodSection:NPTf}Constant-pressure integration with a
431   flexible box (NPTf)}
432  
# Line 491 | Line 441 | method are
441   \dot{{\bf r}} & = & {\bf v} + \overleftrightarrow{\eta} \cdot \left( {\bf r} - {\bf R}_0 \right), \\
442   \dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\overleftrightarrow{\eta} +
443   \chi \cdot \mathsf{1}) {\bf v}, \\
444 < \dot{\mathsf{A}} & = & \mathsf{A} \cdot
444 > \dot{\mathsf{Q}} & = & \mathsf{Q} \cdot
445   \mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) ,\\
446   \dot{{\bf j}} & = & {\bf j} \times \left(
447   \overleftrightarrow{I}^{-1} \cdot {\bf j} \right) - \mbox{
448 < rot}\left(\mathsf{A}^{T} \cdot \frac{\partial
449 < V}{\partial \mathsf{A}} \right) - \chi {\bf j} ,\\
448 > rot}\left(\mathsf{Q}^{T} \cdot \frac{\partial
449 > V}{\partial \mathsf{Q}} \right) - \chi {\bf j} ,\\
450   \dot{\chi} & = & \frac{1}{\tau_{T}^2} \left(
451   \frac{T}{T_{\mathrm{target}}} - 1 \right) ,\\
452   \dot{\overleftrightarrow{\eta}} & = & \frac{1}{\tau_{B}^2 f k_B
# Line 529 | Line 479 | r}(t)\right\},
479      + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
480      \chi(t) \right), \\
481   %
482 < \mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h *
482 > \mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}\left(h *
483      {\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1}
484      \right), \\
485   %
# Line 550 | Line 500 | r}(t)\right\},
500   \mathsf{H}(t + h) &\leftarrow \mathsf{H}(t) \cdot e^{-h
501      \overleftrightarrow{\eta}(t + h / 2)} .
502   \end{align*}
503 < {\sc oopse} uses a power series expansion truncated at second order
504 < for the exponential operation which scales the simulation box.
505 <
506 < The {\tt moveB} portion of the algorithm is largely unchanged from
557 < the NPTi integrator:
503 > Here, a power series expansion truncated at second order for the
504 > exponential operation is used to scale the simulation box. The {\tt
505 > moveB} portion of the algorithm is largely unchanged from the NPTi
506 > integrator:
507  
508   {\tt moveB:}
509   \begin{align*}
# Line 585 | Line 534 | T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
534      + h / 2 \right) + \frac{h}{2} \left( {\bf \tau}^b(t
535      + h) - {\bf j}(t + h) \chi(t + h) \right) .
536   \end{align*}
588
537   The iterative schemes for both {\tt moveA} and {\tt moveB} are
538 < identical to those described for the NPTi integrator.
539 <
540 < The NPTf integrator is known to conserve the following Hamiltonian:
541 < \begin{equation}
594 < H_{\mathrm{NPTf}} = V + K + f k_B T_{\mathrm{target}} \left(
538 > identical to those described for the NPTi integrator. The NPTf
539 > integrator is known to conserve the following Hamiltonian:
540 > \begin{eqnarray*}
541 > H_{\mathrm{NPTf}} & = & V + K + f k_B T_{\mathrm{target}} \left(
542   \frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime)
543 < dt^\prime \right) + P_{\mathrm{target}} \mathcal{V}(t) + \frac{f k_B
543 > dt^\prime \right) \\
544 > & & + P_{\mathrm{target}} \mathcal{V}(t) + \frac{f k_B
545   T_{\mathrm{target}}}{2}
546   \mathrm{Tr}\left[\overleftrightarrow{\eta}(t)\right]^2 \tau_B^2.
547 < \end{equation}
600 <
547 > \end{eqnarray*}
548   This integrator must be used with care, particularly in liquid
549   simulations.  Liquids have very small restoring forces in the
550   off-diagonal directions, and the simulation box can very quickly
# Line 606 | Line 553 | assume non-orthorhombic geometries.
553   finds most use in simulating crystals or liquid crystals which
554   assume non-orthorhombic geometries.
555  
556 < \subsection{\label{methodSection:otherSpecialEnsembles}Other Special Ensembles}
556 > \subsection{\label{methodSection:NPAT}NPAT Ensemble}
557  
558 < \subsubsection{\label{methodSection:NPAT}NPAT Ensemble}
558 > A comprehensive understanding of relations between structures and
559 > functions in biological membrane system ultimately relies on
560 > structure and dynamics of lipid bilayers, which are strongly
561 > affected by the interfacial interaction between lipid molecules and
562 > surrounding media. One quantity used to describe the interfacial
563 > interaction is the average surface area per lipid.
564 > Constant area and constant lateral pressure simulations can be
565 > achieved by extending the standard NPT ensemble with a different
566 > pressure control strategy
567  
613 A comprehensive understanding of structure¨Cfunction relations of
614 biological membrane system ultimately relies on structure and
615 dynamics of lipid bilayer, which are strongly affected by the
616 interfacial interaction between lipid molecules and surrounding
617 media. One quantity to describe the interfacial interaction is so
618 called the average surface area per lipid. Constat area and constant
619 lateral pressure simulation can be achieved by extending the
620 standard NPT ensemble with a different pressure control strategy
621
568   \begin{equation}
569   \dot {\overleftrightarrow{\eta}} _{\alpha \beta}=\left\{\begin{array}{ll}
570                    \frac{{V(P_{\alpha \beta }  - P_{{\rm{target}}} )}}{{\tau_{\rm{B}}^{\rm{2}} fk_B T_{{\rm{target}}} }}
# Line 627 | Line 573 | standard NPT ensemble with a different pressure contro
573             \end{array}
574      \right.
575   \end{equation}
630
576   Note that the iterative schemes for NPAT are identical to those
577   described for the NPTi integrator.
578  
579 < \subsubsection{\label{methodSection:NPrT}NP$\gamma$T Ensemble}
579 > \subsection{\label{methodSection:NPrT}NP$\gamma$T
580 > Ensemble}
581  
582   Theoretically, the surface tension $\gamma$ of a stress free
583   membrane system should be zero since its surface free energy $G$ is
584 < minimum with respect to surface area $A$
585 < \[
584 > minimum with respect to surface area $A$,
585 > \begin{equation}
586   \gamma  = \frac{{\partial G}}{{\partial A}}.
587 < \]
588 < However, a surface tension of zero is not appropriate for relatively
589 < small patches of membrane. In order to eliminate the edge effect of
590 < the membrane simulation, a special ensemble, NP$\gamma$T, is
591 < proposed to maintain the lateral surface tension and normal
592 < pressure. The equation of motion for cell size control tensor,
593 < $\eta$, in $NP\gamma T$ is
587 > \end{equation}0
588 > However, a surface tension of zero is not
589 > appropriate for relatively small patches of membrane. In order to
590 > eliminate the edge effect of membrane simulations, a special
591 > ensemble, NP$\gamma$T, has been proposed to maintain the lateral
592 > surface tension and normal pressure. The equation of motion for the
593 > cell size control tensor, $\eta$, in $NP\gamma T$ is
594   \begin{equation}
595   \dot {\overleftrightarrow{\eta}} _{\alpha \beta}=\left\{\begin{array}{ll}
596      - A_{xy} (\gamma _\alpha   - \gamma _{{\rm{target}}} ) & \mbox{$\alpha  = \beta  = x$ or $\alpha  = \beta  = y$}\\
# Line 660 | Line 606 | - P_{{\rm{target}}} )
606   - P_{{\rm{target}}} )
607   \label{methodEquation:instantaneousSurfaceTensor}
608   \end{equation}
663
609   There is one additional extended system integrator (NPTxyz), in
610   which each attempt to preserve the target pressure along the box
611   walls perpendicular to that particular axis.  The lengths of the box
612   axes are allowed to fluctuate independently, but the angle between
613   the box axes does not change. It should be noted that the NPTxyz
614   integrator is a special case of $NP\gamma T$ if the surface tension
615 < $\gamma$ is set to zero.
615 > $\gamma$ is set to zero, and if $x$ and $y$ can move independently.
616  
617 < \section{\label{methodSection:zcons}Z-Constraint Method}
617 > \section{\label{methodSection:zcons}The Z-Constraint Method}
618  
619   Based on the fluctuation-dissipation theorem, a force
620   auto-correlation method was developed by Roux and Karplus to
# Line 683 | Line 628 | where%
628   \begin{equation}
629   \delta F(z,t)=F(z,t)-\langle F(z,t)\rangle.
630   \end{equation}
686
631   If the time-dependent friction decays rapidly, the static friction
632   coefficient can be approximated by
633   \begin{equation}
# Line 696 | Line 640 | D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B
640   D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty
641   }\langle\delta F(z,t)\delta F(z,0)\rangle dt}.%
642   \end{equation}
699
643   The Z-Constraint method, which fixes the z coordinates of the
644   molecules with respect to the center of the mass of the system, has
645   been a method suggested to obtain the forces required for the force
# Line 705 | Line 648 | each time step instead of resetting the coordinate.
648   whole system. To avoid this problem, we reset the forces of
649   z-constrained molecules as well as subtract the total constraint
650   forces from the rest of the system after the force calculation at
651 < each time step instead of resetting the coordinate.
652 <
710 < After the force calculation, define $G_\alpha$ as
651 > each time step instead of resetting the coordinate. After the force
652 > calculation, we define $G_\alpha$ as
653   \begin{equation}
654   G_{\alpha} = \sum_i F_{\alpha i}, \label{oopseEq:zc1}
655   \end{equation}
# Line 742 | Line 684 | v_{\beta i} = v_{\beta i} + \sum_{\alpha}
684   v_{\beta i} = v_{\beta i} + \sum_{\alpha}
685      \frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}.
686   \end{equation}
745
687   At the very beginning of the simulation, the molecules may not be at
688   their constrained positions. To move a z-constrained molecule to its
689   specified position, a simple harmonic potential is used
# Line 758 | Line 699 | F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\part
699   F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}=
700      -k_{\text{Harmonic}}(z(t)-z_{\text{cons}}).
701   \end{equation}
761
762
763 \section{\label{methodSection:langevin}Integrators for Langevin Dynamics of Rigid Bodies}
764
765 %\subsection{\label{methodSection:temperature}Temperature Control}
766
767 %\subsection{\label{methodSection:pressureControl}Pressure Control}
768
769 %\section{\label{methodSection:hydrodynamics}Hydrodynamics}
770
771 %applications of langevin dynamics
772 As an excellent alternative to newtonian dynamics, Langevin
773 dynamics, which mimics a simple heat bath with stochastic and
774 dissipative forces, has been applied in a variety of studies. The
775 stochastic treatment of the solvent enables us to carry out
776 substantially longer time simulation. Implicit solvent Langevin
777 dynamics simulation of met-enkephalin not only outperforms explicit
778 solvent simulation on computation efficiency, but also agrees very
779 well with explicit solvent simulation on dynamics
780 properties\cite{Shen2002}. Recently, applying Langevin dynamics with
781 UNRES model, Liow and his coworkers suggest that protein folding
782 pathways can be possibly exploited within a reasonable amount of
783 time\cite{Liwo2005}. The stochastic nature of the Langevin dynamics
784 also enhances the sampling of the system and increases the
785 probability of crossing energy barrier\cite{Banerjee2004, Cui2003}.
786 Combining Langevin dynamics with Kramers's theory, Klimov and
787 Thirumalai identified the free-energy barrier by studying the
788 viscosity dependence of the protein folding rates\cite{Klimov1997}.
789 In order to account for solvent induced interactions missing from
790 implicit solvent model, Kaya incorporated desolvation free energy
791 barrier into implicit coarse-grained solvent model in protein
792 folding/unfolding study and discovered a higher free energy barrier
793 between the native and denatured states. Because of its stability
794 against noise, Langevin dynamics is very suitable for studying
795 remagnetization processes in various
796 systems\cite{Garcia-Palacios1998,Berkov2002,Denisov2003}. For
797 instance, the oscillation power spectrum of nanoparticles from
798 Langevin dynamics simulation has the same peak frequencies for
799 different wave vectors,which recovers the property of magnetic
800 excitations in small finite structures\cite{Berkov2005a}. In an
801 attempt to reduce the computational cost of simulation, multiple
802 time stepping (MTS) methods have been introduced and have been of
803 great interest to macromolecule and protein
804 community\cite{Tuckerman1992}. Relying on the observation that
805 forces between distant atoms generally demonstrate slower
806 fluctuations than forces between close atoms, MTS method are
807 generally implemented by evaluating the slowly fluctuating forces
808 less frequently than the fast ones. Unfortunately, nonlinear
809 instability resulting from increasing timestep in MTS simulation
810 have became a critical obstruction preventing the long time
811 simulation. Due to the coupling to the heat bath, Langevin dynamics
812 has been shown to be able to damp out the resonance artifact more
813 efficiently\cite{Sandu1999}.
814
815 %review rigid body dynamics
816 Rigid bodies are frequently involved in the modeling of different
817 areas, from engineering, physics, to chemistry. For example,
818 missiles and vehicle are usually modeled by rigid bodies.  The
819 movement of the objects in 3D gaming engine or other physics
820 simulator is governed by the rigid body dynamics. In molecular
821 simulation, rigid body is used to simplify the model in
822 protein-protein docking study\cite{Gray2003}.
823
824 It is very important to develop stable and efficient methods to
825 integrate the equations of motion of orientational degrees of
826 freedom. Euler angles are the nature choice to describe the
827 rotational degrees of freedom. However, due to its singularity, the
828 numerical integration of corresponding equations of motion is very
829 inefficient and inaccurate. Although an alternative integrator using
830 different sets of Euler angles can overcome this
831 difficulty\cite{Ryckaert1977, Andersen1983}, the computational
832 penalty and the lost of angular momentum conservation still remain.
833 In 1977, a singularity free representation utilizing quaternions was
834 developed by Evans\cite{Evans1977}. Unfortunately, this approach
835 suffer from the nonseparable Hamiltonian resulted from quaternion
836 representation, which prevents the symplectic algorithm to be
837 utilized. Another different approach is to apply holonomic
838 constraints to the atoms belonging to the rigid
839 body\cite{Barojas1973}. Each atom moves independently under the
840 normal forces deriving from potential energy and constraint forces
841 which are used to guarantee the rigidness. However, due to their
842 iterative nature, SHAKE and Rattle algorithm converge very slowly
843 when the number of constraint increases.
844
845 The break through in geometric literature suggests that, in order to
846 develop a long-term integration scheme, one should preserve the
847 geometric structure of the flow. Matubayasi and Nakahara developed a
848 time-reversible integrator for rigid bodies in quaternion
849 representation. Although it is not symplectic, this integrator still
850 demonstrates a better long-time energy conservation than traditional
851 methods because of the time-reversible nature. Extending
852 Trotter-Suzuki to general system with a flat phase space, Miller and
853 his colleagues devised an novel symplectic, time-reversible and
854 volume-preserving integrator in quaternion representation. However,
855 all of the integrators in quaternion representation suffer from the
856 computational penalty of constructing a rotation matrix from
857 quaternions to evolve coordinates and velocities at every time step.
858 An alternative integration scheme utilizing rotation matrix directly
859 is RSHAKE\cite{Kol1997}, in which a conjugate momentum to rotation
860 matrix is introduced to re-formulate the Hamiltonian's equation and
861 the Hamiltonian is evolved in a constraint manifold by iteratively
862 satisfying the orthogonality constraint. However, RSHAKE is
863 inefficient because of the iterative procedure. An extremely
864 efficient integration scheme in rotation matrix representation,
865 which also preserves the same structural properties of the
866 Hamiltonian flow as Miller's integrator, is proposed by Dullweber,
867 Leimkuhler and McLachlan (DLM)\cite{Dullweber1997}.
868
869 %review langevin/browninan dynamics for arbitrarily shaped rigid body
870 Combining Langevin or Brownian dynamics with rigid body dynamics,
871 one can study the slow processes in biomolecular systems. Modeling
872 the DNA as a chain of rigid spheres beads, which subject to harmonic
873 potentials as well as excluded volume potentials, Mielke and his
874 coworkers discover rapid superhelical stress generations from the
875 stochastic simulation of twin supercoiling DNA with response to
876 induced torques\cite{Mielke2004}. Membrane fusion is another key
877 biological process which controls a variety of physiological
878 functions, such as release of neurotransmitters \textit{etc}. A
879 typical fusion event happens on the time scale of millisecond, which
880 is impracticable to study using all atomistic model with newtonian
881 mechanics. With the help of coarse-grained rigid body model and
882 stochastic dynamics, the fusion pathways were exploited by many
883 researchers\cite{Noguchi2001,Noguchi2002,Shillcock2005}. Due to the
884 difficulty of numerical integration of anisotropy rotation, most of
885 the rigid body models are simply modeled by sphere, cylinder,
886 ellipsoid or other regular shapes in stochastic simulations. In an
887 effort to account for the diffusion anisotropy of the arbitrary
888 particles, Fernandes and de la Torre improved the original Brownian
889 dynamics simulation algorithm\cite{Ermak1978,Allison1991} by
890 incorporating a generalized $6\times6$ diffusion tensor and
891 introducing a simple rotation evolution scheme consisting of three
892 consecutive rotations\cite{Fernandes2002}. Unfortunately, unexpected
893 error and bias are introduced into the system due to the arbitrary
894 order of applying the noncommuting rotation
895 operators\cite{Beard2003}. Based on the observation the momentum
896 relaxation time is much less than the time step, one may ignore the
897 inertia in Brownian dynamics. However, assumption of the zero
898 average acceleration is not always true for cooperative motion which
899 is common in protein motion. An inertial Brownian dynamics (IBD) was
900 proposed to address this issue by adding an inertial correction
901 term\cite{Beard2001}. As a complement to IBD which has a lower bound
902 in time step because of the inertial relaxation time, long-time-step
903 inertial dynamics (LTID) can be used to investigate the inertial
904 behavior of the polymer segments in low friction
905 regime\cite{Beard2001}. LTID can also deal with the rotational
906 dynamics for nonskew bodies without translation-rotation coupling by
907 separating the translation and rotation motion and taking advantage
908 of the analytical solution of hydrodynamics properties. However,
909 typical nonskew bodies like cylinder and ellipsoid are inadequate to
910 represent most complex macromolecule assemblies. These intricate
911 molecules have been represented by a set of beads and their
912 hydrodynamics properties can be calculated using variant
913 hydrodynamic interaction tensors.
914
915 The goal of the present work is to develop a Langevin dynamics
916 algorithm for arbitrary rigid particles by integrating the accurate
917 estimation of friction tensor from hydrodynamics theory into the
918 sophisticated rigid body dynamics.
919
920
921 \subsection{Friction Tensor}
922
923 For an arbitrary rigid body moves in a fluid, it may experience
924 friction force $f_r$ or friction torque $\tau _r$ along the opposite
925 direction of the velocity $v$ or angular velocity $\omega$ at
926 arbitrary origin $P$,
927 \begin{equation}
928 \left( \begin{array}{l}
929 f_r  \\
930 \tau _r  \\
931 \end{array} \right) =  - \left( {\begin{array}{*{20}c}
932   {\Xi _{P,t} } & {\Xi _{P,c}^T }  \\
933   {\Xi _{P,c} } & {\Xi _{P,r} }  \\
934 \end{array}} \right)\left( \begin{array}{l}
935 \nu  \\
936 \omega  \\
937 \end{array} \right)
938 \end{equation}
939 where $\Xi _{P,t}t$ is the translation friction tensor, $\Xi _{P,r}$
940 is the rotational friction tensor and $\Xi _{P,c}$ is the
941 translation-rotation coupling tensor. The procedure of calculating
942 friction tensor using hydrodynamic tensor and comparison between
943 bead model and shell model were elaborated by Carrasco \textit{et
944 al}\cite{Carrasco1999}. An important property of the friction tensor
945 is that the translational friction tensor is independent of origin
946 while the rotational and coupling are sensitive to the choice of the
947 origin \cite{Brenner1967}, which can be described by
948 \begin{equation}
949 \begin{array}{c}
950 \Xi _{P,t}  = \Xi _{O,t}  = \Xi _t  \\
951 \Xi _{P,c}  = \Xi _{O,c}  - r_{OP}  \times \Xi _t  \\
952 \Xi _{P,r}  = \Xi _{O,r}  - r_{OP}  \times \Xi _t  \times r_{OP}  + \Xi _{O,c}  \times r_{OP}  - r_{OP}  \times \Xi _{O,c}^T  \\
953 \end{array}
954 \end{equation}
955 Where $O$ is another origin and $r_{OP}$ is the vector joining $O$
956 and $P$. It is also worthy of mention that both of translational and
957 rotational frictional tensors are always symmetric. In contrast,
958 coupling tensor is only symmetric at center of reaction:
959 \begin{equation}
960 \Xi _{R,c}  = \Xi _{R,c}^T
961 \end{equation}
962 The proper location for applying friction force is the center of
963 reaction, at which the trace of rotational resistance tensor reaches
964 minimum.
965
966 \subsection{Rigid body dynamics}
967
968 The Hamiltonian of rigid body can be separated in terms of potential
969 energy $V(r,A)$ and kinetic energy $T(p,\pi)$,
970 \[
971 H = V(r,A) + T(v,\pi )
972 \]
973 A second-order symplectic method is now obtained by the composition
974 of the flow maps,
975 \[
976 \varphi _{\Delta t}  = \varphi _{\Delta t/2,V}  \circ \varphi
977 _{\Delta t,T}  \circ \varphi _{\Delta t/2,V}.
978 \]
979 Moreover, $\varphi _{\Delta t/2,V}$ can be divided into two
980 sub-flows which corresponding to force and torque respectively,
981 \[
982 \varphi _{\Delta t/2,V}  = \varphi _{\Delta t/2,F}  \circ \varphi
983 _{\Delta t/2,\tau }.
984 \]
985 Since the associated operators of $\varphi _{\Delta t/2,F} $ and
986 $\circ \varphi _{\Delta t/2,\tau }$ are commuted, the composition
987 order inside $\varphi _{\Delta t/2,V}$ does not matter.
988
989 Furthermore, kinetic potential can be separated to translational
990 kinetic term, $T^t (p)$, and rotational kinetic term, $T^r (\pi )$,
991 \begin{equation}
992 T(p,\pi ) =T^t (p) + T^r (\pi ).
993 \end{equation}
994 where $ T^t (p) = \frac{1}{2}v^T m v $ and $T^r (\pi )$ is defined
995 by \ref{introEquation:rotationalKineticRB}. Therefore, the
996 corresponding flow maps are given by
997 \[
998 \varphi _{\Delta t,T}  = \varphi _{\Delta t,T^t }  \circ \varphi
999 _{\Delta t,T^r }.
1000 \]
1001 The free rigid body is an example of Lie-Poisson system with
1002 Hamiltonian function
1003 \begin{equation}
1004 T^r (\pi ) = T_1 ^r (\pi _1 ) + T_2^r (\pi _2 ) + T_3^r (\pi _3 )
1005 \label{introEquation:rotationalKineticRB}
1006 \end{equation}
1007 where $T_i^r (\pi _i ) = \frac{{\pi _i ^2 }}{{2I_i }}$ and
1008 Lie-Poisson structure matrix,
1009 \begin{equation}
1010 J(\pi ) = \left( {\begin{array}{*{20}c}
1011   0 & {\pi _3 } & { - \pi _2 }  \\
1012   { - \pi _3 } & 0 & {\pi _1 }  \\
1013   {\pi _2 } & { - \pi _1 } & 0  \\
1014 \end{array}} \right)
1015 \end{equation}
1016 Thus, the dynamics of free rigid body is governed by
1017 \begin{equation}
1018 \frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi )
1019 \end{equation}
1020 One may notice that each $T_i^r$ in Equation
1021 \ref{introEquation:rotationalKineticRB} can be solved exactly. For
1022 instance, the equations of motion due to $T_1^r$ are given by
1023 \begin{equation}
1024 \frac{d}{{dt}}\pi  = R_1 \pi ,\frac{d}{{dt}}A = AR_1
1025 \label{introEqaution:RBMotionSingleTerm}
1026 \end{equation}
1027 where
1028 \[ R_1  = \left( {\begin{array}{*{20}c}
1029   0 & 0 & 0  \\
1030   0 & 0 & {\pi _1 }  \\
1031   0 & { - \pi _1 } & 0  \\
1032 \end{array}} \right).
1033 \]
1034 The solutions of Equation \ref{introEqaution:RBMotionSingleTerm} is
1035 \[
1036 \pi (\Delta t) = e^{\Delta tR_1 } \pi (0),A(\Delta t) =
1037 A(0)e^{\Delta tR_1 }
1038 \]
1039 with
1040 \[
1041 e^{\Delta tR_1 }  = \left( {\begin{array}{*{20}c}
1042   0 & 0 & 0  \\
1043   0 & {\cos \theta _1 } & {\sin \theta _1 }  \\
1044   0 & { - \sin \theta _1 } & {\cos \theta _1 }  \\
1045 \end{array}} \right),\theta _1  = \frac{{\pi _1 }}{{I_1 }}\Delta t.
1046 \]
1047 To reduce the cost of computing expensive functions in $e^{\Delta
1048 tR_1 }$, we can use Cayley transformation,
1049 \[
1050 e^{\Delta tR_1 }  \approx (1 - \Delta tR_1 )^{ - 1} (1 + \Delta tR_1
1051 )
1052 \]
1053 The flow maps for $T_2^r$ and $T_3^r$ can be found in the same
1054 manner.
1055
1056 In order to construct a second-order symplectic method, we split the
1057 angular kinetic Hamiltonian function into five terms
1058 \[
1059 T^r (\pi ) = \frac{1}{2}T_1 ^r (\pi _1 ) + \frac{1}{2}T_2^r (\pi _2
1060 ) + T_3^r (\pi _3 ) + \frac{1}{2}T_2^r (\pi _2 ) + \frac{1}{2}T_1 ^r
1061 (\pi _1 )
1062 \].
1063 Concatenating flows corresponding to these five terms, we can obtain
1064 the flow map for free rigid body,
1065 \[
1066 \varphi _{\Delta t,T^r }  = \varphi _{\Delta t/2,\pi _1 }  \circ
1067 \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }
1068 \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi
1069 _1 }.
1070 \]
1071
1072 The equations of motion corresponding to potential energy and
1073 kinetic energy are listed in the below table,
1074 \begin{center}
1075 \begin{tabular}{|l|l|}
1076  \hline
1077  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
1078  Potential & Kinetic \\
1079  $\frac{{dq}}{{dt}} = \frac{p}{m}$ & $\frac{d}{{dt}}q = p$ \\
1080  $\frac{d}{{dt}}p =  - \frac{{\partial V}}{{\partial q}}$ & $ \frac{d}{{dt}}p = 0$ \\
1081  $\frac{d}{{dt}}Q = 0$ & $ \frac{d}{{dt}}Q = Qskew(I^{ - 1} j)$ \\
1082  $ \frac{d}{{dt}}\pi  = \sum\limits_i {\left( {Q^T F_i (r,Q)} \right) \times X_i }$ & $\frac{d}{{dt}}\pi  = \pi  \times I^{ - 1} \pi$\\
1083  \hline
1084 \end{tabular}
1085 \end{center}
1086
1087 Finally, we obtain the overall symplectic flow maps for free moving
1088 rigid body
1089 \begin{align*}
1090 \varphi _{\Delta t}  = &\varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau } \circ  \\
1091  &\varphi _{\Delta t,T^t }  \circ \varphi _{\Delta t/2,\pi _1 }  \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }  \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi _1 } \circ  \\
1092  &\varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .\\
1093 \label{introEquation:overallRBFlowMaps}
1094 \end{align*}
1095
1096 \subsection{Langevin dynamics for rigid particles of arbitrary shape}
1097
1098 Consider a Langevin equation of motions in generalized coordinates
1099 \begin{equation}
1100 M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)}  + F_{r,i} (t)
1101 \label{LDGeneralizedForm}
1102 \end{equation}
1103 where $M_i$ is a $6\times6$ generalized diagonal mass (include mass
1104 and moment of inertial) matrix and $V_i$ is a generalized velocity,
1105 $V_i = V_i(v_i,\omega _i)$. The right side of Eq.
1106 (\ref{LDGeneralizedForm}) consists of three generalized forces in
1107 lab-fixed frame, systematic force $F_{s,i}$, dissipative force
1108 $F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the
1109 system in Newtownian mechanics typically refers to lab-fixed frame,
1110 it is also convenient to handle the rotation of rigid body in
1111 body-fixed frame. Thus the friction and random forces are calculated
1112 in body-fixed frame and converted back to lab-fixed frame by:
1113 \[
1114 \begin{array}{l}
1115 F_{f,i}^l (t) = A^T F_{f,i}^b (t), \\
1116 F_{r,i}^l (t) = A^T F_{r,i}^b (t) \\
1117 \end{array}.
1118 \]
1119 Here, the body-fixed friction force $F_{r,i}^b$ is proportional to
1120 the body-fixed velocity at center of resistance $v_{R,i}^b$ and
1121 angular velocity $\omega _i$,
1122 \begin{equation}
1123 F_{r,i}^b (t) = \left( \begin{array}{l}
1124 f_{r,i}^b (t) \\
1125 \tau _{r,i}^b (t) \\
1126 \end{array} \right) =  - \left( {\begin{array}{*{20}c}
1127   {\Xi _{R,t} } & {\Xi _{R,c}^T }  \\
1128   {\Xi _{R,c} } & {\Xi _{R,r} }  \\
1129 \end{array}} \right)\left( \begin{array}{l}
1130 v_{R,i}^b (t) \\
1131 \omega _i (t) \\
1132 \end{array} \right),
1133 \end{equation}
1134 while the random force $F_{r,i}^l$ is a Gaussian stochastic variable
1135 with zero mean and variance
1136 \begin{equation}
1137 \left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle  =
1138 \left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle  =
1139 2k_B T\Xi _R \delta (t - t').
1140 \end{equation}
1141 The equation of motion for $v_i$ can be written as
1142 \begin{equation}
1143 m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) +
1144 f_{r,i}^l (t)
1145 \end{equation}
1146 Since the frictional force is applied at the center of resistance
1147 which generally does not coincide with the center of mass, an extra
1148 torque is exerted at the center of mass. Thus, the net body-fixed
1149 frictional torque at the center of mass, $\tau _{n,i}^b (t)$, is
1150 given by
1151 \begin{equation}
1152 \tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b
1153 \end{equation}
1154 where $r_{MR}$ is the vector from the center of mass to the center
1155 of the resistance. Instead of integrating angular velocity in
1156 lab-fixed frame, we consider the equation of motion of angular
1157 momentum in body-fixed frame
1158 \begin{equation}
1159 \dot \pi _i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b
1160 (t) + \tau _{r,i}^b(t)
1161 \end{equation}
1162
1163 Embedding the friction terms into force and torque, one can
1164 integrate the langevin equations of motion for rigid body of
1165 arbitrary shape in a velocity-Verlet style 2-part algorithm, where
1166 $h= \delta t$:
1167
1168 {\tt part one:}
1169 \begin{align*}
1170 v_i (t + h/2) &\leftarrow v_i (t) + \frac{{hf_{t,i}^l (t)}}{{2m_i }} \\
1171 \pi _i (t + h/2) &\leftarrow \pi _i (t) + \frac{{h\tau _{t,i}^b (t)}}{2} \\
1172 r_i (t + h) &\leftarrow r_i (t) + hv_i (t + h/2) \\
1173 A_i (t + h) &\leftarrow rotate\left( {h\pi _i (t + h/2)I^{ - 1} } \right) \\
1174 \end{align*}
1175 In this context, the $\mathrm{rotate}$ function is the reversible
1176 product of five consecutive body-fixed rotations,
1177 \begin{equation}
1178 \mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
1179 \mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y
1180 / 2) \cdot \mathsf{G}_x(a_x /2),
1181 \end{equation}
1182 where each rotational propagator, $\mathsf{G}_\alpha(\theta)$,
1183 rotates both the rotation matrix ($\mathsf{A}$) and the body-fixed
1184 angular momentum ($\pi$) by an angle $\theta$ around body-fixed axis
1185 $\alpha$,
1186 \begin{equation}
1187 \mathsf{G}_\alpha( \theta ) = \left\{
1188 \begin{array}{lcl}
1189 \mathsf{A}(t) & \leftarrow & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
1190 {\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf
1191 j}(0).
1192 \end{array}
1193 \right.
1194 \end{equation}
1195 $\mathsf{R}_\alpha$ is a quadratic approximation to the single-axis
1196 rotation matrix.  For example, in the small-angle limit, the
1197 rotation matrix around the body-fixed x-axis can be approximated as
1198 \begin{equation}
1199 \mathsf{R}_x(\theta) \approx \left(
1200 \begin{array}{ccc}
1201 1 & 0 & 0 \\
1202 0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}  & -\frac{\theta}{1+
1203 \theta^2 / 4} \\
1204 0 & \frac{\theta}{1+ \theta^2 / 4} & \frac{1-\theta^2 / 4}{1 +
1205 \theta^2 / 4}
1206 \end{array}
1207 \right).
1208 \end{equation}
1209 All other rotations follow in a straightforward manner.
1210
1211 After the first part of the propagation, the friction and random
1212 forces are generated at the center of resistance in body-fixed frame
1213 and converted back into lab-fixed frame
1214 \[
1215 f_{t,i}^l (t + h) =  - \left( {\frac{{\partial V}}{{\partial r_i }}}
1216 \right)_{r_i (t + h)}  + A_i^T (t + h)[F_{f,i}^b (t + h) + F_{r,i}^b
1217 (t + h)],
1218 \]
1219 while the system torque in lab-fixed frame is transformed into
1220 body-fixed frame,
1221 \[
1222 \tau _{t,i}^b (t + h) = A\tau _{s,i}^l (t) + \tau _{n,i}^b (t) +
1223 \tau _{r,i}^b (t).
1224 \]
1225 Once the forces and torques have been obtained at the new time step,
1226 the velocities can be advanced to the same time value.
1227
1228 {\tt part two:}
1229 \begin{align*}
1230 v_i (t) &\leftarrow v_i (t + h/2) + \frac{{hf_{t,i}^l (t + h)}}{{2m_i }} \\
1231 \pi _i (t) &\leftarrow \pi _i (t + h/2) + \frac{{h\tau _{t,i}^b (t + h)}}{2} \\
1232 \end{align*}
1233
1234 \subsection{Results and discussion}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines