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

Comparing trunk/tengDissertation/Introduction.tex (file contents):
Revision 2685 by tim, Mon Apr 3 18:07:54 2006 UTC vs.
Revision 2700 by tim, Tue Apr 11 03:38:09 2006 UTC

# Line 1 | Line 1
1   \chapter{\label{chapt:introduction}INTRODUCTION AND THEORETICAL BACKGROUND}
2  
3 < \section{\label{introSection:classicalMechanics}Classical Mechanics}
3 > \section{\label{introSection:classicalMechanics}Classical
4 > Mechanics}
5  
6 < \section{\label{introSection:rigidBody}Dynamics of Rigid Bodies}
6 > Closely related to Classical Mechanics, Molecular Dynamics
7 > simulations are carried out by integrating the equations of motion
8 > for a given system of particles. There are three fundamental ideas
9 > behind classical mechanics. Firstly, One can determine the state of
10 > a mechanical system at any time of interest; Secondly, all the
11 > mechanical properties of the system at that time can be determined
12 > by combining the knowledge of the properties of the system with the
13 > specification of this state; Finally, the specification of the state
14 > when further combine with the laws of mechanics will also be
15 > sufficient to predict the future behavior of the system.
16  
17 < \section{\label{introSection:statisticalMechanics}Statistical Mechanics}
17 > \subsection{\label{introSection:newtonian}Newtonian Mechanics}
18 > The discovery of Newton's three laws of mechanics which govern the
19 > motion of particles is the foundation of the classical mechanics.
20 > Newton¡¯s first law defines a class of inertial frames. Inertial
21 > frames are reference frames where a particle not interacting with
22 > other bodies will move with constant speed in the same direction.
23 > With respect to inertial frames Newton¡¯s second law has the form
24 > \begin{equation}
25 > F = \frac {dp}{dt} = \frac {mv}{dt}
26 > \label{introEquation:newtonSecondLaw}
27 > \end{equation}
28 > A point mass interacting with other bodies moves with the
29 > acceleration along the direction of the force acting on it. Let
30 > $F_ij$ be the force that particle $i$ exerts on particle $j$, and
31 > $F_ji$ be the force that particle $j$ exerts on particle $i$.
32 > Newton¡¯s third law states that
33 > \begin{equation}
34 > F_ij = -F_ji
35 > \label{introEquation:newtonThirdLaw}
36 > \end{equation}
37  
38 + Conservation laws of Newtonian Mechanics play very important roles
39 + in solving mechanics problems. The linear momentum of a particle is
40 + conserved if it is free or it experiences no force. The second
41 + conservation theorem concerns the angular momentum of a particle.
42 + The angular momentum $L$ of a particle with respect to an origin
43 + from which $r$ is measured is defined to be
44 + \begin{equation}
45 + L \equiv r \times p \label{introEquation:angularMomentumDefinition}
46 + \end{equation}
47 + The torque $\tau$ with respect to the same origin is defined to be
48 + \begin{equation}
49 + N \equiv r \times F \label{introEquation:torqueDefinition}
50 + \end{equation}
51 + Differentiating Eq.~\ref{introEquation:angularMomentumDefinition},
52 + \[
53 + \dot L = \frac{d}{{dt}}(r \times p) = (\dot r \times p) + (r \times
54 + \dot p)
55 + \]
56 + since
57 + \[
58 + \dot r \times p = \dot r \times mv = m\dot r \times \dot r \equiv 0
59 + \]
60 + thus,
61 + \begin{equation}
62 + \dot L = r \times \dot p = N
63 + \end{equation}
64 + If there are no external torques acting on a body, the angular
65 + momentum of it is conserved. The last conservation theorem state
66 + that if all forces are conservative, Energy
67 + \begin{equation}E = T + V \label{introEquation:energyConservation}
68 + \end{equation}
69 + is conserved. All of these conserved quantities are
70 + important factors to determine the quality of numerical integration
71 + scheme for rigid body \cite{Dullweber1997}.
72 +
73 + \subsection{\label{introSection:lagrangian}Lagrangian Mechanics}
74 +
75 + Newtonian Mechanics suffers from two important limitations: it
76 + describes their motion in special cartesian coordinate systems.
77 + Another limitation of Newtonian mechanics becomes obvious when we
78 + try to describe systems with large numbers of particles. It becomes
79 + very difficult to predict the properties of the system by carrying
80 + out calculations involving the each individual interaction between
81 + all the particles, even if we know all of the details of the
82 + interaction. In order to overcome some of the practical difficulties
83 + which arise in attempts to apply Newton's equation to complex
84 + system, alternative procedures may be developed.
85 +
86 + \subsubsection{\label{introSection:halmiltonPrinciple}Hamilton's
87 + Principle}
88 +
89 + Hamilton introduced the dynamical principle upon which it is
90 + possible to base all of mechanics and, indeed, most of classical
91 + physics. Hamilton's Principle may be stated as follow,
92 +
93 + The actual trajectory, along which a dynamical system may move from
94 + one point to another within a specified time, is derived by finding
95 + the path which minimizes the time integral of the difference between
96 + the kinetic, $K$, and potential energies, $U$ \cite{tolman79}.
97 + \begin{equation}
98 + \delta \int_{t_1 }^{t_2 } {(K - U)dt = 0} ,
99 + \label{introEquation:halmitonianPrinciple1}
100 + \end{equation}
101 +
102 + For simple mechanical systems, where the forces acting on the
103 + different part are derivable from a potential and the velocities are
104 + small compared with that of light, the Lagrangian function $L$ can
105 + be define as the difference between the kinetic energy of the system
106 + and its potential energy,
107 + \begin{equation}
108 + L \equiv K - U = L(q_i ,\dot q_i ) ,
109 + \label{introEquation:lagrangianDef}
110 + \end{equation}
111 + then Eq.~\ref{introEquation:halmitonianPrinciple1} becomes
112 + \begin{equation}
113 + \delta \int_{t_1 }^{t_2 } {L dt = 0} ,
114 + \label{introEquation:halmitonianPrinciple2}
115 + \end{equation}
116 +
117 + \subsubsection{\label{introSection:equationOfMotionLagrangian}The
118 + Equations of Motion in Lagrangian Mechanics}
119 +
120 + For a holonomic system of $f$ degrees of freedom, the equations of
121 + motion in the Lagrangian form is
122 + \begin{equation}
123 + \frac{d}{{dt}}\frac{{\partial L}}{{\partial \dot q_i }} -
124 + \frac{{\partial L}}{{\partial q_i }} = 0,{\rm{ }}i = 1, \ldots,f
125 + \label{introEquation:eqMotionLagrangian}
126 + \end{equation}
127 + where $q_{i}$ is generalized coordinate and $\dot{q_{i}}$ is
128 + generalized velocity.
129 +
130 + \subsection{\label{introSection:hamiltonian}Hamiltonian Mechanics}
131 +
132 + Arising from Lagrangian Mechanics, Hamiltonian Mechanics was
133 + introduced by William Rowan Hamilton in 1833 as a re-formulation of
134 + classical mechanics. If the potential energy of a system is
135 + independent of generalized velocities, the generalized momenta can
136 + be defined as
137 + \begin{equation}
138 + p_i = \frac{\partial L}{\partial \dot q_i}
139 + \label{introEquation:generalizedMomenta}
140 + \end{equation}
141 + The Lagrange equations of motion are then expressed by
142 + \begin{equation}
143 + p_i  = \frac{{\partial L}}{{\partial q_i }}
144 + \label{introEquation:generalizedMomentaDot}
145 + \end{equation}
146 +
147 + With the help of the generalized momenta, we may now define a new
148 + quantity $H$ by the equation
149 + \begin{equation}
150 + H = \sum\limits_k {p_k \dot q_k }  - L ,
151 + \label{introEquation:hamiltonianDefByLagrangian}
152 + \end{equation}
153 + where $ \dot q_1  \ldots \dot q_f $ are generalized velocities and
154 + $L$ is the Lagrangian function for the system.
155 +
156 + Differentiating Eq.~\ref{introEquation:hamiltonianDefByLagrangian},
157 + one can obtain
158 + \begin{equation}
159 + dH = \sum\limits_k {\left( {p_k d\dot q_k  + \dot q_k dp_k  -
160 + \frac{{\partial L}}{{\partial q_k }}dq_k  - \frac{{\partial
161 + L}}{{\partial \dot q_k }}d\dot q_k } \right)}  - \frac{{\partial
162 + L}}{{\partial t}}dt \label{introEquation:diffHamiltonian1}
163 + \end{equation}
164 + Making use of  Eq.~\ref{introEquation:generalizedMomenta}, the
165 + second and fourth terms in the parentheses cancel. Therefore,
166 + Eq.~\ref{introEquation:diffHamiltonian1} can be rewritten as
167 + \begin{equation}
168 + dH = \sum\limits_k {\left( {\dot q_k dp_k  - \dot p_k dq_k }
169 + \right)}  - \frac{{\partial L}}{{\partial t}}dt
170 + \label{introEquation:diffHamiltonian2}
171 + \end{equation}
172 + By identifying the coefficients of $dq_k$, $dp_k$ and dt, we can
173 + find
174 + \begin{equation}
175 + \frac{{\partial H}}{{\partial p_k }} = q_k
176 + \label{introEquation:motionHamiltonianCoordinate}
177 + \end{equation}
178 + \begin{equation}
179 + \frac{{\partial H}}{{\partial q_k }} =  - p_k
180 + \label{introEquation:motionHamiltonianMomentum}
181 + \end{equation}
182 + and
183 + \begin{equation}
184 + \frac{{\partial H}}{{\partial t}} =  - \frac{{\partial L}}{{\partial
185 + t}}
186 + \label{introEquation:motionHamiltonianTime}
187 + \end{equation}
188 +
189 + Eq.~\ref{introEquation:motionHamiltonianCoordinate} and
190 + Eq.~\ref{introEquation:motionHamiltonianMomentum} are Hamilton's
191 + equation of motion. Due to their symmetrical formula, they are also
192 + known as the canonical equations of motions \cite{Goldstein01}.
193 +
194 + An important difference between Lagrangian approach and the
195 + Hamiltonian approach is that the Lagrangian is considered to be a
196 + function of the generalized velocities $\dot q_i$ and the
197 + generalized coordinates $q_i$, while the Hamiltonian is considered
198 + to be a function of the generalized momenta $p_i$ and the conjugate
199 + generalized coordinate $q_i$. Hamiltonian Mechanics is more
200 + appropriate for application to statistical mechanics and quantum
201 + mechanics, since it treats the coordinate and its time derivative as
202 + independent variables and it only works with 1st-order differential
203 + equations\cite{Marion90}.
204 +
205 + In Newtonian Mechanics, a system described by conservative forces
206 + conserves the total energy \ref{introEquation:energyConservation}.
207 + It follows that Hamilton's equations of motion conserve the total
208 + Hamiltonian.
209 + \begin{equation}
210 + \frac{{dH}}{{dt}} = \sum\limits_i {\left( {\frac{{\partial
211 + H}}{{\partial q_i }}\dot q_i  + \frac{{\partial H}}{{\partial p_i
212 + }}\dot p_i } \right)}  = \sum\limits_i {\left( {\frac{{\partial
213 + H}}{{\partial q_i }}\frac{{\partial H}}{{\partial p_i }} -
214 + \frac{{\partial H}}{{\partial p_i }}\frac{{\partial H}}{{\partial
215 + q_i }}} \right) = 0} \label{introEquation:conserveHalmitonian}
216 + \end{equation}
217 +
218 + \section{\label{introSection:statisticalMechanics}Statistical
219 + Mechanics}
220 +
221 + The thermodynamic behaviors and properties of Molecular Dynamics
222 + simulation are governed by the principle of Statistical Mechanics.
223 + The following section will give a brief introduction to some of the
224 + Statistical Mechanics concepts and theorem presented in this
225 + dissertation.
226 +
227 + \subsection{\label{introSection:ensemble}Phase Space and Ensemble}
228 +
229 + Mathematically, phase space is the space which represents all
230 + possible states. Each possible state of the system corresponds to
231 + one unique point in the phase space. For mechanical systems, the
232 + phase space usually consists of all possible values of position and
233 + momentum variables. Consider a dynamic system in a cartesian space,
234 + where each of the $6f$ coordinates and momenta is assigned to one of
235 + $6f$ mutually orthogonal axes, the phase space of this system is a
236 + $6f$ dimensional space. A point, $x = (q_1 , \ldots ,q_f ,p_1 ,
237 + \ldots ,p_f )$, with a unique set of values of $6f$ coordinates and
238 + momenta is a phase space vector.
239 +
240 + A microscopic state or microstate of a classical system is
241 + specification of the complete phase space vector of a system at any
242 + instant in time. An ensemble is defined as a collection of systems
243 + sharing one or more macroscopic characteristics but each being in a
244 + unique microstate. The complete ensemble is specified by giving all
245 + systems or microstates consistent with the common macroscopic
246 + characteristics of the ensemble. Although the state of each
247 + individual system in the ensemble could be precisely described at
248 + any instance in time by a suitable phase space vector, when using
249 + ensembles for statistical purposes, there is no need to maintain
250 + distinctions between individual systems, since the numbers of
251 + systems at any time in the different states which correspond to
252 + different regions of the phase space are more interesting. Moreover,
253 + in the point of view of statistical mechanics, one would prefer to
254 + use ensembles containing a large enough population of separate
255 + members so that the numbers of systems in such different states can
256 + be regarded as changing continuously as we traverse different
257 + regions of the phase space. The condition of an ensemble at any time
258 + can be regarded as appropriately specified by the density $\rho$
259 + with which representative points are distributed over the phase
260 + space. The density of distribution for an ensemble with $f$ degrees
261 + of freedom is defined as,
262 + \begin{equation}
263 + \rho  = \rho (q_1 , \ldots ,q_f ,p_1 , \ldots ,p_f ,t).
264 + \label{introEquation:densityDistribution}
265 + \end{equation}
266 + Governed by the principles of mechanics, the phase points change
267 + their value which would change the density at any time at phase
268 + space. Hence, the density of distribution is also to be taken as a
269 + function of the time.
270 +
271 + The number of systems $\delta N$ at time $t$ can be determined by,
272 + \begin{equation}
273 + \delta N = \rho (q,p,t)dq_1  \ldots dq_f dp_1  \ldots dp_f.
274 + \label{introEquation:deltaN}
275 + \end{equation}
276 + Assuming a large enough population of systems are exploited, we can
277 + sufficiently approximate $\delta N$ without introducing
278 + discontinuity when we go from one region in the phase space to
279 + another. By integrating over the whole phase space,
280 + \begin{equation}
281 + N = \int { \ldots \int {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f
282 + \label{introEquation:totalNumberSystem}
283 + \end{equation}
284 + gives us an expression for the total number of the systems. Hence,
285 + the probability per unit in the phase space can be obtained by,
286 + \begin{equation}
287 + \frac{{\rho (q,p,t)}}{N} = \frac{{\rho (q,p,t)}}{{\int { \ldots \int
288 + {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}.
289 + \label{introEquation:unitProbability}
290 + \end{equation}
291 + With the help of Equation(\ref{introEquation:unitProbability}) and
292 + the knowledge of the system, it is possible to calculate the average
293 + value of any desired quantity which depends on the coordinates and
294 + momenta of the system. Even when the dynamics of the real system is
295 + complex, or stochastic, or even discontinuous, the average
296 + properties of the ensemble of possibilities as a whole may still
297 + remain well defined. For a classical system in thermal equilibrium
298 + with its environment, the ensemble average of a mechanical quantity,
299 + $\langle A(q , p) \rangle_t$, takes the form of an integral over the
300 + phase space of the system,
301 + \begin{equation}
302 + \langle  A(q , p) \rangle_t = \frac{{\int { \ldots \int {A(q,p)\rho
303 + (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}{{\int { \ldots \int {\rho
304 + (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}
305 + \label{introEquation:ensembelAverage}
306 + \end{equation}
307 +
308 + There are several different types of ensembles with different
309 + statistical characteristics. As a function of macroscopic
310 + parameters, such as temperature \textit{etc}, partition function can
311 + be used to describe the statistical properties of a system in
312 + thermodynamic equilibrium.
313 +
314 + As an ensemble of systems, each of which is known to be thermally
315 + isolated and conserve energy, Microcanonical ensemble(NVE) has a
316 + partition function like,
317 + \begin{equation}
318 + \Omega (N,V,E) = e^{\beta TS}
319 + \label{introEqaution:NVEPartition}.
320 + \end{equation}
321 + A canonical ensemble(NVT)is an ensemble of systems, each of which
322 + can share its energy with a large heat reservoir. The distribution
323 + of the total energy amongst the possible dynamical states is given
324 + by the partition function,
325 + \begin{equation}
326 + \Omega (N,V,T) = e^{ - \beta A}
327 + \label{introEquation:NVTPartition}
328 + \end{equation}
329 + Here, $A$ is the Helmholtz free energy which is defined as $ A = U -
330 + TS$. Since most experiment are carried out under constant pressure
331 + condition, isothermal-isobaric ensemble(NPT) play a very important
332 + role in molecular simulation. The isothermal-isobaric ensemble allow
333 + the system to exchange energy with a heat bath of temperature $T$
334 + and to change the volume as well. Its partition function is given as
335 + \begin{equation}
336 + \Delta (N,P,T) =  - e^{\beta G}.
337 + \label{introEquation:NPTPartition}
338 + \end{equation}
339 + Here, $G = U - TS + PV$, and $G$ is called Gibbs free energy.
340 +
341 + \subsection{\label{introSection:liouville}Liouville's theorem}
342 +
343 + The Liouville's theorem is the foundation on which statistical
344 + mechanics rests. It describes the time evolution of phase space
345 + distribution function. In order to calculate the rate of change of
346 + $\rho$, we begin from Equation(\ref{introEquation:deltaN}). If we
347 + consider the two faces perpendicular to the $q_1$ axis, which are
348 + located at $q_1$ and $q_1 + \delta q_1$, the number of phase points
349 + leaving the opposite face is given by the expression,
350 + \begin{equation}
351 + \left( {\rho  + \frac{{\partial \rho }}{{\partial q_1 }}\delta q_1 }
352 + \right)\left( {\dot q_1  + \frac{{\partial \dot q_1 }}{{\partial q_1
353 + }}\delta q_1 } \right)\delta q_2  \ldots \delta q_f \delta p_1
354 + \ldots \delta p_f .
355 + \end{equation}
356 + Summing all over the phase space, we obtain
357 + \begin{equation}
358 + \frac{{d(\delta N)}}{{dt}} =  - \sum\limits_{i = 1}^f {\left[ {\rho
359 + \left( {\frac{{\partial \dot q_i }}{{\partial q_i }} +
360 + \frac{{\partial \dot p_i }}{{\partial p_i }}} \right) + \left(
361 + {\frac{{\partial \rho }}{{\partial q_i }}\dot q_i  + \frac{{\partial
362 + \rho }}{{\partial p_i }}\dot p_i } \right)} \right]} \delta q_1
363 + \ldots \delta q_f \delta p_1  \ldots \delta p_f .
364 + \end{equation}
365 + Differentiating the equations of motion in Hamiltonian formalism
366 + (\ref{introEquation:motionHamiltonianCoordinate},
367 + \ref{introEquation:motionHamiltonianMomentum}), we can show,
368 + \begin{equation}
369 + \sum\limits_i {\left( {\frac{{\partial \dot q_i }}{{\partial q_i }}
370 + + \frac{{\partial \dot p_i }}{{\partial p_i }}} \right)}  = 0 ,
371 + \end{equation}
372 + which cancels the first terms of the right hand side. Furthermore,
373 + divining $ \delta q_1  \ldots \delta q_f \delta p_1  \ldots \delta
374 + p_f $ in both sides, we can write out Liouville's theorem in a
375 + simple form,
376 + \begin{equation}
377 + \frac{{\partial \rho }}{{\partial t}} + \sum\limits_{i = 1}^f
378 + {\left( {\frac{{\partial \rho }}{{\partial q_i }}\dot q_i  +
379 + \frac{{\partial \rho }}{{\partial p_i }}\dot p_i } \right)}  = 0 .
380 + \label{introEquation:liouvilleTheorem}
381 + \end{equation}
382 +
383 + Liouville's theorem states that the distribution function is
384 + constant along any trajectory in phase space. In classical
385 + statistical mechanics, since the number of particles in the system
386 + is huge, we may be able to believe the system is stationary,
387 + \begin{equation}
388 + \frac{{\partial \rho }}{{\partial t}} = 0.
389 + \label{introEquation:stationary}
390 + \end{equation}
391 + In such stationary system, the density of distribution $\rho$ can be
392 + connected to the Hamiltonian $H$ through Maxwell-Boltzmann
393 + distribution,
394 + \begin{equation}
395 + \rho  \propto e^{ - \beta H}
396 + \label{introEquation:densityAndHamiltonian}
397 + \end{equation}
398 +
399 + Liouville's theorem can be expresses in a variety of different forms
400 + which are convenient within different contexts. For any two function
401 + $F$ and $G$ of the coordinates and momenta of a system, the Poisson
402 + bracket ${F, G}$ is defined as
403 + \begin{equation}
404 + \left\{ {F,G} \right\} = \sum\limits_i {\left( {\frac{{\partial
405 + F}}{{\partial q_i }}\frac{{\partial G}}{{\partial p_i }} -
406 + \frac{{\partial F}}{{\partial p_i }}\frac{{\partial G}}{{\partial
407 + q_i }}} \right)}.
408 + \label{introEquation:poissonBracket}
409 + \end{equation}
410 + Substituting equations of motion in Hamiltonian formalism(
411 + \ref{introEquation:motionHamiltonianCoordinate} ,
412 + \ref{introEquation:motionHamiltonianMomentum} ) into
413 + (\ref{introEquation:liouvilleTheorem}), we can rewrite Liouville's
414 + theorem using Poisson bracket notion,
415 + \begin{equation}
416 + \left( {\frac{{\partial \rho }}{{\partial t}}} \right) =  - \left\{
417 + {\rho ,H} \right\}.
418 + \label{introEquation:liouvilleTheromInPoissin}
419 + \end{equation}
420 + Moreover, the Liouville operator is defined as
421 + \begin{equation}
422 + iL = \sum\limits_{i = 1}^f {\left( {\frac{{\partial H}}{{\partial
423 + p_i }}\frac{\partial }{{\partial q_i }} - \frac{{\partial
424 + H}}{{\partial q_i }}\frac{\partial }{{\partial p_i }}} \right)}
425 + \label{introEquation:liouvilleOperator}
426 + \end{equation}
427 + In terms of Liouville operator, Liouville's equation can also be
428 + expressed as
429 + \begin{equation}
430 + \left( {\frac{{\partial \rho }}{{\partial t}}} \right) =  - iL\rho
431 + \label{introEquation:liouvilleTheoremInOperator}
432 + \end{equation}
433 +
434 +
435 + \subsection{\label{introSection:ergodic}The Ergodic Hypothesis}
436 +
437 + Various thermodynamic properties can be calculated from Molecular
438 + Dynamics simulation. By comparing experimental values with the
439 + calculated properties, one can determine the accuracy of the
440 + simulation and the quality of the underlying model. However, both of
441 + experiment and computer simulation are usually performed during a
442 + certain time interval and the measurements are averaged over a
443 + period of them which is different from the average behavior of
444 + many-body system in Statistical Mechanics. Fortunately, Ergodic
445 + Hypothesis is proposed to make a connection between time average and
446 + ensemble average. It states that time average and average over the
447 + statistical ensemble are identical \cite{Frenkel1996, leach01:mm}.
448 + \begin{equation}
449 + \langle A(q , p) \rangle_t = \mathop {\lim }\limits_{t \to \infty }
450 + \frac{1}{t}\int\limits_0^t {A(q(t),p(t))dt = \int\limits_\Gamma
451 + {A(q(t),p(t))} } \rho (q(t), p(t)) dqdp
452 + \end{equation}
453 + where $\langle  A(q , p) \rangle_t$ is an equilibrium value of a
454 + physical quantity and $\rho (p(t), q(t))$ is the equilibrium
455 + distribution function. If an observation is averaged over a
456 + sufficiently long time (longer than relaxation time), all accessible
457 + microstates in phase space are assumed to be equally probed, giving
458 + a properly weighted statistical average. This allows the researcher
459 + freedom of choice when deciding how best to measure a given
460 + observable. In case an ensemble averaged approach sounds most
461 + reasonable, the Monte Carlo techniques\cite{metropolis:1949} can be
462 + utilized. Or if the system lends itself to a time averaging
463 + approach, the Molecular Dynamics techniques in
464 + Sec.~\ref{introSection:molecularDynamics} will be the best
465 + choice\cite{Frenkel1996}.
466 +
467 + \section{\label{introSection:geometricIntegratos}Geometric Integrators}
468 + A variety of numerical integrators were proposed to simulate the
469 + motions. They usually begin with an initial conditionals and move
470 + the objects in the direction governed by the differential equations.
471 + However, most of them ignore the hidden physical law contained
472 + within the equations. Since 1990, geometric integrators, which
473 + preserve various phase-flow invariants such as symplectic structure,
474 + volume and time reversal symmetry, are developed to address this
475 + issue. The velocity verlet method, which happens to be a simple
476 + example of symplectic integrator, continues to gain its popularity
477 + in molecular dynamics community. This fact can be partly explained
478 + by its geometric nature.
479 +
480 + \subsection{\label{introSection:symplecticManifold}Symplectic Manifold}
481 + A \emph{manifold} is an abstract mathematical space. It locally
482 + looks like Euclidean space, but when viewed globally, it may have
483 + more complicate structure. A good example of manifold is the surface
484 + of Earth. It seems to be flat locally, but it is round if viewed as
485 + a whole. A \emph{differentiable manifold} (also known as
486 + \emph{smooth manifold}) is a manifold with an open cover in which
487 + the covering neighborhoods are all smoothly isomorphic to one
488 + another. In other words,it is possible to apply calculus on
489 + \emph{differentiable manifold}. A \emph{symplectic manifold} is
490 + defined as a pair $(M, \omega)$ which consisting of a
491 + \emph{differentiable manifold} $M$ and a close, non-degenerated,
492 + bilinear symplectic form, $\omega$. A symplectic form on a vector
493 + space $V$ is a function $\omega(x, y)$ which satisfies
494 + $\omega(\lambda_1x_1+\lambda_2x_2, y) = \lambda_1\omega(x_1, y)+
495 + \lambda_2\omega(x_2, y)$, $\omega(x, y) = - \omega(y, x)$ and
496 + $\omega(x, x) = 0$. Cross product operation in vector field is an
497 + example of symplectic form.
498 +
499 + One of the motivations to study \emph{symplectic manifold} in
500 + Hamiltonian Mechanics is that a symplectic manifold can represent
501 + all possible configurations of the system and the phase space of the
502 + system can be described by it's cotangent bundle. Every symplectic
503 + manifold is even dimensional. For instance, in Hamilton equations,
504 + coordinate and momentum always appear in pairs.
505 +
506 + Let  $(M,\omega)$ and $(N, \eta)$ be symplectic manifolds. A map
507 + \[
508 + f : M \rightarrow N
509 + \]
510 + is a \emph{symplectomorphism} if it is a \emph{diffeomorphims} and
511 + the \emph{pullback} of $\eta$ under f is equal to $\omega$.
512 + Canonical transformation is an example of symplectomorphism in
513 + classical mechanics.
514 +
515 + \subsection{\label{introSection:ODE}Ordinary Differential Equations}
516 +
517 + For a ordinary differential system defined as
518 + \begin{equation}
519 + \dot x = f(x)
520 + \end{equation}
521 + where $x = x(q,p)^T$, this system is canonical Hamiltonian, if
522 + \begin{equation}
523 + f(r) = J\nabla _x H(r).
524 + \end{equation}
525 + $H = H (q, p)$ is Hamiltonian function and $J$ is the skew-symmetric
526 + matrix
527 + \begin{equation}
528 + J = \left( {\begin{array}{*{20}c}
529 +   0 & I  \\
530 +   { - I} & 0  \\
531 + \end{array}} \right)
532 + \label{introEquation:canonicalMatrix}
533 + \end{equation}
534 + where $I$ is an identity matrix. Using this notation, Hamiltonian
535 + system can be rewritten as,
536 + \begin{equation}
537 + \frac{d}{{dt}}x = J\nabla _x H(x)
538 + \label{introEquation:compactHamiltonian}
539 + \end{equation}In this case, $f$ is
540 + called a \emph{Hamiltonian vector field}.
541 +
542 + Another generalization of Hamiltonian dynamics is Poisson Dynamics,
543 + \begin{equation}
544 + \dot x = J(x)\nabla _x H \label{introEquation:poissonHamiltonian}
545 + \end{equation}
546 + The most obvious change being that matrix $J$ now depends on $x$.
547 + The free rigid body is an example of Poisson system (actually a
548 + Lie-Poisson system) with Hamiltonian function of angular kinetic
549 + energy.
550 + \begin{equation}
551 + J(\pi ) = \left( {\begin{array}{*{20}c}
552 +   0 & {\pi _3 } & { - \pi _2 }  \\
553 +   { - \pi _3 } & 0 & {\pi _1 }  \\
554 +   {\pi _2 } & { - \pi _1 } & 0  \\
555 + \end{array}} \right)
556 + \end{equation}
557 +
558 + \begin{equation}
559 + H = \frac{1}{2}\left( {\frac{{\pi _1^2 }}{{I_1 }} + \frac{{\pi _2^2
560 + }}{{I_2 }} + \frac{{\pi _3^2 }}{{I_3 }}} \right)
561 + \end{equation}
562 +
563 + \subsection{\label{introSection:geometricProperties}Geometric Properties}
564 + Let $x(t)$ be the exact solution of the ODE system,
565 + \begin{equation}
566 + \frac{{dx}}{{dt}} = f(x) \label{introEquation:ODE}
567 + \end{equation}
568 + The exact flow(solution) $\varphi_\tau$ is defined by
569 + \[
570 + x(t+\tau) =\varphi_\tau(x(t))
571 + \]
572 + where $\tau$ is a fixed time step and $\varphi$ is a map from phase
573 + space to itself. In most cases, it is not easy to find the exact
574 + flow $\varphi_\tau$. Instead, we use a approximate map, $\psi_\tau$,
575 + which is usually called integrator. The order of an integrator
576 + $\psi_\tau$ is $p$, if the Taylor series of $\psi_\tau$ agree to
577 + order $p$,
578 + \begin{equation}
579 + \psi_tau(x) = x + \tau f(x) + O(\tau^{p+1})
580 + \end{equation}
581 +
582 + The hidden geometric properties of ODE and its flow play important
583 + roles in numerical studies. Let $\varphi$ be the flow of Hamiltonian
584 + vector field, $\varphi$ is a \emph{symplectic} flow if it satisfies,
585 + \begin{equation}
586 + '\varphi^T J '\varphi = J.
587 + \end{equation}
588 + According to Liouville's theorem, the symplectic volume is invariant
589 + under a Hamiltonian flow, which is the basis for classical
590 + statistical mechanics. Furthermore, the flow of a Hamiltonian vector
591 + field on a symplectic manifold can be shown to be a
592 + symplectomorphism. As to the Poisson system,
593 + \begin{equation}
594 + '\varphi ^T J '\varphi  = J \circ \varphi
595 + \end{equation}
596 + is the property must be preserved by the integrator. It is possible
597 + to construct a \emph{volume-preserving} flow for a source free($
598 + \nabla \cdot f = 0 $) ODE, if the flow satisfies $ \det d\varphi  =
599 + 1$. Changing the variables $y = h(x)$ in a
600 + ODE\ref{introEquation:ODE} will result in a new system,
601 + \[
602 + \dot y = \tilde f(y) = ((dh \cdot f)h^{ - 1} )(y).
603 + \]
604 + The vector filed $f$ has reversing symmetry $h$ if $f = - \tilde f$.
605 + In other words, the flow of this vector field is reversible if and
606 + only if $ h \circ \varphi ^{ - 1}  = \varphi  \circ h $. When
607 + designing any numerical methods, one should always try to preserve
608 + the structural properties of the original ODE and its flow.
609 +
610 + \subsection{\label{introSection:constructionSymplectic}Construction of Symplectic Methods}
611 + A lot of well established and very effective numerical methods have
612 + been successful precisely because of their symplecticities even
613 + though this fact was not recognized when they were first
614 + constructed. The most famous example is leapfrog methods in
615 + molecular dynamics. In general, symplectic integrators can be
616 + constructed using one of four different methods.
617 + \begin{enumerate}
618 + \item Generating functions
619 + \item Variational methods
620 + \item Runge-Kutta methods
621 + \item Splitting methods
622 + \end{enumerate}
623 +
624 + Generating function tends to lead to methods which are cumbersome
625 + and difficult to use\cite{}. In dissipative systems, variational
626 + methods can capture the decay of energy accurately\cite{}. Since
627 + their geometrically unstable nature against non-Hamiltonian
628 + perturbations, ordinary implicit Runge-Kutta methods are not
629 + suitable for Hamiltonian system. Recently, various high-order
630 + explicit Runge--Kutta methods have been developed to overcome this
631 + instability \cite{}. However, due to computational penalty involved
632 + in implementing the Runge-Kutta methods, they do not attract too
633 + much attention from Molecular Dynamics community. Instead, splitting
634 + have been widely accepted since they exploit natural decompositions
635 + of the system\cite{Tuckerman92}. The main idea behind splitting
636 + methods is to decompose the discrete $\varphi_h$ as a composition of
637 + simpler flows,
638 + \begin{equation}
639 + \varphi _h  = \varphi _{h_1 }  \circ \varphi _{h_2 }  \ldots  \circ
640 + \varphi _{h_n }
641 + \label{introEquation:FlowDecomposition}
642 + \end{equation}
643 + where each of the sub-flow is chosen such that each represent a
644 + simpler integration of the system. Let $\phi$ and $\psi$ both be
645 + symplectic maps, it is easy to show that any composition of
646 + symplectic flows yields a symplectic map,
647 + \begin{equation}
648 + (\varphi '\phi ')^T J\varphi '\phi ' = \phi '^T \varphi '^T J\varphi
649 + '\phi ' = \phi '^T J\phi ' = J.
650 + \label{introEquation:SymplecticFlowComposition}
651 + \end{equation}
652 + Suppose that a Hamiltonian system has a form with $H = T + V$
653 +
654   \section{\label{introSection:molecularDynamics}Molecular Dynamics}
655  
656 + As a special discipline of molecular modeling, Molecular dynamics
657 + has proven to be a powerful tool for studying the functions of
658 + biological systems, providing structural, thermodynamic and
659 + dynamical information.
660 +
661 + \subsection{\label{introSec:mdInit}Initialization}
662 +
663 + \subsection{\label{introSection:mdIntegration} Integration of the Equations of Motion}
664 +
665 + \section{\label{introSection:rigidBody}Dynamics of Rigid Bodies}
666 +
667 + A rigid body is a body in which the distance between any two given
668 + points of a rigid body remains constant regardless of external
669 + forces exerted on it. A rigid body therefore conserves its shape
670 + during its motion.
671 +
672 + Applications of dynamics of rigid bodies.
673 +
674 + \subsection{\label{introSection:lieAlgebra}Lie Algebra}
675 +
676 + \subsection{\label{introSection:DLMMotionEquation}The Euler Equations of Rigid Body Motion}
677 +
678 + \subsection{\label{introSection:otherRBMotionEquation}Other Formulations for Rigid Body Motion}
679 +
680 + %\subsection{\label{introSection:poissonBrackets}Poisson Brackets}
681 +
682 + \section{\label{introSection:correlationFunctions}Correlation Functions}
683 +
684   \section{\label{introSection:langevinDynamics}Langevin Dynamics}
685  
686 + \subsection{\label{introSection:LDIntroduction}Introduction and application of Langevin Dynamics}
687 +
688 + \subsection{\label{introSection:generalizedLangevinDynamics}Generalized Langevin Dynamics}
689 +
690 + \begin{equation}
691 + H = \frac{{p^2 }}{{2m}} + U(x) + H_B  + \Delta U(x,x_1 , \ldots x_N)
692 + \label{introEquation:bathGLE}
693 + \end{equation}
694 + where $H_B$ is harmonic bath Hamiltonian,
695 + \[
696 + H_B  =\sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
697 + }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  w_\alpha ^2 } \right\}}
698 + \]
699 + and $\Delta U$ is bilinear system-bath coupling,
700 + \[
701 + \Delta U =  - \sum\limits_{\alpha  = 1}^N {g_\alpha  x_\alpha  x}
702 + \]
703 + Completing the square,
704 + \[
705 + H_B  + \Delta U = \sum\limits_{\alpha  = 1}^N {\left\{
706 + {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
707 + w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
708 + w_\alpha ^2 }}x} \right)^2 } \right\}}  - \sum\limits_{\alpha  =
709 + 1}^N {\frac{{g_\alpha ^2 }}{{2m_\alpha  w_\alpha ^2 }}} x^2
710 + \]
711 + and putting it back into Eq.~\ref{introEquation:bathGLE},
712 + \[
713 + H = \frac{{p^2 }}{{2m}} + W(x) + \sum\limits_{\alpha  = 1}^N
714 + {\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
715 + w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
716 + w_\alpha ^2 }}x} \right)^2 } \right\}}
717 + \]
718 + where
719 + \[
720 + W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
721 + }}{{2m_\alpha  w_\alpha ^2 }}} x^2
722 + \]
723 + Since the first two terms of the new Hamiltonian depend only on the
724 + system coordinates, we can get the equations of motion for
725 + Generalized Langevin Dynamics by Hamilton's equations
726 + \ref{introEquation:motionHamiltonianCoordinate,
727 + introEquation:motionHamiltonianMomentum},
728 + \begin{align}
729 + \dot p &=  - \frac{{\partial H}}{{\partial x}}
730 +       &= m\ddot x
731 +       &= - \frac{{\partial W(x)}}{{\partial x}} - \sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right)}
732 + \label{introEq:Lp5}
733 + \end{align}
734 + , and
735 + \begin{align}
736 + \dot p_\alpha   &=  - \frac{{\partial H}}{{\partial x_\alpha  }}
737 +                &= m\ddot x_\alpha
738 +                &= \- m_\alpha  w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha}}{{m_\alpha  w_\alpha ^2 }}x} \right)
739 + \end{align}
740 +
741 + \subsection{\label{introSection:laplaceTransform}The Laplace Transform}
742 +
743 + \[
744 + L(x) = \int_0^\infty  {x(t)e^{ - pt} dt}
745 + \]
746 +
747 + \[
748 + L(x + y) = L(x) + L(y)
749 + \]
750 +
751 + \[
752 + L(ax) = aL(x)
753 + \]
754 +
755 + \[
756 + L(\dot x) = pL(x) - px(0)
757 + \]
758 +
759 + \[
760 + L(\ddot x) = p^2 L(x) - px(0) - \dot x(0)
761 + \]
762 +
763 + \[
764 + L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right) = G(p)H(p)
765 + \]
766 +
767 + Some relatively important transformation,
768 + \[
769 + L(\cos at) = \frac{p}{{p^2  + a^2 }}
770 + \]
771 +
772 + \[
773 + L(\sin at) = \frac{a}{{p^2  + a^2 }}
774 + \]
775 +
776 + \[
777 + L(1) = \frac{1}{p}
778 + \]
779 +
780 + First, the bath coordinates,
781 + \[
782 + p^2 L(x_\alpha  ) - px_\alpha  (0) - \dot x_\alpha  (0) =  - \omega
783 + _\alpha ^2 L(x_\alpha  ) + \frac{{g_\alpha  }}{{\omega _\alpha
784 + }}L(x)
785 + \]
786 + \[
787 + L(x_\alpha  ) = \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) +
788 + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }}
789 + \]
790 + Then, the system coordinates,
791 + \begin{align}
792 + mL(\ddot x) &=  - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
793 + \sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{\frac{{g_\alpha
794 + }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha
795 + (0)}}{{p^2  + \omega _\alpha ^2 }} - \frac{{g_\alpha ^2 }}{{m_\alpha
796 + }}\omega _\alpha ^2 L(x)} \right\}}
797 + %
798 + &= - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
799 + \sum\limits_{\alpha  = 1}^N {\left\{ { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}\frac{p}{{p^2  + \omega _\alpha ^2 }}pL(x)
800 + - \frac{p}{{p^2  + \omega _\alpha ^2 }}g_\alpha  x_\alpha  (0)
801 + - \frac{1}{{p^2  + \omega _\alpha ^2 }}g_\alpha  \dot x_\alpha  (0)} \right\}}
802 + \end{align}
803 + Then, the inverse transform,
804 +
805 + \begin{align}
806 + m\ddot x &=  - \frac{{\partial W(x)}}{{\partial x}} -
807 + \sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
808 + }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\int_0^t {\cos (\omega
809 + _\alpha  t)\dot x(t - \tau )d\tau  - \left[ {g_\alpha  x_\alpha  (0)
810 + - \frac{{g_\alpha  }}{{m_\alpha  \omega _\alpha  }}} \right]\cos
811 + (\omega _\alpha  t) - \frac{{g_\alpha  \dot x_\alpha  (0)}}{{\omega
812 + _\alpha  }}\sin (\omega _\alpha  t)} } \right\}}
813 + %
814 + &= - \frac{{\partial W(x)}}{{\partial x}} - \int_0^t
815 + {\sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
816 + }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha
817 + t)\dot x(t - \tau )d} \tau }  + \sum\limits_{\alpha  = 1}^N {\left\{
818 + {\left[ {g_\alpha  x_\alpha  (0) - \frac{{g_\alpha  }}{{m_\alpha
819 + \omega _\alpha  }}} \right]\cos (\omega _\alpha  t) +
820 + \frac{{g_\alpha  \dot x_\alpha  (0)}}{{\omega _\alpha  }}\sin
821 + (\omega _\alpha  t)} \right\}}
822 + \end{align}
823 +
824 + \begin{equation}
825 + m\ddot x =  - \frac{{\partial W}}{{\partial x}} - \int_0^t {\xi
826 + (t)\dot x(t - \tau )d\tau }  + R(t)
827 + \label{introEuqation:GeneralizedLangevinDynamics}
828 + \end{equation}
829 + %where $ {\xi (t)}$ is friction kernel, $R(t)$ is random force and
830 + %$W$ is the potential of mean force. $W(x) =  - kT\ln p(x)$
831 + \[
832 + \xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
833 + }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
834 + \]
835 + For an infinite harmonic bath, we can use the spectral density and
836 + an integral over frequencies.
837 +
838 + \[
839 + R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
840 + - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
841 + \right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
842 + (0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)
843 + \]
844 + The random forces depend only on initial conditions.
845 +
846 + \subsubsection{\label{introSection:secondFluctuationDissipation}The Second Fluctuation Dissipation Theorem}
847 + So we can define a new set of coordinates,
848 + \[
849 + q_\alpha  (t) = x_\alpha  (t) - \frac{1}{{m_\alpha  \omega _\alpha
850 + ^2 }}x(0)
851 + \]
852 + This makes
853 + \[
854 + R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}
855 + \]
856 + And since the $q$ coordinates are harmonic oscillators,
857 + \[
858 + \begin{array}{l}
859 + \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  = \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
860 + \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle  = \delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
861 + \end{array}
862 + \]
863 +
864 + \begin{align}
865 + \left\langle {R(t)R(0)} \right\rangle  &= \sum\limits_\alpha
866 + {\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha
867 + (t)q_\beta  (0)} \right\rangle } }
868 + %
869 + &= \sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)}
870 + \right\rangle \cos (\omega _\alpha  t)}
871 + %
872 + &= kT\xi (t)
873 + \end{align}
874 +
875 + \begin{equation}
876 + \xi (t) = \left\langle {R(t)R(0)} \right\rangle
877 + \label{introEquation:secondFluctuationDissipation}
878 + \end{equation}
879 +
880   \section{\label{introSection:hydroynamics}Hydrodynamics}
881  
882 < \section{\label{introSection:correlationFunctions}Correlation Functions}
882 > \subsection{\label{introSection:frictionTensor} Friction Tensor}
883 > \subsection{\label{introSection:analyticalApproach}Analytical
884 > Approach}
885 >
886 > \subsection{\label{introSection:approximationApproach}Approximation
887 > Approach}
888 >
889 > \subsection{\label{introSection:centersRigidBody}Centers of Rigid
890 > Body}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines