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 2729 by tim, Mon Apr 24 18:49:04 2006 UTC vs.
Revision 2864 by tim, Fri Jun 16 01:53:48 2006 UTC

# Line 4 | Line 4 | methods have been developed to generate statistical en
4  
5   In order to mimic the experiments, which are usually performed under
6   constant temperature and/or pressure, extended Hamiltonian system
7 < methods have been developed to generate statistical ensemble, such
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
# Line 16 | Line 16 | two decades. Matubayasi and Nakahara developed a time-
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
19 > two decades. Matubayasi developed a time-reversible integrator for
20 > rigid bodies in quaternion representation. Although it is not
21 > symplectic, this integrator still demonstrates a better long-time
22 > energy conservation than traditional methods because of the
23 > time-reversible nature. Extending Trotter-Suzuki to general system
24 > with a flat phase space, Miller and his colleagues devised an novel
25 > symplectic, time-reversible and volume-preserving integrator in
26 > quaternion representation, which was shown to be superior to the
27 > Matubayasi's time-reversible integrator. However, all of the
28 > integrators in quaternion representation suffer from the
29   computational penalty of constructing a rotation matrix from
30   quaternions to evolve coordinates and velocities at every time step.
31   An alternative integration scheme utilizing rotation matrix directly
# Line 117 | Line 117 | torques are calculated at the new positions and orient
117      \cdot {\bf \tau}^s(t + h).
118   \end{align*}
119  
120 < {\sc oopse} automatically updates ${\bf u}$ when the rotation matrix
120 > ${\bf u}$ will be automatically updated when the rotation matrix
121   $\mathsf{A}$ is calculated in {\tt moveA}.  Once the forces and
122   torques have been obtained at the new time step, the velocities can
123   be advanced to the same time value.
# Line 139 | Line 139 | in Fig.~\ref{timestep}.
139   average 7\% increase in computation time using the DLM method in
140   place of quaternions. This cost is more than justified when
141   comparing the energy conservation of the two methods as illustrated
142 < in Fig.~\ref{timestep}.
142 > in Fig.~\ref{methodFig:timestep}.
143  
144   \begin{figure}
145   \centering
# Line 150 | Line 150 | from the true energy baseline for clarity.} \label{tim
150   increasing time step. For each time step, the dotted line is total
151   energy using the DLM integrator, and the solid line comes from the
152   quaternion integrator. The larger time step plots are shifted up
153 < from the true energy baseline for clarity.} \label{timestep}
153 > from the true energy baseline for clarity.}
154 > \label{methodFig:timestep}
155   \end{figure}
156  
157 < In Fig.~\ref{timestep}, the resulting energy drift at various time
158 < steps for both the DLM and quaternion integration schemes is
159 < compared. All of the 1000 molecule water simulations started with
160 < the same configuration, and the only difference was the method for
161 < handling rotational motion. At time steps of 0.1 and 0.5 fs, both
162 < methods for propagating molecule rotation conserve energy fairly
163 < well, with the quaternion method showing a slight energy drift over
164 < time in the 0.5 fs time step simulation. At time steps of 1 and 2
165 < fs, the energy conservation benefits of the DLM method are clearly
166 < demonstrated. Thus, while maintaining the same degree of energy
167 < conservation, one can take considerably longer time steps, leading
168 < to an overall reduction in computation time.
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.
169  
170   \subsection{\label{methodSection:NVT}Nos\'{e}-Hoover Thermostatting}
171  
172 < The Nos\'e-Hoover equations of motion are given by\cite{Hoover85}
172 > The Nos\'e-Hoover equations of motion are given by\cite{Hoover1985}
173   \begin{eqnarray}
174   \dot{{\bf r}} & = & {\bf v}, \\
175   \dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} ,\\
# Line 205 | Line 206 | relaxation of the temperature to the target value.  To
206   \end{equation}
207  
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
211 < {\tt .bass} file.  The units for {\tt tauThermostat} are fs, and the
212 < units for the {\tt targetTemperature} are degrees K.   The
213 < integration of the equations of motion is carried out in a
214 < velocity-Verlet style 2 part algorithm:
209 > relaxation of the temperature to the target value. The integration
210 > of the equations of motion is carried out in a velocity-Verlet style
211 > 2 part algorithm:
212  
213   {\tt moveA:}
214   \begin{align*}
# Line 277 | Line 274 | self-consistent.  The relative tolerance for the self-
274   caclculate $T(t + h)$ as well as $\chi(t + h)$, they indirectly
275   depend on their own values at time $t + h$.  {\tt moveB} is
276   therefore done in an iterative fashion until $\chi(t + h)$ becomes
277 < self-consistent.  The relative tolerance for the self-consistency
281 < check defaults to a value of $\mbox{10}^{-6}$, but {\sc oopse} will
282 < terminate the iteration after 4 loops even if the consistency check
283 < has not been satisfied.
277 > self-consistent.
278  
279   The Nos\'e-Hoover algorithm is known to conserve a Hamiltonian for
280   the extended system that is, to within a constant, identical to the
281 < Helmholtz free energy,\cite{melchionna93}
281 > Helmholtz free energy,\cite{Melchionna1993}
282   \begin{equation}
283   H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left(
284   \frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime)
# Line 295 | Line 289 | Bond constraints are applied at the end of both the {\
289   last column of the {\tt .stat} file to allow checks on the quality
290   of the integration.
291  
298 Bond constraints are applied at the end of both the {\tt moveA} and
299 {\tt moveB} portions of the algorithm.  Details on the constraint
300 algorithms are given in section \ref{oopseSec:rattle}.
301
292   \subsection{\label{methodSection:NPTi}Constant-pressure integration with
293   isotropic box deformations (NPTi)}
294  
295 < To carry out isobaric-isothermal ensemble calculations {\sc oopse}
296 < implements the Melchionna modifications to the
297 < Nos\'e-Hoover-Andersen equations of motion,\cite{melchionna93}
295 > Isobaric-isothermal ensemble integrator is implemented using the
296 > Melchionna modifications to the Nos\'e-Hoover-Andersen equations of
297 > motion,\cite{Melchionna1993}
298  
299   \begin{eqnarray}
300   \dot{{\bf r}} & = & {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right), \\
# Line 359 | Line 349 | relaxation of the pressure to the target value.  To se
349   \end{equation}
350  
351   In eq.(\ref{eq:melchionna1}), $\tau_B$ is the time constant for
352 < relaxation of the pressure to the target value.  To set values for
363 < $\tau_B$ or $P_{\mathrm{target}}$ in a simulation, one would use the
364 < {\tt tauBarostat} and {\tt targetPressure} keywords in the {\tt
365 < .bass} file.  The units for {\tt tauBarostat} are fs, and the units
366 < for the {\tt targetPressure} are atmospheres.  Like in the NVT
352 > relaxation of the pressure to the target value. Like in the NVT
353   integrator, the integration of the equations of motion is carried
354   out in a velocity-Verlet style 2 part algorithm:
355  
# Line 404 | Line 390 | depends on the positions at the same time.  {\sc oopse
390  
391   Most of these equations are identical to their counterparts in the
392   NVT integrator, but the propagation of positions to time $t + h$
393 < depends on the positions at the same time.  {\sc oopse} carries out
394 < this step iteratively (with a limit of 5 passes through the
395 < iterative loop).  Also, the simulation box $\mathsf{H}$ is scaled
396 < uniformly for one full time step by an exponential factor that
397 < depends on the value of $\eta$ at time $t + h / 2$.  Reshaping the
412 < box uniformly also scales the volume of the box by
393 > depends on the positions at the same time. The simulation box
394 > $\mathsf{H}$ is scaled uniformly for one full time step by an
395 > exponential factor that depends on the value of $\eta$ at time $t +
396 > h / 2$.  Reshaping the box uniformly also scales the volume of the
397 > box by
398   \begin{equation}
399   \mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)}.
400   \mathcal{V}(t)
# Line 451 | Line 436 | + h)$ and $\eta(t + h)$ become self-consistent.  The r
436   to caclculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, and $\eta(t +
437   h)$, they indirectly depend on their own values at time $t + h$.
438   {\tt moveB} is therefore done in an iterative fashion until $\chi(t
439 < + h)$ and $\eta(t + h)$ become self-consistent.  The relative
455 < tolerance for the self-consistency check defaults to a value of
456 < $\mbox{10}^{-6}$, but {\sc oopse} will terminate the iteration after
457 < 4 loops even if the consistency check has not been satisfied.
439 > + h)$ and $\eta(t + h)$ become self-consistent.
440  
441   The Melchionna modification of the Nos\'e-Hoover-Andersen algorithm
442   is known to conserve a Hamiltonian for the extended system that is,
# Line 476 | Line 458 | Bond constraints are applied at the end of both the {\
458   P_{\mathrm{target}} \mathcal{V}(t).
459   \end{equation}
460  
479 Bond constraints are applied at the end of both the {\tt moveA} and
480 {\tt moveB} portions of the algorithm.  Details on the constraint
481 algorithms are given in section \ref{oopseSec:rattle}.
482
461   \subsection{\label{methodSection:NPTf}Constant-pressure integration with a
462   flexible box (NPTf)}
463  
# Line 553 | Line 531 | r}(t)\right\},
531   \mathsf{H}(t + h) &\leftarrow \mathsf{H}(t) \cdot e^{-h
532      \overleftrightarrow{\eta}(t + h / 2)} .
533   \end{align*}
534 < {\sc oopse} uses a power series expansion truncated at second order
535 < for the exponential operation which scales the simulation box.
534 > Here, a power series expansion truncated at second order for the
535 > exponential operation is used to scale the simulation box.
536  
537   The {\tt moveB} portion of the algorithm is largely unchanged from
538   the NPTi integrator:
# Line 593 | Line 571 | The NPTf integrator is known to conserve the following
571   identical to those described for the NPTi integrator.
572  
573   The NPTf integrator is known to conserve the following Hamiltonian:
574 < \begin{equation}
575 < H_{\mathrm{NPTf}} = V + K + f k_B T_{\mathrm{target}} \left(
574 > \begin{eqnarray*}
575 > H_{\mathrm{NPTf}} & = & V + K + f k_B T_{\mathrm{target}} \left(
576   \frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime)
577 < dt^\prime \right) + P_{\mathrm{target}} \mathcal{V}(t) + \frac{f k_B
577 > dt^\prime \right) \\
578 > & & + P_{\mathrm{target}} \mathcal{V}(t) + \frac{f k_B
579   T_{\mathrm{target}}}{2}
580   \mathrm{Tr}\left[\overleftrightarrow{\eta}(t)\right]^2 \tau_B^2.
581 < \end{equation}
581 > \end{eqnarray*}
582  
583   This integrator must be used with care, particularly in liquid
584   simulations.  Liquids have very small restoring forces in the
# Line 609 | Line 588 | assume non-orthorhombic geometries.
588   finds most use in simulating crystals or liquid crystals which
589   assume non-orthorhombic geometries.
590  
591 < \subsection{\label{methodSection:NPAT}Constant Lateral Pressure and Constant Surface Area (NPAT)}
591 > \subsection{\label{methodSection:NPAT}NPAT Ensemble}
592  
593 < \subsection{\label{methodSection:NPrT}Constant lateral Pressure and Constant Surface Tension (NP\gamma T) }
593 > A comprehensive understanding of structure¨Cfunction relations of
594 > biological membrane system ultimately relies on structure and
595 > dynamics of lipid bilayer, which are strongly affected by the
596 > interfacial interaction between lipid molecules and surrounding
597 > media. One quantity to describe the interfacial interaction is so
598 > called the average surface area per lipid. Constat area and constant
599 > lateral pressure simulation can be achieved by extending the
600 > standard NPT ensemble with a different pressure control strategy
601  
602 < \subsection{\label{methodSection:NPTxyz}Constant pressure in 3 axes (NPTxyz)}
602 > \begin{equation}
603 > \dot {\overleftrightarrow{\eta}} _{\alpha \beta}=\left\{\begin{array}{ll}
604 >                  \frac{{V(P_{\alpha \beta }  - P_{{\rm{target}}} )}}{{\tau_{\rm{B}}^{\rm{2}} fk_B T_{{\rm{target}}} }}
605 >                  & \mbox{if $ \alpha = \beta  = z$}\\
606 >                  0 & \mbox{otherwise}\\
607 >           \end{array}
608 >    \right.
609 > \end{equation}
610  
611 < There is one additional extended system integrator which is somewhat
612 < simpler than the NPTf method described above.  In this case, the
620 < three axes have independent barostats which each attempt to preserve
621 < the target pressure along the box walls perpendicular to that
622 < particular axis.  The lengths of the box axes are allowed to
623 < fluctuate independently, but the angle between the box axes does not
624 < change. The equations of motion are identical to those described
625 < above, but only the {\it diagonal} elements of
626 < $\overleftrightarrow{\eta}$ are computed.  The off-diagonal elements
627 < are set to zero (even when the pressure tensor has non-zero
628 < off-diagonal elements). It should be noted that the NPTxyz
629 < integrator is a special case of $NP\gamma T$ if the surface tension
630 < $\gamma$ is set to zero.
611 > Note that the iterative schemes for NPAT are identical to those
612 > described for the NPTi integrator.
613  
614 + \subsection{\label{methodSection:NPrT}NP$\gamma$T
615 + Ensemble}
616  
617 < \section{\label{methodSection:constraintMethods}Constraint Methods}
617 > Theoretically, the surface tension $\gamma$ of a stress free
618 > membrane system should be zero since its surface free energy $G$ is
619 > minimum with respect to surface area $A$
620 > \[
621 > \gamma  = \frac{{\partial G}}{{\partial A}}.
622 > \]
623 > However, a surface tension of zero is not appropriate for relatively
624 > small patches of membrane. In order to eliminate the edge effect of
625 > the membrane simulation, a special ensemble, NP$\gamma$T, is
626 > proposed to maintain the lateral surface tension and normal
627 > pressure. The equation of motion for cell size control tensor,
628 > $\eta$, in $NP\gamma T$ is
629 > \begin{equation}
630 > \dot {\overleftrightarrow{\eta}} _{\alpha \beta}=\left\{\begin{array}{ll}
631 >    - A_{xy} (\gamma _\alpha   - \gamma _{{\rm{target}}} ) & \mbox{$\alpha  = \beta  = x$ or $\alpha  = \beta  = y$}\\
632 >    \frac{{V(P_{\alpha \beta }  - P_{{\rm{target}}} )}}{{\tau _{\rm{B}}^{\rm{2}} fk_B T_{{\rm{target}}}}} & \mbox{$\alpha  = \beta  = z$} \\
633 >    0 & \mbox{$\alpha  \ne \beta$} \\
634 >       \end{array}
635 >    \right.
636 > \end{equation}
637 > where $ \gamma _{{\rm{target}}}$ is the external surface tension and
638 > the instantaneous surface tensor $\gamma _\alpha$ is given by
639 > \begin{equation}
640 > \gamma _\alpha   =  - h_z( \overleftrightarrow{P} _{\alpha \alpha }
641 > - P_{{\rm{target}}} )
642 > \label{methodEquation:instantaneousSurfaceTensor}
643 > \end{equation}
644  
645 < \subsection{\label{methodSection:rattle}The {\sc rattle} Method for Bond
646 <    Constraints}
645 > There is one additional extended system integrator (NPTxyz), in
646 > which each attempt to preserve the target pressure along the box
647 > walls perpendicular to that particular axis.  The lengths of the box
648 > axes are allowed to fluctuate independently, but the angle between
649 > the box axes does not change. It should be noted that the NPTxyz
650 > integrator is a special case of $NP\gamma T$ if the surface tension
651 > $\gamma$ is set to zero.
652  
653 < \subsection{\label{methodSection:zcons}Z-Constraint Method}
653 > \section{\label{methodSection:zcons}Z-Constraint Method}
654  
655   Based on the fluctuation-dissipation theorem, a force
656   auto-correlation method was developed by Roux and Karplus to
657 < investigate the dynamics of ions inside ion channels.\cite{Roux91}
657 > investigate the dynamics of ions inside ion channels\cite{Roux1991}.
658   The time-dependent friction coefficient can be calculated from the
659   deviation of the instantaneous force from its mean force.
660   \begin{equation}
# Line 657 | Line 672 | Einstein relation:\cite{Marrink94}
672   F(z,0)\rangle dt.
673   \end{equation}
674   Allowing diffusion constant to then be calculated through the
675 < Einstein relation:\cite{Marrink94}
675 > Einstein relation:\cite{Marrink1994}
676   \begin{equation}
677   D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty
678   }\langle\delta F(z,t)\delta F(z,0)\rangle dt}.%
# Line 666 | Line 681 | auto-correlation calculation.\cite{Marrink94} However,
681   The Z-Constraint method, which fixes the z coordinates of the
682   molecules with respect to the center of the mass of the system, has
683   been a method suggested to obtain the forces required for the force
684 < auto-correlation calculation.\cite{Marrink94} However, simply
684 > auto-correlation calculation.\cite{Marrink1994} However, simply
685   resetting the coordinate will move the center of the mass of the
686 < whole system. To avoid this problem, a new method was used in {\sc
672 < oopse}. Instead of resetting the coordinate, we reset the forces of
686 > whole system. To avoid this problem, we reset the forces of
687   z-constrained molecules as well as subtract the total constraint
688   forces from the rest of the system after the force calculation at
689 < each time step.
689 > each time step instead of resetting the coordinate.
690  
691   After the force calculation, define $G_\alpha$ as
692   \begin{equation}
# Line 725 | Line 739 | F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\part
739   F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}=
740      -k_{\text{Harmonic}}(z(t)-z_{\text{cons}}).
741   \end{equation}
728
729 \section{\label{methodSection:langevin}Integrators for Langevin Dynamics of Rigid Bodies}
730
731 \subsection{\label{methodSection:temperature}Temperature Control}
732
733 \subsection{\label{methodSection:pressureControl}Pressure Control}
734
735 \section{\label{methodSection:hydrodynamics}Hydrodynamics}
736
737 %\section{\label{methodSection:coarseGrained}Coarse-Grained Modeling}
738
739 %\section{\label{methodSection:moleculeScale}Molecular-Scale Modeling}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines