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 2699 by tim, Mon Apr 10 05:35:55 2006 UTC vs.
Revision 2909 by tim, Thu Jun 29 23:00:35 2006 UTC

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

Diff Legend

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