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 2713 by tim, Sat Apr 15 01:08:18 2006 UTC vs.
Revision 2719 by tim, Tue Apr 18 22:38:19 2006 UTC

# Line 822 | Line 822 | q(\Delta t)} \right]. %
822   %
823   q(\Delta t) &= q(0) + \frac{{\Delta t}}{2}\left[ {\dot q(0) + \dot
824   q(\Delta t)} \right]. %
825 < \label{introEquation:positionVerlet1}
825 > \label{introEquation:positionVerlet2}
826   \end{align}
827  
828   \subsubsection{\label{introSection:errorAnalysis}Error Analysis and Higher Order Methods}
# Line 887 | Line 887 | dynamical information.
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}
# Line 927 | Line 941 | rotation matrix $A$ and re-formulating Hamiltonian's e
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 $A$ and re-formulating Hamiltonian's equation, a
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 $A_t A = 1$. An alternative
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
# Line 1136 | Line 1150 | To reduce the cost of computing expensive functions in
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,
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   )
# Line 1213 | Line 1227 | Moreover, \varphi _{\Delta t/2,V} can be divided into
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 sub-flows
1231 < which corresponding to force and torque respectively,
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.
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 )$,
# Line 1247 | Line 1261 | rigid body
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:LDIntroduction}Introduction and application of Langevin Dynamics}
1273 > \subsection{\label{introSection:generalizedLangevinDynamics}Derivation of Generalized Langevin Equation}
1274  
1275 < \subsection{\label{introSection:generalizedLangevinDynamics}Generalized Langevin Dynamics}
1276 <
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}
1285 > \label{introEquation:bathGLE}.
1286   \end{equation}
1287 < where $H_B$ is harmonic bath Hamiltonian,
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  w_\alpha ^2 } \right\}}
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 < and $\Delta U$ is bilinear system-bath coupling,
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 < Completing the square,
1301 > where $g_\alpha$ are the coupling constants between the bath and the
1302 > coordinate $x$. Introducing
1303   \[
1304 < H_B  + \Delta U = \sum\limits_{\alpha  = 1}^N {\left\{
1305 < {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
1306 < w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
1307 < w_\alpha ^2 }}x} \right)^2 } \right\}}  - \sum\limits_{\alpha  =
1308 < 1}^N {\frac{{g_\alpha ^2 }}{{2m_\alpha  w_\alpha ^2 }}} x^2
1275 < \]
1276 < and putting it back into Eq.~\ref{introEquation:bathGLE},
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\}}
1282 \]
1283 where
1284 \[
1285 W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
1286 }}{{2m_\alpha  w_\alpha ^2 }}} x^2
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{align}
1321 < \dot p &=  - \frac{{\partial H}}{{\partial x}}
1322 <       &= m\ddot x
1323 <       &= - \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)}
1324 < \label{introEquation:Lp5}
1325 < \end{align}
1326 < , and
1327 < \begin{align}
1328 < \dot p_\alpha   &=  - \frac{{\partial H}}{{\partial x_\alpha  }}
1329 <                &= m\ddot x_\alpha
1330 <                &= \- m_\alpha  w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha}}{{m_\alpha  w_\alpha ^2 }}x} \right)
1331 < \end{align}
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 < \subsection{\label{introSection:laplaceTransform}The Laplace Transform}
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(x) = \int_0^\infty  {x(t)e^{ - pt} dt}
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 < L(x + y) = L(x) + L(y)
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 <
1367 > By the same way, the system coordinates become
1368   \[
1369 < L(ax) = aL(x)
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 < L(\dot x) = pL(x) - px(0)
1379 < \]
1380 <
1381 < \[
1382 < L(\ddot x) = p^2 L(x) - px(0) - \dot x(0)
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 <
1328 < \[
1329 < L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right) = G(p)H(p)
1330 < \]
1331 <
1332 < Some relatively important transformation,
1333 < \[
1334 < L(\cos at) = \frac{p}{{p^2  + a^2 }}
1335 < \]
1336 <
1337 < \[
1338 < L(\sin at) = \frac{a}{{p^2  + a^2 }}
1339 < \]
1340 <
1341 < \[
1342 < L(1) = \frac{1}{p}
1343 < \]
1344 <
1345 < First, the bath coordinates,
1346 < \[
1347 < p^2 L(x_\alpha  ) - px_\alpha  (0) - \dot x_\alpha  (0) =  - \omega
1348 < _\alpha ^2 L(x_\alpha  ) + \frac{{g_\alpha  }}{{\omega _\alpha
1349 < }}L(x)
1350 < \]
1351 < \[
1352 < L(x_\alpha  ) = \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) +
1353 < px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }}
1354 < \]
1355 < Then, the system coordinates,
1384 > , we obtain
1385   \begin{align}
1357 mL(\ddot x) &=  - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
1358 \sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{\frac{{g_\alpha
1359 }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha
1360 (0)}}{{p^2  + \omega _\alpha ^2 }} - \frac{{g_\alpha ^2 }}{{m_\alpha
1361 }}\omega _\alpha ^2 L(x)} \right\}}
1362 %
1363 &= - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
1364 \sum\limits_{\alpha  = 1}^N {\left\{ { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}\frac{p}{{p^2  + \omega _\alpha ^2 }}pL(x)
1365 - \frac{p}{{p^2  + \omega _\alpha ^2 }}g_\alpha  x_\alpha  (0)
1366 - \frac{1}{{p^2  + \omega _\alpha ^2 }}g_\alpha  \dot x_\alpha  (0)} \right\}}
1367 \end{align}
1368 Then, the inverse transform,
1369
1370 \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
# Line 1386 | Line 1401 | t)\dot x(t - \tau )d} \tau }  + \sum\limits_{\alpha  =
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 < %where $ {\xi (t)}$ is friction kernel, $R(t)$ is random force and
1425 < %$W$ is the potential of mean force. $W(x) =  - kT\ln p(x)$
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 < \xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
1434 < }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
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 < For an infinite harmonic bath, we can use the spectral density and
1439 < an integral over frequencies.
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 < R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
1405 < - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
1406 < \right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
1407 < (0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)
1446 > \int_0^t {\xi (t)\dot x(t - \tau )d\tau }
1447   \]
1448 < The random forces depend only on initial conditions.
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 < So we can define a new set of coordinates,
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 < This makes
1491 > \],
1492 > we can rewrite $R(T)$ as
1493   \[
1494 < R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}
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}{l}
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 <
1429 < \begin{align}
1430 < \left\langle {R(t)R(0)} \right\rangle  &= \sum\limits_\alpha
1431 < {\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha
1432 < (t)q_\beta  (0)} \right\rangle } }
1433 < %
1434 < &= \sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)}
1435 < \right\rangle \cos (\omega _\alpha  t)}
1436 < %
1437 < &= kT\xi (t)
1438 < \end{align}
1439 <
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}
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  
1445 \section{\label{introSection:hydroynamics}Hydrodynamics}
1446
1515   \subsection{\label{introSection:frictionTensor} Friction Tensor}
1516 < \subsection{\label{introSection:analyticalApproach}Analytical
1517 < Approach}
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 < \subsection{\label{introSection:approximationApproach}Approximation
1452 < Approach}
1551 > \subsubsection{\label{introSection:resistanceTensorRegular}The Resistance Tensor for Regular Shape}
1552  
1553 < \subsection{\label{introSection:centersRigidBody}Centers of Rigid
1554 < Body}
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 < \section{\label{introSection:correlationFunctions}Correlation Functions}
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