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

Comparing trunk/tengDissertation/Introduction.tex (file contents):
Revision 2692 by tim, Tue Apr 4 21:32:58 2006 UTC vs.
Revision 2719 by tim, Tue Apr 18 22:38:19 2006 UTC

# Line 1 | Line 1
1   \chapter{\label{chapt:introduction}INTRODUCTION AND THEORETICAL BACKGROUND}
2  
3 < \section{\label{introSection:molecularDynamics}Molecular Dynamics}
3 > \section{\label{introSection:classicalMechanics}Classical
4 > Mechanics}
5  
5 As a special discipline of molecular modeling, Molecular dynamics
6 has proven to be a powerful tool for studying the functions of
7 biological systems, providing structural, thermodynamic and
8 dynamical information.
9
10 \subsection{\label{introSection:classicalMechanics}Classical Mechanics}
11
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
# Line 20 | Line 14 | sufficient to predict the future behavior of the syste
14   when further combine with the laws of mechanics will also be
15   sufficient to predict the future behavior of the system.
16  
17 < \subsubsection{\label{introSection:newtonian}Newtonian Mechanics}
17 > \subsection{\label{introSection:newtonian}Newtonian Mechanics}
18 > The discovery of Newton's three laws of mechanics which govern the
19 > motion of particles is the foundation of the classical mechanics.
20 > Newton¡¯s first law defines a class of inertial frames. Inertial
21 > frames are reference frames where a particle not interacting with
22 > other bodies will move with constant speed in the same direction.
23 > With respect to inertial frames Newton¡¯s second law has the form
24 > \begin{equation}
25 > F = \frac {dp}{dt} = \frac {mv}{dt}
26 > \label{introEquation:newtonSecondLaw}
27 > \end{equation}
28 > A point mass interacting with other bodies moves with the
29 > acceleration along the direction of the force acting on it. Let
30 > $F_{ij}$ be the force that particle $i$ exerts on particle $j$, and
31 > $F_{ji}$ be the force that particle $j$ exerts on particle $i$.
32 > Newton¡¯s third law states that
33 > \begin{equation}
34 > F_{ij} = -F_{ji}
35 > \label{introEquation:newtonThirdLaw}
36 > \end{equation}
37  
38 < \subsubsection{\label{introSection:lagrangian}Lagrangian Mechanics}
38 > Conservation laws of Newtonian Mechanics play very important roles
39 > in solving mechanics problems. The linear momentum of a particle is
40 > conserved if it is free or it experiences no force. The second
41 > conservation theorem concerns the angular momentum of a particle.
42 > The angular momentum $L$ of a particle with respect to an origin
43 > from which $r$ is measured is defined to be
44 > \begin{equation}
45 > L \equiv r \times p \label{introEquation:angularMomentumDefinition}
46 > \end{equation}
47 > The torque $\tau$ with respect to the same origin is defined to be
48 > \begin{equation}
49 > N \equiv r \times F \label{introEquation:torqueDefinition}
50 > \end{equation}
51 > Differentiating Eq.~\ref{introEquation:angularMomentumDefinition},
52 > \[
53 > \dot L = \frac{d}{{dt}}(r \times p) = (\dot r \times p) + (r \times
54 > \dot p)
55 > \]
56 > since
57 > \[
58 > \dot r \times p = \dot r \times mv = m\dot r \times \dot r \equiv 0
59 > \]
60 > thus,
61 > \begin{equation}
62 > \dot L = r \times \dot p = N
63 > \end{equation}
64 > If there are no external torques acting on a body, the angular
65 > momentum of it is conserved. The last conservation theorem state
66 > that if all forces are conservative, Energy
67 > \begin{equation}E = T + V \label{introEquation:energyConservation}
68 > \end{equation}
69 > is conserved. All of these conserved quantities are
70 > important factors to determine the quality of numerical integration
71 > scheme for rigid body \cite{Dullweber1997}.
72  
73 + \subsection{\label{introSection:lagrangian}Lagrangian Mechanics}
74 +
75   Newtonian Mechanics suffers from two important limitations: it
76   describes their motion in special cartesian coordinate systems.
77   Another limitation of Newtonian mechanics becomes obvious when we
# Line 35 | Line 83 | system, alternative procedures may be developed.
83   which arise in attempts to apply Newton's equation to complex
84   system, alternative procedures may be developed.
85  
86 < \subsubsubsection{\label{introSection:halmiltonPrinciple}Hamilton's
86 > \subsubsection{\label{introSection:halmiltonPrinciple}Hamilton's
87   Principle}
88  
89   Hamilton introduced the dynamical principle upon which it is
# Line 45 | Line 93 | the kinetic, $K$, and potential energies, $U$.
93   The actual trajectory, along which a dynamical system may move from
94   one point to another within a specified time, is derived by finding
95   the path which minimizes the time integral of the difference between
96 < the kinetic, $K$, and potential energies, $U$.
96 > the kinetic, $K$, and potential energies, $U$ \cite{tolman79}.
97   \begin{equation}
98   \delta \int_{t_1 }^{t_2 } {(K - U)dt = 0} ,
99 < \lable{introEquation:halmitonianPrinciple1}
99 > \label{introEquation:halmitonianPrinciple1}
100   \end{equation}
101  
102   For simple mechanical systems, where the forces acting on the
# Line 62 | Line 110 | then Eq.~\ref{introEquation:halmitonianPrinciple1} bec
110   \end{equation}
111   then Eq.~\ref{introEquation:halmitonianPrinciple1} becomes
112   \begin{equation}
113 < \delta \int_{t_1 }^{t_2 } {K dt = 0} ,
114 < \lable{introEquation:halmitonianPrinciple2}
113 > \delta \int_{t_1 }^{t_2 } {L dt = 0} ,
114 > \label{introEquation:halmitonianPrinciple2}
115   \end{equation}
116  
117 < \subsubsubsection{\label{introSection:equationOfMotionLagrangian}The
117 > \subsubsection{\label{introSection:equationOfMotionLagrangian}The
118   Equations of Motion in Lagrangian Mechanics}
119  
120 < for a holonomic system of $f$ degrees of freedom, the equations of
120 > For a holonomic system of $f$ degrees of freedom, the equations of
121   motion in the Lagrangian form is
122   \begin{equation}
123   \frac{d}{{dt}}\frac{{\partial L}}{{\partial \dot q_i }} -
124   \frac{{\partial L}}{{\partial q_i }} = 0,{\rm{ }}i = 1, \ldots,f
125 < \lable{introEquation:eqMotionLagrangian}
125 > \label{introEquation:eqMotionLagrangian}
126   \end{equation}
127   where $q_{i}$ is generalized coordinate and $\dot{q_{i}}$ is
128   generalized velocity.
129  
130 < \subsubsection{\label{introSection:hamiltonian}Hamiltonian Mechanics}
130 > \subsection{\label{introSection:hamiltonian}Hamiltonian Mechanics}
131  
132   Arising from Lagrangian Mechanics, Hamiltonian Mechanics was
133   introduced by William Rowan Hamilton in 1833 as a re-formulation of
# Line 90 | Line 138 | With the help of these momenta, we may now define a ne
138   p_i = \frac{\partial L}{\partial \dot q_i}
139   \label{introEquation:generalizedMomenta}
140   \end{equation}
141 < With the help of these momenta, we may now define a new quantity $H$
94 < by the equation
141 > The Lagrange equations of motion are then expressed by
142   \begin{equation}
143 < H = p_1 \dot q_1  +  \ldots  + p_f \dot q_f  - L,
143 > p_i  = \frac{{\partial L}}{{\partial q_i }}
144 > \label{introEquation:generalizedMomentaDot}
145 > \end{equation}
146 >
147 > With the help of the generalized momenta, we may now define a new
148 > quantity $H$ by the equation
149 > \begin{equation}
150 > H = \sum\limits_k {p_k \dot q_k }  - L ,
151   \label{introEquation:hamiltonianDefByLagrangian}
152   \end{equation}
153   where $ \dot q_1  \ldots \dot q_f $ are generalized velocities and
154   $L$ is the Lagrangian function for the system.
155  
156 + Differentiating Eq.~\ref{introEquation:hamiltonianDefByLagrangian},
157 + one can obtain
158 + \begin{equation}
159 + dH = \sum\limits_k {\left( {p_k d\dot q_k  + \dot q_k dp_k  -
160 + \frac{{\partial L}}{{\partial q_k }}dq_k  - \frac{{\partial
161 + L}}{{\partial \dot q_k }}d\dot q_k } \right)}  - \frac{{\partial
162 + L}}{{\partial t}}dt \label{introEquation:diffHamiltonian1}
163 + \end{equation}
164 + Making use of  Eq.~\ref{introEquation:generalizedMomenta}, the
165 + second and fourth terms in the parentheses cancel. Therefore,
166 + Eq.~\ref{introEquation:diffHamiltonian1} can be rewritten as
167 + \begin{equation}
168 + dH = \sum\limits_k {\left( {\dot q_k dp_k  - \dot p_k dq_k }
169 + \right)}  - \frac{{\partial L}}{{\partial t}}dt
170 + \label{introEquation:diffHamiltonian2}
171 + \end{equation}
172 + By identifying the coefficients of $dq_k$, $dp_k$ and dt, we can
173 + find
174 + \begin{equation}
175 + \frac{{\partial H}}{{\partial p_k }} = q_k
176 + \label{introEquation:motionHamiltonianCoordinate}
177 + \end{equation}
178 + \begin{equation}
179 + \frac{{\partial H}}{{\partial q_k }} =  - p_k
180 + \label{introEquation:motionHamiltonianMomentum}
181 + \end{equation}
182 + and
183 + \begin{equation}
184 + \frac{{\partial H}}{{\partial t}} =  - \frac{{\partial L}}{{\partial
185 + t}}
186 + \label{introEquation:motionHamiltonianTime}
187 + \end{equation}
188 +
189 + Eq.~\ref{introEquation:motionHamiltonianCoordinate} and
190 + Eq.~\ref{introEquation:motionHamiltonianMomentum} are Hamilton's
191 + equation of motion. Due to their symmetrical formula, they are also
192 + known as the canonical equations of motions \cite{Goldstein01}.
193 +
194   An important difference between Lagrangian approach and the
195   Hamiltonian approach is that the Lagrangian is considered to be a
196   function of the generalized velocities $\dot q_i$ and the
# Line 108 | Line 200 | equations.
200   appropriate for application to statistical mechanics and quantum
201   mechanics, since it treats the coordinate and its time derivative as
202   independent variables and it only works with 1st-order differential
203 < equations.
203 > equations\cite{Marion90}.
204  
205 + In Newtonian Mechanics, a system described by conservative forces
206 + conserves the total energy \ref{introEquation:energyConservation}.
207 + It follows that Hamilton's equations of motion conserve the total
208 + Hamiltonian.
209 + \begin{equation}
210 + \frac{{dH}}{{dt}} = \sum\limits_i {\left( {\frac{{\partial
211 + H}}{{\partial q_i }}\dot q_i  + \frac{{\partial H}}{{\partial p_i
212 + }}\dot p_i } \right)}  = \sum\limits_i {\left( {\frac{{\partial
213 + H}}{{\partial q_i }}\frac{{\partial H}}{{\partial p_i }} -
214 + \frac{{\partial H}}{{\partial p_i }}\frac{{\partial H}}{{\partial
215 + q_i }}} \right) = 0} \label{introEquation:conserveHalmitonian}
216 + \end{equation}
217  
218 < \subsubsection{\label{introSection:canonicalTransformation}Canonical Transformation}
218 > \section{\label{introSection:statisticalMechanics}Statistical
219 > Mechanics}
220  
221 < \subsection{\label{introSection:statisticalMechanics}Statistical Mechanics}
117 <
118 < The thermodynamic behaviors and properties  of Molecular Dynamics
221 > The thermodynamic behaviors and properties of Molecular Dynamics
222   simulation are governed by the principle of Statistical Mechanics.
223   The following section will give a brief introduction to some of the
224 < Statistical Mechanics concepts presented in this dissertation.
224 > Statistical Mechanics concepts and theorem presented in this
225 > dissertation.
226  
227 < \subsubsection{\label{introSection::ensemble}Ensemble}
227 > \subsection{\label{introSection:ensemble}Phase Space and Ensemble}
228  
229 < \subsubsection{\label{introSection:ergodic}The Ergodic Hypothesis}
229 > Mathematically, phase space is the space which represents all
230 > possible states. Each possible state of the system corresponds to
231 > one unique point in the phase space. For mechanical systems, the
232 > phase space usually consists of all possible values of position and
233 > momentum variables. Consider a dynamic system in a cartesian space,
234 > where each of the $6f$ coordinates and momenta is assigned to one of
235 > $6f$ mutually orthogonal axes, the phase space of this system is a
236 > $6f$ dimensional space. A point, $x = (q_1 , \ldots ,q_f ,p_1 ,
237 > \ldots ,p_f )$, with a unique set of values of $6f$ coordinates and
238 > momenta is a phase space vector.
239  
240 < \subsection{\label{introSection:rigidBody}Dynamics of Rigid Bodies}
240 > A microscopic state or microstate of a classical system is
241 > specification of the complete phase space vector of a system at any
242 > instant in time. An ensemble is defined as a collection of systems
243 > sharing one or more macroscopic characteristics but each being in a
244 > unique microstate. The complete ensemble is specified by giving all
245 > systems or microstates consistent with the common macroscopic
246 > characteristics of the ensemble. Although the state of each
247 > individual system in the ensemble could be precisely described at
248 > any instance in time by a suitable phase space vector, when using
249 > ensembles for statistical purposes, there is no need to maintain
250 > distinctions between individual systems, since the numbers of
251 > systems at any time in the different states which correspond to
252 > different regions of the phase space are more interesting. Moreover,
253 > in the point of view of statistical mechanics, one would prefer to
254 > use ensembles containing a large enough population of separate
255 > members so that the numbers of systems in such different states can
256 > be regarded as changing continuously as we traverse different
257 > regions of the phase space. The condition of an ensemble at any time
258 > can be regarded as appropriately specified by the density $\rho$
259 > with which representative points are distributed over the phase
260 > space. The density of distribution for an ensemble with $f$ degrees
261 > of freedom is defined as,
262 > \begin{equation}
263 > \rho  = \rho (q_1 , \ldots ,q_f ,p_1 , \ldots ,p_f ,t).
264 > \label{introEquation:densityDistribution}
265 > \end{equation}
266 > Governed by the principles of mechanics, the phase points change
267 > their value which would change the density at any time at phase
268 > space. Hence, the density of distribution is also to be taken as a
269 > function of the time.
270  
271 < \subsection{\label{introSection:correlationFunctions}Correlation Functions}
272 <
273 < \section{\label{introSection:langevinDynamics}Langevin Dynamics}
271 > The number of systems $\delta N$ at time $t$ can be determined by,
272 > \begin{equation}
273 > \delta N = \rho (q,p,t)dq_1  \ldots dq_f dp_1  \ldots dp_f.
274 > \label{introEquation:deltaN}
275 > \end{equation}
276 > Assuming a large enough population of systems are exploited, we can
277 > sufficiently approximate $\delta N$ without introducing
278 > discontinuity when we go from one region in the phase space to
279 > another. By integrating over the whole phase space,
280 > \begin{equation}
281 > N = \int { \ldots \int {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f
282 > \label{introEquation:totalNumberSystem}
283 > \end{equation}
284 > gives us an expression for the total number of the systems. Hence,
285 > the probability per unit in the phase space can be obtained by,
286 > \begin{equation}
287 > \frac{{\rho (q,p,t)}}{N} = \frac{{\rho (q,p,t)}}{{\int { \ldots \int
288 > {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}.
289 > \label{introEquation:unitProbability}
290 > \end{equation}
291 > With the help of Equation(\ref{introEquation:unitProbability}) and
292 > the knowledge of the system, it is possible to calculate the average
293 > value of any desired quantity which depends on the coordinates and
294 > momenta of the system. Even when the dynamics of the real system is
295 > complex, or stochastic, or even discontinuous, the average
296 > properties of the ensemble of possibilities as a whole may still
297 > remain well defined. For a classical system in thermal equilibrium
298 > with its environment, the ensemble average of a mechanical quantity,
299 > $\langle A(q , p) \rangle_t$, takes the form of an integral over the
300 > phase space of the system,
301 > \begin{equation}
302 > \langle  A(q , p) \rangle_t = \frac{{\int { \ldots \int {A(q,p)\rho
303 > (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}{{\int { \ldots \int {\rho
304 > (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}
305 > \label{introEquation:ensembelAverage}
306 > \end{equation}
307  
308 < \subsection{\label{introSection:generalizedLangevinDynamics}Generalized Langevin Dynamics}
308 > There are several different types of ensembles with different
309 > statistical characteristics. As a function of macroscopic
310 > parameters, such as temperature \textit{etc}, partition function can
311 > be used to describe the statistical properties of a system in
312 > thermodynamic equilibrium.
313  
314 < \subsection{\label{introSection:hydroynamics}Hydrodynamics}
314 > As an ensemble of systems, each of which is known to be thermally
315 > isolated and conserve energy, Microcanonical ensemble(NVE) has a
316 > partition function like,
317 > \begin{equation}
318 > \Omega (N,V,E) = e^{\beta TS} \label{introEquation:NVEPartition}.
319 > \end{equation}
320 > A canonical ensemble(NVT)is an ensemble of systems, each of which
321 > can share its energy with a large heat reservoir. The distribution
322 > of the total energy amongst the possible dynamical states is given
323 > by the partition function,
324 > \begin{equation}
325 > \Omega (N,V,T) = e^{ - \beta A}
326 > \label{introEquation:NVTPartition}
327 > \end{equation}
328 > Here, $A$ is the Helmholtz free energy which is defined as $ A = U -
329 > TS$. Since most experiment are carried out under constant pressure
330 > condition, isothermal-isobaric ensemble(NPT) play a very important
331 > role in molecular simulation. The isothermal-isobaric ensemble allow
332 > the system to exchange energy with a heat bath of temperature $T$
333 > and to change the volume as well. Its partition function is given as
334 > \begin{equation}
335 > \Delta (N,P,T) =  - e^{\beta G}.
336 > \label{introEquation:NPTPartition}
337 > \end{equation}
338 > Here, $G = U - TS + PV$, and $G$ is called Gibbs free energy.
339 >
340 > \subsection{\label{introSection:liouville}Liouville's theorem}
341 >
342 > The Liouville's theorem is the foundation on which statistical
343 > mechanics rests. It describes the time evolution of phase space
344 > distribution function. In order to calculate the rate of change of
345 > $\rho$, we begin from Equation(\ref{introEquation:deltaN}). If we
346 > consider the two faces perpendicular to the $q_1$ axis, which are
347 > located at $q_1$ and $q_1 + \delta q_1$, the number of phase points
348 > leaving the opposite face is given by the expression,
349 > \begin{equation}
350 > \left( {\rho  + \frac{{\partial \rho }}{{\partial q_1 }}\delta q_1 }
351 > \right)\left( {\dot q_1  + \frac{{\partial \dot q_1 }}{{\partial q_1
352 > }}\delta q_1 } \right)\delta q_2  \ldots \delta q_f \delta p_1
353 > \ldots \delta p_f .
354 > \end{equation}
355 > Summing all over the phase space, we obtain
356 > \begin{equation}
357 > \frac{{d(\delta N)}}{{dt}} =  - \sum\limits_{i = 1}^f {\left[ {\rho
358 > \left( {\frac{{\partial \dot q_i }}{{\partial q_i }} +
359 > \frac{{\partial \dot p_i }}{{\partial p_i }}} \right) + \left(
360 > {\frac{{\partial \rho }}{{\partial q_i }}\dot q_i  + \frac{{\partial
361 > \rho }}{{\partial p_i }}\dot p_i } \right)} \right]} \delta q_1
362 > \ldots \delta q_f \delta p_1  \ldots \delta p_f .
363 > \end{equation}
364 > Differentiating the equations of motion in Hamiltonian formalism
365 > (\ref{introEquation:motionHamiltonianCoordinate},
366 > \ref{introEquation:motionHamiltonianMomentum}), we can show,
367 > \begin{equation}
368 > \sum\limits_i {\left( {\frac{{\partial \dot q_i }}{{\partial q_i }}
369 > + \frac{{\partial \dot p_i }}{{\partial p_i }}} \right)}  = 0 ,
370 > \end{equation}
371 > which cancels the first terms of the right hand side. Furthermore,
372 > divining $ \delta q_1  \ldots \delta q_f \delta p_1  \ldots \delta
373 > p_f $ in both sides, we can write out Liouville's theorem in a
374 > simple form,
375 > \begin{equation}
376 > \frac{{\partial \rho }}{{\partial t}} + \sum\limits_{i = 1}^f
377 > {\left( {\frac{{\partial \rho }}{{\partial q_i }}\dot q_i  +
378 > \frac{{\partial \rho }}{{\partial p_i }}\dot p_i } \right)}  = 0 .
379 > \label{introEquation:liouvilleTheorem}
380 > \end{equation}
381 >
382 > Liouville's theorem states that the distribution function is
383 > constant along any trajectory in phase space. In classical
384 > statistical mechanics, since the number of particles in the system
385 > is huge, we may be able to believe the system is stationary,
386 > \begin{equation}
387 > \frac{{\partial \rho }}{{\partial t}} = 0.
388 > \label{introEquation:stationary}
389 > \end{equation}
390 > In such stationary system, the density of distribution $\rho$ can be
391 > connected to the Hamiltonian $H$ through Maxwell-Boltzmann
392 > distribution,
393 > \begin{equation}
394 > \rho  \propto e^{ - \beta H}
395 > \label{introEquation:densityAndHamiltonian}
396 > \end{equation}
397 >
398 > \subsubsection{\label{introSection:phaseSpaceConservation}Conservation of Phase Space}
399 > Lets consider a region in the phase space,
400 > \begin{equation}
401 > \delta v = \int { \ldots \int {dq_1 } ...dq_f dp_1 } ..dp_f .
402 > \end{equation}
403 > If this region is small enough, the density $\rho$ can be regarded
404 > as uniform over the whole phase space. Thus, the number of phase
405 > points inside this region is given by,
406 > \begin{equation}
407 > \delta N = \rho \delta v = \rho \int { \ldots \int {dq_1 } ...dq_f
408 > dp_1 } ..dp_f.
409 > \end{equation}
410 >
411 > \begin{equation}
412 > \frac{{d(\delta N)}}{{dt}} = \frac{{d\rho }}{{dt}}\delta v + \rho
413 > \frac{d}{{dt}}(\delta v) = 0.
414 > \end{equation}
415 > With the help of stationary assumption
416 > (\ref{introEquation:stationary}), we obtain the principle of the
417 > \emph{conservation of extension in phase space},
418 > \begin{equation}
419 > \frac{d}{{dt}}(\delta v) = \frac{d}{{dt}}\int { \ldots \int {dq_1 }
420 > ...dq_f dp_1 } ..dp_f  = 0.
421 > \label{introEquation:volumePreserving}
422 > \end{equation}
423 >
424 > \subsubsection{\label{introSection:liouvilleInOtherForms}Liouville's Theorem in Other Forms}
425 >
426 > Liouville's theorem can be expresses in a variety of different forms
427 > which are convenient within different contexts. For any two function
428 > $F$ and $G$ of the coordinates and momenta of a system, the Poisson
429 > bracket ${F, G}$ is defined as
430 > \begin{equation}
431 > \left\{ {F,G} \right\} = \sum\limits_i {\left( {\frac{{\partial
432 > F}}{{\partial q_i }}\frac{{\partial G}}{{\partial p_i }} -
433 > \frac{{\partial F}}{{\partial p_i }}\frac{{\partial G}}{{\partial
434 > q_i }}} \right)}.
435 > \label{introEquation:poissonBracket}
436 > \end{equation}
437 > Substituting equations of motion in Hamiltonian formalism(
438 > \ref{introEquation:motionHamiltonianCoordinate} ,
439 > \ref{introEquation:motionHamiltonianMomentum} ) into
440 > (\ref{introEquation:liouvilleTheorem}), we can rewrite Liouville's
441 > theorem using Poisson bracket notion,
442 > \begin{equation}
443 > \left( {\frac{{\partial \rho }}{{\partial t}}} \right) =  - \left\{
444 > {\rho ,H} \right\}.
445 > \label{introEquation:liouvilleTheromInPoissin}
446 > \end{equation}
447 > Moreover, the Liouville operator is defined as
448 > \begin{equation}
449 > iL = \sum\limits_{i = 1}^f {\left( {\frac{{\partial H}}{{\partial
450 > p_i }}\frac{\partial }{{\partial q_i }} - \frac{{\partial
451 > H}}{{\partial q_i }}\frac{\partial }{{\partial p_i }}} \right)}
452 > \label{introEquation:liouvilleOperator}
453 > \end{equation}
454 > In terms of Liouville operator, Liouville's equation can also be
455 > expressed as
456 > \begin{equation}
457 > \left( {\frac{{\partial \rho }}{{\partial t}}} \right) =  - iL\rho
458 > \label{introEquation:liouvilleTheoremInOperator}
459 > \end{equation}
460 >
461 > \subsection{\label{introSection:ergodic}The Ergodic Hypothesis}
462 >
463 > Various thermodynamic properties can be calculated from Molecular
464 > Dynamics simulation. By comparing experimental values with the
465 > calculated properties, one can determine the accuracy of the
466 > simulation and the quality of the underlying model. However, both of
467 > experiment and computer simulation are usually performed during a
468 > certain time interval and the measurements are averaged over a
469 > period of them which is different from the average behavior of
470 > many-body system in Statistical Mechanics. Fortunately, Ergodic
471 > Hypothesis is proposed to make a connection between time average and
472 > ensemble average. It states that time average and average over the
473 > statistical ensemble are identical \cite{Frenkel1996, leach01:mm}.
474 > \begin{equation}
475 > \langle A(q , p) \rangle_t = \mathop {\lim }\limits_{t \to \infty }
476 > \frac{1}{t}\int\limits_0^t {A(q(t),p(t))dt = \int\limits_\Gamma
477 > {A(q(t),p(t))} } \rho (q(t), p(t)) dqdp
478 > \end{equation}
479 > where $\langle  A(q , p) \rangle_t$ is an equilibrium value of a
480 > physical quantity and $\rho (p(t), q(t))$ is the equilibrium
481 > distribution function. If an observation is averaged over a
482 > sufficiently long time (longer than relaxation time), all accessible
483 > microstates in phase space are assumed to be equally probed, giving
484 > a properly weighted statistical average. This allows the researcher
485 > freedom of choice when deciding how best to measure a given
486 > observable. In case an ensemble averaged approach sounds most
487 > reasonable, the Monte Carlo techniques\cite{metropolis:1949} can be
488 > utilized. Or if the system lends itself to a time averaging
489 > approach, the Molecular Dynamics techniques in
490 > Sec.~\ref{introSection:molecularDynamics} will be the best
491 > choice\cite{Frenkel1996}.
492 >
493 > \section{\label{introSection:geometricIntegratos}Geometric Integrators}
494 > A variety of numerical integrators were proposed to simulate the
495 > motions. They usually begin with an initial conditionals and move
496 > the objects in the direction governed by the differential equations.
497 > However, most of them ignore the hidden physical law contained
498 > within the equations. Since 1990, geometric integrators, which
499 > preserve various phase-flow invariants such as symplectic structure,
500 > volume and time reversal symmetry, are developed to address this
501 > issue. The velocity verlet method, which happens to be a simple
502 > example of symplectic integrator, continues to gain its popularity
503 > in molecular dynamics community. This fact can be partly explained
504 > by its geometric nature.
505 >
506 > \subsection{\label{introSection:symplecticManifold}Symplectic Manifold}
507 > A \emph{manifold} is an abstract mathematical space. It locally
508 > looks like Euclidean space, but when viewed globally, it may have
509 > more complicate structure. A good example of manifold is the surface
510 > of Earth. It seems to be flat locally, but it is round if viewed as
511 > a whole. A \emph{differentiable manifold} (also known as
512 > \emph{smooth manifold}) is a manifold with an open cover in which
513 > the covering neighborhoods are all smoothly isomorphic to one
514 > another. In other words,it is possible to apply calculus on
515 > \emph{differentiable manifold}. A \emph{symplectic manifold} is
516 > defined as a pair $(M, \omega)$ which consisting of a
517 > \emph{differentiable manifold} $M$ and a close, non-degenerated,
518 > bilinear symplectic form, $\omega$. A symplectic form on a vector
519 > space $V$ is a function $\omega(x, y)$ which satisfies
520 > $\omega(\lambda_1x_1+\lambda_2x_2, y) = \lambda_1\omega(x_1, y)+
521 > \lambda_2\omega(x_2, y)$, $\omega(x, y) = - \omega(y, x)$ and
522 > $\omega(x, x) = 0$. Cross product operation in vector field is an
523 > example of symplectic form.
524 >
525 > One of the motivations to study \emph{symplectic manifold} in
526 > Hamiltonian Mechanics is that a symplectic manifold can represent
527 > all possible configurations of the system and the phase space of the
528 > system can be described by it's cotangent bundle. Every symplectic
529 > manifold is even dimensional. For instance, in Hamilton equations,
530 > coordinate and momentum always appear in pairs.
531 >
532 > Let  $(M,\omega)$ and $(N, \eta)$ be symplectic manifolds. A map
533 > \[
534 > f : M \rightarrow N
535 > \]
536 > is a \emph{symplectomorphism} if it is a \emph{diffeomorphims} and
537 > the \emph{pullback} of $\eta$ under f is equal to $\omega$.
538 > Canonical transformation is an example of symplectomorphism in
539 > classical mechanics.
540 >
541 > \subsection{\label{introSection:ODE}Ordinary Differential Equations}
542 >
543 > For a ordinary differential system defined as
544 > \begin{equation}
545 > \dot x = f(x)
546 > \end{equation}
547 > where $x = x(q,p)^T$, this system is canonical Hamiltonian, if
548 > \begin{equation}
549 > f(r) = J\nabla _x H(r).
550 > \end{equation}
551 > $H = H (q, p)$ is Hamiltonian function and $J$ is the skew-symmetric
552 > matrix
553 > \begin{equation}
554 > J = \left( {\begin{array}{*{20}c}
555 >   0 & I  \\
556 >   { - I} & 0  \\
557 > \end{array}} \right)
558 > \label{introEquation:canonicalMatrix}
559 > \end{equation}
560 > where $I$ is an identity matrix. Using this notation, Hamiltonian
561 > system can be rewritten as,
562 > \begin{equation}
563 > \frac{d}{{dt}}x = J\nabla _x H(x)
564 > \label{introEquation:compactHamiltonian}
565 > \end{equation}In this case, $f$ is
566 > called a \emph{Hamiltonian vector field}.
567 >
568 > Another generalization of Hamiltonian dynamics is Poisson Dynamics,
569 > \begin{equation}
570 > \dot x = J(x)\nabla _x H \label{introEquation:poissonHamiltonian}
571 > \end{equation}
572 > The most obvious change being that matrix $J$ now depends on $x$.
573 >
574 > \subsection{\label{introSection:exactFlow}Exact Flow}
575 >
576 > Let $x(t)$ be the exact solution of the ODE system,
577 > \begin{equation}
578 > \frac{{dx}}{{dt}} = f(x) \label{introEquation:ODE}
579 > \end{equation}
580 > The exact flow(solution) $\varphi_\tau$ is defined by
581 > \[
582 > x(t+\tau) =\varphi_\tau(x(t))
583 > \]
584 > where $\tau$ is a fixed time step and $\varphi$ is a map from phase
585 > space to itself. The flow has the continuous group property,
586 > \begin{equation}
587 > \varphi _{\tau _1 }  \circ \varphi _{\tau _2 }  = \varphi _{\tau _1
588 > + \tau _2 } .
589 > \end{equation}
590 > In particular,
591 > \begin{equation}
592 > \varphi _\tau   \circ \varphi _{ - \tau }  = I
593 > \end{equation}
594 > Therefore, the exact flow is self-adjoint,
595 > \begin{equation}
596 > \varphi _\tau   = \varphi _{ - \tau }^{ - 1}.
597 > \end{equation}
598 > The exact flow can also be written in terms of the of an operator,
599 > \begin{equation}
600 > \varphi _\tau  (x) = e^{\tau \sum\limits_i {f_i (x)\frac{\partial
601 > }{{\partial x_i }}} } (x) \equiv \exp (\tau f)(x).
602 > \label{introEquation:exponentialOperator}
603 > \end{equation}
604 >
605 > In most cases, it is not easy to find the exact flow $\varphi_\tau$.
606 > Instead, we use a approximate map, $\psi_\tau$, which is usually
607 > called integrator. The order of an integrator $\psi_\tau$ is $p$, if
608 > the Taylor series of $\psi_\tau$ agree to order $p$,
609 > \begin{equation}
610 > \psi_tau(x) = x + \tau f(x) + O(\tau^{p+1})
611 > \end{equation}
612 >
613 > \subsection{\label{introSection:geometricProperties}Geometric Properties}
614 >
615 > The hidden geometric properties of ODE and its flow play important
616 > roles in numerical studies. Many of them can be found in systems
617 > which occur naturally in applications.
618 >
619 > Let $\varphi$ be the flow of Hamiltonian vector field, $\varphi$ is
620 > a \emph{symplectic} flow if it satisfies,
621 > \begin{equation}
622 > {\varphi '}^T J \varphi ' = J.
623 > \end{equation}
624 > According to Liouville's theorem, the symplectic volume is invariant
625 > under a Hamiltonian flow, which is the basis for classical
626 > statistical mechanics. Furthermore, the flow of a Hamiltonian vector
627 > field on a symplectic manifold can be shown to be a
628 > symplectomorphism. As to the Poisson system,
629 > \begin{equation}
630 > {\varphi '}^T J \varphi ' = J \circ \varphi
631 > \end{equation}
632 > is the property must be preserved by the integrator.
633 >
634 > It is possible to construct a \emph{volume-preserving} flow for a
635 > source free($ \nabla \cdot f = 0 $) ODE, if the flow satisfies $
636 > \det d\varphi  = 1$. One can show easily that a symplectic flow will
637 > be volume-preserving.
638 >
639 > Changing the variables $y = h(x)$ in a ODE\ref{introEquation:ODE}
640 > will result in a new system,
641 > \[
642 > \dot y = \tilde f(y) = ((dh \cdot f)h^{ - 1} )(y).
643 > \]
644 > The vector filed $f$ has reversing symmetry $h$ if $f = - \tilde f$.
645 > In other words, the flow of this vector field is reversible if and
646 > only if $ h \circ \varphi ^{ - 1}  = \varphi  \circ h $.
647 >
648 > A \emph{first integral}, or conserved quantity of a general
649 > differential function is a function $ G:R^{2d}  \to R^d $ which is
650 > constant for all solutions of the ODE $\frac{{dx}}{{dt}} = f(x)$ ,
651 > \[
652 > \frac{{dG(x(t))}}{{dt}} = 0.
653 > \]
654 > Using chain rule, one may obtain,
655 > \[
656 > \sum\limits_i {\frac{{dG}}{{dx_i }}} f_i (x) = f \bullet \nabla G,
657 > \]
658 > which is the condition for conserving \emph{first integral}. For a
659 > canonical Hamiltonian system, the time evolution of an arbitrary
660 > smooth function $G$ is given by,
661 > \begin{equation}
662 > \begin{array}{c}
663 > \frac{{dG(x(t))}}{{dt}} = [\nabla _x G(x(t))]^T \dot x(t) \\
664 >  = [\nabla _x G(x(t))]^T J\nabla _x H(x(t)). \\
665 > \end{array}
666 > \label{introEquation:firstIntegral1}
667 > \end{equation}
668 > Using poisson bracket notion, Equation
669 > \ref{introEquation:firstIntegral1} can be rewritten as
670 > \[
671 > \frac{d}{{dt}}G(x(t)) = \left\{ {G,H} \right\}(x(t)).
672 > \]
673 > Therefore, the sufficient condition for $G$ to be the \emph{first
674 > integral} of a Hamiltonian system is
675 > \[
676 > \left\{ {G,H} \right\} = 0.
677 > \]
678 > As well known, the Hamiltonian (or energy) H of a Hamiltonian system
679 > is a \emph{first integral}, which is due to the fact $\{ H,H\}  =
680 > 0$.
681 >
682 >
683 > When designing any numerical methods, one should always try to
684 > preserve the structural properties of the original ODE and its flow.
685 >
686 > \subsection{\label{introSection:constructionSymplectic}Construction of Symplectic Methods}
687 > A lot of well established and very effective numerical methods have
688 > been successful precisely because of their symplecticities even
689 > though this fact was not recognized when they were first
690 > constructed. The most famous example is leapfrog methods in
691 > molecular dynamics. In general, symplectic integrators can be
692 > constructed using one of four different methods.
693 > \begin{enumerate}
694 > \item Generating functions
695 > \item Variational methods
696 > \item Runge-Kutta methods
697 > \item Splitting methods
698 > \end{enumerate}
699 >
700 > Generating function tends to lead to methods which are cumbersome
701 > and difficult to use. In dissipative systems, variational methods
702 > can capture the decay of energy accurately. Since their
703 > geometrically unstable nature against non-Hamiltonian perturbations,
704 > ordinary implicit Runge-Kutta methods are not suitable for
705 > Hamiltonian system. Recently, various high-order explicit
706 > Runge--Kutta methods have been developed to overcome this
707 > instability. However, due to computational penalty involved in
708 > implementing the Runge-Kutta methods, they do not attract too much
709 > attention from Molecular Dynamics community. Instead, splitting have
710 > been widely accepted since they exploit natural decompositions of
711 > the system\cite{Tuckerman92}.
712 >
713 > \subsubsection{\label{introSection:splittingMethod}Splitting Method}
714 >
715 > The main idea behind splitting methods is to decompose the discrete
716 > $\varphi_h$ as a composition of simpler flows,
717 > \begin{equation}
718 > \varphi _h  = \varphi _{h_1 }  \circ \varphi _{h_2 }  \ldots  \circ
719 > \varphi _{h_n }
720 > \label{introEquation:FlowDecomposition}
721 > \end{equation}
722 > where each of the sub-flow is chosen such that each represent a
723 > simpler integration of the system.
724 >
725 > Suppose that a Hamiltonian system takes the form,
726 > \[
727 > H = H_1 + H_2.
728 > \]
729 > Here, $H_1$ and $H_2$ may represent different physical processes of
730 > the system. For instance, they may relate to kinetic and potential
731 > energy respectively, which is a natural decomposition of the
732 > problem. If $H_1$ and $H_2$ can be integrated using exact flows
733 > $\varphi_1(t)$ and $\varphi_2(t)$, respectively, a simple first
734 > order is then given by the Lie-Trotter formula
735 > \begin{equation}
736 > \varphi _h  = \varphi _{1,h}  \circ \varphi _{2,h},
737 > \label{introEquation:firstOrderSplitting}
738 > \end{equation}
739 > where $\varphi _h$ is the result of applying the corresponding
740 > continuous $\varphi _i$ over a time $h$. By definition, as
741 > $\varphi_i(t)$ is the exact solution of a Hamiltonian system, it
742 > must follow that each operator $\varphi_i(t)$ is a symplectic map.
743 > It is easy to show that any composition of symplectic flows yields a
744 > symplectic map,
745 > \begin{equation}
746 > (\varphi '\phi ')^T J\varphi '\phi ' = \phi '^T \varphi '^T J\varphi
747 > '\phi ' = \phi '^T J\phi ' = J,
748 > \label{introEquation:SymplecticFlowComposition}
749 > \end{equation}
750 > where $\phi$ and $\psi$ both are symplectic maps. Thus operator
751 > splitting in this context automatically generates a symplectic map.
752 >
753 > The Lie-Trotter splitting(\ref{introEquation:firstOrderSplitting})
754 > introduces local errors proportional to $h^2$, while Strang
755 > splitting gives a second-order decomposition,
756 > \begin{equation}
757 > \varphi _h  = \varphi _{1,h/2}  \circ \varphi _{2,h}  \circ \varphi
758 > _{1,h/2} , \label{introEquation:secondOrderSplitting}
759 > \end{equation}
760 > which has a local error proportional to $h^3$. Sprang splitting's
761 > popularity in molecular simulation community attribute to its
762 > symmetric property,
763 > \begin{equation}
764 > \varphi _h^{ - 1} = \varphi _{ - h}.
765 > \label{introEquation:timeReversible}
766 > \end{equation}
767 >
768 > \subsubsection{\label{introSection:exampleSplittingMethod}Example of Splitting Method}
769 > The classical equation for a system consisting of interacting
770 > particles can be written in Hamiltonian form,
771 > \[
772 > H = T + V
773 > \]
774 > where $T$ is the kinetic energy and $V$ is the potential energy.
775 > Setting $H_1 = T, H_2 = V$ and applying Strang splitting, one
776 > obtains the following:
777 > \begin{align}
778 > q(\Delta t) &= q(0) + \dot{q}(0)\Delta t +
779 >    \frac{F[q(0)]}{m}\frac{\Delta t^2}{2}, %
780 > \label{introEquation:Lp10a} \\%
781 > %
782 > \dot{q}(\Delta t) &= \dot{q}(0) + \frac{\Delta t}{2m}
783 >    \biggl [F[q(0)] + F[q(\Delta t)] \biggr]. %
784 > \label{introEquation:Lp10b}
785 > \end{align}
786 > where $F(t)$ is the force at time $t$. This integration scheme is
787 > known as \emph{velocity verlet} which is
788 > symplectic(\ref{introEquation:SymplecticFlowComposition}),
789 > time-reversible(\ref{introEquation:timeReversible}) and
790 > volume-preserving (\ref{introEquation:volumePreserving}). These
791 > geometric properties attribute to its long-time stability and its
792 > popularity in the community. However, the most commonly used
793 > velocity verlet integration scheme is written as below,
794 > \begin{align}
795 > \dot{q}\biggl (\frac{\Delta t}{2}\biggr ) &=
796 >    \dot{q}(0) + \frac{\Delta t}{2m}\, F[q(0)], \label{introEquation:Lp9a}\\%
797 > %
798 > q(\Delta t) &= q(0) + \Delta t\, \dot{q}\biggl (\frac{\Delta t}{2}\biggr ),%
799 >    \label{introEquation:Lp9b}\\%
800 > %
801 > \dot{q}(\Delta t) &= \dot{q}\biggl (\frac{\Delta t}{2}\biggr ) +
802 >    \frac{\Delta t}{2m}\, F[q(0)]. \label{introEquation:Lp9c}
803 > \end{align}
804 > From the preceding splitting, one can see that the integration of
805 > the equations of motion would follow:
806 > \begin{enumerate}
807 > \item calculate the velocities at the half step, $\frac{\Delta t}{2}$, from the forces calculated at the initial position.
808 >
809 > \item Use the half step velocities to move positions one whole step, $\Delta t$.
810 >
811 > \item Evaluate the forces at the new positions, $\mathbf{r}(\Delta t)$, and use the new forces to complete the velocity move.
812 >
813 > \item Repeat from step 1 with the new position, velocities, and forces assuming the roles of the initial values.
814 > \end{enumerate}
815 >
816 > Simply switching the order of splitting and composing, a new
817 > integrator, the \emph{position verlet} integrator, can be generated,
818 > \begin{align}
819 > \dot q(\Delta t) &= \dot q(0) + \Delta tF(q(0))\left[ {q(0) +
820 > \frac{{\Delta t}}{{2m}}\dot q(0)} \right], %
821 > \label{introEquation:positionVerlet1} \\%
822 > %
823 > q(\Delta t) &= q(0) + \frac{{\Delta t}}{2}\left[ {\dot q(0) + \dot
824 > q(\Delta t)} \right]. %
825 > \label{introEquation:positionVerlet2}
826 > \end{align}
827 >
828 > \subsubsection{\label{introSection:errorAnalysis}Error Analysis and Higher Order Methods}
829 >
830 > Baker-Campbell-Hausdorff formula can be used to determine the local
831 > error of splitting method in terms of commutator of the
832 > operators(\ref{introEquation:exponentialOperator}) associated with
833 > the sub-flow. For operators $hX$ and $hY$ which are associate to
834 > $\varphi_1(t)$ and $\varphi_2(t$ respectively , we have
835 > \begin{equation}
836 > \exp (hX + hY) = \exp (hZ)
837 > \end{equation}
838 > where
839 > \begin{equation}
840 > hZ = hX + hY + \frac{{h^2 }}{2}[X,Y] + \frac{{h^3 }}{2}\left(
841 > {[X,[X,Y]] + [Y,[Y,X]]} \right) +  \ldots .
842 > \end{equation}
843 > Here, $[X,Y]$ is the commutators of operator $X$ and $Y$ given by
844 > \[
845 > [X,Y] = XY - YX .
846 > \]
847 > Applying Baker-Campbell-Hausdorff formula to Sprang splitting, we
848 > can obtain
849 > \begin{eqnarray*}
850 > \exp (h X/2)\exp (h Y)\exp (h X/2) & = & \exp (h X + h Y + h^2
851 > [X,Y]/4 + h^2 [Y,X]/4 \\ & & \mbox{} + h^2 [X,X]/8 + h^2 [Y,Y]/8 \\
852 > & & \mbox{} + h^3 [Y,[Y,X]]/12 - h^3 [X,[X,Y]]/24 & & \mbox{} +
853 > \ldots )
854 > \end{eqnarray*}
855 > Since \[ [X,Y] + [Y,X] = 0\] and \[ [X,X] = 0\], the dominant local
856 > error of Spring splitting is proportional to $h^3$. The same
857 > procedure can be applied to general splitting,  of the form
858 > \begin{equation}
859 > \varphi _{b_m h}^2  \circ \varphi _{a_m h}^1  \circ \varphi _{b_{m -
860 > 1} h}^2  \circ  \ldots  \circ \varphi _{a_1 h}^1 .
861 > \end{equation}
862 > Careful choice of coefficient $a_1 ,\ldot , b_m$ will lead to higher
863 > order method. Yoshida proposed an elegant way to compose higher
864 > order methods based on symmetric splitting. Given a symmetric second
865 > order base method $ \varphi _h^{(2)} $, a fourth-order symmetric
866 > method can be constructed by composing,
867 > \[
868 > \varphi _h^{(4)}  = \varphi _{\alpha h}^{(2)}  \circ \varphi _{\beta
869 > h}^{(2)}  \circ \varphi _{\alpha h}^{(2)}
870 > \]
871 > where $ \alpha  =  - \frac{{2^{1/3} }}{{2 - 2^{1/3} }}$ and $ \beta
872 > = \frac{{2^{1/3} }}{{2 - 2^{1/3} }}$. Moreover, a symmetric
873 > integrator $ \varphi _h^{(2n + 2)}$ can be composed by
874 > \begin{equation}
875 > \varphi _h^{(2n + 2)}  = \varphi _{\alpha h}^{(2n)}  \circ \varphi
876 > _{\beta h}^{(2n)}  \circ \varphi _{\alpha h}^{(2n)}
877 > \end{equation}
878 > , if the weights are chosen as
879 > \[
880 > \alpha  =  - \frac{{2^{1/(2n + 1)} }}{{2 - 2^{1/(2n + 1)} }},\beta =
881 > \frac{{2^{1/(2n + 1)} }}{{2 - 2^{1/(2n + 1)} }} .
882 > \]
883 >
884 > \section{\label{introSection:molecularDynamics}Molecular Dynamics}
885 >
886 > As a special discipline of molecular modeling, Molecular dynamics
887 > has proven to be a powerful tool for studying the functions of
888 > biological systems, providing structural, thermodynamic and
889 > dynamical information.
890 >
891 > One of the principal tools for modeling proteins, nucleic acids and
892 > their complexes. Stability of proteins Folding of proteins.
893 > Molecular recognition by:proteins, DNA, RNA, lipids, hormones STP,
894 > etc. Enzyme reactions Rational design of biologically active
895 > molecules (drug design) Small and large-scale conformational
896 > changes. determination and construction of 3D structures (homology,
897 > Xray diffraction, NMR) Dynamic processes such as ion transport in
898 > biological systems.
899 >
900 > Macroscopic properties are related to microscopic behavior.
901 >
902 > Time dependent (and independent) microscopic behavior of a molecule
903 > can be calculated by molecular dynamics simulations.
904 >
905 > \subsection{\label{introSec:mdInit}Initialization}
906 >
907 > \subsection{\label{introSec:forceEvaluation}Force Evaluation}
908 >
909 > \subsection{\label{introSection:mdIntegration} Integration of the Equations of Motion}
910 >
911 > \section{\label{introSection:rigidBody}Dynamics of Rigid Bodies}
912 >
913 > Rigid bodies are frequently involved in the modeling of different
914 > areas, from engineering, physics, to chemistry. For example,
915 > missiles and vehicle are usually modeled by rigid bodies.  The
916 > movement of the objects in 3D gaming engine or other physics
917 > simulator is governed by the rigid body dynamics. In molecular
918 > simulation, rigid body is used to simplify the model in
919 > protein-protein docking study{\cite{Gray03}}.
920 >
921 > It is very important to develop stable and efficient methods to
922 > integrate the equations of motion of orientational degrees of
923 > freedom. Euler angles are the nature choice to describe the
924 > rotational degrees of freedom. However, due to its singularity, the
925 > numerical integration of corresponding equations of motion is very
926 > inefficient and inaccurate. Although an alternative integrator using
927 > different sets of Euler angles can overcome this difficulty\cite{},
928 > the computational penalty and the lost of angular momentum
929 > conservation still remain. A singularity free representation
930 > utilizing quaternions was developed by Evans in 1977. Unfortunately,
931 > this approach suffer from the nonseparable Hamiltonian resulted from
932 > quaternion representation, which prevents the symplectic algorithm
933 > to be utilized. Another different approach is to apply holonomic
934 > constraints to the atoms belonging to the rigid body. Each atom
935 > moves independently under the normal forces deriving from potential
936 > energy and constraint forces which are used to guarantee the
937 > rigidness. However, due to their iterative nature, SHAKE and Rattle
938 > algorithm converge very slowly when the number of constraint
939 > increases.
940 >
941 > The break through in geometric literature suggests that, in order to
942 > develop a long-term integration scheme, one should preserve the
943 > symplectic structure of the flow. Introducing conjugate momentum to
944 > rotation matrix $Q$ and re-formulating Hamiltonian's equation, a
945 > symplectic integrator, RSHAKE, was proposed to evolve the
946 > Hamiltonian system in a constraint manifold by iteratively
947 > satisfying the orthogonality constraint $Q_T Q = 1$. An alternative
948 > method using quaternion representation was developed by Omelyan.
949 > However, both of these methods are iterative and inefficient. In
950 > this section, we will present a symplectic Lie-Poisson integrator
951 > for rigid body developed by Dullweber and his
952 > coworkers\cite{Dullweber1997} in depth.
953 >
954 > \subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Body}
955 > The motion of the rigid body is Hamiltonian with the Hamiltonian
956 > function
957 > \begin{equation}
958 > H = \frac{1}{2}(p^T m^{ - 1} p) + \frac{1}{2}tr(PJ^{ - 1} P) +
959 > V(q,Q) + \frac{1}{2}tr[(QQ^T  - 1)\Lambda ].
960 > \label{introEquation:RBHamiltonian}
961 > \end{equation}
962 > Here, $q$ and $Q$  are the position and rotation matrix for the
963 > rigid-body, $p$ and $P$  are conjugate momenta to $q$  and $Q$ , and
964 > $J$, a diagonal matrix, is defined by
965 > \[
966 > I_{ii}^{ - 1}  = \frac{1}{2}\sum\limits_{i \ne j} {J_{jj}^{ - 1} }
967 > \]
968 > where $I_{ii}$ is the diagonal element of the inertia tensor. This
969 > constrained Hamiltonian equation subjects to a holonomic constraint,
970 > \begin{equation}
971 > Q^T Q = 1$, \label{introEquation:orthogonalConstraint}
972 > \end{equation}
973 > which is used to ensure rotation matrix's orthogonality.
974 > Differentiating \ref{introEquation:orthogonalConstraint} and using
975 > Equation \ref{introEquation:RBMotionMomentum}, one may obtain,
976 > \begin{equation}
977 > Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0 . \\
978 > \label{introEquation:RBFirstOrderConstraint}
979 > \end{equation}
980 >
981 > Using Equation (\ref{introEquation:motionHamiltonianCoordinate},
982 > \ref{introEquation:motionHamiltonianMomentum}), one can write down
983 > the equations of motion,
984 > \[
985 > \begin{array}{c}
986 > \frac{{dq}}{{dt}} = \frac{p}{m} \label{introEquation:RBMotionPosition}\\
987 > \frac{{dp}}{{dt}} =  - \nabla _q V(q,Q) \label{introEquation:RBMotionMomentum}\\
988 > \frac{{dQ}}{{dt}} = PJ^{ - 1}  \label{introEquation:RBMotionRotation}\\
989 > \frac{{dP}}{{dt}} =  - \nabla _Q V(q,Q) - 2Q\Lambda . \label{introEquation:RBMotionP}\\
990 > \end{array}
991 > \]
992 >
993 > In general, there are two ways to satisfy the holonomic constraints.
994 > We can use constraint force provided by lagrange multiplier on the
995 > normal manifold to keep the motion on constraint space. Or we can
996 > simply evolve the system in constraint manifold. The two method are
997 > proved to be equivalent. The holonomic constraint and equations of
998 > motions define a constraint manifold for rigid body
999 > \[
1000 > M = \left\{ {(Q,P):Q^T Q = 1,Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0}
1001 > \right\}.
1002 > \]
1003 >
1004 > Unfortunately, this constraint manifold is not the cotangent bundle
1005 > $T_{\star}SO(3)$. However, it turns out that under symplectic
1006 > transformation, the cotangent space and the phase space are
1007 > diffeomorphic. Introducing
1008 > \[
1009 > \tilde Q = Q,\tilde P = \frac{1}{2}\left( {P - QP^T Q} \right),
1010 > \]
1011 > the mechanical system subject to a holonomic constraint manifold $M$
1012 > can be re-formulated as a Hamiltonian system on the cotangent space
1013 > \[
1014 > T^* SO(3) = \left\{ {(\tilde Q,\tilde P):\tilde Q^T \tilde Q =
1015 > 1,\tilde Q^T \tilde PJ^{ - 1}  + J^{ - 1} P^T \tilde Q = 0} \right\}
1016 > \]
1017 >
1018 > For a body fixed vector $X_i$ with respect to the center of mass of
1019 > the rigid body, its corresponding lab fixed vector $X_0^{lab}$  is
1020 > given as
1021 > \begin{equation}
1022 > X_i^{lab} = Q X_i + q.
1023 > \end{equation}
1024 > Therefore, potential energy $V(q,Q)$ is defined by
1025 > \[
1026 > V(q,Q) = V(Q X_0 + q).
1027 > \]
1028 > Hence, the force and torque are given by
1029 > \[
1030 > \nabla _q V(q,Q) = F(q,Q) = \sum\limits_i {F_i (q,Q)},
1031 > \]
1032 > and
1033 > \[
1034 > \nabla _Q V(q,Q) = F(q,Q)X_i^t
1035 > \]
1036 > respectively.
1037 >
1038 > As a common choice to describe the rotation dynamics of the rigid
1039 > body, angular momentum on body frame $\Pi  = Q^t P$ is introduced to
1040 > rewrite the equations of motion,
1041 > \begin{equation}
1042 > \begin{array}{l}
1043 > \mathop \Pi \limits^ \bullet   = J^{ - 1} \Pi ^T \Pi  + Q^T \sum\limits_i {F_i (q,Q)X_i^T }  - \Lambda  \\
1044 > \mathop Q\limits^{{\rm{   }} \bullet }  = Q\Pi {\rm{ }}J^{ - 1}  \\
1045 > \end{array}
1046 > \label{introEqaution:RBMotionPI}
1047 > \end{equation}
1048 > , as well as holonomic constraints,
1049 > \[
1050 > \begin{array}{l}
1051 > \Pi J^{ - 1}  + J^{ - 1} \Pi ^t  = 0 \\
1052 > Q^T Q = 1 \\
1053 > \end{array}
1054 > \]
1055 >
1056 > For a vector $v(v_1 ,v_2 ,v_3 ) \in R^3$ and a matrix $\hat v \in
1057 > so(3)^ \star$, the hat-map isomorphism,
1058 > \begin{equation}
1059 > v(v_1 ,v_2 ,v_3 ) \Leftrightarrow \hat v = \left(
1060 > {\begin{array}{*{20}c}
1061 >   0 & { - v_3 } & {v_2 }  \\
1062 >   {v_3 } & 0 & { - v_1 }  \\
1063 >   { - v_2 } & {v_1 } & 0  \\
1064 > \end{array}} \right),
1065 > \label{introEquation:hatmapIsomorphism}
1066 > \end{equation}
1067 > will let us associate the matrix products with traditional vector
1068 > operations
1069 > \[
1070 > \hat vu = v \times u
1071 > \]
1072 >
1073 > Using \ref{introEqaution:RBMotionPI}, one can construct a skew
1074 > matrix,
1075 > \begin{equation}
1076 > (\mathop \Pi \limits^ \bullet   - \mathop \Pi \limits^ \bullet  ^T
1077 > ){\rm{ }} = {\rm{ }}(\Pi  - \Pi ^T ){\rm{ }}(J^{ - 1} \Pi  + \Pi J^{
1078 > - 1} ) + \sum\limits_i {[Q^T F_i (r,Q)X_i^T  - X_i F_i (r,Q)^T Q]} -
1079 > (\Lambda  - \Lambda ^T ) . \label{introEquation:skewMatrixPI}
1080 > \end{equation}
1081 > Since $\Lambda$ is symmetric, the last term of Equation
1082 > \ref{introEquation:skewMatrixPI} is zero, which implies the Lagrange
1083 > multiplier $\Lambda$ is absent from the equations of motion. This
1084 > unique property eliminate the requirement of iterations which can
1085 > not be avoided in other methods\cite{}.
1086 >
1087 > Applying hat-map isomorphism, we obtain the equation of motion for
1088 > angular momentum on body frame
1089 > \begin{equation}
1090 > \dot \pi  = \pi  \times I^{ - 1} \pi  + \sum\limits_i {\left( {Q^T
1091 > F_i (r,Q)} \right) \times X_i }.
1092 > \label{introEquation:bodyAngularMotion}
1093 > \end{equation}
1094 > In the same manner, the equation of motion for rotation matrix is
1095 > given by
1096 > \[
1097 > \dot Q = Qskew(I^{ - 1} \pi )
1098 > \]
1099 >
1100 > \subsection{\label{introSection:SymplecticFreeRB}Symplectic
1101 > Lie-Poisson Integrator for Free Rigid Body}
1102 >
1103 > If there is not external forces exerted on the rigid body, the only
1104 > contribution to the rotational is from the kinetic potential (the
1105 > first term of \ref{ introEquation:bodyAngularMotion}). The free
1106 > rigid body is an example of Lie-Poisson system with Hamiltonian
1107 > function
1108 > \begin{equation}
1109 > T^r (\pi ) = T_1 ^r (\pi _1 ) + T_2^r (\pi _2 ) + T_3^r (\pi _3 )
1110 > \label{introEquation:rotationalKineticRB}
1111 > \end{equation}
1112 > where $T_i^r (\pi _i ) = \frac{{\pi _i ^2 }}{{2I_i }}$ and
1113 > Lie-Poisson structure matrix,
1114 > \begin{equation}
1115 > J(\pi ) = \left( {\begin{array}{*{20}c}
1116 >   0 & {\pi _3 } & { - \pi _2 }  \\
1117 >   { - \pi _3 } & 0 & {\pi _1 }  \\
1118 >   {\pi _2 } & { - \pi _1 } & 0  \\
1119 > \end{array}} \right)
1120 > \end{equation}
1121 > Thus, the dynamics of free rigid body is governed by
1122 > \begin{equation}
1123 > \frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi )
1124 > \end{equation}
1125 >
1126 > One may notice that each $T_i^r$ in Equation
1127 > \ref{introEquation:rotationalKineticRB} can be solved exactly. For
1128 > instance, the equations of motion due to $T_1^r$ are given by
1129 > \begin{equation}
1130 > \frac{d}{{dt}}\pi  = R_1 \pi ,\frac{d}{{dt}}Q = QR_1
1131 > \label{introEqaution:RBMotionSingleTerm}
1132 > \end{equation}
1133 > where
1134 > \[ R_1  = \left( {\begin{array}{*{20}c}
1135 >   0 & 0 & 0  \\
1136 >   0 & 0 & {\pi _1 }  \\
1137 >   0 & { - \pi _1 } & 0  \\
1138 > \end{array}} \right).
1139 > \]
1140 > The solutions of Equation \ref{introEqaution:RBMotionSingleTerm} is
1141 > \[
1142 > \pi (\Delta t) = e^{\Delta tR_1 } \pi (0),Q(\Delta t) =
1143 > Q(0)e^{\Delta tR_1 }
1144 > \]
1145 > with
1146 > \[
1147 > e^{\Delta tR_1 }  = \left( {\begin{array}{*{20}c}
1148 >   0 & 0 & 0  \\
1149 >   0 & {\cos \theta _1 } & {\sin \theta _1 }  \\
1150 >   0 & { - \sin \theta _1 } & {\cos \theta _1 }  \\
1151 > \end{array}} \right),\theta _1  = \frac{{\pi _1 }}{{I_1 }}\Delta t.
1152 > \]
1153 > To reduce the cost of computing expensive functions in $e^{\Delta
1154 > tR_1 }$, we can use Cayley transformation,
1155 > \[
1156 > e^{\Delta tR_1 }  \approx (1 - \Delta tR_1 )^{ - 1} (1 + \Delta tR_1
1157 > )
1158 > \]
1159 >
1160 > The flow maps for $T_2^r$ and $T_2^r$ can be found in the same
1161 > manner.
1162 >
1163 > In order to construct a second-order symplectic method, we split the
1164 > angular kinetic Hamiltonian function can into five terms
1165 > \[
1166 > T^r (\pi ) = \frac{1}{2}T_1 ^r (\pi _1 ) + \frac{1}{2}T_2^r (\pi _2
1167 > ) + T_3^r (\pi _3 ) + \frac{1}{2}T_2^r (\pi _2 ) + \frac{1}{2}T_1 ^r
1168 > (\pi _1 )
1169 > \].
1170 > Concatenating flows corresponding to these five terms, we can obtain
1171 > an symplectic integrator,
1172 > \[
1173 > \varphi _{\Delta t,T^r }  = \varphi _{\Delta t/2,\pi _1 }  \circ
1174 > \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }
1175 > \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi
1176 > _1 }.
1177 > \]
1178 >
1179 > The non-canonical Lie-Poisson bracket ${F, G}$ of two function
1180 > $F(\pi )$ and $G(\pi )$ is defined by
1181 > \[
1182 > \{ F,G\} (\pi ) = [\nabla _\pi  F(\pi )]^T J(\pi )\nabla _\pi  G(\pi
1183 > )
1184 > \]
1185 > If the Poisson bracket of a function $F$ with an arbitrary smooth
1186 > function $G$ is zero, $F$ is a \emph{Casimir}, which is the
1187 > conserved quantity in Poisson system. We can easily verify that the
1188 > norm of the angular momentum, $\parallel \pi
1189 > \parallel$, is a \emph{Casimir}. Let$ F(\pi ) = S(\frac{{\parallel
1190 > \pi \parallel ^2 }}{2})$ for an arbitrary function $ S:R \to R$ ,
1191 > then by the chain rule
1192 > \[
1193 > \nabla _\pi  F(\pi ) = S'(\frac{{\parallel \pi \parallel ^2
1194 > }}{2})\pi
1195 > \]
1196 > Thus $ [\nabla _\pi  F(\pi )]^T J(\pi ) =  - S'(\frac{{\parallel \pi
1197 > \parallel ^2 }}{2})\pi  \times \pi  = 0 $. This explicit
1198 > Lie-Poisson integrator is found to be extremely efficient and stable
1199 > which can be explained by the fact the small angle approximation is
1200 > used and the norm of the angular momentum is conserved.
1201 >
1202 > \subsection{\label{introSection:RBHamiltonianSplitting} Hamiltonian
1203 > Splitting for Rigid Body}
1204 >
1205 > The Hamiltonian of rigid body can be separated in terms of kinetic
1206 > energy and potential energy,
1207 > \[
1208 > H = T(p,\pi ) + V(q,Q)
1209 > \]
1210 > The equations of motion corresponding to potential energy and
1211 > kinetic energy are listed in the below table,
1212 > \begin{center}
1213 > \begin{tabular}{|l|l|}
1214 >  \hline
1215 >  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
1216 >  Potential & Kinetic \\
1217 >  $\frac{{dq}}{{dt}} = \frac{p}{m}$ & $\frac{d}{{dt}}q = p$ \\
1218 >  $\frac{d}{{dt}}p =  - \frac{{\partial V}}{{\partial q}}$ & $ \frac{d}{{dt}}p = 0$ \\
1219 >  $\frac{d}{{dt}}Q = 0$ & $ \frac{d}{{dt}}Q = Qskew(I^{ - 1} j)$ \\
1220 >  $ \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$\\
1221 >  \hline
1222 > \end{tabular}
1223 > \end{center}
1224 > A second-order symplectic method is now obtained by the composition
1225 > of the flow maps,
1226 > \[
1227 > \varphi _{\Delta t}  = \varphi _{\Delta t/2,V}  \circ \varphi
1228 > _{\Delta t,T}  \circ \varphi _{\Delta t/2,V}.
1229 > \]
1230 > Moreover, $\varphi _{\Delta t/2,V}$ can be divided into two
1231 > sub-flows which corresponding to force and torque respectively,
1232 > \[
1233 > \varphi _{\Delta t/2,V}  = \varphi _{\Delta t/2,F}  \circ \varphi
1234 > _{\Delta t/2,\tau }.
1235 > \]
1236 > Since the associated operators of $\varphi _{\Delta t/2,F} $ and
1237 > $\circ \varphi _{\Delta t/2,\tau }$ are commuted, the composition
1238 > order inside $\varphi _{\Delta t/2,V}$ does not matter.
1239 >
1240 > Furthermore, kinetic potential can be separated to translational
1241 > kinetic term, $T^t (p)$, and rotational kinetic term, $T^r (\pi )$,
1242 > \begin{equation}
1243 > T(p,\pi ) =T^t (p) + T^r (\pi ).
1244 > \end{equation}
1245 > where $ T^t (p) = \frac{1}{2}p^T m^{ - 1} p $ and $T^r (\pi )$ is
1246 > defined by \ref{introEquation:rotationalKineticRB}. Therefore, the
1247 > corresponding flow maps are given by
1248 > \[
1249 > \varphi _{\Delta t,T}  = \varphi _{\Delta t,T^t }  \circ \varphi
1250 > _{\Delta t,T^r }.
1251 > \]
1252 > Finally, we obtain the overall symplectic flow maps for free moving
1253 > rigid body
1254 > \begin{equation}
1255 > \begin{array}{c}
1256 > \varphi _{\Delta t}  = \varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau }  \\
1257 >  \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 }  \\
1258 >  \circ \varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .\\
1259 > \end{array}
1260 > \label{introEquation:overallRBFlowMaps}
1261 > \end{equation}
1262 >
1263 > \section{\label{introSection:langevinDynamics}Langevin Dynamics}
1264 > As an alternative to newtonian dynamics, Langevin dynamics, which
1265 > mimics a simple heat bath with stochastic and dissipative forces,
1266 > has been applied in a variety of studies. This section will review
1267 > the theory of Langevin dynamics simulation. A brief derivation of
1268 > generalized Langevin equation will be given first. Follow that, we
1269 > will discuss the physical meaning of the terms appearing in the
1270 > equation as well as the calculation of friction tensor from
1271 > hydrodynamics theory.
1272 >
1273 > \subsection{\label{introSection:generalizedLangevinDynamics}Derivation of Generalized Langevin Equation}
1274 >
1275 > Harmonic bath model, in which an effective set of harmonic
1276 > oscillators are used to mimic the effect of a linearly responding
1277 > environment, has been widely used in quantum chemistry and
1278 > statistical mechanics. One of the successful applications of
1279 > Harmonic bath model is the derivation of Deriving Generalized
1280 > Langevin Dynamics. Lets consider a system, in which the degree of
1281 > freedom $x$ is assumed to couple to the bath linearly, giving a
1282 > Hamiltonian of the form
1283 > \begin{equation}
1284 > H = \frac{{p^2 }}{{2m}} + U(x) + H_B  + \Delta U(x,x_1 , \ldots x_N)
1285 > \label{introEquation:bathGLE}.
1286 > \end{equation}
1287 > Here $p$ is a momentum conjugate to $q$, $m$ is the mass associated
1288 > with this degree of freedom, $H_B$ is harmonic bath Hamiltonian,
1289 > \[
1290 > H_B  = \sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
1291 > }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  \omega _\alpha ^2 }
1292 > \right\}}
1293 > \]
1294 > where the index $\alpha$ runs over all the bath degrees of freedom,
1295 > $\omega _\alpha$ are the harmonic bath frequencies, $m_\alpha$ are
1296 > the harmonic bath masses, and $\Delta U$ is bilinear system-bath
1297 > coupling,
1298 > \[
1299 > \Delta U =  - \sum\limits_{\alpha  = 1}^N {g_\alpha  x_\alpha  x}
1300 > \]
1301 > where $g_\alpha$ are the coupling constants between the bath and the
1302 > coordinate $x$. Introducing
1303 > \[
1304 > W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
1305 > }}{{2m_\alpha  w_\alpha ^2 }}} x^2
1306 > \] and combining the last two terms in Equation
1307 > \ref{introEquation:bathGLE}, we may rewrite the Harmonic bath
1308 > Hamiltonian as
1309 > \[
1310 > H = \frac{{p^2 }}{{2m}} + W(x) + \sum\limits_{\alpha  = 1}^N
1311 > {\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
1312 > w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
1313 > w_\alpha ^2 }}x} \right)^2 } \right\}}
1314 > \]
1315 > Since the first two terms of the new Hamiltonian depend only on the
1316 > system coordinates, we can get the equations of motion for
1317 > Generalized Langevin Dynamics by Hamilton's equations
1318 > \ref{introEquation:motionHamiltonianCoordinate,
1319 > introEquation:motionHamiltonianMomentum},
1320 > \begin{equation}
1321 > m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} -
1322 > \sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   -
1323 > \frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right)},
1324 > \label{introEquation:coorMotionGLE}
1325 > \end{equation}
1326 > and
1327 > \begin{equation}
1328 > m\ddot x_\alpha   =  - m_\alpha  w_\alpha ^2 \left( {x_\alpha   -
1329 > \frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right).
1330 > \label{introEquation:bathMotionGLE}
1331 > \end{equation}
1332 >
1333 > In order to derive an equation for $x$, the dynamics of the bath
1334 > variables $x_\alpha$ must be solved exactly first. As an integral
1335 > transform which is particularly useful in solving linear ordinary
1336 > differential equations, Laplace transform is the appropriate tool to
1337 > solve this problem. The basic idea is to transform the difficult
1338 > differential equations into simple algebra problems which can be
1339 > solved easily. Then applying inverse Laplace transform, also known
1340 > as the Bromwich integral, we can retrieve the solutions of the
1341 > original problems.
1342 >
1343 > Let $f(t)$ be a function defined on $ [0,\infty ) $. The Laplace
1344 > transform of f(t) is a new function defined as
1345 > \[
1346 > L(f(t)) \equiv F(p) = \int_0^\infty  {f(t)e^{ - pt} dt}
1347 > \]
1348 > where  $p$ is real and  $L$ is called the Laplace Transform
1349 > Operator. Below are some important properties of Laplace transform
1350 > \begin{equation}
1351 > \begin{array}{c}
1352 > L(x + y) = L(x) + L(y) \\
1353 > L(ax) = aL(x) \\
1354 > L(\dot x) = pL(x) - px(0) \\
1355 > L(\ddot x) = p^2 L(x) - px(0) - \dot x(0) \\
1356 > L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right) = G(p)H(p) \\
1357 > \end{array}
1358 > \end{equation}
1359 >
1360 > Applying Laplace transform to the bath coordinates, we obtain
1361 > \[
1362 > \begin{array}{c}
1363 > 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) \\
1364 > L(x_\alpha  ) = \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }} \\
1365 > \end{array}
1366 > \]
1367 > By the same way, the system coordinates become
1368 > \[
1369 > \begin{array}{c}
1370 > mL(\ddot x) =  - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} \\
1371 >  - \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\}}  \\
1372 > \end{array}
1373 > \]
1374 >
1375 > With the help of some relatively important inverse Laplace
1376 > transformations:
1377 > \[
1378 > \begin{array}{c}
1379 > L(\cos at) = \frac{p}{{p^2  + a^2 }} \\
1380 > L(\sin at) = \frac{a}{{p^2  + a^2 }} \\
1381 > L(1) = \frac{1}{p} \\
1382 > \end{array}
1383 > \]
1384 > , we obtain
1385 > \begin{align}
1386 > m\ddot x &=  - \frac{{\partial W(x)}}{{\partial x}} -
1387 > \sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
1388 > }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\int_0^t {\cos (\omega
1389 > _\alpha  t)\dot x(t - \tau )d\tau  - \left[ {g_\alpha  x_\alpha  (0)
1390 > - \frac{{g_\alpha  }}{{m_\alpha  \omega _\alpha  }}} \right]\cos
1391 > (\omega _\alpha  t) - \frac{{g_\alpha  \dot x_\alpha  (0)}}{{\omega
1392 > _\alpha  }}\sin (\omega _\alpha  t)} } \right\}}
1393 > %
1394 > &= - \frac{{\partial W(x)}}{{\partial x}} - \int_0^t
1395 > {\sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
1396 > }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha
1397 > t)\dot x(t - \tau )d} \tau }  + \sum\limits_{\alpha  = 1}^N {\left\{
1398 > {\left[ {g_\alpha  x_\alpha  (0) - \frac{{g_\alpha  }}{{m_\alpha
1399 > \omega _\alpha  }}} \right]\cos (\omega _\alpha  t) +
1400 > \frac{{g_\alpha  \dot x_\alpha  (0)}}{{\omega _\alpha  }}\sin
1401 > (\omega _\alpha  t)} \right\}}
1402 > \end{align}
1403 >
1404 > Introducing a \emph{dynamic friction kernel}
1405 > \begin{equation}
1406 > \xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
1407 > }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
1408 > \label{introEquation:dynamicFrictionKernelDefinition}
1409 > \end{equation}
1410 > and \emph{a random force}
1411 > \begin{equation}
1412 > R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
1413 > - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
1414 > \right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
1415 > (0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t),
1416 > \label{introEquation:randomForceDefinition}
1417 > \end{equation}
1418 > the equation of motion can be rewritten as
1419 > \begin{equation}
1420 > m\ddot x =  - \frac{{\partial W}}{{\partial x}} - \int_0^t {\xi
1421 > (t)\dot x(t - \tau )d\tau }  + R(t)
1422 > \label{introEuqation:GeneralizedLangevinDynamics}
1423 > \end{equation}
1424 > which is known as the \emph{generalized Langevin equation}.
1425 >
1426 > \subsubsection{\label{introSection:randomForceDynamicFrictionKernel}Random Force and Dynamic Friction Kernel}
1427 >
1428 > One may notice that $R(t)$ depends only on initial conditions, which
1429 > implies it is completely deterministic within the context of a
1430 > harmonic bath. However, it is easy to verify that $R(t)$ is totally
1431 > uncorrelated to $x$ and $\dot x$,
1432 > \[
1433 > \begin{array}{l}
1434 > \left\langle {x(t)R(t)} \right\rangle  = 0, \\
1435 > \left\langle {\dot x(t)R(t)} \right\rangle  = 0. \\
1436 > \end{array}
1437 > \]
1438 > This property is what we expect from a truly random process. As long
1439 > as the model, which is gaussian distribution in general, chosen for
1440 > $R(t)$ is a truly random process, the stochastic nature of the GLE
1441 > still remains.
1442 >
1443 > %dynamic friction kernel
1444 > The convolution integral
1445 > \[
1446 > \int_0^t {\xi (t)\dot x(t - \tau )d\tau }
1447 > \]
1448 > depends on the entire history of the evolution of $x$, which implies
1449 > that the bath retains memory of previous motions. In other words,
1450 > the bath requires a finite time to respond to change in the motion
1451 > of the system. For a sluggish bath which responds slowly to changes
1452 > in the system coordinate, we may regard $\xi(t)$ as a constant
1453 > $\xi(t) = \Xi_0$. Hence, the convolution integral becomes
1454 > \[
1455 > \int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = \xi _0 (x(t) - x(0))
1456 > \]
1457 > and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
1458 > \[
1459 > m\ddot x =  - \frac{\partial }{{\partial x}}\left( {W(x) +
1460 > \frac{1}{2}\xi _0 (x - x_0 )^2 } \right) + R(t),
1461 > \]
1462 > which can be used to describe dynamic caging effect. The other
1463 > extreme is the bath that responds infinitely quickly to motions in
1464 > the system. Thus, $\xi (t)$ can be taken as a $delta$ function in
1465 > time:
1466 > \[
1467 > \xi (t) = 2\xi _0 \delta (t)
1468 > \]
1469 > Hence, the convolution integral becomes
1470 > \[
1471 > \int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = 2\xi _0 \int_0^t
1472 > {\delta (t)\dot x(t - \tau )d\tau }  = \xi _0 \dot x(t),
1473 > \]
1474 > and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
1475 > \begin{equation}
1476 > m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} - \xi _0 \dot
1477 > x(t) + R(t) \label{introEquation:LangevinEquation}
1478 > \end{equation}
1479 > which is known as the Langevin equation. The static friction
1480 > coefficient $\xi _0$ can either be calculated from spectral density
1481 > or be determined by Stokes' law for regular shaped particles.A
1482 > briefly review on calculating friction tensor for arbitrary shaped
1483 > particles is given in section \ref{introSection:frictionTensor}.
1484 >
1485 > \subsubsection{\label{introSection:secondFluctuationDissipation}The Second Fluctuation Dissipation Theorem}
1486 >
1487 > Defining a new set of coordinates,
1488 > \[
1489 > q_\alpha  (t) = x_\alpha  (t) - \frac{1}{{m_\alpha  \omega _\alpha
1490 > ^2 }}x(0)
1491 > \],
1492 > we can rewrite $R(T)$ as
1493 > \[
1494 > R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}.
1495 > \]
1496 > And since the $q$ coordinates are harmonic oscillators,
1497 > \[
1498 > \begin{array}{c}
1499 > \left\langle {q_\alpha ^2 } \right\rangle  = \frac{{kT}}{{m_\alpha  \omega _\alpha ^2 }} \\
1500 > \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  = \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
1501 > \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle  = \delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
1502 > \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 } }  \\
1503 >  = \sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t)}  \\
1504 >  = kT\xi (t) \\
1505 > \end{array}
1506 > \]
1507 > Thus, we recover the \emph{second fluctuation dissipation theorem}
1508 > \begin{equation}
1509 > \xi (t) = \left\langle {R(t)R(0)} \right\rangle
1510 > \label{introEquation:secondFluctuationDissipation}.
1511 > \end{equation}
1512 > In effect, it acts as a constraint on the possible ways in which one
1513 > can model the random force and friction kernel.
1514 >
1515 > \subsection{\label{introSection:frictionTensor} Friction Tensor}
1516 > Theoretically, the friction kernel can be determined using velocity
1517 > autocorrelation function. However, this approach become impractical
1518 > when the system become more and more complicate. Instead, various
1519 > approaches based on hydrodynamics have been developed to calculate
1520 > the friction coefficients. The friction effect is isotropic in
1521 > Equation, \zeta can be taken as a scalar. In general, friction
1522 > tensor \Xi is a $6\times 6$ matrix given by
1523 > \[
1524 > \Xi  = \left( {\begin{array}{*{20}c}
1525 >   {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
1526 >   {\Xi _{}^{tr} } & {\Xi _{}^{rr} }  \\
1527 > \end{array}} \right).
1528 > \]
1529 > Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are translational friction
1530 > tensor and rotational resistance (friction) tensor respectively,
1531 > while ${\Xi^{tr} }$ is translation-rotation coupling tensor and $
1532 > {\Xi^{rt} }$ is rotation-translation coupling tensor. When a
1533 > particle moves in a fluid, it may experience friction force or
1534 > torque along the opposite direction of the velocity or angular
1535 > velocity,
1536 > \[
1537 > \left( \begin{array}{l}
1538 > F_R  \\
1539 > \tau _R  \\
1540 > \end{array} \right) =  - \left( {\begin{array}{*{20}c}
1541 >   {\Xi ^{tt} } & {\Xi ^{rt} }  \\
1542 >   {\Xi ^{tr} } & {\Xi ^{rr} }  \\
1543 > \end{array}} \right)\left( \begin{array}{l}
1544 > v \\
1545 > w \\
1546 > \end{array} \right)
1547 > \]
1548 > where $F_r$ is the friction force and $\tau _R$ is the friction
1549 > toque.
1550 >
1551 > \subsubsection{\label{introSection:resistanceTensorRegular}The Resistance Tensor for Regular Shape}
1552 >
1553 > For a spherical particle, the translational and rotational friction
1554 > constant can be calculated from Stoke's law,
1555 > \[
1556 > \Xi ^{tt}  = \left( {\begin{array}{*{20}c}
1557 >   {6\pi \eta R} & 0 & 0  \\
1558 >   0 & {6\pi \eta R} & 0  \\
1559 >   0 & 0 & {6\pi \eta R}  \\
1560 > \end{array}} \right)
1561 > \]
1562 > and
1563 > \[
1564 > \Xi ^{rr}  = \left( {\begin{array}{*{20}c}
1565 >   {8\pi \eta R^3 } & 0 & 0  \\
1566 >   0 & {8\pi \eta R^3 } & 0  \\
1567 >   0 & 0 & {8\pi \eta R^3 }  \\
1568 > \end{array}} \right)
1569 > \]
1570 > where $\eta$ is the viscosity of the solvent and $R$ is the
1571 > hydrodynamics radius.
1572 >
1573 > Other non-spherical shape, such as cylinder and ellipsoid
1574 > \textit{etc}, are widely used as reference for developing new
1575 > hydrodynamics theory, because their properties can be calculated
1576 > exactly. In 1936, Perrin extended Stokes's law to general ellipsoid,
1577 > also called a triaxial ellipsoid, which is given in Cartesian
1578 > coordinates by
1579 > \[
1580 > \frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2
1581 > }} = 1
1582 > \]
1583 > where the semi-axes are of lengths $a$, $b$, and $c$. Unfortunately,
1584 > due to the complexity of the elliptic integral, only the ellipsoid
1585 > with the restriction of two axes having to be equal, \textit{i.e.}
1586 > prolate($ a \ge b = c$) and oblate ($ a < b = c $), can be solved
1587 > exactly. Introducing an elliptic integral parameter $S$ for prolate,
1588 > \[
1589 > S = \frac{2}{{\sqrt {a^2  - b^2 } }}\ln \frac{{a + \sqrt {a^2  - b^2
1590 > } }}{b},
1591 > \]
1592 > and oblate,
1593 > \[
1594 > S = \frac{2}{{\sqrt {b^2  - a^2 } }}arctg\frac{{\sqrt {b^2  - a^2 }
1595 > }}{a}
1596 > \],
1597 > one can write down the translational and rotational resistance
1598 > tensors
1599 > \[
1600 > \begin{array}{l}
1601 > \Xi _a^{tt}  = 16\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - b^2 )S - 2a}} \\
1602 > \Xi _b^{tt}  = \Xi _c^{tt}  = 32\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - 3b^2 )S + 2a}} \\
1603 > \end{array},
1604 > \]
1605 > and
1606 > \[
1607 > \begin{array}{l}
1608 > \Xi _a^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^2  - b^2 )b^2 }}{{2a - b^2 S}} \\
1609 > \Xi _b^{rr}  = \Xi _c^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^4  - b^4 )}}{{(2a^2  - b^2 )S - 2a}} \\
1610 > \end{array}.
1611 > \]
1612 >
1613 > \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}The Resistance Tensor for Arbitrary Shape}
1614 >
1615 > Unlike spherical and other regular shaped molecules, there is not
1616 > analytical solution for friction tensor of any arbitrary shaped
1617 > rigid molecules. The ellipsoid of revolution model and general
1618 > triaxial ellipsoid model have been used to approximate the
1619 > hydrodynamic properties of rigid bodies. However, since the mapping
1620 > from all possible ellipsoidal space, $r$-space, to all possible
1621 > combination of rotational diffusion coefficients, $D$-space is not
1622 > unique\cite{Wegener79} as well as the intrinsic coupling between
1623 > translational and rotational motion of rigid body\cite{}, general
1624 > ellipsoid is not always suitable for modeling arbitrarily shaped
1625 > rigid molecule. A number of studies have been devoted to determine
1626 > the friction tensor for irregularly shaped rigid bodies using more
1627 > advanced method\cite{} where the molecule of interest was modeled by
1628 > combinations of spheres(beads)\cite{} and the hydrodynamics
1629 > properties of the molecule can be calculated using the hydrodynamic
1630 > interaction tensor. Let us consider a rigid assembly of $N$ beads
1631 > immersed in a continuous medium. Due to hydrodynamics interaction,
1632 > the ``net'' velocity of $i$th bead, $v'_i$ is different than its
1633 > unperturbed velocity $v_i$,
1634 > \[
1635 > v'_i  = v_i  - \sum\limits_{j \ne i} {T_{ij} F_j }
1636 > \]
1637 > where $F_i$ is the frictional force, and $T_{ij}$ is the
1638 > hydrodynamic interaction tensor. The friction force of $i$th bead is
1639 > proportional to its ``net'' velocity
1640 > \begin{equation}
1641 > F_i  = \zeta _i v_i  - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }.
1642 > \label{introEquation:tensorExpression}
1643 > \end{equation}
1644 > This equation is the basis for deriving the hydrodynamic tensor. In
1645 > 1930, Oseen and Burgers gave a simple solution to Equation
1646 > \ref{introEquation:tensorExpression}
1647 > \begin{equation}
1648 > T_{ij}  = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij}
1649 > R_{ij}^T }}{{R_{ij}^2 }}} \right).
1650 > \label{introEquation:oseenTensor}
1651 > \end{equation}
1652 > Here $R_{ij}$ is the distance vector between bead $i$ and bead $j$.
1653 > A second order expression for element of different size was
1654 > introduced by Rotne and Prager\cite{} and improved by Garc\'{i}a de
1655 > la Torre and Bloomfield,
1656 > \begin{equation}
1657 > T_{ij}  = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {I +
1658 > \frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right) + R\frac{{\sigma
1659 > _i^2  + \sigma _j^2 }}{{r_{ij}^2 }}\left( {\frac{I}{3} -
1660 > \frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right].
1661 > \label{introEquation:RPTensorNonOverlapped}
1662 > \end{equation}
1663 > Both of the Equation \ref{introEquation:oseenTensor} and Equation
1664 > \ref{introEquation:RPTensorNonOverlapped} have an assumption $R_{ij}
1665 > \ge \sigma _i  + \sigma _j$. An alternative expression for
1666 > overlapping beads with the same radius, $\sigma$, is given by
1667 > \begin{equation}
1668 > T_{ij}  = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 -
1669 > \frac{2}{{32}}\frac{{R_{ij} }}{\sigma }} \right)I +
1670 > \frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right]
1671 > \label{introEquation:RPTensorOverlapped}
1672 > \end{equation}
1673 >
1674 > To calculate the resistance tensor at an arbitrary origin $O$, we
1675 > construct a $3N \times 3N$ matrix consisting of $N \times N$
1676 > $B_{ij}$ blocks
1677 > \begin{equation}
1678 > B = \left( {\begin{array}{*{20}c}
1679 >   {B_{11} } &  \ldots  & {B_{1N} }  \\
1680 >    \vdots  &  \ddots  &  \vdots   \\
1681 >   {B_{N1} } &  \cdots  & {B_{NN} }  \\
1682 > \end{array}} \right),
1683 > \end{equation}
1684 > where $B_{ij}$ is given by
1685 > \[
1686 > B_{ij}  = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
1687 > )T_{ij}
1688 > \]
1689 > where $\delta _{ij}$ is Kronecker delta function. Inverting matrix
1690 > $B$, we obtain
1691 >
1692 > \[
1693 > C = B^{ - 1}  = \left( {\begin{array}{*{20}c}
1694 >   {C_{11} } &  \ldots  & {C_{1N} }  \\
1695 >    \vdots  &  \ddots  &  \vdots   \\
1696 >   {C_{N1} } &  \cdots  & {C_{NN} }  \\
1697 > \end{array}} \right)
1698 > \]
1699 > , which can be partitioned into $N \times N$ $3 \times 3$ block
1700 > $C_{ij}$. With the help of $C_{ij}$ and skew matrix $U_i$
1701 > \[
1702 > U_i  = \left( {\begin{array}{*{20}c}
1703 >   0 & { - z_i } & {y_i }  \\
1704 >   {z_i } & 0 & { - x_i }  \\
1705 >   { - y_i } & {x_i } & 0  \\
1706 > \end{array}} \right)
1707 > \]
1708 > where $x_i$, $y_i$, $z_i$ are the components of the vector joining
1709 > bead $i$ and origin $O$. Hence, the elements of resistance tensor at
1710 > arbitrary origin $O$ can be written as
1711 > \begin{equation}
1712 > \begin{array}{l}
1713 > \Xi _{}^{tt}  = \sum\limits_i {\sum\limits_j {C_{ij} } } , \\
1714 > \Xi _{}^{tr}  = \Xi _{}^{rt}  = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\
1715 > \Xi _{}^{rr}  =  - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j  \\
1716 > \end{array}
1717 > \label{introEquation:ResistanceTensorArbitraryOrigin}
1718 > \end{equation}
1719 >
1720 > The resistance tensor depends on the origin to which they refer. The
1721 > proper location for applying friction force is the center of
1722 > resistance (reaction), at which the trace of rotational resistance
1723 > tensor, $ \Xi ^{rr}$ reaches minimum. Mathematically, the center of
1724 > resistance is defined as an unique point of the rigid body at which
1725 > the translation-rotation coupling tensor are symmetric,
1726 > \begin{equation}
1727 > \Xi^{tr}  = \left( {\Xi^{tr} } \right)^T
1728 > \label{introEquation:definitionCR}
1729 > \end{equation}
1730 > Form Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
1731 > we can easily find out that the translational resistance tensor is
1732 > origin independent, while the rotational resistance tensor and
1733 > translation-rotation coupling resistance tensor depend on the
1734 > origin. Given resistance tensor at an arbitrary origin $O$, and a
1735 > vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
1736 > obtain the resistance tensor at $P$ by
1737 > \begin{equation}
1738 > \begin{array}{l}
1739 > \Xi _P^{tt}  = \Xi _O^{tt}  \\
1740 > \Xi _P^{tr}  = \Xi _P^{rt}  = \Xi _O^{tr}  - U_{OP} \Xi _O^{tt}  \\
1741 > \Xi _P^{rr}  = \Xi _O^{rr}  - U_{OP} \Xi _O^{tt} U_{OP}  + \Xi _O^{tr} U_{OP}  - U_{OP} \Xi _O^{tr} ^{^T }  \\
1742 > \end{array}
1743 > \label{introEquation:resistanceTensorTransformation}
1744 > \end{equation}
1745 > where
1746 > \[
1747 > U_{OP}  = \left( {\begin{array}{*{20}c}
1748 >   0 & { - z_{OP} } & {y_{OP} }  \\
1749 >   {z_i } & 0 & { - x_{OP} }  \\
1750 >   { - y_{OP} } & {x_{OP} } & 0  \\
1751 > \end{array}} \right)
1752 > \]
1753 > Using Equations \ref{introEquation:definitionCR} and
1754 > \ref{introEquation:resistanceTensorTransformation}, one can locate
1755 > the position of center of resistance,
1756 > \[
1757 > \left( \begin{array}{l}
1758 > x_{OR}  \\
1759 > y_{OR}  \\
1760 > z_{OR}  \\
1761 > \end{array} \right) = \left( {\begin{array}{*{20}c}
1762 >   {(\Xi _O^{rr} )_{yy}  + (\Xi _O^{rr} )_{zz} } & { - (\Xi _O^{rr} )_{xy} } & { - (\Xi _O^{rr} )_{xz} }  \\
1763 >   { - (\Xi _O^{rr} )_{xy} } & {(\Xi _O^{rr} )_{zz}  + (\Xi _O^{rr} )_{xx} } & { - (\Xi _O^{rr} )_{yz} }  \\
1764 >   { - (\Xi _O^{rr} )_{xz} } & { - (\Xi _O^{rr} )_{yz} } & {(\Xi _O^{rr} )_{xx}  + (\Xi _O^{rr} )_{yy} }  \\
1765 > \end{array}} \right)^{ - 1} \left( \begin{array}{l}
1766 > (\Xi _O^{tr} )_{yz}  - (\Xi _O^{tr} )_{zy}  \\
1767 > (\Xi _O^{tr} )_{zx}  - (\Xi _O^{tr} )_{xz}  \\
1768 > (\Xi _O^{tr} )_{xy}  - (\Xi _O^{tr} )_{yx}  \\
1769 > \end{array} \right).
1770 > \]
1771 > where $x_OR$, $y_OR$, $z_OR$ are the components of the vector
1772 > joining center of resistance $R$ and origin $O$.
1773 >
1774 > %\section{\label{introSection:correlationFunctions}Correlation Functions}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines