ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/tengDissertation/Introduction.tex

Comparing trunk/tengDissertation/Introduction.tex (file contents):
Revision 2697 by tim, Fri Apr 7 05:03:54 2006 UTC vs.
Revision 2850 by tim, Sun Jun 11 01:55:48 2006 UTC

#	Line 6 | Line 6 | behind classical mechanics. Firstly, One can determine
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
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
#	Line 17 | Line 17 | Newton¡¯s first law defines a class of inertial frames
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
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
23	>	With respect to inertial frames, Newton's second law has the form
24		\begin{equation}
25	<	F = \frac {dp}{dt} = \frac {mv}{dt}
25	>	F = \frac {dp}{dt} = \frac {mdv}{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
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
34	>	F_{ij} = -F_{ji}
35		\label{introEquation:newtonThirdLaw}
36		\end{equation}
37
#	Line 46 | Line 46 | N \equiv r \times F \label{introEquation:torqueDefinit
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}
49	>	\tau \equiv r \times F \label{introEquation:torqueDefinition}
50		\end{equation}
51		Differentiating Eq.~\ref{introEquation:angularMomentumDefinition},
52		\[
#	Line 59 | Line 59 | thus,
59		\]
60		thus,
61		\begin{equation}
62	<	\dot L = r \times \dot p = N
62	>	\dot L = r \times \dot p = \tau
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
#	Line 68 | Line 68 | scheme for rigid body \cite{Dullweber1997}.
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}.
71	>	schemes for rigid bodies \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.
75	>	Newtonian Mechanics suffers from two important limitations: motions
76	>	can only be described in cartesian coordinate systems. Moreover, It
77	>	become impossible to predict analytically the properties of the
78	>	system even if we know all of the details of the interaction. In
79	>	order to overcome some of the practical difficulties which arise in
80	>	attempts to apply Newton's equation to complex system, approximate
81	>	numerical procedures may be developed.
82
83	<	\subsubsection{\label{introSection:halmiltonPrinciple}Hamilton's
84	<	Principle}
83	>	\subsubsection{\label{introSection:halmiltonPrinciple}\textbf{Hamilton's
84	>	Principle}}
85
86		Hamilton introduced the dynamical principle upon which it is
87	<	possible to base all of mechanics and, indeed, most of classical
88	<	physics. Hamilton's Principle may be stated as follow,
87	>	possible to base all of mechanics and most of classical physics.
88	>	Hamilton's Principle may be stated as follows,
89
90		The actual trajectory, along which a dynamical system may move from
91		one point to another within a specified time, is derived by finding
92		the path which minimizes the time integral of the difference between
93	<	the kinetic, $K$, and potential energies, $U$ \cite{tolman79}.
93	>	the kinetic, $K$, and potential energies, $U$.
94		\begin{equation}
95		\delta \int_{t_1 }^{t_2 } {(K - U)dt = 0} ,
96		\label{introEquation:halmitonianPrinciple1}
97		\end{equation}
98
99		For simple mechanical systems, where the forces acting on the
100	<	different part are derivable from a potential and the velocities are
101	<	small compared with that of light, the Lagrangian function $L$ can
102	<	be define as the difference between the kinetic energy of the system
106	<	and its potential energy,
100	>	different parts are derivable from a potential, the Lagrangian
101	>	function $L$ can be defined as the difference between the kinetic
102	>	energy of the system and its potential energy,
103		\begin{equation}
104		L \equiv K - U = L(q_i ,\dot q_i ) ,
105		\label{introEquation:lagrangianDef}
#	Line 114 | Line 110 | then Eq.~\ref{introEquation:halmitonianPrinciple1} bec
110		\label{introEquation:halmitonianPrinciple2}
111		\end{equation}
112
113	<	\subsubsection{\label{introSection:equationOfMotionLagrangian}The
114	<	Equations of Motion in Lagrangian Mechanics}
113	>	\subsubsection{\label{introSection:equationOfMotionLagrangian}\textbf{The
114	>	Equations of Motion in Lagrangian Mechanics}}
115
116	<	for a holonomic system of $f$ degrees of freedom, the equations of
117	<	motion in the Lagrangian form is
116	>	For a system of $f$ degrees of freedom, the equations of motion in
117	>	the Lagrangian form is
118		\begin{equation}
119		\frac{d}{{dt}}\frac{{\partial L}}{{\partial \dot q_i }} -
120		\frac{{\partial L}}{{\partial q_i }} = 0,{\rm{ }}i = 1, \ldots,f
#	Line 132 | Line 128 | independent of generalized velocities, the generalized
128		Arising from Lagrangian Mechanics, Hamiltonian Mechanics was
129		introduced by William Rowan Hamilton in 1833 as a re-formulation of
130		classical mechanics. If the potential energy of a system is
131	<	independent of generalized velocities, the generalized momenta can
136	<	be defined as
131	>	independent of velocities, the momenta can be defined as
132		\begin{equation}
133		p_i = \frac{\partial L}{\partial \dot q_i}
134		\label{introEquation:generalizedMomenta}
#	Line 172 | Line 167 | find
167		By identifying the coefficients of $dq_k$, $dp_k$ and dt, we can
168		find
169		\begin{equation}
170	<	\frac{{\partial H}}{{\partial p_k }} = q_k
170	>	\frac{{\partial H}}{{\partial p_k }} = \dot {q_k}
171		\label{introEquation:motionHamiltonianCoordinate}
172		\end{equation}
173		\begin{equation}
174	<	\frac{{\partial H}}{{\partial q_k }} =  - p_k
174	>	\frac{{\partial H}}{{\partial q_k }} =  - \dot {p_k}
175		\label{introEquation:motionHamiltonianMomentum}
176		\end{equation}
177		and
#	Line 189 | Line 184 | known as the canonical equations of motions \cite{Gold
184		Eq.~\ref{introEquation:motionHamiltonianCoordinate} and
185		Eq.~\ref{introEquation:motionHamiltonianMomentum} are Hamilton's
186		equation of motion. Due to their symmetrical formula, they are also
187	<	known as the canonical equations of motions \cite{Goldstein01}.
187	>	known as the canonical equations of motions \cite{Goldstein2001}.
188
189		An important difference between Lagrangian approach and the
190		Hamiltonian approach is that the Lagrangian is considered to be a
191	<	function of the generalized velocities $\dot q_i$ and the
192	<	generalized coordinates $q_i$, while the Hamiltonian is considered
193	<	to be a function of the generalized momenta $p_i$ and the conjugate
194	<	generalized coordinate $q_i$. Hamiltonian Mechanics is more
195	<	appropriate for application to statistical mechanics and quantum
196	<	mechanics, since it treats the coordinate and its time derivative as
197	<	independent variables and it only works with 1st-order differential
203	<	equations\cite{Marion90}.
191	>	function of the generalized velocities $\dot q_i$ and coordinates
192	>	$q_i$, while the Hamiltonian is considered to be a function of the
193	>	generalized momenta $p_i$ and the conjugate coordinates $q_i$.
194	>	Hamiltonian Mechanics is more appropriate for application to
195	>	statistical mechanics and quantum mechanics, since it treats the
196	>	coordinate and its time derivative as independent variables and it
197	>	only works with 1st-order differential equations\cite{Marion1990}.
198
199		In Newtonian Mechanics, a system described by conservative forces
200		conserves the total energy \ref{introEquation:energyConservation}.
#	Line 212 | Line 206 | q_i }}} \right) = 0}
206		}}\dot p_i } \right)}  = \sum\limits_i {\left( {\frac{{\partial
207		H}}{{\partial q_i }}\frac{{\partial H}}{{\partial p_i }} -
208		\frac{{\partial H}}{{\partial p_i }}\frac{{\partial H}}{{\partial
209	<	q_i }}} \right) = 0}
216	<	\label{introEquation:conserveHalmitonian}
209	>	q_i }}} \right) = 0} \label{introEquation:conserveHalmitonian}
210		\end{equation}
211
219	–	When studying Hamiltonian system, it is more convenient to use
220	–	notation
221	–	\begin{equation}
222	–	r = r(q,p)^T
223	–	\end{equation}
224	–	and to introduce a $2n \times 2n$ canonical structure matrix $J$,
225	–	\begin{equation}
226	–	J = \left( {\begin{array}{*{20}c}
227	–	0 & I  \\
228	–	{ - I} & 0  \\
229	–	\end{array}} \right)
230	–	\label{introEquation:canonicalMatrix}
231	–	\end{equation}
232	–	where $I$ is a $n \times n$ identity matrix and $J$ is a
233	–	skew-symmetric matrix ($ J^T  =  - J $). Thus, Hamiltonian system
234	–	can be rewritten as,
235	–	\begin{equation}
236	–	\frac{d}{{dt}}r = J\nabla _r H(r)
237	–	\label{introEquation:compactHamiltonian}
238	–	\end{equation}
239	–
212		\section{\label{introSection:statisticalMechanics}Statistical
213		Mechanics}
214
215		The thermodynamic behaviors and properties of Molecular Dynamics
216		simulation are governed by the principle of Statistical Mechanics.
217		The following section will give a brief introduction to some of the
218	<	Statistical Mechanics concepts presented in this dissertation.
218	>	Statistical Mechanics concepts and theorem presented in this
219	>	dissertation.
220
221	<	\subsection{\label{introSection:ensemble}Ensemble and Phase Space}
221	>	\subsection{\label{introSection:ensemble}Phase Space and Ensemble}
222
223	+	Mathematically, phase space is the space which represents all
224	+	possible states. Each possible state of the system corresponds to
225	+	one unique point in the phase space. For mechanical systems, the
226	+	phase space usually consists of all possible values of position and
227	+	momentum variables. Consider a dynamic system of $f$ particles in a
228	+	cartesian space, where each of the $6f$ coordinates and momenta is
229	+	assigned to one of $6f$ mutually orthogonal axes, the phase space of
230	+	this system is a $6f$ dimensional space. A point, $x = (q_1 , \ldots
231	+	,q_f ,p_1 , \ldots ,p_f )$, with a unique set of values of $6f$
232	+	coordinates and momenta is a phase space vector.
233	+
234	+	%%%fix me
235	+	A microscopic state or microstate of a classical system is
236	+	specification of the complete phase space vector of a system at any
237	+	instant in time. An ensemble is defined as a collection of systems
238	+	sharing one or more macroscopic characteristics but each being in a
239	+	unique microstate. The complete ensemble is specified by giving all
240	+	systems or microstates consistent with the common macroscopic
241	+	characteristics of the ensemble. Although the state of each
242	+	individual system in the ensemble could be precisely described at
243	+	any instance in time by a suitable phase space vector, when using
244	+	ensembles for statistical purposes, there is no need to maintain
245	+	distinctions between individual systems, since the numbers of
246	+	systems at any time in the different states which correspond to
247	+	different regions of the phase space are more interesting. Moreover,
248	+	in the point of view of statistical mechanics, one would prefer to
249	+	use ensembles containing a large enough population of separate
250	+	members so that the numbers of systems in such different states can
251	+	be regarded as changing continuously as we traverse different
252	+	regions of the phase space. The condition of an ensemble at any time
253	+	can be regarded as appropriately specified by the density $\rho$
254	+	with which representative points are distributed over the phase
255	+	space. The density distribution for an ensemble with $f$ degrees of
256	+	freedom is defined as,
257	+	\begin{equation}
258	+	\rho  = \rho (q_1 , \ldots ,q_f ,p_1 , \ldots ,p_f ,t).
259	+	\label{introEquation:densityDistribution}
260	+	\end{equation}
261	+	Governed by the principles of mechanics, the phase points change
262	+	their locations which would change the density at any time at phase
263	+	space. Hence, the density distribution is also to be taken as a
264	+	function of the time.
265	+
266	+	The number of systems $\delta N$ at time $t$ can be determined by,
267	+	\begin{equation}
268	+	\delta N = \rho (q,p,t)dq_1  \ldots dq_f dp_1  \ldots dp_f.
269	+	\label{introEquation:deltaN}
270	+	\end{equation}
271	+	Assuming a large enough population of systems, we can sufficiently
272	+	approximate $\delta N$ without introducing discontinuity when we go
273	+	from one region in the phase space to another. By integrating over
274	+	the whole phase space,
275	+	\begin{equation}
276	+	N = \int { \ldots \int {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f
277	+	\label{introEquation:totalNumberSystem}
278	+	\end{equation}
279	+	gives us an expression for the total number of the systems. Hence,
280	+	the probability per unit in the phase space can be obtained by,
281	+	\begin{equation}
282	+	\frac{{\rho (q,p,t)}}{N} = \frac{{\rho (q,p,t)}}{{\int { \ldots \int
283	+	{\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}.
284	+	\label{introEquation:unitProbability}
285	+	\end{equation}
286	+	With the help of Eq.~\ref{introEquation:unitProbability} and the
287	+	knowledge of the system, it is possible to calculate the average
288	+	value of any desired quantity which depends on the coordinates and
289	+	momenta of the system. Even when the dynamics of the real system is
290	+	complex, or stochastic, or even discontinuous, the average
291	+	properties of the ensemble of possibilities as a whole remaining
292	+	well defined. For a classical system in thermal equilibrium with its
293	+	environment, the ensemble average of a mechanical quantity, $\langle
294	+	A(q , p) \rangle_t$, takes the form of an integral over the phase
295	+	space of the system,
296	+	\begin{equation}
297	+	\langle  A(q , p) \rangle_t = \frac{{\int { \ldots \int {A(q,p)\rho
298	+	(q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}{{\int { \ldots \int {\rho
299	+	(q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}
300	+	\label{introEquation:ensembelAverage}
301	+	\end{equation}
302	+
303	+	There are several different types of ensembles with different
304	+	statistical characteristics. As a function of macroscopic
305	+	parameters, such as temperature \textit{etc}, the partition function
306	+	can be used to describe the statistical properties of a system in
307	+	thermodynamic equilibrium.
308	+
309	+	As an ensemble of systems, each of which is known to be thermally
310	+	isolated and conserve energy, the Microcanonical ensemble (NVE) has
311	+	a partition function like,
312	+	\begin{equation}
313	+	\Omega (N,V,E) = e^{\beta TS} \label{introEquation:NVEPartition}.
314	+	\end{equation}
315	+	A canonical ensemble (NVT)is an ensemble of systems, each of which
316	+	can share its energy with a large heat reservoir. The distribution
317	+	of the total energy amongst the possible dynamical states is given
318	+	by the partition function,
319	+	\begin{equation}
320	+	\Omega (N,V,T) = e^{ - \beta A}
321	+	\label{introEquation:NVTPartition}
322	+	\end{equation}
323	+	Here, $A$ is the Helmholtz free energy which is defined as $ A = U -
324	+	TS$. Since most experiments are carried out under constant pressure
325	+	condition, the isothermal-isobaric ensemble (NPT) plays a very
326	+	important role in molecular simulations. The isothermal-isobaric
327	+	ensemble allow the system to exchange energy with a heat bath of
328	+	temperature $T$ and to change the volume as well. Its partition
329	+	function is given as
330	+	\begin{equation}
331	+	\Delta (N,P,T) =  - e^{\beta G}.
332	+	\label{introEquation:NPTPartition}
333	+	\end{equation}
334	+	Here, $G = U - TS + PV$, and $G$ is called Gibbs free energy.
335	+
336	+	\subsection{\label{introSection:liouville}Liouville's theorem}
337	+
338	+	Liouville's theorem is the foundation on which statistical mechanics
339	+	rests. It describes the time evolution of the phase space
340	+	distribution function. In order to calculate the rate of change of
341	+	$\rho$, we begin from Eq.~\ref{introEquation:deltaN}. If we consider
342	+	the two faces perpendicular to the $q_1$ axis, which are located at
343	+	$q_1$ and $q_1 + \delta q_1$, the number of phase points leaving the
344	+	opposite face is given by the expression,
345	+	\begin{equation}
346	+	\left( {\rho  + \frac{{\partial \rho }}{{\partial q_1 }}\delta q_1 }
347	+	\right)\left( {\dot q_1  + \frac{{\partial \dot q_1 }}{{\partial q_1
348	+	}}\delta q_1 } \right)\delta q_2  \ldots \delta q_f \delta p_1
349	+	\ldots \delta p_f .
350	+	\end{equation}
351	+	Summing all over the phase space, we obtain
352	+	\begin{equation}
353	+	\frac{{d(\delta N)}}{{dt}} =  - \sum\limits_{i = 1}^f {\left[ {\rho
354	+	\left( {\frac{{\partial \dot q_i }}{{\partial q_i }} +
355	+	\frac{{\partial \dot p_i }}{{\partial p_i }}} \right) + \left(
356	+	{\frac{{\partial \rho }}{{\partial q_i }}\dot q_i  + \frac{{\partial
357	+	\rho }}{{\partial p_i }}\dot p_i } \right)} \right]} \delta q_1
358	+	\ldots \delta q_f \delta p_1  \ldots \delta p_f .
359	+	\end{equation}
360	+	Differentiating the equations of motion in Hamiltonian formalism
361	+	(\ref{introEquation:motionHamiltonianCoordinate},
362	+	\ref{introEquation:motionHamiltonianMomentum}), we can show,
363	+	\begin{equation}
364	+	\sum\limits_i {\left( {\frac{{\partial \dot q_i }}{{\partial q_i }}
365	+	+ \frac{{\partial \dot p_i }}{{\partial p_i }}} \right)}  = 0 ,
366	+	\end{equation}
367	+	which cancels the first terms of the right hand side. Furthermore,
368	+	dividing $ \delta q_1  \ldots \delta q_f \delta p_1  \ldots \delta
369	+	p_f $ in both sides, we can write out Liouville's theorem in a
370	+	simple form,
371	+	\begin{equation}
372	+	\frac{{\partial \rho }}{{\partial t}} + \sum\limits_{i = 1}^f
373	+	{\left( {\frac{{\partial \rho }}{{\partial q_i }}\dot q_i  +
374	+	\frac{{\partial \rho }}{{\partial p_i }}\dot p_i } \right)}  = 0 .
375	+	\label{introEquation:liouvilleTheorem}
376	+	\end{equation}
377	+
378	+	Liouville's theorem states that the distribution function is
379	+	constant along any trajectory in phase space. In classical
380	+	statistical mechanics, since the number of members in an ensemble is
381	+	huge and constant, we can assume the local density has no reason
382	+	(other than classical mechanics) to change,
383	+	\begin{equation}
384	+	\frac{{\partial \rho }}{{\partial t}} = 0.
385	+	\label{introEquation:stationary}
386	+	\end{equation}
387	+	In such stationary system, the density of distribution $\rho$ can be
388	+	connected to the Hamiltonian $H$ through Maxwell-Boltzmann
389	+	distribution,
390	+	\begin{equation}
391	+	\rho  \propto e^{ - \beta H}
392	+	\label{introEquation:densityAndHamiltonian}
393	+	\end{equation}
394	+
395	+	\subsubsection{\label{introSection:phaseSpaceConservation}\textbf{Conservation of Phase Space}}
396	+	Lets consider a region in the phase space,
397	+	\begin{equation}
398	+	\delta v = \int { \ldots \int {dq_1 } ...dq_f dp_1 } ..dp_f .
399	+	\end{equation}
400	+	If this region is small enough, the density $\rho$ can be regarded
401	+	as uniform over the whole integral. Thus, the number of phase points
402	+	inside this region is given by,
403	+	\begin{equation}
404	+	\delta N = \rho \delta v = \rho \int { \ldots \int {dq_1 } ...dq_f
405	+	dp_1 } ..dp_f.
406	+	\end{equation}
407	+
408	+	\begin{equation}
409	+	\frac{{d(\delta N)}}{{dt}} = \frac{{d\rho }}{{dt}}\delta v + \rho
410	+	\frac{d}{{dt}}(\delta v) = 0.
411	+	\end{equation}
412	+	With the help of stationary assumption
413	+	(\ref{introEquation:stationary}), we obtain the principle of the
414	+	\emph{conservation of volume in phase space},
415	+	\begin{equation}
416	+	\frac{d}{{dt}}(\delta v) = \frac{d}{{dt}}\int { \ldots \int {dq_1 }
417	+	...dq_f dp_1 } ..dp_f  = 0.
418	+	\label{introEquation:volumePreserving}
419	+	\end{equation}
420	+
421	+	\subsubsection{\label{introSection:liouvilleInOtherForms}\textbf{Liouville's Theorem in Other Forms}}
422	+
423	+	Liouville's theorem can be expresses in a variety of different forms
424	+	which are convenient within different contexts. For any two function
425	+	$F$ and $G$ of the coordinates and momenta of a system, the Poisson
426	+	bracket ${F, G}$ is defined as
427	+	\begin{equation}
428	+	\left\{ {F,G} \right\} = \sum\limits_i {\left( {\frac{{\partial
429	+	F}}{{\partial q_i }}\frac{{\partial G}}{{\partial p_i }} -
430	+	\frac{{\partial F}}{{\partial p_i }}\frac{{\partial G}}{{\partial
431	+	q_i }}} \right)}.
432	+	\label{introEquation:poissonBracket}
433	+	\end{equation}
434	+	Substituting equations of motion in Hamiltonian formalism(
435	+	Eq.~\ref{introEquation:motionHamiltonianCoordinate} ,
436	+	Eq.~\ref{introEquation:motionHamiltonianMomentum} ) into
437	+	(Eq.~\ref{introEquation:liouvilleTheorem}), we can rewrite
438	+	Liouville's theorem using Poisson bracket notion,
439	+	\begin{equation}
440	+	\left( {\frac{{\partial \rho }}{{\partial t}}} \right) =  - \left\{
441	+	{\rho ,H} \right\}.
442	+	\label{introEquation:liouvilleTheromInPoissin}
443	+	\end{equation}
444	+	Moreover, the Liouville operator is defined as
445	+	\begin{equation}
446	+	iL = \sum\limits_{i = 1}^f {\left( {\frac{{\partial H}}{{\partial
447	+	p_i }}\frac{\partial }{{\partial q_i }} - \frac{{\partial
448	+	H}}{{\partial q_i }}\frac{\partial }{{\partial p_i }}} \right)}
449	+	\label{introEquation:liouvilleOperator}
450	+	\end{equation}
451	+	In terms of Liouville operator, Liouville's equation can also be
452	+	expressed as
453	+	\begin{equation}
454	+	\left( {\frac{{\partial \rho }}{{\partial t}}} \right) =  - iL\rho
455	+	\label{introEquation:liouvilleTheoremInOperator}
456	+	\end{equation}
457	+
458		\subsection{\label{introSection:ergodic}The Ergodic Hypothesis}
459
460		Various thermodynamic properties can be calculated from Molecular
461		Dynamics simulation. By comparing experimental values with the
462		calculated properties, one can determine the accuracy of the
463	<	simulation and the quality of the underlying model. However, both of
464	<	experiment and computer simulation are usually performed during a
463	>	simulation and the quality of the underlying model. However, both
464	>	experiments and computer simulations are usually performed during a
465		certain time interval and the measurements are averaged over a
466		period of them which is different from the average behavior of
467	<	many-body system in Statistical Mechanics. Fortunately, Ergodic
468	<	Hypothesis is proposed to make a connection between time average and
469	<	ensemble average. It states that time average and average over the
470	<	statistical ensemble are identical \cite{Frenkel1996, leach01:mm}.
467	>	many-body system in Statistical Mechanics. Fortunately, the Ergodic
468	>	Hypothesis makes a connection between time average and the ensemble
469	>	average. It states that the time average and average over the
470	>	statistical ensemble are identical \cite{Frenkel1996, Leach2001}.
471		\begin{equation}
472	<	\langle A \rangle_t = \mathop {\lim }\limits_{t \to \infty }
473	<	\frac{1}{t}\int\limits_0^t {A(p(t),q(t))dt = \int\limits_\Gamma
474	<	{A(p(t),q(t))} } \rho (p(t), q(t)) dpdq
472	>	\langle A(q , p) \rangle_t = \mathop {\lim }\limits_{t \to \infty }
473	>	\frac{1}{t}\int\limits_0^t {A(q(t),p(t))dt = \int\limits_\Gamma
474	>	{A(q(t),p(t))} } \rho (q(t), p(t)) dqdp
475		\end{equation}
476	<	where $\langle A \rangle_t$ is an equilibrium value of a physical
477	<	quantity and $\rho (p(t), q(t))$ is the equilibrium distribution
478	<	function. If an observation is averaged over a sufficiently long
479	<	time (longer than relaxation time), all accessible microstates in
480	<	phase space are assumed to be equally probed, giving a properly
481	<	weighted statistical average. This allows the researcher freedom of
482	<	choice when deciding how best to measure a given observable. In case
483	<	an ensemble averaged approach sounds most reasonable, the Monte
484	<	Carlo techniques\cite{metropolis:1949} can be utilized. Or if the
485	<	system lends itself to a time averaging approach, the Molecular
486	<	Dynamics techniques in Sec.~\ref{introSection:molecularDynamics}
487	<	will be the best choice\cite{Frenkel1996}.
476	>	where $\langle  A(q , p) \rangle_t$ is an equilibrium value of a
477	>	physical quantity and $\rho (p(t), q(t))$ is the equilibrium
478	>	distribution function. If an observation is averaged over a
479	>	sufficiently long time (longer than relaxation time), all accessible
480	>	microstates in phase space are assumed to be equally probed, giving
481	>	a properly weighted statistical average. This allows the researcher
482	>	freedom of choice when deciding how best to measure a given
483	>	observable. In case an ensemble averaged approach sounds most
484	>	reasonable, the Monte Carlo techniques\cite{Metropolis1949} can be
485	>	utilized. Or if the system lends itself to a time averaging
486	>	approach, the Molecular Dynamics techniques in
487	>	Sec.~\ref{introSection:molecularDynamics} will be the best
488	>	choice\cite{Frenkel1996}.
489
490		\section{\label{introSection:geometricIntegratos}Geometric Integrators}
491	<	A variety of numerical integrators were proposed to simulate the
492	<	motions. They usually begin with an initial conditionals and move
493	<	the objects in the direction governed by the differential equations.
494	<	However, most of them ignore the hidden physical law contained
495	<	within the equations. Since 1990, geometric integrators, which
496	<	preserve various phase-flow invariants such as symplectic structure,
497	<	volume and time reversal symmetry, are developed to address this
498	<	issue. The velocity verlet method, which happens to be a simple
499	<	example of symplectic integrator, continues to gain its popularity
500	<	in molecular dynamics community. This fact can be partly explained
501	<	by its geometric nature.
491	>	A variety of numerical integrators have been proposed to simulate
492	>	the motions of atoms in MD simulation. They usually begin with
493	>	initial conditionals and move the objects in the direction governed
494	>	by the differential equations. However, most of them ignore the
495	>	hidden physical laws contained within the equations. Since 1990,
496	>	geometric integrators, which preserve various phase-flow invariants
497	>	such as symplectic structure, volume and time reversal symmetry, are
498	>	developed to address this issue\cite{Dullweber1997, McLachlan1998,
499	>	Leimkuhler1999}. The velocity verlet method, which happens to be a
500	>	simple example of symplectic integrator, continues to gain
501	>	popularity in the molecular dynamics community. This fact can be
502	>	partly explained by its geometric nature.
503
504	<	\subsection{\label{introSection:symplecticManifold}Symplectic Manifold}
505	<	A \emph{manifold} is an abstract mathematical space. It locally
506	<	looks like Euclidean space, but when viewed globally, it may have
507	<	more complicate structure. A good example of manifold is the surface
508	<	of Earth. It seems to be flat locally, but it is round if viewed as
509	<	a whole. A \emph{differentiable manifold} (also known as
510	<	\emph{smooth manifold}) is a manifold with an open cover in which
511	<	the covering neighborhoods are all smoothly isomorphic to one
512	<	another. In other words,it is possible to apply calculus on
303	<	\emph{differentiable manifold}. A \emph{symplectic manifold} is
304	<	defined as a pair $(M, \omega)$ which consisting of a
504	>	\subsection{\label{introSection:symplecticManifold}Symplectic Manifolds}
505	>	A \emph{manifold} is an abstract mathematical space. It looks
506	>	locally like Euclidean space, but when viewed globally, it may have
507	>	more complicated structure. A good example of manifold is the
508	>	surface of Earth. It seems to be flat locally, but it is round if
509	>	viewed as a whole. A \emph{differentiable manifold} (also known as
510	>	\emph{smooth manifold}) is a manifold on which it is possible to
511	>	apply calculus on \emph{differentiable manifold}. A \emph{symplectic
512	>	manifold} is defined as a pair $(M, \omega)$ which consists of a
513		\emph{differentiable manifold} $M$ and a close, non-degenerated,
514		bilinear symplectic form, $\omega$. A symplectic form on a vector
515		space $V$ is a function $\omega(x, y)$ which satisfies
516		$\omega(\lambda_1x_1+\lambda_2x_2, y) = \lambda_1\omega(x_1, y)+
517		\lambda_2\omega(x_2, y)$, $\omega(x, y) = - \omega(y, x)$ and
518	<	$\omega(x, x) = 0$. Cross product operation in vector field is an
519	<	example of symplectic form.
518	>	$\omega(x, x) = 0$. The cross product operation in vector field is
519	>	an example of symplectic form.
520
521	<	One of the motivations to study \emph{symplectic manifold} in
521	>	One of the motivations to study \emph{symplectic manifolds} in
522		Hamiltonian Mechanics is that a symplectic manifold can represent
523		all possible configurations of the system and the phase space of the
524		system can be described by it's cotangent bundle. Every symplectic
525		manifold is even dimensional. For instance, in Hamilton equations,
526		coordinate and momentum always appear in pairs.
527
528	<	Let  $(M,\omega)$ and $(N, \eta)$ be symplectic manifolds. A map
321	<	\[
322	<	f : M \rightarrow N
323	<	\]
324	<	is a \emph{symplectomorphism} if it is a \emph{diffeomorphims} and
325	<	the \emph{pullback} of $\eta$ under f is equal to $\omega$.
326	<	Canonical transformation is an example of symplectomorphism in
327	<	classical mechanics. According to Liouville's theorem, the
328	<	Hamiltonian \emph{flow} or \emph{symplectomorphism} generated by the
329	<	Hamiltonian vector filed preserves the volume form on the phase
330	<	space, which is the basis of classical statistical mechanics.
528	>	\subsection{\label{introSection:ODE}Ordinary Differential Equations}
529
530	<	\subsection{\label{introSection:exactFlow}The Exact Flow of ODE}
530	>	For an ordinary differential system defined as
531	>	\begin{equation}
532	>	\dot x = f(x)
533	>	\end{equation}
534	>	where $x = x(q,p)^T$, this system is a canonical Hamiltonian, if
535	>	\begin{equation}
536	>	f(r) = J\nabla _x H(r).
537	>	\end{equation}
538	>	$H = H (q, p)$ is Hamiltonian function and $J$ is the skew-symmetric
539	>	matrix
540	>	\begin{equation}
541	>	J = \left( {\begin{array}{*{20}c}
542	>	0 & I  \\
543	>	{ - I} & 0  \\
544	>	\end{array}} \right)
545	>	\label{introEquation:canonicalMatrix}
546	>	\end{equation}
547	>	where $I$ is an identity matrix. Using this notation, Hamiltonian
548	>	system can be rewritten as,
549	>	\begin{equation}
550	>	\frac{d}{{dt}}x = J\nabla _x H(x)
551	>	\label{introEquation:compactHamiltonian}
552	>	\end{equation}In this case, $f$ is
553	>	called a \emph{Hamiltonian vector field}.
554
555	<	\subsection{\label{introSection:hamiltonianSplitting}Hamiltonian Splitting}
556	<
557	<	\section{\label{introSection:molecularDynamics}Molecular Dynamics}
558	<
559	<	As a special discipline of molecular modeling, Molecular dynamics
560	<	has proven to be a powerful tool for studying the functions of
340	<	biological systems, providing structural, thermodynamic and
341	<	dynamical information.
555	>	Another generalization of Hamiltonian dynamics is Poisson
556	>	Dynamics\cite{Olver1986},
557	>	\begin{equation}
558	>	\dot x = J(x)\nabla _x H \label{introEquation:poissonHamiltonian}
559	>	\end{equation}
560	>	The most obvious change being that matrix $J$ now depends on $x$.
561
562	<	\subsection{\label{introSec:mdInit}Initialization}
562	>	\subsection{\label{introSection:exactFlow}Exact Flow}
563
564	<	\subsection{\label{introSection:mdIntegration} Integration of the Equations of Motion}
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. The flow has the continuous group property,
574	>	\begin{equation}
575	>	\varphi _{\tau _1 }  \circ \varphi _{\tau _2 }  = \varphi _{\tau _1
576	>	+ \tau _2 } .
577	>	\end{equation}
578	>	In particular,
579	>	\begin{equation}
580	>	\varphi _\tau   \circ \varphi _{ - \tau }  = I
581	>	\end{equation}
582	>	Therefore, the exact flow is self-adjoint,
583	>	\begin{equation}
584	>	\varphi _\tau   = \varphi _{ - \tau }^{ - 1}.
585	>	\end{equation}
586	>	The exact flow can also be written in terms of the of an operator,
587	>	\begin{equation}
588	>	\varphi _\tau  (x) = e^{\tau \sum\limits_i {f_i (x)\frac{\partial
589	>	}{{\partial x_i }}} } (x) \equiv \exp (\tau f)(x).
590	>	\label{introEquation:exponentialOperator}
591	>	\end{equation}
592
593	<	\section{\label{introSection:rigidBody}Dynamics of Rigid Bodies}
593	>	In most cases, it is not easy to find the exact flow $\varphi_\tau$.
594	>	Instead, we use a approximate map, $\psi_\tau$, which is usually
595	>	called integrator. The order of an integrator $\psi_\tau$ is $p$, if
596	>	the Taylor series of $\psi_\tau$ agree to order $p$,
597	>	\begin{equation}
598	>	\psi_tau(x) = x + \tau f(x) + O(\tau^{p+1})
599	>	\end{equation}
600
601	<	A rigid body is a body in which the distance between any two given
350	<	points of a rigid body remains constant regardless of external
351	<	forces exerted on it. A rigid body therefore conserves its shape
352	<	during its motion.
601	>	\subsection{\label{introSection:geometricProperties}Geometric Properties}
602
603	<	Applications of dynamics of rigid bodies.
603	>	The hidden geometric properties\cite{Budd1999, Marsden1998} of ODE
604	>	and its flow play important roles in numerical studies. Many of them
605	>	can be found in systems which occur naturally in applications.
606
607	<	\subsection{\label{introSection:lieAlgebra}Lie Algebra}
607	>	Let $\varphi$ be the flow of Hamiltonian vector field, $\varphi$ is
608	>	a \emph{symplectic} flow if it satisfies,
609	>	\begin{equation}
610	>	{\varphi '}^T J \varphi ' = J.
611	>	\end{equation}
612	>	According to Liouville's theorem, the symplectic volume is invariant
613	>	under a Hamiltonian flow, which is the basis for classical
614	>	statistical mechanics. Furthermore, the flow of a Hamiltonian vector
615	>	field on a symplectic manifold can be shown to be a
616	>	symplectomorphism. As to the Poisson system,
617	>	\begin{equation}
618	>	{\varphi '}^T J \varphi ' = J \circ \varphi
619	>	\end{equation}
620	>	is the property must be preserved by the integrator.
621
622	<	\subsection{\label{introSection:DLMMotionEquation}The Euler Equations of Rigid Body Motion}
622	>	It is possible to construct a \emph{volume-preserving} flow for a
623	>	source free($ \nabla \cdot f = 0 $) ODE, if the flow satisfies $
624	>	\det d\varphi  = 1$. One can show easily that a symplectic flow will
625	>	be volume-preserving.
626
627	<	\subsection{\label{introSection:otherRBMotionEquation}Other Formulations for Rigid Body Motion}
627	>	Changing the variables $y = h(x)$ in a ODE\ref{introEquation:ODE}
628	>	will result in a new system,
629	>	\[
630	>	\dot y = \tilde f(y) = ((dh \cdot f)h^{ - 1} )(y).
631	>	\]
632	>	The vector filed $f$ has reversing symmetry $h$ if $f = - \tilde f$.
633	>	In other words, the flow of this vector field is reversible if and
634	>	only if $ h \circ \varphi ^{ - 1}  = \varphi  \circ h $.
635
636	<	%\subsection{\label{introSection:poissonBrackets}Poisson Brackets}
636	>	A \emph{first integral}, or conserved quantity of a general
637	>	differential function is a function $ G:R^{2d}  \to R^d $ which is
638	>	constant for all solutions of the ODE $\frac{{dx}}{{dt}} = f(x)$ ,
639	>	\[
640	>	\frac{{dG(x(t))}}{{dt}} = 0.
641	>	\]
642	>	Using chain rule, one may obtain,
643	>	\[
644	>	\sum\limits_i {\frac{{dG}}{{dx_i }}} f_i (x) = f \bullet \nabla G,
645	>	\]
646	>	which is the condition for conserving \emph{first integral}. For a
647	>	canonical Hamiltonian system, the time evolution of an arbitrary
648	>	smooth function $G$ is given by,
649
650	<	\section{\label{introSection:correlationFunctions}Correlation Functions}
650	>	\begin{eqnarray}
651	>	\frac{{dG(x(t))}}{{dt}} & = & [\nabla _x G(x(t))]^T \dot x(t) \\
652	>	& = & [\nabla _x G(x(t))]^T J\nabla _x H(x(t)). \\
653	>	\label{introEquation:firstIntegral1}
654	>	\end{eqnarray}
655
366	–	\section{\label{introSection:langevinDynamics}Langevin Dynamics}
656
657	<	\subsection{\label{introSection:LDIntroduction}Introduction and application of Langevin Dynamics}
657	>	Using poisson bracket notion, Equation
658	>	\ref{introEquation:firstIntegral1} can be rewritten as
659	>	\[
660	>	\frac{d}{{dt}}G(x(t)) = \left\{ {G,H} \right\}(x(t)).
661	>	\]
662	>	Therefore, the sufficient condition for $G$ to be the \emph{first
663	>	integral} of a Hamiltonian system is
664	>	\[
665	>	\left\{ {G,H} \right\} = 0.
666	>	\]
667	>	As well known, the Hamiltonian (or energy) H of a Hamiltonian system
668	>	is a \emph{first integral}, which is due to the fact $\{ H,H\}  =
669	>	0$.
670
671	<	\subsection{\label{introSection:generalizedLangevinDynamics}Generalized Langevin Dynamics}
671	>	When designing any numerical methods, one should always try to
672	>	preserve the structural properties of the original ODE and its flow.
673
674	+	\subsection{\label{introSection:constructionSymplectic}Construction of Symplectic Methods}
675	+	A lot of well established and very effective numerical methods have
676	+	been successful precisely because of their symplecticities even
677	+	though this fact was not recognized when they were first
678	+	constructed. The most famous example is the Verlet-leapfrog methods
679	+	in molecular dynamics. In general, symplectic integrators can be
680	+	constructed using one of four different methods.
681	+	\begin{enumerate}
682	+	\item Generating functions
683	+	\item Variational methods
684	+	\item Runge-Kutta methods
685	+	\item Splitting methods
686	+	\end{enumerate}
687	+
688	+	Generating function\cite{Channell1990} tends to lead to methods
689	+	which are cumbersome and difficult to use. In dissipative systems,
690	+	variational methods can capture the decay of energy
691	+	accurately\cite{Kane2000}. Since their geometrically unstable nature
692	+	against non-Hamiltonian perturbations, ordinary implicit Runge-Kutta
693	+	methods are not suitable for Hamiltonian system. Recently, various
694	+	high-order explicit Runge-Kutta methods
695	+	\cite{Owren1992,Chen2003}have been developed to overcome this
696	+	instability. However, due to computational penalty involved in
697	+	implementing the Runge-Kutta methods, they have not attracted much
698	+	attention from the Molecular Dynamics community. Instead, splitting
699	+	methods have been widely accepted since they exploit natural
700	+	decompositions of the system\cite{Tuckerman1992, McLachlan1998}.
701	+
702	+	\subsubsection{\label{introSection:splittingMethod}\textbf{Splitting Methods}}
703	+
704	+	The main idea behind splitting methods is to decompose the discrete
705	+	$\varphi_h$ as a composition of simpler flows,
706		\begin{equation}
707	<	H = \frac{{p^2 }}{{2m}} + U(x) + H_B  + \Delta U(x,x_1 , \ldots x_N)
708	<	\label{introEquation:bathGLE}
707	>	\varphi _h  = \varphi _{h_1 }  \circ \varphi _{h_2 }  \ldots  \circ
708	>	\varphi _{h_n }
709	>	\label{introEquation:FlowDecomposition}
710		\end{equation}
711	<	where $H_B$ is harmonic bath Hamiltonian,
711	>	where each of the sub-flow is chosen such that each represent a
712	>	simpler integration of the system.
713	>
714	>	Suppose that a Hamiltonian system takes the form,
715		\[
716	<	H_B  =\sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
379	<	}}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  w_\alpha ^2 } \right\}}
716	>	H = H_1 + H_2.
717		\]
718	<	and $\Delta U$ is bilinear system-bath coupling,
718	>	Here, $H_1$ and $H_2$ may represent different physical processes of
719	>	the system. For instance, they may relate to kinetic and potential
720	>	energy respectively, which is a natural decomposition of the
721	>	problem. If $H_1$ and $H_2$ can be integrated using exact flows
722	>	$\varphi_1(t)$ and $\varphi_2(t)$, respectively, a simple first
723	>	order expression is then given by the Lie-Trotter formula
724	>	\begin{equation}
725	>	\varphi _h  = \varphi _{1,h}  \circ \varphi _{2,h},
726	>	\label{introEquation:firstOrderSplitting}
727	>	\end{equation}
728	>	where $\varphi _h$ is the result of applying the corresponding
729	>	continuous $\varphi _i$ over a time $h$. By definition, as
730	>	$\varphi_i(t)$ is the exact solution of a Hamiltonian system, it
731	>	must follow that each operator $\varphi_i(t)$ is a symplectic map.
732	>	It is easy to show that any composition of symplectic flows yields a
733	>	symplectic map,
734	>	\begin{equation}
735	>	(\varphi '\phi ')^T J\varphi '\phi ' = \phi '^T \varphi '^T J\varphi
736	>	'\phi ' = \phi '^T J\phi ' = J,
737	>	\label{introEquation:SymplecticFlowComposition}
738	>	\end{equation}
739	>	where $\phi$ and $\psi$ both are symplectic maps. Thus operator
740	>	splitting in this context automatically generates a symplectic map.
741	>
742	>	The Lie-Trotter splitting(\ref{introEquation:firstOrderSplitting})
743	>	introduces local errors proportional to $h^2$, while Strang
744	>	splitting gives a second-order decomposition,
745	>	\begin{equation}
746	>	\varphi _h  = \varphi _{1,h/2}  \circ \varphi _{2,h}  \circ \varphi
747	>	_{1,h/2} , \label{introEquation:secondOrderSplitting}
748	>	\end{equation}
749	>	which has a local error proportional to $h^3$. The Sprang
750	>	splitting's popularity in molecular simulation community attribute
751	>	to its symmetric property,
752	>	\begin{equation}
753	>	\varphi _h^{ - 1} = \varphi _{ - h}.
754	>	\label{introEquation:timeReversible}
755	>	\end{equation},appendixFig:architecture
756	>
757	>	\subsubsection{\label{introSection:exampleSplittingMethod}\textbf{Example of Splitting Method}}
758	>	The classical equation for a system consisting of interacting
759	>	particles can be written in Hamiltonian form,
760		\[
761	<	\Delta U =  - \sum\limits_{\alpha  = 1}^N {g_\alpha  x_\alpha  x}
761	>	H = T + V
762		\]
763	<	Completing the square,
764	<	\[
765	<	H_B  + \Delta U = \sum\limits_{\alpha  = 1}^N {\left\{
388	<	{\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
389	<	w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
390	<	w_\alpha ^2 }}x} \right)^2 } \right\}}  - \sum\limits_{\alpha  =
391	<	1}^N {\frac{{g_\alpha ^2 }}{{2m_\alpha  w_\alpha ^2 }}} x^2
392	<	\]
393	<	and putting it back into Eq.~\ref{introEquation:bathGLE},
394	<	\[
395	<	H = \frac{{p^2 }}{{2m}} + W(x) + \sum\limits_{\alpha  = 1}^N
396	<	{\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
397	<	w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
398	<	w_\alpha ^2 }}x} \right)^2 } \right\}}
399	<	\]
400	<	where
401	<	\[
402	<	W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
403	<	}}{{2m_\alpha  w_\alpha ^2 }}} x^2
404	<	\]
405	<	Since the first two terms of the new Hamiltonian depend only on the
406	<	system coordinates, we can get the equations of motion for
407	<	Generalized Langevin Dynamics by Hamilton's equations
408	<	\ref{introEquation:motionHamiltonianCoordinate,
409	<	introEquation:motionHamiltonianMomentum},
763	>	where $T$ is the kinetic energy and $V$ is the potential energy.
764	>	Setting $H_1 = T, H_2 = V$ and applying Strang splitting, one
765	>	obtains the following:
766		\begin{align}
767	<	\dot p &=  - \frac{{\partial H}}{{\partial x}}
768	<	&= m\ddot x
769	<	&= - \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)}
770	<	\label{introEq:Lp5}
767	>	q(\Delta t) &= q(0) + \dot{q}(0)\Delta t +
768	>	\frac{F[q(0)]}{m}\frac{\Delta t^2}{2}, %
769	>	\label{introEquation:Lp10a} \\%
770	>	%
771	>	\dot{q}(\Delta t) &= \dot{q}(0) + \frac{\Delta t}{2m}
772	>	\biggl [F[q(0)] + F[q(\Delta t)] \biggr]. %
773	>	\label{introEquation:Lp10b}
774		\end{align}
775	<	, and
775	>	where $F(t)$ is the force at time $t$. This integration scheme is
776	>	known as \emph{velocity verlet} which is
777	>	symplectic(\ref{introEquation:SymplecticFlowComposition}),
778	>	time-reversible(\ref{introEquation:timeReversible}) and
779	>	volume-preserving (\ref{introEquation:volumePreserving}). These
780	>	geometric properties attribute to its long-time stability and its
781	>	popularity in the community. However, the most commonly used
782	>	velocity verlet integration scheme is written as below,
783		\begin{align}
784	<	\dot p_\alpha   &=  - \frac{{\partial H}}{{\partial x_\alpha  }}
785	<	&= m\ddot x_\alpha
786	<	&= \- m_\alpha  w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha}}{{m_\alpha  w_\alpha ^2 }}x} \right)
784	>	\dot{q}\biggl (\frac{\Delta t}{2}\biggr ) &=
785	>	\dot{q}(0) + \frac{\Delta t}{2m}\, F[q(0)], \label{introEquation:Lp9a}\\%
786	>	%
787	>	q(\Delta t) &= q(0) + \Delta t\, \dot{q}\biggl (\frac{\Delta t}{2}\biggr ),%
788	>	\label{introEquation:Lp9b}\\%
789	>	%
790	>	\dot{q}(\Delta t) &= \dot{q}\biggl (\frac{\Delta t}{2}\biggr ) +
791	>	\frac{\Delta t}{2m}\, F[q(0)]. \label{introEquation:Lp9c}
792		\end{align}
793	+	From the preceding splitting, one can see that the integration of
794	+	the equations of motion would follow:
795	+	\begin{enumerate}
796	+	\item calculate the velocities at the half step, $\frac{\Delta t}{2}$, from the forces calculated at the initial position.
797
798	<	\subsection{\label{introSection:laplaceTransform}The Laplace Transform}
798	>	\item Use the half step velocities to move positions one whole step, $\Delta t$.
799
800	<	\[
426	<	L(x) = \int_0^\infty  {x(t)e^{ - pt} dt}
427	<	\]
800	>	\item Evaluate the forces at the new positions, $\mathbf{r}(\Delta t)$, and use the new forces to complete the velocity move.
801
802	<	\[
803	<	L(x + y) = L(x) + L(y)
431	<	\]
802	>	\item Repeat from step 1 with the new position, velocities, and forces assuming the roles of the initial values.
803	>	\end{enumerate}
804
805	<	\[
806	<	L(ax) = aL(x)
807	<	\]
805	>	Simply switching the order of splitting and composing, a new
806	>	integrator, the \emph{position verlet} integrator, can be generated,
807	>	\begin{align}
808	>	\dot q(\Delta t) &= \dot q(0) + \Delta tF(q(0))\left[ {q(0) +
809	>	\frac{{\Delta t}}{{2m}}\dot q(0)} \right], %
810	>	\label{introEquation:positionVerlet1} \\%
811	>	%
812	>	q(\Delta t) &= q(0) + \frac{{\Delta t}}{2}\left[ {\dot q(0) + \dot
813	>	q(\Delta t)} \right]. %
814	>	\label{introEquation:positionVerlet2}
815	>	\end{align}
816
817	<	\[
438	<	L(\dot x) = pL(x) - px(0)
439	<	\]
817	>	\subsubsection{\label{introSection:errorAnalysis}\textbf{Error Analysis and Higher Order Methods}}
818
819	+	Baker-Campbell-Hausdorff formula can be used to determine the local
820	+	error of splitting method in terms of commutator of the
821	+	operators(\ref{introEquation:exponentialOperator}) associated with
822	+	the sub-flow. For operators $hX$ and $hY$ which are associate to
823	+	$\varphi_1(t)$ and $\varphi_2(t)$ respectively , we have
824	+	\begin{equation}
825	+	\exp (hX + hY) = \exp (hZ)
826	+	\end{equation}
827	+	where
828	+	\begin{equation}
829	+	hZ = hX + hY + \frac{{h^2 }}{2}[X,Y] + \frac{{h^3 }}{2}\left(
830	+	{[X,[X,Y]] + [Y,[Y,X]]} \right) +  \ldots .
831	+	\end{equation}
832	+	Here, $[X,Y]$ is the commutators of operator $X$ and $Y$ given by
833		\[
834	<	L(\ddot x) = p^2 L(x) - px(0) - \dot x(0)
834	>	[X,Y] = XY - YX .
835		\]
836	<
836	>	Applying Baker-Campbell-Hausdorff formula\cite{Varadarajan1974} to
837	>	Sprang splitting, we can obtain
838	>	\begin{eqnarray*}
839	>	\exp (h X/2)\exp (h Y)\exp (h X/2) & = & \exp (h X + h Y + h^2 [X,Y]/4 + h^2 [Y,X]/4 \\
840	>	&   & \mbox{} + h^2 [X,X]/8 + h^2 [Y,Y]/8 \\
841	>	&   & \mbox{} + h^3 [Y,[Y,X]]/12 - h^3[X,[X,Y]]/24 + \ldots )
842	>	\end{eqnarray*}
843	>	Since \[ [X,Y] + [Y,X] = 0\] and \[ [X,X] = 0\], the dominant local
844	>	error of Spring splitting is proportional to $h^3$. The same
845	>	procedure can be applied to general splitting,  of the form
846	>	\begin{equation}
847	>	\varphi _{b_m h}^2  \circ \varphi _{a_m h}^1  \circ \varphi _{b_{m -
848	>	1} h}^2  \circ  \ldots  \circ \varphi _{a_1 h}^1 .
849	>	\end{equation}
850	>	Careful choice of coefficient $a_1 \ldots b_m$ will lead to higher
851	>	order method. Yoshida proposed an elegant way to compose higher
852	>	order methods based on symmetric splitting\cite{Yoshida1990}. Given
853	>	a symmetric second order base method $ \varphi _h^{(2)} $, a
854	>	fourth-order symmetric method can be constructed by composing,
855		\[
856	<	L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right) = G(p)H(p)
856	>	\varphi _h^{(4)}  = \varphi _{\alpha h}^{(2)}  \circ \varphi _{\beta
857	>	h}^{(2)}  \circ \varphi _{\alpha h}^{(2)}
858		\]
859	<
860	<	Some relatively important transformation,
859	>	where $ \alpha  =  - \frac{{2^{1/3} }}{{2 - 2^{1/3} }}$ and $ \beta
860	>	= \frac{{2^{1/3} }}{{2 - 2^{1/3} }}$. Moreover, a symmetric
861	>	integrator $ \varphi _h^{(2n + 2)}$ can be composed by
862	>	\begin{equation}
863	>	\varphi _h^{(2n + 2)}  = \varphi _{\alpha h}^{(2n)}  \circ \varphi
864	>	_{\beta h}^{(2n)}  \circ \varphi _{\alpha h}^{(2n)}
865	>	\end{equation}
866	>	, if the weights are chosen as
867		\[
868	<	L(\cos at) = \frac{p}{{p^2  + a^2 }}
868	>	\alpha  =  - \frac{{2^{1/(2n + 1)} }}{{2 - 2^{1/(2n + 1)} }},\beta =
869	>	\frac{{2^{1/(2n + 1)} }}{{2 - 2^{1/(2n + 1)} }} .
870		\]
871
872	<	\[
455	<	L(\sin at) = \frac{a}{{p^2  + a^2 }}
456	<	\]
872	>	\section{\label{introSection:molecularDynamics}Molecular Dynamics}
873
874	+	As one of the principal tools of molecular modeling, Molecular
875	+	dynamics has proven to be a powerful tool for studying the functions
876	+	of biological systems, providing structural, thermodynamic and
877	+	dynamical information. The basic idea of molecular dynamics is that
878	+	macroscopic properties are related to microscopic behavior and
879	+	microscopic behavior can be calculated from the trajectories in
880	+	simulations. For instance, instantaneous temperature of an
881	+	Hamiltonian system of $N$ particle can be measured by
882		\[
883	<	L(1) = \frac{1}{p}
883	>	T = \sum\limits_{i = 1}^N {\frac{{m_i v_i^2 }}{{fk_B }}}
884		\]
885	+	where $m_i$ and $v_i$ are the mass and velocity of $i$th particle
886	+	respectively, $f$ is the number of degrees of freedom, and $k_B$ is
887	+	the boltzman constant.
888
889	<	First, the bath coordinates,
890	<	\[
891	<	p^2 L(x_\alpha  ) - px_\alpha  (0) - \dot x_\alpha  (0) =  - \omega
892	<	_\alpha ^2 L(x_\alpha  ) + \frac{{g_\alpha  }}{{\omega _\alpha
893	<	}}L(x)
889	>	A typical molecular dynamics run consists of three essential steps:
890	>	\begin{enumerate}
891	>	\item Initialization
892	>	\begin{enumerate}
893	>	\item Preliminary preparation
894	>	\item Minimization
895	>	\item Heating
896	>	\item Equilibration
897	>	\end{enumerate}
898	>	\item Production
899	>	\item Analysis
900	>	\end{enumerate}
901	>	These three individual steps will be covered in the following
902	>	sections. Sec.~\ref{introSec:initialSystemSettings} deals with the
903	>	initialization of a simulation. Sec.~\ref{introSection:production}
904	>	will discusses issues in production run.
905	>	Sec.~\ref{introSection:Analysis} provides the theoretical tools for
906	>	trajectory analysis.
907	>
908	>	\subsection{\label{introSec:initialSystemSettings}Initialization}
909	>
910	>	\subsubsection{\textbf{Preliminary preparation}}
911	>
912	>	When selecting the starting structure of a molecule for molecular
913	>	simulation, one may retrieve its Cartesian coordinates from public
914	>	databases, such as RCSB Protein Data Bank \textit{etc}. Although
915	>	thousands of crystal structures of molecules are discovered every
916	>	year, many more remain unknown due to the difficulties of
917	>	purification and crystallization. Even for the molecule with known
918	>	structure, some important information is missing. For example, the
919	>	missing hydrogen atom which acts as donor in hydrogen bonding must
920	>	be added. Moreover, in order to include electrostatic interaction,
921	>	one may need to specify the partial charges for individual atoms.
922	>	Under some circumstances, we may even need to prepare the system in
923	>	a special setup. For instance, when studying transport phenomenon in
924	>	membrane system, we may prepare the lipids in bilayer structure
925	>	instead of placing lipids randomly in solvent, since we are not
926	>	interested in self-aggregation and it takes a long time to happen.
927	>
928	>	\subsubsection{\textbf{Minimization}}
929	>
930	>	It is quite possible that some of molecules in the system from
931	>	preliminary preparation may be overlapped with each other. This
932	>	close proximity leads to high potential energy which consequently
933	>	jeopardizes any molecular dynamics simulations. To remove these
934	>	steric overlaps, one typically performs energy minimization to find
935	>	a more reasonable conformation. Several energy minimization methods
936	>	have been developed to exploit the energy surface and to locate the
937	>	local minimum. While converging slowly near the minimum, steepest
938	>	descent method is extremely robust when systems are far from
939	>	harmonic. Thus, it is often used to refine structure from
940	>	crystallographic data. Relied on the gradient or hessian, advanced
941	>	methods like conjugate gradient and Newton-Raphson converge rapidly
942	>	to a local minimum, while become unstable if the energy surface is
943	>	far from quadratic. Another factor must be taken into account, when
944	>	choosing energy minimization method, is the size of the system.
945	>	Steepest descent and conjugate gradient can deal with models of any
946	>	size. Because of the limit of computation power to calculate hessian
947	>	matrix and insufficient storage capacity to store them, most
948	>	Newton-Raphson methods can not be used with very large models.
949	>
950	>	\subsubsection{\textbf{Heating}}
951	>
952	>	Typically, Heating is performed by assigning random velocities
953	>	according to a Gaussian distribution for a temperature. Beginning at
954	>	a lower temperature and gradually increasing the temperature by
955	>	assigning greater random velocities, we end up with setting the
956	>	temperature of the system to a final temperature at which the
957	>	simulation will be conducted. In heating phase, we should also keep
958	>	the system from drifting or rotating as a whole. Equivalently, the
959	>	net linear momentum and angular momentum of the system should be
960	>	shifted to zero.
961	>
962	>	\subsubsection{\textbf{Equilibration}}
963	>
964	>	The purpose of equilibration is to allow the system to evolve
965	>	spontaneously for a period of time and reach equilibrium. The
966	>	procedure is continued until various statistical properties, such as
967	>	temperature, pressure, energy, volume and other structural
968	>	properties \textit{etc}, become independent of time. Strictly
969	>	speaking, minimization and heating are not necessary, provided the
970	>	equilibration process is long enough. However, these steps can serve
971	>	as a means to arrive at an equilibrated structure in an effective
972	>	way.
973	>
974	>	\subsection{\label{introSection:production}Production}
975	>
976	>	Production run is the most important step of the simulation, in
977	>	which the equilibrated structure is used as a starting point and the
978	>	motions of the molecules are collected for later analysis. In order
979	>	to capture the macroscopic properties of the system, the molecular
980	>	dynamics simulation must be performed in correct and efficient way.
981	>
982	>	The most expensive part of a molecular dynamics simulation is the
983	>	calculation of non-bonded forces, such as van der Waals force and
984	>	Coulombic forces \textit{etc}. For a system of $N$ particles, the
985	>	complexity of the algorithm for pair-wise interactions is $O(N^2 )$,
986	>	which making large simulations prohibitive in the absence of any
987	>	computation saving techniques.
988	>
989	>	A natural approach to avoid system size issue is to represent the
990	>	bulk behavior by a finite number of the particles. However, this
991	>	approach will suffer from the surface effect. To offset this,
992	>	\textit{Periodic boundary condition} (see Fig.~\ref{introFig:pbc})
993	>	is developed to simulate bulk properties with a relatively small
994	>	number of particles. In this method, the simulation box is
995	>	replicated throughout space to form an infinite lattice. During the
996	>	simulation, when a particle moves in the primary cell, its image in
997	>	other cells move in exactly the same direction with exactly the same
998	>	orientation. Thus, as a particle leaves the primary cell, one of its
999	>	images will enter through the opposite face.
1000	>	\begin{figure}
1001	>	\centering
1002	>	\includegraphics[width=\linewidth]{pbc.eps}
1003	>	\caption[An illustration of periodic boundary conditions]{A 2-D
1004	>	illustration of periodic boundary conditions. As one particle leaves
1005	>	the left of the simulation box, an image of it enters the right.}
1006	>	\label{introFig:pbc}
1007	>	\end{figure}
1008	>
1009	>	%cutoff and minimum image convention
1010	>	Another important technique to improve the efficiency of force
1011	>	evaluation is to apply cutoff where particles farther than a
1012	>	predetermined distance, are not included in the calculation
1013	>	\cite{Frenkel1996}. The use of a cutoff radius will cause a
1014	>	discontinuity in the potential energy curve. Fortunately, one can
1015	>	shift the potential to ensure the potential curve go smoothly to
1016	>	zero at the cutoff radius. Cutoff strategy works pretty well for
1017	>	Lennard-Jones interaction because of its short range nature.
1018	>	However, simply truncating the electrostatic interaction with the
1019	>	use of cutoff has been shown to lead to severe artifacts in
1020	>	simulations. Ewald summation, in which the slowly conditionally
1021	>	convergent Coulomb potential is transformed into direct and
1022	>	reciprocal sums with rapid and absolute convergence, has proved to
1023	>	minimize the periodicity artifacts in liquid simulations. Taking the
1024	>	advantages of the fast Fourier transform (FFT) for calculating
1025	>	discrete Fourier transforms, the particle mesh-based
1026	>	methods\cite{Hockney1981,Shimada1993, Luty1994} are accelerated from
1027	>	$O(N^{3/2})$ to $O(N logN)$. An alternative approach is \emph{fast
1028	>	multipole method}\cite{Greengard1987, Greengard1994}, which treats
1029	>	Coulombic interaction exactly at short range, and approximate the
1030	>	potential at long range through multipolar expansion. In spite of
1031	>	their wide acceptances at the molecular simulation community, these
1032	>	two methods are hard to be implemented correctly and efficiently.
1033	>	Instead, we use a damped and charge-neutralized Coulomb potential
1034	>	method developed by Wolf and his coworkers\cite{Wolf1999}. The
1035	>	shifted Coulomb potential for particle $i$ and particle $j$ at
1036	>	distance $r_{rj}$ is given by:
1037	>	\begin{equation}
1038	>	V(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha
1039	>	r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow
1040	>	R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha
1041	>	r_{ij})}{r_{ij}}\right\}. \label{introEquation:shiftedCoulomb}
1042	>	\end{equation}
1043	>	where $\alpha$ is the convergence parameter. Due to the lack of
1044	>	inherent periodicity and rapid convergence,this method is extremely
1045	>	efficient and easy to implement.
1046	>	\begin{figure}
1047	>	\centering
1048	>	\includegraphics[width=\linewidth]{shifted_coulomb.eps}
1049	>	\caption[An illustration of shifted Coulomb potential]{An
1050	>	illustration of shifted Coulomb potential.}
1051	>	\label{introFigure:shiftedCoulomb}
1052	>	\end{figure}
1053	>
1054	>	%multiple time step
1055	>
1056	>	\subsection{\label{introSection:Analysis} Analysis}
1057	>
1058	>	Recently, advanced visualization technique are widely applied to
1059	>	monitor the motions of molecules. Although the dynamics of the
1060	>	system can be described qualitatively from animation, quantitative
1061	>	trajectory analysis are more appreciable. According to the
1062	>	principles of Statistical Mechanics,
1063	>	Sec.~\ref{introSection:statisticalMechanics}, one can compute
1064	>	thermodynamics properties, analyze fluctuations of structural
1065	>	parameters, and investigate time-dependent processes of the molecule
1066	>	from the trajectories.
1067	>
1068	>	\subsubsection{\label{introSection:thermodynamicsProperties}\textbf{Thermodynamics Properties}}
1069	>
1070	>	Thermodynamics properties, which can be expressed in terms of some
1071	>	function of the coordinates and momenta of all particles in the
1072	>	system, can be directly computed from molecular dynamics. The usual
1073	>	way to measure the pressure is based on virial theorem of Clausius
1074	>	which states that the virial is equal to $-3Nk_BT$. For a system
1075	>	with forces between particles, the total virial, $W$, contains the
1076	>	contribution from external pressure and interaction between the
1077	>	particles:
1078	>	\[
1079	>	W =  - 3PV + \left\langle {\sum\limits_{i < j} {r{}_{ij} \cdot
1080	>	f_{ij} } } \right\rangle
1081		\]
1082	+	where $f_{ij}$ is the force between particle $i$ and $j$ at a
1083	+	distance $r_{ij}$. Thus, the expression for the pressure is given
1084	+	by:
1085	+	\begin{equation}
1086	+	P = \frac{{Nk_B T}}{V} - \frac{1}{{3V}}\left\langle {\sum\limits_{i
1087	+	< j} {r{}_{ij} \cdot f_{ij} } } \right\rangle
1088	+	\end{equation}
1089	+
1090	+	\subsubsection{\label{introSection:structuralProperties}\textbf{Structural Properties}}
1091	+
1092	+	Structural Properties of a simple fluid can be described by a set of
1093	+	distribution functions. Among these functions,\emph{pair
1094	+	distribution function}, also known as \emph{radial distribution
1095	+	function}, is of most fundamental importance to liquid-state theory.
1096	+	Pair distribution function can be gathered by Fourier transforming
1097	+	raw data from a series of neutron diffraction experiments and
1098	+	integrating over the surface factor \cite{Powles1973}. The
1099	+	experiment result can serve as a criterion to justify the
1100	+	correctness of the theory. Moreover, various equilibrium
1101	+	thermodynamic and structural properties can also be expressed in
1102	+	terms of radial distribution function \cite{Allen1987}.
1103	+
1104	+	A pair distribution functions $g(r)$ gives the probability that a
1105	+	particle $i$ will be located at a distance $r$ from a another
1106	+	particle $j$ in the system
1107		\[
1108	<	L(x_\alpha  ) = \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) +
1109	<	px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }}
1108	>	g(r) = \frac{V}{{N^2 }}\left\langle {\sum\limits_i {\sum\limits_{j
1109	>	\ne i} {\delta (r - r_{ij} )} } } \right\rangle.
1110		\]
1111	<	Then, the system coordinates,
1112	<	\begin{align}
1113	<	mL(\ddot x) &=  - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
1114	<	\sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{\frac{{g_\alpha
1115	<	}}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha
1116	<	(0)}}{{p^2  + \omega _\alpha ^2 }} - \frac{{g_\alpha ^2 }}{{m_\alpha
1117	<	}}\omega _\alpha ^2 L(x)} \right\}}
1118	<	%
1119	<	&= - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
1120	<	\sum\limits_{\alpha  = 1}^N {\left\{ { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}\frac{p}{{p^2  + \omega _\alpha ^2 }}pL(x)
482	<	- \frac{p}{{p^2  + \omega _\alpha ^2 }}g_\alpha  x_\alpha  (0)
483	<	- \frac{1}{{p^2  + \omega _\alpha ^2 }}g_\alpha  \dot x_\alpha  (0)} \right\}}
484	<	\end{align}
485	<	Then, the inverse transform,
1111	>	Note that the delta function can be replaced by a histogram in
1112	>	computer simulation. Figure
1113	>	\ref{introFigure:pairDistributionFunction} shows a typical pair
1114	>	distribution function for the liquid argon system. The occurrence of
1115	>	several peaks in the plot of $g(r)$ suggests that it is more likely
1116	>	to find particles at certain radial values than at others. This is a
1117	>	result of the attractive interaction at such distances. Because of
1118	>	the strong repulsive forces at short distance, the probability of
1119	>	locating particles at distances less than about 2.5{\AA} from each
1120	>	other is essentially zero.
1121
1122	<	\begin{align}
1123	<	m\ddot x &=  - \frac{{\partial W(x)}}{{\partial x}} -
1124	<	\sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
1125	<	}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\int_0^t {\cos (\omega
1126	<	_\alpha  t)\dot x(t - \tau )d\tau  - \left[ {g_\alpha  x_\alpha  (0)
1127	<	- \frac{{g_\alpha  }}{{m_\alpha  \omega _\alpha  }}} \right]\cos
1128	<	(\omega _\alpha  t) - \frac{{g_\alpha  \dot x_\alpha  (0)}}{{\omega
494	<	_\alpha  }}\sin (\omega _\alpha  t)} } \right\}}
495	<	%
496	<	&= - \frac{{\partial W(x)}}{{\partial x}} - \int_0^t
497	<	{\sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
498	<	}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha
499	<	t)\dot x(t - \tau )d} \tau }  + \sum\limits_{\alpha  = 1}^N {\left\{
500	<	{\left[ {g_\alpha  x_\alpha  (0) - \frac{{g_\alpha  }}{{m_\alpha
501	<	\omega _\alpha  }}} \right]\cos (\omega _\alpha  t) +
502	<	\frac{{g_\alpha  \dot x_\alpha  (0)}}{{\omega _\alpha  }}\sin
503	<	(\omega _\alpha  t)} \right\}}
504	<	\end{align}
1122	>	%\begin{figure}
1123	>	%\centering
1124	>	%\includegraphics[width=\linewidth]{pdf.eps}
1125	>	%\caption[Pair distribution function for the liquid argon
1126	>	%]{Pair distribution function for the liquid argon}
1127	>	%\label{introFigure:pairDistributionFunction}
1128	>	%\end{figure}
1129
1130	+	\subsubsection{\label{introSection:timeDependentProperties}\textbf{Time-dependent
1131	+	Properties}}
1132	+
1133	+	Time-dependent properties are usually calculated using \emph{time
1134	+	correlation function}, which correlates random variables $A$ and $B$
1135	+	at two different time
1136		\begin{equation}
1137	<	m\ddot x =  - \frac{{\partial W}}{{\partial x}} - \int_0^t {\xi
1138	<	(t)\dot x(t - \tau )d\tau }  + R(t)
509	<	\label{introEuqation:GeneralizedLangevinDynamics}
1137	>	C_{AB} (t) = \left\langle {A(t)B(0)} \right\rangle.
1138	>	\label{introEquation:timeCorrelationFunction}
1139		\end{equation}
1140	<	%where $ {\xi (t)}$ is friction kernel, $R(t)$ is random force and
1141	<	%$W$ is the potential of mean force. $W(x) =  - kT\ln p(x)$
1140	>	If $A$ and $B$ refer to same variable, this kind of correlation
1141	>	function is called \emph{auto correlation function}. One example of
1142	>	auto correlation function is velocity auto-correlation function
1143	>	which is directly related to transport properties of molecular
1144	>	liquids:
1145		\[
1146	<	\xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
1147	<	}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
1146	>	D = \frac{1}{3}\int\limits_0^\infty  {\left\langle {v(t) \cdot v(0)}
1147	>	\right\rangle } dt
1148		\]
1149	<	For an infinite harmonic bath, we can use the spectral density and
1150	<	an integral over frequencies.
1149	>	where $D$ is diffusion constant. Unlike velocity autocorrelation
1150	>	function which is averaging over time origins and over all the
1151	>	atoms, dipole autocorrelation are calculated for the entire system.
1152	>	The dipole autocorrelation function is given by:
1153	>	\[
1154	>	c_{dipole}  = \left\langle {u_{tot} (t) \cdot u_{tot} (t)}
1155	>	\right\rangle
1156	>	\]
1157	>	Here $u_{tot}$ is the net dipole of the entire system and is given
1158	>	by
1159	>	\[
1160	>	u_{tot} (t) = \sum\limits_i {u_i (t)}
1161	>	\]
1162	>	In principle, many time correlation functions can be related with
1163	>	Fourier transforms of the infrared, Raman, and inelastic neutron
1164	>	scattering spectra of molecular liquids. In practice, one can
1165	>	extract the IR spectrum from the intensity of dipole fluctuation at
1166	>	each frequency using the following relationship:
1167	>	\[
1168	>	\hat c_{dipole} (v) = \int_{ - \infty }^\infty  {c_{dipole} (t)e^{ -
1169	>	i2\pi vt} dt}
1170	>	\]
1171
1172	+	\section{\label{introSection:rigidBody}Dynamics of Rigid Bodies}
1173	+
1174	+	Rigid bodies are frequently involved in the modeling of different
1175	+	areas, from engineering, physics, to chemistry. For example,
1176	+	missiles and vehicle are usually modeled by rigid bodies.  The
1177	+	movement of the objects in 3D gaming engine or other physics
1178	+	simulator is governed by the rigid body dynamics. In molecular
1179	+	simulation, rigid body is used to simplify the model in
1180	+	protein-protein docking study\cite{Gray2003}.
1181	+
1182	+	It is very important to develop stable and efficient methods to
1183	+	integrate the equations of motion of orientational degrees of
1184	+	freedom. Euler angles are the nature choice to describe the
1185	+	rotational degrees of freedom. However, due to its singularity, the
1186	+	numerical integration of corresponding equations of motion is very
1187	+	inefficient and inaccurate. Although an alternative integrator using
1188	+	different sets of Euler angles can overcome this
1189	+	difficulty\cite{Barojas1973}, the computational penalty and the lost
1190	+	of angular momentum conservation still remain. A singularity free
1191	+	representation utilizing quaternions was developed by Evans in
1192	+	1977\cite{Evans1977}. Unfortunately, this approach suffer from the
1193	+	nonseparable Hamiltonian resulted from quaternion representation,
1194	+	which prevents the symplectic algorithm to be utilized. Another
1195	+	different approach is to apply holonomic constraints to the atoms
1196	+	belonging to the rigid body. Each atom moves independently under the
1197	+	normal forces deriving from potential energy and constraint forces
1198	+	which are used to guarantee the rigidness. However, due to their
1199	+	iterative nature, SHAKE and Rattle algorithm converge very slowly
1200	+	when the number of constraint increases\cite{Ryckaert1977,
1201	+	Andersen1983}.
1202	+
1203	+	The break through in geometric literature suggests that, in order to
1204	+	develop a long-term integration scheme, one should preserve the
1205	+	symplectic structure of the flow. Introducing conjugate momentum to
1206	+	rotation matrix $Q$ and re-formulating Hamiltonian's equation, a
1207	+	symplectic integrator, RSHAKE\cite{Kol1997}, was proposed to evolve
1208	+	the Hamiltonian system in a constraint manifold by iteratively
1209	+	satisfying the orthogonality constraint $Q_T Q = 1$. An alternative
1210	+	method using quaternion representation was developed by
1211	+	Omelyan\cite{Omelyan1998}. However, both of these methods are
1212	+	iterative and inefficient. In this section, we will present a
1213	+	symplectic Lie-Poisson integrator for rigid body developed by
1214	+	Dullweber and his coworkers\cite{Dullweber1997} in depth.
1215	+
1216	+	\subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Body}
1217	+	The motion of the rigid body is Hamiltonian with the Hamiltonian
1218	+	function
1219	+	\begin{equation}
1220	+	H = \frac{1}{2}(p^T m^{ - 1} p) + \frac{1}{2}tr(PJ^{ - 1} P) +
1221	+	V(q,Q) + \frac{1}{2}tr[(QQ^T  - 1)\Lambda ].
1222	+	\label{introEquation:RBHamiltonian}
1223	+	\end{equation}
1224	+	Here, $q$ and $Q$  are the position and rotation matrix for the
1225	+	rigid-body, $p$ and $P$  are conjugate momenta to $q$  and $Q$ , and
1226	+	$J$, a diagonal matrix, is defined by
1227		\[
1228	<	R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
522	<	- \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
523	<	\right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
524	<	(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)
1228	>	I_{ii}^{ - 1}  = \frac{1}{2}\sum\limits_{i \ne j} {J_{jj}^{ - 1} }
1229		\]
1230	<	The random forces depend only on initial conditions.
1230	>	where $I_{ii}$ is the diagonal element of the inertia tensor. This
1231	>	constrained Hamiltonian equation subjects to a holonomic constraint,
1232	>	\begin{equation}
1233	>	Q^T Q = 1, \label{introEquation:orthogonalConstraint}
1234	>	\end{equation}
1235	>	which is used to ensure rotation matrix's orthogonality.
1236	>	Differentiating \ref{introEquation:orthogonalConstraint} and using
1237	>	Equation \ref{introEquation:RBMotionMomentum}, one may obtain,
1238	>	\begin{equation}
1239	>	Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0 . \\
1240	>	\label{introEquation:RBFirstOrderConstraint}
1241	>	\end{equation}
1242
1243	<	\subsubsection{\label{introSection:secondFluctuationDissipation}The Second Fluctuation Dissipation Theorem}
1244	<	So we can define a new set of coordinates,
1243	>	Using Equation (\ref{introEquation:motionHamiltonianCoordinate},
1244	>	\ref{introEquation:motionHamiltonianMomentum}), one can write down
1245	>	the equations of motion,
1246	>
1247	>	\begin{eqnarray}
1248	>	\frac{{dq}}{{dt}} & = & \frac{p}{m} \label{introEquation:RBMotionPosition}\\
1249	>	\frac{{dp}}{{dt}} & = & - \nabla _q V(q,Q) \label{introEquation:RBMotionMomentum}\\
1250	>	\frac{{dQ}}{{dt}} & = & PJ^{ - 1}  \label{introEquation:RBMotionRotation}\\
1251	>	\frac{{dP}}{{dt}} & = & - \nabla _Q V(q,Q) - 2Q\Lambda . \label{introEquation:RBMotionP}
1252	>	\end{eqnarray}
1253	>
1254	>	In general, there are two ways to satisfy the holonomic constraints.
1255	>	We can use constraint force provided by lagrange multiplier on the
1256	>	normal manifold to keep the motion on constraint space. Or we can
1257	>	simply evolve the system in constraint manifold. These two methods
1258	>	are proved to be equivalent. The holonomic constraint and equations
1259	>	of motions define a constraint manifold for rigid body
1260		\[
1261	<	q_\alpha  (t) = x_\alpha  (t) - \frac{1}{{m_\alpha  \omega _\alpha
1262	<	^2 }}x(0)
1261	>	M = \left\{ {(Q,P):Q^T Q = 1,Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0}
1262	>	\right\}.
1263		\]
1264	<	This makes
1264	>
1265	>	Unfortunately, this constraint manifold is not the cotangent bundle
1266	>	$T_{\star}SO(3)$. However, it turns out that under symplectic
1267	>	transformation, the cotangent space and the phase space are
1268	>	diffeomorphic. Introducing
1269		\[
1270	<	R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}
1270	>	\tilde Q = Q,\tilde P = \frac{1}{2}\left( {P - QP^T Q} \right),
1271		\]
1272	<	And since the $q$ coordinates are harmonic oscillators,
1272	>	the mechanical system subject to a holonomic constraint manifold $M$
1273	>	can be re-formulated as a Hamiltonian system on the cotangent space
1274		\[
1275	<	\begin{array}{l}
1276	<	\left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  = \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
542	<	\left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle  = \delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
543	<	\end{array}
1275	>	T^* SO(3) = \left\{ {(\tilde Q,\tilde P):\tilde Q^T \tilde Q =
1276	>	1,\tilde Q^T \tilde PJ^{ - 1}  + J^{ - 1} P^T \tilde Q = 0} \right\}
1277		\]
1278
1279	<	\begin{align}
1280	<	\left\langle {R(t)R(0)} \right\rangle  &= \sum\limits_\alpha
1281	<	{\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha
1282	<	(t)q_\beta  (0)} \right\rangle } }
1283	<	%
1284	<	&= \sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)}
1285	<	\right\rangle \cos (\omega _\alpha  t)}
1286	<	%
1287	<	&= kT\xi (t)
1288	<	\end{align}
1279	>	For a body fixed vector $X_i$ with respect to the center of mass of
1280	>	the rigid body, its corresponding lab fixed vector $X_0^{lab}$  is
1281	>	given as
1282	>	\begin{equation}
1283	>	X_i^{lab} = Q X_i + q.
1284	>	\end{equation}
1285	>	Therefore, potential energy $V(q,Q)$ is defined by
1286	>	\[
1287	>	V(q,Q) = V(Q X_0 + q).
1288	>	\]
1289	>	Hence, the force and torque are given by
1290	>	\[
1291	>	\nabla _q V(q,Q) = F(q,Q) = \sum\limits_i {F_i (q,Q)},
1292	>	\]
1293	>	and
1294	>	\[
1295	>	\nabla _Q V(q,Q) = F(q,Q)X_i^t
1296	>	\]
1297	>	respectively.
1298
1299	+	As a common choice to describe the rotation dynamics of the rigid
1300	+	body, angular momentum on body frame $\Pi  = Q^t P$ is introduced to
1301	+	rewrite the equations of motion,
1302		\begin{equation}
1303	<	\xi (t) = \left\langle {R(t)R(0)} \right\rangle
1304	<	\label{introEquation:secondFluctuationDissipation}
1303	>	\begin{array}{l}
1304	>	\mathop \Pi \limits^ \bullet   = J^{ - 1} \Pi ^T \Pi  + Q^T \sum\limits_i {F_i (q,Q)X_i^T }  - \Lambda  \\
1305	>	\mathop Q\limits^{{\rm{   }} \bullet }  = Q\Pi {\rm{ }}J^{ - 1}  \\
1306	>	\end{array}
1307	>	\label{introEqaution:RBMotionPI}
1308		\end{equation}
1309	+	, as well as holonomic constraints,
1310	+	\[
1311	+	\begin{array}{l}
1312	+	\Pi J^{ - 1}  + J^{ - 1} \Pi ^t  = 0 \\
1313	+	Q^T Q = 1 \\
1314	+	\end{array}
1315	+	\]
1316
1317	<	\section{\label{introSection:hydroynamics}Hydrodynamics}
1318	<
1319	<	\subsection{\label{introSection:frictionTensor} Friction Tensor}
1320	<	\subsection{\label{introSection:analyticalApproach}Analytical
1321	<	Approach}
1322	<
1323	<	\subsection{\label{introSection:approximationApproach}Approximation
1324	<	Approach}
1317	>	For a vector $v(v_1 ,v_2 ,v_3 ) \in R^3$ and a matrix $\hat v \in
1318	>	so(3)^ \star$, the hat-map isomorphism,
1319	>	\begin{equation}
1320	>	v(v_1 ,v_2 ,v_3 ) \Leftrightarrow \hat v = \left(
1321	>	{\begin{array}{*{20}c}
1322	>	0 & { - v_3 } & {v_2 }  \\
1323	>	{v_3 } & 0 & { - v_1 }  \\
1324	>	{ - v_2 } & {v_1 } & 0  \\
1325	>	\end{array}} \right),
1326	>	\label{introEquation:hatmapIsomorphism}
1327	>	\end{equation}
1328	>	will let us associate the matrix products with traditional vector
1329	>	operations
1330	>	\[
1331	>	\hat vu = v \times u
1332	>	\]
1333	>	Using \ref{introEqaution:RBMotionPI}, one can construct a skew
1334	>	matrix,
1335	>	\begin{equation}
1336	>	(\mathop \Pi \limits^ \bullet   - \mathop \Pi \limits^ {\bullet  ^T}
1337	>	){\rm{ }} = {\rm{ }}(\Pi  - \Pi ^T ){\rm{ }}(J^{ - 1} \Pi  + \Pi J^{
1338	>	- 1} ) + \sum\limits_i {[Q^T F_i (r,Q)X_i^T  - X_i F_i (r,Q)^T Q]} -
1339	>	(\Lambda  - \Lambda ^T ) . \label{introEquation:skewMatrixPI}
1340	>	\end{equation}
1341	>	Since $\Lambda$ is symmetric, the last term of Equation
1342	>	\ref{introEquation:skewMatrixPI} is zero, which implies the Lagrange
1343	>	multiplier $\Lambda$ is absent from the equations of motion. This
1344	>	unique property eliminate the requirement of iterations which can
1345	>	not be avoided in other methods\cite{Kol1997, Omelyan1998}.
1346
1347	<	\subsection{\label{introSection:centersRigidBody}Centers of Rigid
1348	<	Body}
1347	>	Applying hat-map isomorphism, we obtain the equation of motion for
1348	>	angular momentum on body frame
1349	>	\begin{equation}
1350	>	\dot \pi  = \pi  \times I^{ - 1} \pi  + \sum\limits_i {\left( {Q^T
1351	>	F_i (r,Q)} \right) \times X_i }.
1352	>	\label{introEquation:bodyAngularMotion}
1353	>	\end{equation}
1354	>	In the same manner, the equation of motion for rotation matrix is
1355	>	given by
1356	>	\[
1357	>	\dot Q = Qskew(I^{ - 1} \pi )
1358	>	\]
1359	>
1360	>	\subsection{\label{introSection:SymplecticFreeRB}Symplectic
1361	>	Lie-Poisson Integrator for Free Rigid Body}
1362	>
1363	>	If there is not external forces exerted on the rigid body, the only
1364	>	contribution to the rotational is from the kinetic potential (the
1365	>	first term of \ref{introEquation:bodyAngularMotion}). The free rigid
1366	>	body is an example of Lie-Poisson system with Hamiltonian function
1367	>	\begin{equation}
1368	>	T^r (\pi ) = T_1 ^r (\pi _1 ) + T_2^r (\pi _2 ) + T_3^r (\pi _3 )
1369	>	\label{introEquation:rotationalKineticRB}
1370	>	\end{equation}
1371	>	where $T_i^r (\pi _i ) = \frac{{\pi _i ^2 }}{{2I_i }}$ and
1372	>	Lie-Poisson structure matrix,
1373	>	\begin{equation}
1374	>	J(\pi ) = \left( {\begin{array}{*{20}c}
1375	>	0 & {\pi _3 } & { - \pi _2 }  \\
1376	>	{ - \pi _3 } & 0 & {\pi _1 }  \\
1377	>	{\pi _2 } & { - \pi _1 } & 0  \\
1378	>	\end{array}} \right)
1379	>	\end{equation}
1380	>	Thus, the dynamics of free rigid body is governed by
1381	>	\begin{equation}
1382	>	\frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi )
1383	>	\end{equation}
1384	>
1385	>	One may notice that each $T_i^r$ in Equation
1386	>	\ref{introEquation:rotationalKineticRB} can be solved exactly. For
1387	>	instance, the equations of motion due to $T_1^r$ are given by
1388	>	\begin{equation}
1389	>	\frac{d}{{dt}}\pi  = R_1 \pi ,\frac{d}{{dt}}Q = QR_1
1390	>	\label{introEqaution:RBMotionSingleTerm}
1391	>	\end{equation}
1392	>	where
1393	>	\[ R_1  = \left( {\begin{array}{*{20}c}
1394	>	0 & 0 & 0  \\
1395	>	0 & 0 & {\pi _1 }  \\
1396	>	0 & { - \pi _1 } & 0  \\
1397	>	\end{array}} \right).
1398	>	\]
1399	>	The solutions of Equation \ref{introEqaution:RBMotionSingleTerm} is
1400	>	\[
1401	>	\pi (\Delta t) = e^{\Delta tR_1 } \pi (0),Q(\Delta t) =
1402	>	Q(0)e^{\Delta tR_1 }
1403	>	\]
1404	>	with
1405	>	\[
1406	>	e^{\Delta tR_1 }  = \left( {\begin{array}{*{20}c}
1407	>	0 & 0 & 0  \\
1408	>	0 & {\cos \theta _1 } & {\sin \theta _1 }  \\
1409	>	0 & { - \sin \theta _1 } & {\cos \theta _1 }  \\
1410	>	\end{array}} \right),\theta _1  = \frac{{\pi _1 }}{{I_1 }}\Delta t.
1411	>	\]
1412	>	To reduce the cost of computing expensive functions in $e^{\Delta
1413	>	tR_1 }$, we can use Cayley transformation,
1414	>	\[
1415	>	e^{\Delta tR_1 }  \approx (1 - \Delta tR_1 )^{ - 1} (1 + \Delta tR_1
1416	>	)
1417	>	\]
1418	>	The flow maps for $T_2^r$ and $T_3^r$ can be found in the same
1419	>	manner.
1420	>
1421	>	In order to construct a second-order symplectic method, we split the
1422	>	angular kinetic Hamiltonian function can into five terms
1423	>	\[
1424	>	T^r (\pi ) = \frac{1}{2}T_1 ^r (\pi _1 ) + \frac{1}{2}T_2^r (\pi _2
1425	>	) + T_3^r (\pi _3 ) + \frac{1}{2}T_2^r (\pi _2 ) + \frac{1}{2}T_1 ^r
1426	>	(\pi _1 )
1427	>	\].
1428	>	Concatenating flows corresponding to these five terms, we can obtain
1429	>	an symplectic integrator,
1430	>	\[
1431	>	\varphi _{\Delta t,T^r }  = \varphi _{\Delta t/2,\pi _1 }  \circ
1432	>	\varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }
1433	>	\circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi
1434	>	_1 }.
1435	>	\]
1436	>
1437	>	The non-canonical Lie-Poisson bracket ${F, G}$ of two function
1438	>	$F(\pi )$ and $G(\pi )$ is defined by
1439	>	\[
1440	>	\{ F,G\} (\pi ) = [\nabla _\pi  F(\pi )]^T J(\pi )\nabla _\pi  G(\pi
1441	>	)
1442	>	\]
1443	>	If the Poisson bracket of a function $F$ with an arbitrary smooth
1444	>	function $G$ is zero, $F$ is a \emph{Casimir}, which is the
1445	>	conserved quantity in Poisson system. We can easily verify that the
1446	>	norm of the angular momentum, $\parallel \pi
1447	>	\parallel$, is a \emph{Casimir}. Let$ F(\pi ) = S(\frac{{\parallel
1448	>	\pi \parallel ^2 }}{2})$ for an arbitrary function $ S:R \to R$ ,
1449	>	then by the chain rule
1450	>	\[
1451	>	\nabla _\pi  F(\pi ) = S'(\frac{{\parallel \pi \parallel ^2
1452	>	}}{2})\pi
1453	>	\]
1454	>	Thus $ [\nabla _\pi  F(\pi )]^T J(\pi ) =  - S'(\frac{{\parallel \pi
1455	>	\parallel ^2 }}{2})\pi  \times \pi  = 0 $. This explicit
1456	>	Lie-Poisson integrator is found to be extremely efficient and stable
1457	>	which can be explained by the fact the small angle approximation is
1458	>	used and the norm of the angular momentum is conserved.
1459	>
1460	>	\subsection{\label{introSection:RBHamiltonianSplitting} Hamiltonian
1461	>	Splitting for Rigid Body}
1462	>
1463	>	The Hamiltonian of rigid body can be separated in terms of kinetic
1464	>	energy and potential energy,
1465	>	\[
1466	>	H = T(p,\pi ) + V(q,Q)
1467	>	\]
1468	>	The equations of motion corresponding to potential energy and
1469	>	kinetic energy are listed in the below table,
1470	>	\begin{table}
1471	>	\caption{Equations of motion due to Potential and Kinetic Energies}
1472	>	\begin{center}
1473	>	\begin{tabular}{|l|l|}
1474	>	\hline
1475	>	% after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
1476	>	Potential & Kinetic \\
1477	>	$\frac{{dq}}{{dt}} = \frac{p}{m}$ & $\frac{d}{{dt}}q = p$ \\
1478	>	$\frac{d}{{dt}}p =  - \frac{{\partial V}}{{\partial q}}$ & $ \frac{d}{{dt}}p = 0$ \\
1479	>	$\frac{d}{{dt}}Q = 0$ & $ \frac{d}{{dt}}Q = Qskew(I^{ - 1} j)$ \\
1480	>	$ \frac{d}{{dt}}\pi  = \sum\limits_i {\left( {Q^T F_i (r,Q)} \right) \times X_i }$ & $\frac{d}{{dt}}\pi  = \pi  \times I^{ - 1} \pi$\\
1481	>	\hline
1482	>	\end{tabular}
1483	>	\end{center}
1484	>	\end{table}
1485	>	A second-order symplectic method is now obtained by the
1486	>	composition of the flow maps,
1487	>	\[
1488	>	\varphi _{\Delta t}  = \varphi _{\Delta t/2,V}  \circ \varphi
1489	>	_{\Delta t,T}  \circ \varphi _{\Delta t/2,V}.
1490	>	\]
1491	>	Moreover, $\varphi _{\Delta t/2,V}$ can be divided into two
1492	>	sub-flows which corresponding to force and torque respectively,
1493	>	\[
1494	>	\varphi _{\Delta t/2,V}  = \varphi _{\Delta t/2,F}  \circ \varphi
1495	>	_{\Delta t/2,\tau }.
1496	>	\]
1497	>	Since the associated operators of $\varphi _{\Delta t/2,F} $ and
1498	>	$\circ \varphi _{\Delta t/2,\tau }$ are commuted, the composition
1499	>	order inside $\varphi _{\Delta t/2,V}$ does not matter.
1500	>
1501	>	Furthermore, kinetic potential can be separated to translational
1502	>	kinetic term, $T^t (p)$, and rotational kinetic term, $T^r (\pi )$,
1503	>	\begin{equation}
1504	>	T(p,\pi ) =T^t (p) + T^r (\pi ).
1505	>	\end{equation}
1506	>	where $ T^t (p) = \frac{1}{2}p^T m^{ - 1} p $ and $T^r (\pi )$ is
1507	>	defined by \ref{introEquation:rotationalKineticRB}. Therefore, the
1508	>	corresponding flow maps are given by
1509	>	\[
1510	>	\varphi _{\Delta t,T}  = \varphi _{\Delta t,T^t }  \circ \varphi
1511	>	_{\Delta t,T^r }.
1512	>	\]
1513	>	Finally, we obtain the overall symplectic flow maps for free moving
1514	>	rigid body
1515	>	\begin{equation}
1516	>	\begin{array}{c}
1517	>	\varphi _{\Delta t}  = \varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau }  \\
1518	>	\circ \varphi _{\Delta t,T^t }  \circ \varphi _{\Delta t/2,\pi _1 }  \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }  \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi _1 }  \\
1519	>	\circ \varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .\\
1520	>	\end{array}
1521	>	\label{introEquation:overallRBFlowMaps}
1522	>	\end{equation}
1523	>
1524	>	\section{\label{introSection:langevinDynamics}Langevin Dynamics}
1525	>	As an alternative to newtonian dynamics, Langevin dynamics, which
1526	>	mimics a simple heat bath with stochastic and dissipative forces,
1527	>	has been applied in a variety of studies. This section will review
1528	>	the theory of Langevin dynamics simulation. A brief derivation of
1529	>	generalized Langevin equation will be given first. Follow that, we
1530	>	will discuss the physical meaning of the terms appearing in the
1531	>	equation as well as the calculation of friction tensor from
1532	>	hydrodynamics theory.
1533	>
1534	>	\subsection{\label{introSection:generalizedLangevinDynamics}Derivation of Generalized Langevin Equation}
1535	>
1536	>	Harmonic bath model, in which an effective set of harmonic
1537	>	oscillators are used to mimic the effect of a linearly responding
1538	>	environment, has been widely used in quantum chemistry and
1539	>	statistical mechanics. One of the successful applications of
1540	>	Harmonic bath model is the derivation of Deriving Generalized
1541	>	Langevin Dynamics. Lets consider a system, in which the degree of
1542	>	freedom $x$ is assumed to couple to the bath linearly, giving a
1543	>	Hamiltonian of the form
1544	>	\begin{equation}
1545	>	H = \frac{{p^2 }}{{2m}} + U(x) + H_B  + \Delta U(x,x_1 , \ldots x_N)
1546	>	\label{introEquation:bathGLE}.
1547	>	\end{equation}
1548	>	Here $p$ is a momentum conjugate to $q$, $m$ is the mass associated
1549	>	with this degree of freedom, $H_B$ is harmonic bath Hamiltonian,
1550	>	\[
1551	>	H_B  = \sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
1552	>	}}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  \omega _\alpha ^2 }
1553	>	\right\}}
1554	>	\]
1555	>	where the index $\alpha$ runs over all the bath degrees of freedom,
1556	>	$\omega _\alpha$ are the harmonic bath frequencies, $m_\alpha$ are
1557	>	the harmonic bath masses, and $\Delta U$ is bilinear system-bath
1558	>	coupling,
1559	>	\[
1560	>	\Delta U =  - \sum\limits_{\alpha  = 1}^N {g_\alpha  x_\alpha  x}
1561	>	\]
1562	>	where $g_\alpha$ are the coupling constants between the bath and the
1563	>	coordinate $x$. Introducing
1564	>	\[
1565	>	W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
1566	>	}}{{2m_\alpha  w_\alpha ^2 }}} x^2
1567	>	\] and combining the last two terms in Equation
1568	>	\ref{introEquation:bathGLE}, we may rewrite the Harmonic bath
1569	>	Hamiltonian as
1570	>	\[
1571	>	H = \frac{{p^2 }}{{2m}} + W(x) + \sum\limits_{\alpha  = 1}^N
1572	>	{\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
1573	>	w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
1574	>	w_\alpha ^2 }}x} \right)^2 } \right\}}
1575	>	\]
1576	>	Since the first two terms of the new Hamiltonian depend only on the
1577	>	system coordinates, we can get the equations of motion for
1578	>	Generalized Langevin Dynamics by Hamilton's equations
1579	>	\ref{introEquation:motionHamiltonianCoordinate,
1580	>	introEquation:motionHamiltonianMomentum},
1581	>	\begin{equation}
1582	>	m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} -
1583	>	\sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   -
1584	>	\frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right)},
1585	>	\label{introEquation:coorMotionGLE}
1586	>	\end{equation}
1587	>	and
1588	>	\begin{equation}
1589	>	m\ddot x_\alpha   =  - m_\alpha  w_\alpha ^2 \left( {x_\alpha   -
1590	>	\frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right).
1591	>	\label{introEquation:bathMotionGLE}
1592	>	\end{equation}
1593	>
1594	>	In order to derive an equation for $x$, the dynamics of the bath
1595	>	variables $x_\alpha$ must be solved exactly first. As an integral
1596	>	transform which is particularly useful in solving linear ordinary
1597	>	differential equations, Laplace transform is the appropriate tool to
1598	>	solve this problem. The basic idea is to transform the difficult
1599	>	differential equations into simple algebra problems which can be
1600	>	solved easily. Then applying inverse Laplace transform, also known
1601	>	as the Bromwich integral, we can retrieve the solutions of the
1602	>	original problems.
1603	>
1604	>	Let $f(t)$ be a function defined on $ [0,\infty ) $. The Laplace
1605	>	transform of f(t) is a new function defined as
1606	>	\[
1607	>	L(f(t)) \equiv F(p) = \int_0^\infty  {f(t)e^{ - pt} dt}
1608	>	\]
1609	>	where  $p$ is real and  $L$ is called the Laplace Transform
1610	>	Operator. Below are some important properties of Laplace transform
1611	>
1612	>	\begin{eqnarray*}
1613	>	L(x + y)  & = & L(x) + L(y) \\
1614	>	L(ax)     & = & aL(x) \\
1615	>	L(\dot x) & = & pL(x) - px(0) \\
1616	>	L(\ddot x)& = & p^2 L(x) - px(0) - \dot x(0) \\
1617	>	L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right)& = & G(p)H(p) \\
1618	>	\end{eqnarray*}
1619	>
1620	>
1621	>	Applying Laplace transform to the bath coordinates, we obtain
1622	>	\begin{eqnarray*}
1623	>	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) \\
1624	>	L(x_\alpha  ) & = & \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }} \\
1625	>	\end{eqnarray*}
1626	>
1627	>	By the same way, the system coordinates become
1628	>	\begin{eqnarray*}
1629	>	mL(\ddot x) & = & - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} \\
1630	>	& & \mbox{} - \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\}}  \\
1631	>	\end{eqnarray*}
1632	>
1633	>	With the help of some relatively important inverse Laplace
1634	>	transformations:
1635	>	\[
1636	>	\begin{array}{c}
1637	>	L(\cos at) = \frac{p}{{p^2  + a^2 }} \\
1638	>	L(\sin at) = \frac{a}{{p^2  + a^2 }} \\
1639	>	L(1) = \frac{1}{p} \\
1640	>	\end{array}
1641	>	\]
1642	>	, we obtain
1643	>	\begin{eqnarray*}
1644	>	m\ddot x & =  & - \frac{{\partial W(x)}}{{\partial x}} -
1645	>	\sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
1646	>	}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\int_0^t {\cos (\omega
1647	>	_\alpha  t)\dot x(t - \tau )d\tau } } \right\}}  \\
1648	>	& & + \sum\limits_{\alpha  = 1}^N {\left\{ {\left[ {g_\alpha
1649	>	x_\alpha (0) - \frac{{g_\alpha  }}{{m_\alpha  \omega _\alpha  }}}
1650	>	\right]\cos (\omega _\alpha  t) + \frac{{g_\alpha  \dot x_\alpha
1651	>	(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)} \right\}}
1652	>	\end{eqnarray*}
1653	>	\begin{eqnarray*}
1654	>	m\ddot x & = & - \frac{{\partial W(x)}}{{\partial x}} - \int_0^t
1655	>	{\sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
1656	>	}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha
1657	>	t)\dot x(t - \tau )d} \tau }  \\
1658	>	& & + \sum\limits_{\alpha  = 1}^N {\left\{ {\left[ {g_\alpha
1659	>	x_\alpha (0) - \frac{{g_\alpha }}{{m_\alpha \omega _\alpha  }}}
1660	>	\right]\cos (\omega _\alpha  t) + \frac{{g_\alpha  \dot x_\alpha
1661	>	(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)} \right\}}
1662	>	\end{eqnarray*}
1663	>	Introducing a \emph{dynamic friction kernel}
1664	>	\begin{equation}
1665	>	\xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
1666	>	}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
1667	>	\label{introEquation:dynamicFrictionKernelDefinition}
1668	>	\end{equation}
1669	>	and \emph{a random force}
1670	>	\begin{equation}
1671	>	R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
1672	>	- \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
1673	>	\right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
1674	>	(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t),
1675	>	\label{introEquation:randomForceDefinition}
1676	>	\end{equation}
1677	>	the equation of motion can be rewritten as
1678	>	\begin{equation}
1679	>	m\ddot x =  - \frac{{\partial W}}{{\partial x}} - \int_0^t {\xi
1680	>	(t)\dot x(t - \tau )d\tau }  + R(t)
1681	>	\label{introEuqation:GeneralizedLangevinDynamics}
1682	>	\end{equation}
1683	>	which is known as the \emph{generalized Langevin equation}.
1684	>
1685	>	\subsubsection{\label{introSection:randomForceDynamicFrictionKernel}\textbf{Random Force and Dynamic Friction Kernel}}
1686	>
1687	>	One may notice that $R(t)$ depends only on initial conditions, which
1688	>	implies it is completely deterministic within the context of a
1689	>	harmonic bath. However, it is easy to verify that $R(t)$ is totally
1690	>	uncorrelated to $x$ and $\dot x$,
1691	>	\[
1692	>	\begin{array}{l}
1693	>	\left\langle {x(t)R(t)} \right\rangle  = 0, \\
1694	>	\left\langle {\dot x(t)R(t)} \right\rangle  = 0. \\
1695	>	\end{array}
1696	>	\]
1697	>	This property is what we expect from a truly random process. As long
1698	>	as the model, which is gaussian distribution in general, chosen for
1699	>	$R(t)$ is a truly random process, the stochastic nature of the GLE
1700	>	still remains.
1701	>
1702	>	%dynamic friction kernel
1703	>	The convolution integral
1704	>	\[
1705	>	\int_0^t {\xi (t)\dot x(t - \tau )d\tau }
1706	>	\]
1707	>	depends on the entire history of the evolution of $x$, which implies
1708	>	that the bath retains memory of previous motions. In other words,
1709	>	the bath requires a finite time to respond to change in the motion
1710	>	of the system. For a sluggish bath which responds slowly to changes
1711	>	in the system coordinate, we may regard $\xi(t)$ as a constant
1712	>	$\xi(t) = \Xi_0$. Hence, the convolution integral becomes
1713	>	\[
1714	>	\int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = \xi _0 (x(t) - x(0))
1715	>	\]
1716	>	and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
1717	>	\[
1718	>	m\ddot x =  - \frac{\partial }{{\partial x}}\left( {W(x) +
1719	>	\frac{1}{2}\xi _0 (x - x_0 )^2 } \right) + R(t),
1720	>	\]
1721	>	which can be used to describe dynamic caging effect. The other
1722	>	extreme is the bath that responds infinitely quickly to motions in
1723	>	the system. Thus, $\xi (t)$ can be taken as a $delta$ function in
1724	>	time:
1725	>	\[
1726	>	\xi (t) = 2\xi _0 \delta (t)
1727	>	\]
1728	>	Hence, the convolution integral becomes
1729	>	\[
1730	>	\int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = 2\xi _0 \int_0^t
1731	>	{\delta (t)\dot x(t - \tau )d\tau }  = \xi _0 \dot x(t),
1732	>	\]
1733	>	and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
1734	>	\begin{equation}
1735	>	m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} - \xi _0 \dot
1736	>	x(t) + R(t) \label{introEquation:LangevinEquation}
1737	>	\end{equation}
1738	>	which is known as the Langevin equation. The static friction
1739	>	coefficient $\xi _0$ can either be calculated from spectral density
1740	>	or be determined by Stokes' law for regular shaped particles. A
1741	>	briefly review on calculating friction tensor for arbitrary shaped
1742	>	particles is given in Sec.~\ref{introSection:frictionTensor}.
1743	>
1744	>	\subsubsection{\label{introSection:secondFluctuationDissipation}\textbf{The Second Fluctuation Dissipation Theorem}}
1745	>
1746	>	Defining a new set of coordinates,
1747	>	\[
1748	>	q_\alpha  (t) = x_\alpha  (t) - \frac{1}{{m_\alpha  \omega _\alpha
1749	>	^2 }}x(0)
1750	>	\],
1751	>	we can rewrite $R(T)$ as
1752	>	\[
1753	>	R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}.
1754	>	\]
1755	>	And since the $q$ coordinates are harmonic oscillators,
1756	>
1757	>	\begin{eqnarray*}
1758	>	\left\langle {q_\alpha ^2 } \right\rangle  & = & \frac{{kT}}{{m_\alpha  \omega _\alpha ^2 }} \\
1759	>	\left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle & = & \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
1760	>	\left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle & = &\delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
1761	>	\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 } }  \\
1762	>	& = &\sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t)}  \\
1763	>	& = &kT\xi (t) \\
1764	>	\end{eqnarray*}
1765	>
1766	>	Thus, we recover the \emph{second fluctuation dissipation theorem}
1767	>	\begin{equation}
1768	>	\xi (t) = \left\langle {R(t)R(0)} \right\rangle
1769	>	\label{introEquation:secondFluctuationDissipation}.
1770	>	\end{equation}
1771	>	In effect, it acts as a constraint on the possible ways in which one
1772	>	can model the random force and friction kernel.

Diff Legend

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