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 2854 by tim, Sun Jun 11 02:53:15 2006 UTC vs.
Revision 2905 by tim, Wed Jun 28 19:09:25 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 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
32 < proposed by Dullweber, Leimkuhler and McLachlan (DLM) also preserved
33 < the same structural properties of the Hamiltonian flow. In this
34 < section, the integration scheme of DLM method will be reviewed and
35 < extended to other ensembles.
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 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 Euler angle methods because of
23 > the time-reversible nature. Extending the Trotter-Suzuki
24 > factorization to general system with a flat phase space, Miller and
25 > his colleagues devised a novel symplectic, time-reversible and
26 > volume-preserving integrator in the quaternion representation, which
27 > was shown to be superior to the Matubayasi's time-reversible
28 > integrator. However, all of the integrators in the quaternion
29 > representation suffer from the computational penalty of constructing
30 > a rotation matrix from quaternions to evolve coordinates and
31 > velocities at every time step. An alternative integration scheme
32 > utilizing the rotation matrix directly proposed by Dullweber,
33 > Leimkuhler and McLachlan (DLM) also preserved the same structural
34 > properties of the Hamiltonian flow. In this section, the integration
35 > scheme of DLM method will be reviewed and extended to other
36 > ensembles.
37  
38   \subsection{\label{methodSection:DLM}Integrating the Equations of Motion: the
39   DLM method}
# Line 47 | Line 48 | for timesteps of length $h$.
48   \item the error for a single time step is of order $\mathcal{O}\left(h^4\right)$
49   for timesteps of length $h$.
50   \end{enumerate}
50 –
51   The integration of the equations of motion is carried out in a
52   velocity-Verlet style 2-part algorithm, where $h= \delta t$:
53  
# Line 65 | Line 65 | velocity-Verlet style 2-part algorithm, where $h= \del
65   \mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
66      (t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right).
67   \end{align*}
68 –
68   In this context, the $\mathrm{rotate}$ function is the reversible
69   product of the three body-fixed rotations,
70   \begin{equation}
# Line 100 | Line 99 | All other rotations follow in a straightforward manner
99   \end{array}
100   \right).
101   \end{equation}
102 < All other rotations follow in a straightforward manner.
102 > All other rotations follow in a straightforward manner. After the
103 > first part of the propagation, the forces and body-fixed torques are
104 > calculated at the new positions and orientations
105  
105 – After the first part of the propagation, the forces and body-fixed
106 – torques are calculated at the new positions and orientations
107 –
106   {\tt doForces:}
107   \begin{align*}
108   {\bf f}(t + h) &\leftarrow
109      - \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)}, \\
110   %
111   {\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h)
112 <    \times \frac{\partial V}{\partial {\bf u}}, \\
112 >    \times (\frac{\partial V}{\partial {\bf u}})_{u(t+h)}, \\
113   %
114   {\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h)
115      \cdot {\bf \tau}^s(t + h).
116   \end{align*}
117 <
120 < ${\bf u}$ will be automatically updated when the rotation matrix
117 > ${\bf u}$ is automatically updated when the rotation matrix
118   $\mathsf{A}$ is calculated in {\tt moveA}.  Once the forces and
119   torques have been obtained at the new time step, the velocities can
120   be advanced to the same time value.
# Line 132 | Line 129 | be advanced to the same time value.
129   \right)
130      + \frac{h}{2} {\bf \tau}^b(t + h) .
131   \end{align*}
135 –
132   The matrix rotations used in the DLM method end up being more costly
133   computationally than the simpler arithmetic quaternion propagation.
134   With the same time step, a 1000-molecule water simulation shows an
135   average 7\% increase in computation time using the DLM method in
136   place of quaternions. This cost is more than justified when
137   comparing the energy conservation of the two methods as illustrated
138 < in Fig.~\ref{methodFig:timestep}.
138 > in Fig.~\ref{methodFig:timestep} where the resulting energy drift at
139 > various time steps for both the DLM and quaternion integration
140 > schemes is compared. All of the 1000 molecule water simulations
141 > started with the same configuration, and the only difference was the
142 > method for handling rotational motion. At time steps of 0.1 and 0.5
143 > fs, both methods for propagating molecule rotation conserve energy
144 > fairly well, with the quaternion method showing a slight energy
145 > drift over time in the 0.5 fs time step simulation. At time steps of
146 > 1 and 2 fs, the energy conservation benefits of the DLM method are
147 > clearly demonstrated. Thus, while maintaining the same degree of
148 > energy conservation, one can take considerably longer time steps,
149 > leading to an overall reduction in computation time.
150  
151   \begin{figure}
152   \centering
# Line 154 | Line 161 | In Fig.~\ref{methodFig:timestep}, the resulting energy
161   \label{methodFig:timestep}
162   \end{figure}
163  
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 –
164   \subsection{\label{methodSection:NVT}Nos\'{e}-Hoover Thermostatting}
165  
166   The Nos\'e-Hoover equations of motion are given by\cite{Hoover1985}
# Line 180 | Line 174 | rot}\left(\mathsf{A}^{T} \cdot \frac{\partial V}{\part
174   rot}\left(\mathsf{A}^{T} \cdot \frac{\partial V}{\partial
175   \mathsf{A}} \right) - \chi {\bf j}. \label{eq:nosehoovereom}
176   \end{eqnarray}
183 –
177   $\chi$ is an ``extra'' variable included in the extended system, and
178   it is propagated using the first order equation of motion
179   \begin{equation}
180   \dot{\chi} = \frac{1}{\tau_{T}^2} \left(
181   \frac{T}{T_{\mathrm{target}}} - 1 \right). \label{eq:nosehooverext}
182   \end{equation}
190 –
183   The instantaneous temperature $T$ is proportional to the total
184   kinetic energy (both translational and orientational) and is given
185   by
186   \begin{equation}
187 < T = \frac{2 K}{f k_B}
187 > T = \frac{2 K}{f k_B}.
188   \end{equation}
189   Here, $f$ is the total number of degrees of freedom in the system,
190   \begin{equation}
191   f = 3 N + 3 N_{\mathrm{orient}} - N_{\mathrm{constraints}},
192   \end{equation}
193 < and $K$ is the total kinetic energy,
193 > where $N_{\mathrm{orient}}$ is the number of molecules with
194 > orientational degrees of freedom, and $K$ is the total kinetic
195 > energy,
196   \begin{equation}
197   K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i +
198   \sum_{i=1}^{N_{\mathrm{orient}}}  \frac{1}{2} {\bf j}_i^T \cdot
199   \overleftrightarrow{\mathsf{I}}_i^{-1} \cdot {\bf j}_i.
200   \end{equation}
201 <
202 < In eq.(\ref{eq:nosehooverext}), $\tau_T$ is the time constant for
203 < relaxation of the temperature to the target value.  To set values
204 < 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:
201 > In Eq.~\ref{eq:nosehooverext}, $\tau_T$ is the time constant for
202 > relaxation of the temperature to the target value. The integration
203 > of the equations of motion is carried out in a velocity-Verlet style
204 > 2 part algorithm:
205  
206   {\tt moveA:}
207   \begin{align*}
# Line 237 | Line 226 | T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\b
226      + \frac{h}{2 \tau_T^2} \left( \frac{T(t)}
227      {T_{\mathrm{target}}} - 1 \right) .
228   \end{align*}
240 –
229   Here $\mathrm{rotate}(h * {\bf j}
230   \overleftrightarrow{\mathsf{I}}^{-1})$ is the same symplectic
231   Trotter factorization of the three rotation operations that was
# Line 248 | Line 236 | integrator.
236   step.  The new positions (and orientations) are then used to
237   calculate a new set of forces and torques in exactly the same way
238   they are calculated in the {\tt doForces} portion of the DLM
239 < integrator.
239 > integrator. Once the forces and torques have been obtained at the
240 > new time step, the temperature, velocities, and the extended system
241 > variable can be advanced to the same time value.
242  
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
255 – advanced to the same time value.
256 –
243   {\tt moveB:}
244   \begin{align*}
245   T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
# Line 273 | Line 259 | T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
259      \left( {\bf \tau}^b(t + h) - {\bf j}(t + h)
260      \chi(t + h) \right) .
261   \end{align*}
276 –
262   Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to
263   caclculate $T(t + h)$ as well as $\chi(t + h)$, they indirectly
264   depend on their own values at time $t + h$.  {\tt moveB} is
265   therefore done in an iterative fashion until $\chi(t + h)$ becomes
266 < self-consistent.
267 <
268 < The Nos\'e-Hoover algorithm is known to conserve a Hamiltonian for
284 < the extended system that is, to within a constant, identical to the
285 < Helmholtz free energy,\cite{Melchionna1993}
266 > self-consistent. The Nos\'e-Hoover algorithm is known to conserve a
267 > Hamiltonian for the extended system that is, to within a constant,
268 > identical to the Helmholtz free energy,\cite{Melchionna1993}
269   \begin{equation}
270   H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left(
271   \frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime)
# Line 294 | Line 277 | isotropic box deformations (NPTi)}
277   of the integration.
278  
279   \subsection{\label{methodSection:NPTi}Constant-pressure integration with
280 < isotropic box deformations (NPTi)}
280 > isotropic box (NPTi)}
281  
282 < Isobaric-isothermal ensemble integrator is implemented using the
283 < Melchionna modifications to the Nos\'e-Hoover-Andersen equations of
284 < motion,\cite{Melchionna1993}
302 <
282 > We can used an isobaric-isothermal ensemble integrator which is
283 > implemented using the Melchionna modifications to the
284 > Nos\'e-Hoover-Andersen equations of motion\cite{Melchionna1993}
285   \begin{eqnarray}
286   \dot{{\bf r}} & = & {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right), \\
287   \dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\eta + \chi) {\bf v}, \\
# Line 316 | Line 298 | P_{\mathrm{target}} \right), \\
298   P_{\mathrm{target}} \right), \\
299   \dot{\mathcal{V}} & = & 3 \mathcal{V} \eta . \label{eq:melchionna1}
300   \end{eqnarray}
319 –
301   $\chi$ and $\eta$ are the ``extra'' degrees of freedom in the
302   extended system.  $\chi$ is a thermostat, and it has the same
303   function as it does in the Nos\'e-Hoover NVT integrator.  $\eta$ is
# Line 349 | Line 330 | P(t) = \frac{1}{3} \mathrm{Tr} \left(
330   the Pressure tensor,
331   \begin{equation}
332   P(t) = \frac{1}{3} \mathrm{Tr} \left(
333 < \overleftrightarrow{\mathsf{P}}(t). \right)
333 > \overleftrightarrow{\mathsf{P}}(t) \right) .
334   \end{equation}
335 <
355 < In eq.(\ref{eq:melchionna1}), $\tau_B$ is the time constant for
335 > In Eq.~\ref{eq:melchionna1}, $\tau_B$ is the time constant for
336   relaxation of the pressure to the target value. Like in the NVT
337   integrator, the integration of the equations of motion is carried
338   out in a velocity-Verlet style 2 part algorithm:
# Line 391 | Line 371 | P(t) &\leftarrow \left\{{\bf r}(t)\right\}, \left\{{\b
371   \mathsf{H}(t + h) &\leftarrow e^{-h \eta(t + h / 2)}
372      \mathsf{H}(t).
373   \end{align*}
394 –
374   Most of these equations are identical to their counterparts in the
375   NVT integrator, but the propagation of positions to time $t + h$
376   depends on the positions at the same time. The simulation box
# Line 403 | Line 382 | box by
382   \mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)}.
383   \mathcal{V}(t)
384   \end{equation}
406 –
385   The {\tt doForces} step for the NPTi integrator is exactly the same
386   as in both the DLM and NVT integrators.  Once the forces and torques
387   have been obtained at the new time step, the velocities can be
# Line 435 | Line 413 | P(t + h) &\leftarrow  \left\{{\bf r}(t + h)\right\},
413      \tau}^b(t + h) - {\bf j}(t + h)
414      \chi(t + h) \right) .
415   \end{align*}
438 –
416   Once again, since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required
417   to caclculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, and $\eta(t +
418   h)$, they indirectly depend on their own values at time $t + h$.
# Line 536 | Line 513 | exponential operation is used to scale the simulation
513      \overleftrightarrow{\eta}(t + h / 2)} .
514   \end{align*}
515   Here, a power series expansion truncated at second order for the
516 < exponential operation is used to scale the simulation box.
516 > exponential operation is used to scale the simulation box. The {\tt
517 > moveB} portion of the algorithm is largely unchanged from the NPTi
518 > integrator:
519  
541 – The {\tt moveB} portion of the algorithm is largely unchanged from
542 – the NPTi integrator:
543 –
520   {\tt moveB:}
521   \begin{align*}
522   T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
# Line 570 | Line 546 | T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
546      + h / 2 \right) + \frac{h}{2} \left( {\bf \tau}^b(t
547      + h) - {\bf j}(t + h) \chi(t + h) \right) .
548   \end{align*}
573 –
549   The iterative schemes for both {\tt moveA} and {\tt moveB} are
550 < identical to those described for the NPTi integrator.
551 <
577 < The NPTf integrator is known to conserve the following Hamiltonian:
550 > identical to those described for the NPTi integrator. The NPTf
551 > integrator is known to conserve the following Hamiltonian:
552   \begin{eqnarray*}
553   H_{\mathrm{NPTf}} & = & V + K + f k_B T_{\mathrm{target}} \left(
554   \frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime)
555   dt^\prime \right) \\
556 < + P_{\mathrm{target}} \mathcal{V}(t) + \frac{f
557 < k_B T_{\mathrm{target}}}{2}
556 > & & + P_{\mathrm{target}} \mathcal{V}(t) + \frac{f k_B
557 > T_{\mathrm{target}}}{2}
558   \mathrm{Tr}\left[\overleftrightarrow{\eta}(t)\right]^2 \tau_B^2.
559   \end{eqnarray*}
586 –
560   This integrator must be used with care, particularly in liquid
561   simulations.  Liquids have very small restoring forces in the
562   off-diagonal directions, and the simulation box can very quickly
# Line 592 | Line 565 | assume non-orthorhombic geometries.
565   finds most use in simulating crystals or liquid crystals which
566   assume non-orthorhombic geometries.
567  
568 < \subsubsection{\label{methodSection:NPAT}NPAT Ensemble}
568 > \subsection{\label{methodSection:NPAT}NPAT Ensemble}
569  
570 < A comprehensive understanding of structure¨Cfunction relations of
571 < biological membrane system ultimately relies on structure and
572 < dynamics of lipid bilayer, which are strongly affected by the
573 < interfacial interaction between lipid molecules and surrounding
574 < media. One quantity to describe the interfacial interaction is so
575 < called the average surface area per lipid. Constat area and constant
576 < lateral pressure simulation can be achieved by extending the
577 < standard NPT ensemble with a different pressure control strategy
570 > A comprehensive understanding of relations between structures and
571 > functions in biological membrane system ultimately relies on
572 > structure and dynamics of lipid bilayers, which are strongly
573 > affected by the interfacial interaction between lipid molecules and
574 > surrounding media. One quantity to describe the interfacial
575 > interaction is so called the average surface area per lipid.
576 > Constant area and constant lateral pressure simulations can be
577 > achieved by extending the standard NPT ensemble with a different
578 > pressure control strategy
579  
580   \begin{equation}
581   \dot {\overleftrightarrow{\eta}} _{\alpha \beta}=\left\{\begin{array}{ll}
# Line 620 | Line 594 | minimum with respect to surface area $A$
594  
595   Theoretically, the surface tension $\gamma$ of a stress free
596   membrane system should be zero since its surface free energy $G$ is
597 < minimum with respect to surface area $A$
598 < \[
599 < \gamma  = \frac{{\partial G}}{{\partial A}}.
600 < \]
601 < However, a surface tension of zero is not appropriate for relatively
602 < small patches of membrane. In order to eliminate the edge effect of
603 < the membrane simulation, a special ensemble, NP$\gamma$T, is
630 < proposed to maintain the lateral surface tension and normal
631 < pressure. The equation of motion for cell size control tensor,
632 < $\eta$, in $NP\gamma T$ is
597 > minimum with respect to surface area $A$, $\gamma  = \frac{{\partial
598 > G}}{{\partial A}}.$ However, a surface tension of zero is not
599 > appropriate for relatively small patches of membrane. In order to
600 > eliminate the edge effect of the membrane simulation, a special
601 > ensemble, NP$\gamma$T, has been proposed to maintain the lateral
602 > surface tension and normal pressure. The equation of motion for the
603 > cell size control tensor, $\eta$, in $NP\gamma T$ is
604   \begin{equation}
605   \dot {\overleftrightarrow{\eta}} _{\alpha \beta}=\left\{\begin{array}{ll}
606      - A_{xy} (\gamma _\alpha   - \gamma _{{\rm{target}}} ) & \mbox{$\alpha  = \beta  = x$ or $\alpha  = \beta  = y$}\\
# Line 645 | Line 616 | - P_{{\rm{target}}} )
616   - P_{{\rm{target}}} )
617   \label{methodEquation:instantaneousSurfaceTensor}
618   \end{equation}
648 –
619   There is one additional extended system integrator (NPTxyz), in
620   which each attempt to preserve the target pressure along the box
621   walls perpendicular to that particular axis.  The lengths of the box
622   axes are allowed to fluctuate independently, but the angle between
623   the box axes does not change. It should be noted that the NPTxyz
624   integrator is a special case of $NP\gamma T$ if the surface tension
625 < $\gamma$ is set to zero.
625 > $\gamma$ is set to zero, and if $x$ and $y$ can move independently.
626  
627 < \section{\label{methodSection:zcons}Z-Constraint Method}
627 > \section{\label{methodSection:zcons}The Z-Constraint Method}
628  
629   Based on the fluctuation-dissipation theorem, a force
630   auto-correlation method was developed by Roux and Karplus to
# Line 668 | Line 638 | where%
638   \begin{equation}
639   \delta F(z,t)=F(z,t)-\langle F(z,t)\rangle.
640   \end{equation}
671 –
641   If the time-dependent friction decays rapidly, the static friction
642   coefficient can be approximated by
643   \begin{equation}
# Line 681 | Line 650 | D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B
650   D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty
651   }\langle\delta F(z,t)\delta F(z,0)\rangle dt}.%
652   \end{equation}
684 –
653   The Z-Constraint method, which fixes the z coordinates of the
654   molecules with respect to the center of the mass of the system, has
655   been a method suggested to obtain the forces required for the force
# Line 690 | Line 658 | each time step instead of resetting the coordinate.
658   whole system. To avoid this problem, we reset the forces of
659   z-constrained molecules as well as subtract the total constraint
660   forces from the rest of the system after the force calculation at
661 < each time step instead of resetting the coordinate.
662 <
695 < After the force calculation, define $G_\alpha$ as
661 > each time step instead of resetting the coordinate. After the force
662 > calculation, we define $G_\alpha$ as
663   \begin{equation}
664   G_{\alpha} = \sum_i F_{\alpha i}, \label{oopseEq:zc1}
665   \end{equation}
# Line 727 | Line 694 | v_{\beta i} = v_{\beta i} + \sum_{\alpha}
694   v_{\beta i} = v_{\beta i} + \sum_{\alpha}
695      \frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}.
696   \end{equation}
730 –
697   At the very beginning of the simulation, the molecules may not be at
698   their constrained positions. To move a z-constrained molecule to its
699   specified position, a simple harmonic potential is used

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines