trunk/tengDissertation/Introduction.tex

\chapter{\label{chapt:introduction}INTRODUCTION AND THEORETICAL BACKGROUND}

\section{\label{introSection:classicalMechanics}Classical
Mechanics}

Closely related to Classical Mechanics, Molecular Dynamics
simulations are carried out by integrating the equations of motion
for a given system of particles. There are three fundamental ideas
behind classical mechanics. Firstly, One can determine the state of
a mechanical system at any time of interest; Secondly, all the
mechanical properties of the system at that time can be determined
by combining the knowledge of the properties of the system with the
specification of this state; Finally, the specification of the state
when further combine with the laws of mechanics will also be
sufficient to predict the future behavior of the system.

\subsection{\label{introSection:newtonian}Newtonian Mechanics}
The discovery of Newton's three laws of mechanics which govern the
motion of particles is the foundation of the classical mechanics.
Newton’s first law defines a class of inertial frames. Inertial
frames are reference frames where a particle not interacting with
other bodies will move with constant speed in the same direction.
With respect to inertial frames Newton’s second law has the form
\begin{equation}
F = \frac {dp}{dt} = \frac {mv}{dt}
\label{introEquation:newtonSecondLaw}
\end{equation}
A point mass interacting with other bodies moves with the
acceleration along the direction of the force acting on it. Let
$F_ij$ be the force that particle $i$ exerts on particle $j$, and
$F_ji$ be the force that particle $j$ exerts on particle $i$.
Newton’s third law states that
\begin{equation}
F_ij = -F_ji
 \label{introEquation:newtonThirdLaw}
\end{equation}

Conservation laws of Newtonian Mechanics play very important roles
in solving mechanics problems. The linear momentum of a particle is
conserved if it is free or it experiences no force. The second
conservation theorem concerns the angular momentum of a particle.
The angular momentum $L$ of a particle with respect to an origin
from which $r$ is measured is defined to be
\begin{equation}
L \equiv r \times p \label{introEquation:angularMomentumDefinition}
\end{equation}
The torque $\tau$ with respect to the same origin is defined to be
\begin{equation}
N \equiv r \times F \label{introEquation:torqueDefinition}
\end{equation}
Differentiating Eq.~\ref{introEquation:angularMomentumDefinition},
\[
\dot L = \frac{d}{{dt}}(r \times p) = (\dot r \times p) + (r \times
\dot p)
\]
since
\[
\dot r \times p = \dot r \times mv = m\dot r \times \dot r \equiv 0
\]
thus,
\begin{equation}
\dot L = r \times \dot p = N
\end{equation}
If there are no external torques acting on a body, the angular
momentum of it is conserved. The last conservation theorem state
that if all forces are conservative, Energy
\begin{equation}E = T + V \label{introEquation:energyConservation}
\end{equation}
 is conserved. All of these conserved quantities are
important factors to determine the quality of numerical integration
scheme for rigid body \cite{Dullweber1997}.

\subsection{\label{introSection:lagrangian}Lagrangian Mechanics}

Newtonian Mechanics suffers from two important limitations: it
describes their motion in special cartesian coordinate systems.
Another limitation of Newtonian mechanics becomes obvious when we
try to describe systems with large numbers of particles. It becomes
very difficult to predict the properties of the system by carrying
out calculations involving the each individual interaction between
all the particles, even if we know all of the details of the
interaction. In order to overcome some of the practical difficulties
which arise in attempts to apply Newton's equation to complex
system, alternative procedures may be developed.

\subsubsection{\label{introSection:halmiltonPrinciple}Hamilton's
Principle}

Hamilton introduced the dynamical principle upon which it is
possible to base all of mechanics and, indeed, most of classical
physics. Hamilton's Principle may be stated as follow,

The actual trajectory, along which a dynamical system may move from
one point to another within a specified time, is derived by finding
the path which minimizes the time integral of the difference between
the kinetic, $K$, and potential energies, $U$ \cite{tolman79}.
\begin{equation}
\delta \int_{t_1 }^{t_2 } {(K - U)dt = 0} ,
\label{introEquation:halmitonianPrinciple1}
\end{equation}

For simple mechanical systems, where the forces acting on the
different part are derivable from a potential and the velocities are
small compared with that of light, the Lagrangian function $L$ can
be define as the difference between the kinetic energy of the system
and its potential energy,
\begin{equation}
L \equiv K - U = L(q_i ,\dot q_i ) ,
\label{introEquation:lagrangianDef}
\end{equation}
then Eq.~\ref{introEquation:halmitonianPrinciple1} becomes
\begin{equation}
\delta \int_{t_1 }^{t_2 } {L dt = 0} ,
\label{introEquation:halmitonianPrinciple2}
\end{equation}

\subsubsection{\label{introSection:equationOfMotionLagrangian}The
Equations of Motion in Lagrangian Mechanics}

For a holonomic system of $f$ degrees of freedom, the equations of
motion in the Lagrangian form is
\begin{equation}
\frac{d}{{dt}}\frac{{\partial L}}{{\partial \dot q_i }} -
\frac{{\partial L}}{{\partial q_i }} = 0,{\rm{ }}i = 1, \ldots,f
\label{introEquation:eqMotionLagrangian}
\end{equation}
where $q_{i}$ is generalized coordinate and $\dot{q_{i}}$ is
generalized velocity.

\subsection{\label{introSection:hamiltonian}Hamiltonian Mechanics}

Arising from Lagrangian Mechanics, Hamiltonian Mechanics was
introduced by William Rowan Hamilton in 1833 as a re-formulation of
classical mechanics. If the potential energy of a system is
independent of generalized velocities, the generalized momenta can
be defined as
\begin{equation}
p_i = \frac{\partial L}{\partial \dot q_i}
\label{introEquation:generalizedMomenta}
\end{equation}
The Lagrange equations of motion are then expressed by
\begin{equation}
p_i  = \frac{{\partial L}}{{\partial q_i }}
\label{introEquation:generalizedMomentaDot}
\end{equation}

With the help of the generalized momenta, we may now define a new
quantity $H$ by the equation
\begin{equation}
H = \sum\limits_k {p_k \dot q_k }  - L ,
\label{introEquation:hamiltonianDefByLagrangian}
\end{equation}
where $ \dot q_1  \ldots \dot q_f $ are generalized velocities and
$L$ is the Lagrangian function for the system.

Differentiating Eq.~\ref{introEquation:hamiltonianDefByLagrangian},
one can obtain
\begin{equation}
dH = \sum\limits_k {\left( {p_k d\dot q_k  + \dot q_k dp_k  -
\frac{{\partial L}}{{\partial q_k }}dq_k  - \frac{{\partial
L}}{{\partial \dot q_k }}d\dot q_k } \right)}  - \frac{{\partial
L}}{{\partial t}}dt \label{introEquation:diffHamiltonian1}
\end{equation}
Making use of  Eq.~\ref{introEquation:generalizedMomenta}, the
second and fourth terms in the parentheses cancel. Therefore,
Eq.~\ref{introEquation:diffHamiltonian1} can be rewritten as
\begin{equation}
dH = \sum\limits_k {\left( {\dot q_k dp_k  - \dot p_k dq_k }
\right)}  - \frac{{\partial L}}{{\partial t}}dt
\label{introEquation:diffHamiltonian2}
\end{equation}
By identifying the coefficients of $dq_k$, $dp_k$ and dt, we can
find
\begin{equation}
\frac{{\partial H}}{{\partial p_k }} = q_k
\label{introEquation:motionHamiltonianCoordinate}
\end{equation}
\begin{equation}
\frac{{\partial H}}{{\partial q_k }} =  - p_k
\label{introEquation:motionHamiltonianMomentum}
\end{equation}
and
\begin{equation}
\frac{{\partial H}}{{\partial t}} =  - \frac{{\partial L}}{{\partial
t}}
\label{introEquation:motionHamiltonianTime}
\end{equation}

Eq.~\ref{introEquation:motionHamiltonianCoordinate} and
Eq.~\ref{introEquation:motionHamiltonianMomentum} are Hamilton's
equation of motion. Due to their symmetrical formula, they are also
known as the canonical equations of motions \cite{Goldstein01}.

An important difference between Lagrangian approach and the
Hamiltonian approach is that the Lagrangian is considered to be a
function of the generalized velocities $\dot q_i$ and the
generalized coordinates $q_i$, while the Hamiltonian is considered
to be a function of the generalized momenta $p_i$ and the conjugate
generalized coordinate $q_i$. Hamiltonian Mechanics is more
appropriate for application to statistical mechanics and quantum
mechanics, since it treats the coordinate and its time derivative as
independent variables and it only works with 1st-order differential
equations\cite{Marion90}.

In Newtonian Mechanics, a system described by conservative forces
conserves the total energy \ref{introEquation:energyConservation}.
It follows that Hamilton's equations of motion conserve the total
Hamiltonian.
\begin{equation}
\frac{{dH}}{{dt}} = \sum\limits_i {\left( {\frac{{\partial
H}}{{\partial q_i }}\dot q_i  + \frac{{\partial H}}{{\partial p_i
}}\dot p_i } \right)}  = \sum\limits_i {\left( {\frac{{\partial
H}}{{\partial q_i }}\frac{{\partial H}}{{\partial p_i }} -
\frac{{\partial H}}{{\partial p_i }}\frac{{\partial H}}{{\partial
q_i }}} \right) = 0} \label{introEquation:conserveHalmitonian}
\end{equation}

\section{\label{introSection:statisticalMechanics}Statistical
Mechanics}

The thermodynamic behaviors and properties of Molecular Dynamics
simulation are governed by the principle of Statistical Mechanics.
The following section will give a brief introduction to some of the
Statistical Mechanics concepts and theorem presented in this
dissertation.

\subsection{\label{introSection:ensemble}Phase Space and Ensemble}

Mathematically, phase space is the space which represents all
possible states. Each possible state of the system corresponds to
one unique point in the phase space. For mechanical systems, the
phase space usually consists of all possible values of position and
momentum variables. Consider a dynamic system in a cartesian space,
where each of the $6f$ coordinates and momenta is assigned to one of
$6f$ mutually orthogonal axes, the phase space of this system is a
$6f$ dimensional space. A point, $x = (q_1 , \ldots ,q_f ,p_1 ,
\ldots ,p_f )$, with a unique set of values of $6f$ coordinates and
momenta is a phase space vector.

A microscopic state or microstate of a classical system is
specification of the complete phase space vector of a system at any
instant in time. An ensemble is defined as a collection of systems
sharing one or more macroscopic characteristics but each being in a
unique microstate. The complete ensemble is specified by giving all
systems or microstates consistent with the common macroscopic
characteristics of the ensemble. Although the state of each
individual system in the ensemble could be precisely described at
any instance in time by a suitable phase space vector, when using
ensembles for statistical purposes, there is no need to maintain
distinctions between individual systems, since the numbers of
systems at any time in the different states which correspond to
different regions of the phase space are more interesting. Moreover,
in the point of view of statistical mechanics, one would prefer to
use ensembles containing a large enough population of separate
members so that the numbers of systems in such different states can
be regarded as changing continuously as we traverse different
regions of the phase space. The condition of an ensemble at any time
can be regarded as appropriately specified by the density $\rho$
with which representative points are distributed over the phase
space. The density of distribution for an ensemble with $f$ degrees
of freedom is defined as,
\begin{equation}
\rho  = \rho (q_1 , \ldots ,q_f ,p_1 , \ldots ,p_f ,t).
\label{introEquation:densityDistribution}
\end{equation}
Governed by the principles of mechanics, the phase points change
their value which would change the density at any time at phase
space. Hence, the density of distribution is also to be taken as a
function of the time.

The number of systems $\delta N$ at time $t$ can be determined by,
\begin{equation}
\delta N = \rho (q,p,t)dq_1  \ldots dq_f dp_1  \ldots dp_f.
\label{introEquation:deltaN}
\end{equation}
Assuming a large enough population of systems are exploited, we can
sufficiently approximate $\delta N$ without introducing
discontinuity when we go from one region in the phase space to
another. By integrating over the whole phase space,
\begin{equation}
N = \int { \ldots \int {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f
\label{introEquation:totalNumberSystem}
\end{equation}
gives us an expression for the total number of the systems. Hence,
the probability per unit in the phase space can be obtained by,
\begin{equation}
\frac{{\rho (q,p,t)}}{N} = \frac{{\rho (q,p,t)}}{{\int { \ldots \int
{\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}.
\label{introEquation:unitProbability}
\end{equation}
With the help of Equation(\ref{introEquation:unitProbability}) and
the knowledge of the system, it is possible to calculate the average
value of any desired quantity which depends on the coordinates and
momenta of the system. Even when the dynamics of the real system is
complex, or stochastic, or even discontinuous, the average
properties of the ensemble of possibilities as a whole may still
remain well defined. For a classical system in thermal equilibrium
with its environment, the ensemble average of a mechanical quantity,
$\langle A(q , p) \rangle_t$, takes the form of an integral over the
phase space of the system,
\begin{equation}
\langle  A(q , p) \rangle_t = \frac{{\int { \ldots \int {A(q,p)\rho
(q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}{{\int { \ldots \int {\rho
(q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}
\label{introEquation:ensembelAverage}
\end{equation}

There are several different types of ensembles with different
statistical characteristics. As a function of macroscopic
parameters, such as temperature \textit{etc}, partition function can
be used to describe the statistical properties of a system in
thermodynamic equilibrium.

As an ensemble of systems, each of which is known to be thermally
isolated and conserve energy, Microcanonical ensemble(NVE) has a
partition function like,
\begin{equation}
\Omega (N,V,E) = e^{\beta TS}
\label{introEqaution:NVEPartition}.
\end{equation}
A canonical ensemble(NVT)is an ensemble of systems, each of which
can share its energy with a large heat reservoir. The distribution
of the total energy amongst the possible dynamical states is given
by the partition function,
\begin{equation}
\Omega (N,V,T) = e^{ - \beta A}
\label{introEquation:NVTPartition}
\end{equation}
Here, $A$ is the Helmholtz free energy which is defined as $ A = U -
TS$. Since most experiment are carried out under constant pressure
condition, isothermal-isobaric ensemble(NPT) play a very important
role in molecular simulation. The isothermal-isobaric ensemble allow
the system to exchange energy with a heat bath of temperature $T$
and to change the volume as well. Its partition function is given as
\begin{equation}
\Delta (N,P,T) =  - e^{\beta G}.
 \label{introEquation:NPTPartition}
\end{equation}
Here, $G = U - TS + PV$, and $G$ is called Gibbs free energy.

\subsection{\label{introSection:liouville}Liouville's theorem}

The Liouville's theorem is the foundation on which statistical
mechanics rests. It describes the time evolution of phase space
distribution function. In order to calculate the rate of change of
$\rho$, we begin from Equation(\ref{introEquation:deltaN}). If we
consider the two faces perpendicular to the $q_1$ axis, which are
located at $q_1$ and $q_1 + \delta q_1$, the number of phase points
leaving the opposite face is given by the expression,
\begin{equation}
\left( {\rho  + \frac{{\partial \rho }}{{\partial q_1 }}\delta q_1 }
\right)\left( {\dot q_1  + \frac{{\partial \dot q_1 }}{{\partial q_1
}}\delta q_1 } \right)\delta q_2  \ldots \delta q_f \delta p_1
\ldots \delta p_f .
\end{equation}
Summing all over the phase space, we obtain
\begin{equation}
\frac{{d(\delta N)}}{{dt}} =  - \sum\limits_{i = 1}^f {\left[ {\rho
\left( {\frac{{\partial \dot q_i }}{{\partial q_i }} +
\frac{{\partial \dot p_i }}{{\partial p_i }}} \right) + \left(
{\frac{{\partial \rho }}{{\partial q_i }}\dot q_i  + \frac{{\partial
\rho }}{{\partial p_i }}\dot p_i } \right)} \right]} \delta q_1
\ldots \delta q_f \delta p_1  \ldots \delta p_f .
\end{equation}
Differentiating the equations of motion in Hamiltonian formalism
(\ref{introEquation:motionHamiltonianCoordinate},
\ref{introEquation:motionHamiltonianMomentum}), we can show,
\begin{equation}
\sum\limits_i {\left( {\frac{{\partial \dot q_i }}{{\partial q_i }}
+ \frac{{\partial \dot p_i }}{{\partial p_i }}} \right)}  = 0 ,
\end{equation}
which cancels the first terms of the right hand side. Furthermore,
divining $ \delta q_1  \ldots \delta q_f \delta p_1  \ldots \delta
p_f $ in both sides, we can write out Liouville's theorem in a
simple form,
\begin{equation}
\frac{{\partial \rho }}{{\partial t}} + \sum\limits_{i = 1}^f
{\left( {\frac{{\partial \rho }}{{\partial q_i }}\dot q_i  +
\frac{{\partial \rho }}{{\partial p_i }}\dot p_i } \right)}  = 0 .
\label{introEquation:liouvilleTheorem}
\end{equation}

Liouville's theorem states that the distribution function is
constant along any trajectory in phase space. In classical
statistical mechanics, since the number of particles in the system
is huge, we may be able to believe the system is stationary,
\begin{equation}
\frac{{\partial \rho }}{{\partial t}} = 0.
\label{introEquation:stationary}
\end{equation}
In such stationary system, the density of distribution $\rho$ can be
connected to the Hamiltonian $H$ through Maxwell-Boltzmann
distribution,
\begin{equation}
\rho  \propto e^{ - \beta H}
\label{introEquation:densityAndHamiltonian}
\end{equation}

Liouville's theorem can be expresses in a variety of different forms
which are convenient within different contexts. For any two function
$F$ and $G$ of the coordinates and momenta of a system, the Poisson
bracket ${F, G}$ is defined as
\begin{equation}
\left\{ {F,G} \right\} = \sum\limits_i {\left( {\frac{{\partial
F}}{{\partial q_i }}\frac{{\partial G}}{{\partial p_i }} -
\frac{{\partial F}}{{\partial p_i }}\frac{{\partial G}}{{\partial
q_i }}} \right)}.
\label{introEquation:poissonBracket}
\end{equation}
Substituting equations of motion in Hamiltonian formalism(
\ref{introEquation:motionHamiltonianCoordinate} ,
\ref{introEquation:motionHamiltonianMomentum} ) into
(\ref{introEquation:liouvilleTheorem}), we can rewrite Liouville's
theorem using Poisson bracket notion,
\begin{equation}
\left( {\frac{{\partial \rho }}{{\partial t}}} \right) =  - \left\{
{\rho ,H} \right\}.
\label{introEquation:liouvilleTheromInPoissin}
\end{equation}
Moreover, the Liouville operator is defined as
\begin{equation}
iL = \sum\limits_{i = 1}^f {\left( {\frac{{\partial H}}{{\partial
p_i }}\frac{\partial }{{\partial q_i }} - \frac{{\partial
H}}{{\partial q_i }}\frac{\partial }{{\partial p_i }}} \right)}
\label{introEquation:liouvilleOperator}
\end{equation}
In terms of Liouville operator, Liouville's equation can also be
expressed as
\begin{equation}
\left( {\frac{{\partial \rho }}{{\partial t}}} \right) =  - iL\rho
\label{introEquation:liouvilleTheoremInOperator}
\end{equation}


\subsection{\label{introSection:ergodic}The Ergodic Hypothesis}

Various thermodynamic properties can be calculated from Molecular
Dynamics simulation. By comparing experimental values with the
calculated properties, one can determine the accuracy of the
simulation and the quality of the underlying model. However, both of
experiment and computer simulation are usually performed during a
certain time interval and the measurements are averaged over a
period of them which is different from the average behavior of
many-body system in Statistical Mechanics. Fortunately, Ergodic
Hypothesis is proposed to make a connection between time average and
ensemble average. It states that time average and average over the
statistical ensemble are identical \cite{Frenkel1996, leach01:mm}.
\begin{equation}
\langle A(q , p) \rangle_t = \mathop {\lim }\limits_{t \to \infty }
\frac{1}{t}\int\limits_0^t {A(q(t),p(t))dt = \int\limits_\Gamma
{A(q(t),p(t))} } \rho (q(t), p(t)) dqdp
\end{equation}
where $\langle  A(q , p) \rangle_t$ is an equilibrium value of a
physical quantity and $\rho (p(t), q(t))$ is the equilibrium
distribution function. If an observation is averaged over a
sufficiently long time (longer than relaxation time), all accessible
microstates in phase space are assumed to be equally probed, giving
a properly weighted statistical average. This allows the researcher
freedom of choice when deciding how best to measure a given
observable. In case an ensemble averaged approach sounds most
reasonable, the Monte Carlo techniques\cite{metropolis:1949} can be
utilized. Or if the system lends itself to a time averaging
approach, the Molecular Dynamics techniques in
Sec.~\ref{introSection:molecularDynamics} will be the best
choice\cite{Frenkel1996}.

\section{\label{introSection:geometricIntegratos}Geometric Integrators}
A variety of numerical integrators were proposed to simulate the
motions. They usually begin with an initial conditionals and move
the objects in the direction governed by the differential equations.
However, most of them ignore the hidden physical law contained
within the equations. Since 1990, geometric integrators, which
preserve various phase-flow invariants such as symplectic structure,
volume and time reversal symmetry, are developed to address this
issue. The velocity verlet method, which happens to be a simple
example of symplectic integrator, continues to gain its popularity
in molecular dynamics community. This fact can be partly explained
by its geometric nature.

\subsection{\label{introSection:symplecticManifold}Symplectic Manifold}
A \emph{manifold} is an abstract mathematical space. It locally
looks like Euclidean space, but when viewed globally, it may have
more complicate structure. A good example of manifold is the surface
of Earth. It seems to be flat locally, but it is round if viewed as
a whole. A \emph{differentiable manifold} (also known as
\emph{smooth manifold}) is a manifold with an open cover in which
the covering neighborhoods are all smoothly isomorphic to one
another. In other words,it is possible to apply calculus on
\emph{differentiable manifold}. A \emph{symplectic manifold} is
defined as a pair $(M, \omega)$ which consisting of a
\emph{differentiable manifold} $M$ and a close, non-degenerated,
bilinear symplectic form, $\omega$. A symplectic form on a vector
space $V$ is a function $\omega(x, y)$ which satisfies
$\omega(\lambda_1x_1+\lambda_2x_2, y) = \lambda_1\omega(x_1, y)+
\lambda_2\omega(x_2, y)$, $\omega(x, y) = - \omega(y, x)$ and
$\omega(x, x) = 0$. Cross product operation in vector field is an
example of symplectic form.

One of the motivations to study \emph{symplectic manifold} in
Hamiltonian Mechanics is that a symplectic manifold can represent
all possible configurations of the system and the phase space of the
system can be described by it's cotangent bundle. Every symplectic
manifold is even dimensional. For instance, in Hamilton equations,
coordinate and momentum always appear in pairs.

Let  $(M,\omega)$ and $(N, \eta)$ be symplectic manifolds. A map
\[
f : M \rightarrow N
\]
is a \emph{symplectomorphism} if it is a \emph{diffeomorphims} and
the \emph{pullback} of $\eta$ under f is equal to $\omega$.
Canonical transformation is an example of symplectomorphism in
classical mechanics.

\subsection{\label{introSection:ODE}Ordinary Differential Equations}

For a ordinary differential system defined as
\begin{equation}
\dot x = f(x)
\end{equation}
where $x = x(q,p)^T$, this system is canonical Hamiltonian, if
\begin{equation}
f(r) = J\nabla _x H(r).
\end{equation}
$H = H (q, p)$ is Hamiltonian function and $J$ is the skew-symmetric
matrix
\begin{equation}
J = \left( {\begin{array}{*{20}c}
   0 & I  \\
   { - I} & 0  \\
\end{array}} \right)
\label{introEquation:canonicalMatrix}
\end{equation}
where $I$ is an identity matrix. Using this notation, Hamiltonian
system can be rewritten as,
\begin{equation}
\frac{d}{{dt}}x = J\nabla _x H(x)
\label{introEquation:compactHamiltonian}
\end{equation}In this case, $f$ is
called a \emph{Hamiltonian vector field}.

Another generalization of Hamiltonian dynamics is Poisson Dynamics,
\begin{equation}
\dot x = J(x)\nabla _x H \label{introEquation:poissonHamiltonian}
\end{equation}
The most obvious change being that matrix $J$ now depends on $x$.
The free rigid body is an example of Poisson system (actually a
Lie-Poisson system) with Hamiltonian function of angular kinetic
energy.
\begin{equation}
J(\pi ) = \left( {\begin{array}{*{20}c}
   0 & {\pi _3 } & { - \pi _2 }  \\
   { - \pi _3 } & 0 & {\pi _1 }  \\
   {\pi _2 } & { - \pi _1 } & 0  \\
\end{array}} \right)
\end{equation}

\begin{equation}
H = \frac{1}{2}\left( {\frac{{\pi _1^2 }}{{I_1 }} + \frac{{\pi _2^2
}}{{I_2 }} + \frac{{\pi _3^2 }}{{I_3 }}} \right)
\end{equation}

\subsection{\label{introSection:geometricProperties}Geometric Properties}
Let $x(t)$ be the exact solution of the ODE system,
\begin{equation}
\frac{{dx}}{{dt}} = f(x) \label{introEquation:ODE}
\end{equation}
The exact flow(solution) $\varphi_\tau$ is defined by
\[
x(t+\tau) =\varphi_\tau(x(t))
\]
where $\tau$ is a fixed time step and $\varphi$ is a map from phase
space to itself. In most cases, it is not easy to find the exact
flow $\varphi_\tau$. Instead, we use a approximate map, $\psi_\tau$,
which is usually called integrator. The order of an integrator
$\psi_\tau$ is $p$, if the Taylor series of $\psi_\tau$ agree to
order $p$,
\begin{equation}
\psi_tau(x) = x + \tau f(x) + O(\tau^{p+1})
\end{equation}

The hidden geometric properties of ODE and its flow play important
roles in numerical studies. Let $\varphi$ be the flow of Hamiltonian
vector field, $\varphi$ is a \emph{symplectic} flow if it satisfies,
\begin{equation}
'\varphi^T J '\varphi = J.
\end{equation}
According to Liouville's theorem, the symplectic volume is invariant
under a Hamiltonian flow, which is the basis for classical
statistical mechanics. Furthermore, the flow of a Hamiltonian vector
field on a symplectic manifold can be shown to be a
symplectomorphism. As to the Poisson system,
\begin{equation}
'\varphi ^T J '\varphi  = J \circ \varphi
\end{equation}
is the property must be preserved by the integrator. It is possible
to construct a \emph{volume-preserving} flow for a source free($
\nabla \cdot f = 0 $) ODE, if the flow satisfies $ \det d\varphi  =
1$. Changing the variables $y = h(x)$ in a
ODE\ref{introEquation:ODE} will result in a new system,
\[
\dot y = \tilde f(y) = ((dh \cdot f)h^{ - 1} )(y).
\]
The vector filed $f$ has reversing symmetry $h$ if $f = - \tilde f$.
In other words, the flow of this vector field is reversible if and
only if $ h \circ \varphi ^{ - 1}  = \varphi  \circ h $. When
designing any numerical methods, one should always try to preserve
the structural properties of the original ODE and its flow.

\subsection{\label{introSection:constructionSymplectic}Construction of Symplectic Methods}
A lot of well established and very effective numerical methods have
been successful precisely because of their symplecticities even
though this fact was not recognized when they were first
constructed. The most famous example is leapfrog methods in
molecular dynamics. In general, symplectic integrators can be
constructed using one of four different methods.
\begin{enumerate}
\item Generating functions
\item Variational methods
\item Runge-Kutta methods
\item Splitting methods
\end{enumerate}

Generating function tends to lead to methods which are cumbersome
and difficult to use\cite{}. In dissipative systems, variational
methods can capture the decay of energy accurately\cite{}. Since
their geometrically unstable nature against non-Hamiltonian
perturbations, ordinary implicit Runge-Kutta methods are not
suitable for Hamiltonian system. Recently, various high-order
explicit Runge--Kutta methods have been developed to overcome this
instability \cite{}. However, due to computational penalty involved
in implementing the Runge-Kutta methods, they do not attract too
much attention from Molecular Dynamics community. Instead, splitting
have been widely accepted since they exploit natural decompositions
of the system\cite{Tuckerman92}. The main idea behind splitting
methods is to decompose the discrete $\varphi_h$ as a composition of
simpler flows,
\begin{equation}
\varphi _h  = \varphi _{h_1 }  \circ \varphi _{h_2 }  \ldots  \circ
\varphi _{h_n }
\label{introEquation:FlowDecomposition}
\end{equation}
where each of the sub-flow is chosen such that each represent a
simpler integration of the system. Let $\phi$ and $\psi$ both be
symplectic maps, it is easy to show that any composition of
symplectic flows yields a symplectic map,
\begin{equation}
(\varphi '\phi ')^T J\varphi '\phi ' = \phi '^T \varphi '^T J\varphi
'\phi ' = \phi '^T J\phi ' = J.
 \label{introEquation:SymplecticFlowComposition}
\end{equation}
Suppose that a Hamiltonian system has a form with $H = T + V$

\section{\label{introSection:molecularDynamics}Molecular Dynamics}

As a special discipline of molecular modeling, Molecular dynamics
has proven to be a powerful tool for studying the functions of
biological systems, providing structural, thermodynamic and
dynamical information.

\subsection{\label{introSec:mdInit}Initialization}

\subsection{\label{introSection:mdIntegration} Integration of the Equations of Motion}

\section{\label{introSection:rigidBody}Dynamics of Rigid Bodies}

A rigid body is a body in which the distance between any two given
points of a rigid body remains constant regardless of external
forces exerted on it. A rigid body therefore conserves its shape
during its motion.

Applications of dynamics of rigid bodies.

\subsection{\label{introSection:lieAlgebra}Lie Algebra}

\subsection{\label{introSection:DLMMotionEquation}The Euler Equations of Rigid Body Motion}

\subsection{\label{introSection:otherRBMotionEquation}Other Formulations for Rigid Body Motion}

%\subsection{\label{introSection:poissonBrackets}Poisson Brackets}

\section{\label{introSection:correlationFunctions}Correlation Functions}

\section{\label{introSection:langevinDynamics}Langevin Dynamics}

\subsection{\label{introSection:LDIntroduction}Introduction and application of Langevin Dynamics}

\subsection{\label{introSection:generalizedLangevinDynamics}Generalized Langevin Dynamics}

\begin{equation}
H = \frac{{p^2 }}{{2m}} + U(x) + H_B  + \Delta U(x,x_1 , \ldots x_N)
\label{introEquation:bathGLE}
\end{equation}
where $H_B$ is harmonic bath Hamiltonian,
\[
H_B  =\sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
}}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  w_\alpha ^2 } \right\}}
\]
and $\Delta U$ is bilinear system-bath coupling,
\[
\Delta U =  - \sum\limits_{\alpha  = 1}^N {g_\alpha  x_\alpha  x}
\]
Completing the square,
\[
H_B  + \Delta U = \sum\limits_{\alpha  = 1}^N {\left\{
{\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
w_\alpha ^2 }}x} \right)^2 } \right\}}  - \sum\limits_{\alpha  =
1}^N {\frac{{g_\alpha ^2 }}{{2m_\alpha  w_\alpha ^2 }}} x^2
\]
and putting it back into Eq.~\ref{introEquation:bathGLE},
\[
H = \frac{{p^2 }}{{2m}} + W(x) + \sum\limits_{\alpha  = 1}^N
{\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
w_\alpha ^2 }}x} \right)^2 } \right\}}
\]
where
\[
W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
}}{{2m_\alpha  w_\alpha ^2 }}} x^2
\]
Since the first two terms of the new Hamiltonian depend only on the
system coordinates, we can get the equations of motion for
Generalized Langevin Dynamics by Hamilton's equations
\ref{introEquation:motionHamiltonianCoordinate,
introEquation:motionHamiltonianMomentum},
\begin{align}
\dot p &=  - \frac{{\partial H}}{{\partial x}}
       &= m\ddot x
       &= - \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)}
\label{introEq:Lp5}
\end{align}
, and
\begin{align}
\dot p_\alpha   &=  - \frac{{\partial H}}{{\partial x_\alpha  }}
                &= m\ddot x_\alpha
                &= \- m_\alpha  w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha}}{{m_\alpha  w_\alpha ^2 }}x} \right)
\end{align}

\subsection{\label{introSection:laplaceTransform}The Laplace Transform}

\[
L(x) = \int_0^\infty  {x(t)e^{ - pt} dt}
\]

\[
L(x + y) = L(x) + L(y)
\]

\[
L(ax) = aL(x)
\]

\[
L(\dot x) = pL(x) - px(0)
\]

\[
L(\ddot x) = p^2 L(x) - px(0) - \dot x(0)
\]

\[
L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right) = G(p)H(p)
\]

Some relatively important transformation,
\[
L(\cos at) = \frac{p}{{p^2  + a^2 }}
\]

\[
L(\sin at) = \frac{a}{{p^2  + a^2 }}
\]

\[
L(1) = \frac{1}{p}
\]

First, the bath coordinates,
\[
p^2 L(x_\alpha  ) - px_\alpha  (0) - \dot x_\alpha  (0) =  - \omega
_\alpha ^2 L(x_\alpha  ) + \frac{{g_\alpha  }}{{\omega _\alpha
}}L(x)
\]
\[
L(x_\alpha  ) = \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) +
px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }}
\]
Then, the system coordinates,
\begin{align}
mL(\ddot x) &=  - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
\sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{\frac{{g_\alpha
}}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha
(0)}}{{p^2  + \omega _\alpha ^2 }} - \frac{{g_\alpha ^2 }}{{m_\alpha
}}\omega _\alpha ^2 L(x)} \right\}}
%
 &= - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
 \sum\limits_{\alpha  = 1}^N {\left\{ { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}\frac{p}{{p^2  + \omega _\alpha ^2 }}pL(x)
 - \frac{p}{{p^2  + \omega _\alpha ^2 }}g_\alpha  x_\alpha  (0)
 - \frac{1}{{p^2  + \omega _\alpha ^2 }}g_\alpha  \dot x_\alpha  (0)} \right\}}
\end{align}
Then, the inverse transform,

\begin{align}
m\ddot x &=  - \frac{{\partial W(x)}}{{\partial x}} -
\sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\int_0^t {\cos (\omega
_\alpha  t)\dot x(t - \tau )d\tau  - \left[ {g_\alpha  x_\alpha  (0)
- \frac{{g_\alpha  }}{{m_\alpha  \omega _\alpha  }}} \right]\cos
(\omega _\alpha  t) - \frac{{g_\alpha  \dot x_\alpha  (0)}}{{\omega
_\alpha  }}\sin (\omega _\alpha  t)} } \right\}}
%
&= - \frac{{\partial W(x)}}{{\partial x}} - \int_0^t
{\sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha
t)\dot x(t - \tau )d} \tau }  + \sum\limits_{\alpha  = 1}^N {\left\{
{\left[ {g_\alpha  x_\alpha  (0) - \frac{{g_\alpha  }}{{m_\alpha
\omega _\alpha  }}} \right]\cos (\omega _\alpha  t) +
\frac{{g_\alpha  \dot x_\alpha  (0)}}{{\omega _\alpha  }}\sin
(\omega _\alpha  t)} \right\}}
\end{align}

\begin{equation}
m\ddot x =  - \frac{{\partial W}}{{\partial x}} - \int_0^t {\xi
(t)\dot x(t - \tau )d\tau }  + R(t)
\label{introEuqation:GeneralizedLangevinDynamics}
\end{equation}
%where $ {\xi (t)}$ is friction kernel, $R(t)$ is random force and
%$W$ is the potential of mean force. $W(x) =  - kT\ln p(x)$
\[
\xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
\]
For an infinite harmonic bath, we can use the spectral density and
an integral over frequencies.

\[
R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
- \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
\right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)
\]
The random forces depend only on initial conditions.

\subsubsection{\label{introSection:secondFluctuationDissipation}The Second Fluctuation Dissipation Theorem}
So we can define a new set of coordinates,
\[
q_\alpha  (t) = x_\alpha  (t) - \frac{1}{{m_\alpha  \omega _\alpha
^2 }}x(0)
\]
This makes
\[
R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}
\]
And since the $q$ coordinates are harmonic oscillators,
\[
\begin{array}{l}
 \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  = \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
 \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle  = \delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
 \end{array}
\]

\begin{align}
\left\langle {R(t)R(0)} \right\rangle  &= \sum\limits_\alpha
{\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha
(t)q_\beta  (0)} \right\rangle } }
%
&= \sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)}
\right\rangle \cos (\omega _\alpha  t)}
%
&= kT\xi (t)
\end{align}

\begin{equation}
\xi (t) = \left\langle {R(t)R(0)} \right\rangle
\label{introEquation:secondFluctuationDissipation}
\end{equation}

\section{\label{introSection:hydroynamics}Hydrodynamics}

\subsection{\label{introSection:frictionTensor} Friction Tensor}
\subsection{\label{introSection:analyticalApproach}Analytical
Approach}

\subsection{\label{introSection:approximationApproach}Approximation
Approach}

\subsection{\label{introSection:centersRigidBody}Centers of Rigid
Body}
Revision:	2700
Committed:	Tue Apr 11 03:38:09 2006 UTC (19 years, 7 months ago) by tim
Content type:	application/x-tex
File size:	37269 byte(s)
Log Message:	adding Liouville's theorem section and Phase Space and Ensemble section