--- trunk/tengDissertation/Methodology.tex 2006/04/26 01:27:56 2739 +++ trunk/tengDissertation/Methodology.tex 2006/06/23 21:33:52 2882 @@ -2,37 +2,38 @@ In order to mimic the experiments, which are usually p \section{\label{methodSection:rigidBodyIntegrators}Integrators for Rigid Body Motion in Molecular Dynamics} -In order to mimic the experiments, which are usually performed under +In order to mimic experiments which are usually performed under constant temperature and/or pressure, extended Hamiltonian system methods have been developed to generate statistical ensembles, such -as canonical ensemble and isobaric-isothermal ensemble \textit{etc}. -In addition to the standard ensemble, specific ensembles have been -developed to account for the anisotropy between the lateral and -normal directions of membranes. The $NPAT$ ensemble, in which the -normal pressure and the lateral surface area of the membrane are -kept constant, and the $NP\gamma T$ ensemble, in which the normal -pressure and the lateral surface tension are kept constant were -proposed to address this issue. +as the canonical and isobaric-isothermal ensembles. In addition to +the standard ensemble, specific ensembles have been developed to +account for the anisotropy between the lateral and normal directions +of membranes. The $NPAT$ ensemble, in which the normal pressure and +the lateral surface area of the membrane are kept constant, and the +$NP\gamma T$ ensemble, in which the normal pressure and the lateral +surface tension are kept constant were proposed to address the +issues. -Integration schemes for rotational motion of the rigid molecules in -microcanonical ensemble have been extensively studied in the last -two decades. Matubayasi and Nakahara developed a time-reversible +Integration schemes for the rotational motion of the rigid molecules +in the microcanonical ensemble have been extensively studied over +the last two decades. Matubayasi developed a time-reversible integrator for rigid bodies in quaternion representation. Although it is not symplectic, this integrator still demonstrates a better -long-time energy conservation than traditional methods because of -the time-reversible nature. Extending Trotter-Suzuki to general -system with a flat phase space, Miller and his colleagues devised an -novel symplectic, time-reversible and volume-preserving integrator -in quaternion representation, which was shown to be superior to the -time-reversible integrator of Matubayasi and Nakahara. However, all -of the integrators in quaternion representation suffer from the -computational penalty of constructing a rotation matrix from -quaternions to evolve coordinates and velocities at every time step. -An alternative integration scheme utilizing rotation matrix directly -proposed by Dullweber, Leimkuhler and McLachlan (DLM) also preserved -the same structural properties of the Hamiltonian flow. In this -section, the integration scheme of DLM method will be reviewed and -extended to other ensembles. +long-time energy conservation than Euler angle methods because of +the time-reversible nature. Extending the Trotter-Suzuki +factorization to general system with a flat phase space, Miller and +his colleagues devised a novel symplectic, time-reversible and +volume-preserving integrator in the quaternion representation, which +was shown to be superior to the Matubayasi's time-reversible +integrator. However, all of the integrators in the quaternion +representation suffer from the computational penalty of constructing +a rotation matrix from quaternions to evolve coordinates and +velocities at every time step. An alternative integration scheme +utilizing the rotation matrix directly proposed by Dullweber, +Leimkuhler and McLachlan (DLM) also preserved the same structural +properties of the Hamiltonian flow. In this section, the integration +scheme of DLM method will be reviewed and extended to other +ensembles. \subsection{\label{methodSection:DLM}Integrating the Equations of Motion: the DLM method} @@ -111,13 +112,13 @@ torques are calculated at the new positions and orient - \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)}, \\ % {\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h) - \times \frac{\partial V}{\partial {\bf u}}, \\ + \times (\frac{\partial V}{\partial {\bf u}})_{u(t+h)}, \\ % {\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h) \cdot {\bf \tau}^s(t + h). \end{align*} -{\sc oopse} automatically updates ${\bf u}$ when the rotation matrix +${\bf u}$ is automatically updated when the rotation matrix $\mathsf{A}$ is calculated in {\tt moveA}. Once the forces and torques have been obtained at the new time step, the velocities can be advanced to the same time value. @@ -139,7 +140,7 @@ in Fig.~\ref{timestep}. average 7\% increase in computation time using the DLM method in place of quaternions. This cost is more than justified when comparing the energy conservation of the two methods as illustrated -in Fig.~\ref{timestep}. +in Fig.~\ref{methodFig:timestep}. \begin{figure} \centering @@ -150,25 +151,26 @@ from the true energy baseline for clarity.} \label{tim increasing time step. For each time step, the dotted line is total energy using the DLM integrator, and the solid line comes from the quaternion integrator. The larger time step plots are shifted up -from the true energy baseline for clarity.} \label{timestep} +from the true energy baseline for clarity.} +\label{methodFig:timestep} \end{figure} -In Fig.~\ref{timestep}, the resulting energy drift at various time -steps for both the DLM and quaternion integration schemes is -compared. All of the 1000 molecule water simulations started with -the same configuration, and the only difference was the method for -handling rotational motion. At time steps of 0.1 and 0.5 fs, both -methods for propagating molecule rotation conserve energy fairly -well, with the quaternion method showing a slight energy drift over -time in the 0.5 fs time step simulation. At time steps of 1 and 2 -fs, the energy conservation benefits of the DLM method are clearly -demonstrated. Thus, while maintaining the same degree of energy -conservation, one can take considerably longer time steps, leading -to an overall reduction in computation time. +In Fig.~\ref{methodFig:timestep}, the resulting energy drift at +various time steps for both the DLM and quaternion integration +schemes is compared. All of the 1000 molecule water simulations +started with the same configuration, and the only difference was the +method for handling rotational motion. At time steps of 0.1 and 0.5 +fs, both methods for propagating molecule rotation conserve energy +fairly well, with the quaternion method showing a slight energy +drift over time in the 0.5 fs time step simulation. At time steps of +1 and 2 fs, the energy conservation benefits of the DLM method are +clearly demonstrated. Thus, while maintaining the same degree of +energy conservation, one can take considerably longer time steps, +leading to an overall reduction in computation time. \subsection{\label{methodSection:NVT}Nos\'{e}-Hoover Thermostatting} -The Nos\'e-Hoover equations of motion are given by\cite{Hoover85} +The Nos\'e-Hoover equations of motion are given by\cite{Hoover1985} \begin{eqnarray} \dot{{\bf r}} & = & {\bf v}, \\ \dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} ,\\ @@ -197,7 +199,9 @@ and $K$ is the total kinetic energy, \begin{equation} f = 3 N + 3 N_{\mathrm{orient}} - N_{\mathrm{constraints}}, \end{equation} -and $K$ is the total kinetic energy, +where $N_{\mathrm{orient}}$ is the number of molecules with +orientational degrees of freedom, and $K$ is the total kinetic +energy, \begin{equation} K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i + \sum_{i=1}^{N_{\mathrm{orient}}} \frac{1}{2} {\bf j}_i^T \cdot @@ -205,13 +209,9 @@ relaxation of the temperature to the target value. To \end{equation} In eq.(\ref{eq:nosehooverext}), $\tau_T$ is the time constant for -relaxation of the temperature to the target value. To set values -for $\tau_T$ or $T_{\mathrm{target}}$ in a simulation, one would use -the {\tt tauThermostat} and {\tt targetTemperature} keywords in the -{\tt .bass} file. The units for {\tt tauThermostat} are fs, and the -units for the {\tt targetTemperature} are degrees K. The -integration of the equations of motion is carried out in a -velocity-Verlet style 2 part algorithm: +relaxation of the temperature to the target value. The integration +of the equations of motion is carried out in a velocity-Verlet style +2 part algorithm: {\tt moveA:} \begin{align*} @@ -277,14 +277,11 @@ self-consistent. The relative tolerance for the self- caclculate $T(t + h)$ as well as $\chi(t + h)$, they indirectly depend on their own values at time $t + h$. {\tt moveB} is therefore done in an iterative fashion until $\chi(t + h)$ becomes -self-consistent. The relative tolerance for the self-consistency -check defaults to a value of $\mbox{10}^{-6}$, but {\sc oopse} will -terminate the iteration after 4 loops even if the consistency check -has not been satisfied. +self-consistent. The Nos\'e-Hoover algorithm is known to conserve a Hamiltonian for the extended system that is, to within a constant, identical to the -Helmholtz free energy,\cite{melchionna93} +Helmholtz free energy,\cite{Melchionna1993} \begin{equation} H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left( \frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) @@ -295,16 +292,12 @@ Bond constraints are applied at the end of both the {\ last column of the {\tt .stat} file to allow checks on the quality of the integration. -Bond constraints are applied at the end of both the {\tt moveA} and -{\tt moveB} portions of the algorithm. Details on the constraint -algorithms are given in section \ref{oopseSec:rattle}. - \subsection{\label{methodSection:NPTi}Constant-pressure integration with isotropic box deformations (NPTi)} -To carry out isobaric-isothermal ensemble calculations {\sc oopse} -implements the Melchionna modifications to the -Nos\'e-Hoover-Andersen equations of motion,\cite{melchionna93} +We can used an isobaric-isothermal ensemble integrator which is +implemented using the Melchionna modifications to the +Nos\'e-Hoover-Andersen equations of motion,\cite{Melchionna1993} \begin{eqnarray} \dot{{\bf r}} & = & {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right), \\ @@ -359,11 +352,7 @@ relaxation of the pressure to the target value. To se \end{equation} In eq.(\ref{eq:melchionna1}), $\tau_B$ is the time constant for -relaxation of the pressure to the target value. To set values for -$\tau_B$ or $P_{\mathrm{target}}$ in a simulation, one would use the -{\tt tauBarostat} and {\tt targetPressure} keywords in the {\tt -.bass} file. The units for {\tt tauBarostat} are fs, and the units -for the {\tt targetPressure} are atmospheres. Like in the NVT +relaxation of the pressure to the target value. Like in the NVT integrator, the integration of the equations of motion is carried out in a velocity-Verlet style 2 part algorithm: @@ -404,12 +393,11 @@ depends on the positions at the same time. {\sc oopse Most of these equations are identical to their counterparts in the NVT integrator, but the propagation of positions to time $t + h$ -depends on the positions at the same time. {\sc oopse} carries out -this step iteratively (with a limit of 5 passes through the -iterative loop). Also, the simulation box $\mathsf{H}$ is scaled -uniformly for one full time step by an exponential factor that -depends on the value of $\eta$ at time $t + h / 2$. Reshaping the -box uniformly also scales the volume of the box by +depends on the positions at the same time. The simulation box +$\mathsf{H}$ is scaled uniformly for one full time step by an +exponential factor that depends on the value of $\eta$ at time $t + +h / 2$. Reshaping the box uniformly also scales the volume of the +box by \begin{equation} \mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)}. \mathcal{V}(t) @@ -451,10 +439,7 @@ + h)$ and $\eta(t + h)$ become self-consistent. The r to caclculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, and $\eta(t + h)$, they indirectly depend on their own values at time $t + h$. {\tt moveB} is therefore done in an iterative fashion until $\chi(t -+ h)$ and $\eta(t + h)$ become self-consistent. The relative -tolerance for the self-consistency check defaults to a value of -$\mbox{10}^{-6}$, but {\sc oopse} will terminate the iteration after -4 loops even if the consistency check has not been satisfied. ++ h)$ and $\eta(t + h)$ become self-consistent. The Melchionna modification of the Nos\'e-Hoover-Andersen algorithm is known to conserve a Hamiltonian for the extended system that is, @@ -476,10 +461,6 @@ Bond constraints are applied at the end of both the {\ P_{\mathrm{target}} \mathcal{V}(t). \end{equation} -Bond constraints are applied at the end of both the {\tt moveA} and -{\tt moveB} portions of the algorithm. Details on the constraint -algorithms are given in section \ref{oopseSec:rattle}. - \subsection{\label{methodSection:NPTf}Constant-pressure integration with a flexible box (NPTf)} @@ -553,8 +534,8 @@ r}(t)\right\}, \mathsf{H}(t + h) &\leftarrow \mathsf{H}(t) \cdot e^{-h \overleftrightarrow{\eta}(t + h / 2)} . \end{align*} -{\sc oopse} uses a power series expansion truncated at second order -for the exponential operation which scales the simulation box. +Here, a power series expansion truncated at second order for the +exponential operation is used to scale the simulation box. The {\tt moveB} portion of the algorithm is largely unchanged from the NPTi integrator: @@ -593,13 +574,14 @@ The NPTf integrator is known to conserve the following identical to those described for the NPTi integrator. The NPTf integrator is known to conserve the following Hamiltonian: -\begin{equation} -H_{\mathrm{NPTf}} = V + K + f k_B T_{\mathrm{target}} \left( +\begin{eqnarray*} +H_{\mathrm{NPTf}} & = & V + K + f k_B T_{\mathrm{target}} \left( \frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) -dt^\prime \right) + P_{\mathrm{target}} \mathcal{V}(t) + \frac{f k_B +dt^\prime \right) \\ + & & + P_{\mathrm{target}} \mathcal{V}(t) + \frac{f k_B T_{\mathrm{target}}}{2} \mathrm{Tr}\left[\overleftrightarrow{\eta}(t)\right]^2 \tau_B^2. -\end{equation} +\end{eqnarray*} This integrator must be used with care, particularly in liquid simulations. Liquids have very small restoring forces in the @@ -609,31 +591,32 @@ assume non-orthorhombic geometries. finds most use in simulating crystals or liquid crystals which assume non-orthorhombic geometries. -\subsection{\label{methodSection:otherSpecialEnsembles}Other Special Ensembles} +\subsection{\label{methodSection:NPAT}NPAT Ensemble} -\subsubsection{\label{methodSection:NPAT}Constant Normal Pressure, Constant Lateral Surface Area and Constant Temperature (NPAT) Ensemble} +A comprehensive understanding of relations between structures and +functions in biological membrane system ultimately relies on +structure and dynamics of lipid bilayers, which are strongly +affected by the interfacial interaction between lipid molecules and +surrounding media. One quantity to describe the interfacial +interaction is so called the average surface area per lipid. +Constant area and constant lateral pressure simulations can be +achieved by extending the standard NPT ensemble with a different +pressure control strategy -A comprehensive understanding of structure¨Cfunction relations of -biological membrane system ultimately relies on structure and -dynamics of lipid bilayer, which are strongly affected by the -interfacial interaction between lipid molecules and surrounding -media. One quantity to describe the interfacial interaction is so -called the average surface area per lipid. Constat area and constant -lateral pressure simulation can be achieved by extending the -standard NPT ensemble with a different pressure control strategy \begin{equation} -\dot -\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} -\over \eta } _{\alpha \beta } = \left\{ \begin{array}{l} - \frac{{V(P_{\alpha \beta } - P_{{\rm{target}}} )}}{{\tau _{\rm{B}}^{\rm{2}} fk_B T_{{\rm{target}}} }}{\rm{ }}(\alpha = \beta = z), \\ - 0{\rm{ }}(\alpha \ne z{\rm{ }}or{\rm{ }}\beta \ne z) \\ - \end{array} \right. -\label{methodEquation:NPATeta} + \dot {\overleftrightarrow{\eta}} _{\alpha \beta}=\left\{\begin{array}{ll} + \frac{{V(P_{\alpha \beta } - P_{{\rm{target}}} )}}{{\tau_{\rm{B}}^{\rm{2}} fk_B T_{{\rm{target}}} }} + & \mbox{if $ \alpha = \beta = z$}\\ + 0 & \mbox{otherwise}\\ + \end{array} + \right. \end{equation} + Note that the iterative schemes for NPAT are identical to those described for the NPTi integrator. -\subsubsection{\label{methodSection:NPrT}Constant Normal Pressure, Constant Lateral Surface Tension and Constant Temperature (NP\gamma T) Ensemble } +\subsection{\label{methodSection:NPrT}NP$\gamma$T +Ensemble} Theoretically, the surface tension $\gamma$ of a stress free membrane system should be zero since its surface free energy $G$ is @@ -643,26 +626,23 @@ the membrane simulation, a special ensemble, NP\gamma \] However, a surface tension of zero is not appropriate for relatively small patches of membrane. In order to eliminate the edge effect of -the membrane simulation, a special ensemble, NP\gamma T, is proposed -to maintain the lateral surface tension and normal pressure. The -equation of motion for cell size control tensor, $\eta$, in NP\gamma -T is +the membrane simulation, a special ensemble, NP$\gamma$T, has been +proposed to maintain the lateral surface tension and normal +pressure. The equation of motion for the cell size control tensor, +$\eta$, in $NP\gamma T$ is \begin{equation} -\dot -\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} -\over \eta } _{\alpha \beta } = \left\{ \begin{array}{l} - - A_{xy} (\gamma _\alpha - \gamma _{{\rm{target}}} ){\rm{ (}}\alpha = \beta = x{\rm{ or }} = y{\rm{)}} \\ - \frac{{V(P_{\alpha \beta } - P_{{\rm{target}}} )}}{{\tau _{\rm{B}}^{\rm{2}} fk_B T_{{\rm{target}}} }}{\rm{ }}(\alpha = \beta = z) \\ - 0{\rm{ }}(\alpha \ne \beta ) \\ - \end{array} \right. -\label{methodEquation:NPrTeta} + \dot {\overleftrightarrow{\eta}} _{\alpha \beta}=\left\{\begin{array}{ll} + - A_{xy} (\gamma _\alpha - \gamma _{{\rm{target}}} ) & \mbox{$\alpha = \beta = x$ or $\alpha = \beta = y$}\\ + \frac{{V(P_{\alpha \beta } - P_{{\rm{target}}} )}}{{\tau _{\rm{B}}^{\rm{2}} fk_B T_{{\rm{target}}}}} & \mbox{$\alpha = \beta = z$} \\ + 0 & \mbox{$\alpha \ne \beta$} \\ + \end{array} + \right. \end{equation} where $ \gamma _{{\rm{target}}}$ is the external surface tension and the instantaneous surface tensor $\gamma _\alpha$ is given by \begin{equation} -\gamma _\alpha = - h_z -(\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} -\over P} _{\alpha \alpha } - P_{{\rm{target}}} ) +\gamma _\alpha = - h_z( \overleftrightarrow{P} _{\alpha \alpha } +- P_{{\rm{target}}} ) \label{methodEquation:instantaneousSurfaceTensor} \end{equation} @@ -672,16 +652,94 @@ $\gamma$ is set to zero. axes are allowed to fluctuate independently, but the angle between the box axes does not change. It should be noted that the NPTxyz integrator is a special case of $NP\gamma T$ if the surface tension -$\gamma$ is set to zero. +$\gamma$ is set to zero, and if $x$ and $y$ can move independently. -\section{\label{methodSection:langevin}Integrators for Langevin Dynamics of Rigid Bodies} +\section{\label{methodSection:zcons}The Z-Constraint Method} -\subsection{\label{methodSection:temperature}Temperature Control} +Based on the fluctuation-dissipation theorem, a force +auto-correlation method was developed by Roux and Karplus to +investigate the dynamics of ions inside ion channels\cite{Roux1991}. +The time-dependent friction coefficient can be calculated from the +deviation of the instantaneous force from its mean force. +\begin{equation} +\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T, +\end{equation} +where% +\begin{equation} +\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle. +\end{equation} -\subsection{\label{methodSection:pressureControl}Pressure Control} +If the time-dependent friction decays rapidly, the static friction +coefficient can be approximated by +\begin{equation} +\xi_{\text{static}}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta +F(z,0)\rangle dt. +\end{equation} +Allowing diffusion constant to then be calculated through the +Einstein relation:\cite{Marrink1994} +\begin{equation} +D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty +}\langle\delta F(z,t)\delta F(z,0)\rangle dt}.% +\end{equation} -\section{\label{methodSection:hydrodynamics}Hydrodynamics} +The Z-Constraint method, which fixes the z coordinates of the +molecules with respect to the center of the mass of the system, has +been a method suggested to obtain the forces required for the force +auto-correlation calculation.\cite{Marrink1994} However, simply +resetting the coordinate will move the center of the mass of the +whole system. To avoid this problem, we reset the forces of +z-constrained molecules as well as subtract the total constraint +forces from the rest of the system after the force calculation at +each time step instead of resetting the coordinate. -%\section{\label{methodSection:coarseGrained}Coarse-Grained Modeling} +After the force calculation, we define $G_\alpha$ as +\begin{equation} +G_{\alpha} = \sum_i F_{\alpha i}, \label{oopseEq:zc1} +\end{equation} +where $F_{\alpha i}$ is the force in the z direction of atom $i$ in +z-constrained molecule $\alpha$. The forces of the z constrained +molecule are then set to: +\begin{equation} +F_{\alpha i} = F_{\alpha i} - + \frac{m_{\alpha i} G_{\alpha}}{\sum_i m_{\alpha i}}. +\end{equation} +Here, $m_{\alpha i}$ is the mass of atom $i$ in the z-constrained +molecule. Having rescaled the forces, the velocities must also be +rescaled to subtract out any center of mass velocity in the z +direction. +\begin{equation} +v_{\alpha i} = v_{\alpha i} - + \frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}, +\end{equation} +where $v_{\alpha i}$ is the velocity of atom $i$ in the z direction. +Lastly, all of the accumulated z constrained forces must be +subtracted from the system to keep the system center of mass from +drifting. +\begin{equation} +F_{\beta i} = F_{\beta i} - \frac{m_{\beta i} \sum_{\alpha} +G_{\alpha}} + {\sum_{\beta}\sum_i m_{\beta i}}, +\end{equation} +where $\beta$ are all of the unconstrained molecules in the system. +Similarly, the velocities of the unconstrained molecules must also +be scaled. +\begin{equation} +v_{\beta i} = v_{\beta i} + \sum_{\alpha} + \frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}. +\end{equation} -%\section{\label{methodSection:moleculeScale}Molecular-Scale Modeling} +At the very beginning of the simulation, the molecules may not be at +their constrained positions. To move a z-constrained molecule to its +specified position, a simple harmonic potential is used +\begin{equation} +U(t)=\frac{1}{2}k_{\text{Harmonic}}(z(t)-z_{\text{cons}})^{2},% +\end{equation} +where $k_{\text{Harmonic}}$ is the harmonic force constant, $z(t)$ +is the current $z$ coordinate of the center of mass of the +constrained molecule, and $z_{\text{cons}}$ is the constrained +position. The harmonic force operating on the z-constrained molecule +at time $t$ can be calculated by +\begin{equation} +F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}= + -k_{\text{Harmonic}}(z(t)-z_{\text{cons}}). +\end{equation}