--- trunk/tengDissertation/Methodology.tex 2006/06/06 14:56:36 2801 +++ trunk/tengDissertation/Methodology.tex 2006/06/11 02:53:15 2854 @@ -16,16 +16,16 @@ two decades. Matubayasi and Nakahara developed a time- 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 -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 +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 +Matubayasi's time-reversible integrator. 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 @@ -117,7 +117,7 @@ torques are calculated at the new positions and orient \cdot {\bf \tau}^s(t + h). \end{align*} -{\sc oopse} automatically updates ${\bf u}$ when the rotation matrix +${\bf u}$ will be 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 +139,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 @@ -278,10 +278,7 @@ 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 @@ -299,9 +296,9 @@ To carry out isobaric-isothermal ensemble calculations \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{Melchionna1993} +Isobaric-isothermal ensemble integrator 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), \\ @@ -356,11 +353,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: @@ -401,12 +394,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) @@ -448,10 +440,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, @@ -473,10 +462,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)} @@ -550,8 +535,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: @@ -590,13 +575,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( -\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 -T_{\mathrm{target}}}{2} +\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 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 @@ -606,8 +592,6 @@ 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} - \subsubsection{\label{methodSection:NPAT}NPAT Ensemble} A comprehensive understanding of structure¨Cfunction relations of @@ -631,7 +615,8 @@ described for the NPTi integrator. Note that the iterative schemes for NPAT are identical to those described for the NPTi integrator. -\subsubsection{\label{methodSection:NPrT}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 @@ -669,20 +654,92 @@ $\gamma$ is set to zero. integrator is a special case of $NP\gamma T$ if the surface tension $\gamma$ is set to zero. -%\section{\label{methodSection:constraintMethod}Constraint Method} +\section{\label{methodSection:zcons}Z-Constraint Method} -%\subsection{\label{methodSection:bondConstraint}Bond Constraint for Rigid Body} - -%\subsection{\label{methodSection:zcons}Z-constraint Method} - -\section{\label{methodSection:langevin}Integrators for Langevin Dynamics of Rigid Bodies} - -\subsection{\label{methodSection:temperature}Temperature Control} - -\subsection{\label{methodSection:pressureControl}Pressure 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} -\section{\label{methodSection:hydrodynamics}Hydrodynamics} +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:coarseGrained}Coarse-Grained Modeling} +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:moleculeScale}Molecular-Scale Modeling} +After the force calculation, 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} + +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}