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

Comparing trunk/oopseDocs/oopseDoc.tex (file contents):
Revision 3395 by gezelter, Fri May 16 14:25:26 2008 UTC vs.
Revision 3402 by xsun, Fri May 30 15:58:01 2008 UTC

# Line 2369 | Line 2369 | simulations).
2369   for use with liquid simulations, or in cases where there is
2370   orientational anisotropy in the system (i.e. in lipid bilayer
2371   simulations).
2372 +
2373 + \section{Langevin Dynamics (LANGEVINDYNAMICS)\label{LDRB}}
2374 +
2375 + {\sc oopse} implements a Langevin dynamics integrator in order to
2376 + perform molecular dynamics simulations in implicit solvent
2377 + environments.  This results in substantial performance gains when the
2378 + detailed dynamics of the solvent is not important. Since {\sc oopse}
2379 + also handles rigid bodies of arbitrary composition and shape, the
2380 + Langevin integrator is somewhat more complex than in other simulation
2381 + packages.
2382 +
2383 + Consider the Langevin equations of motion in generalized coordinates
2384 + \begin{equation}
2385 + {\bf M} \dot{{\bf V}}(t) = {\bf F}_{s}(t) +
2386 + {\bf F}_{f}(t)  + {\bf F}_{r}(t)
2387 + \label{LDGeneralizedForm}
2388 + \end{equation}
2389 + where ${\bf M}$ is a $6 \times 6$ diagonal mass matrix (which
2390 + includes the mass of the rigid body as well as the moments of inertia
2391 + in the body-fixed frame) and ${\bf V}$ is a generalized velocity,
2392 + ${\bf V} =
2393 + \left\{{\bf v},{\bf \omega}\right\}$. The right side of
2394 + Eq.~\ref{LDGeneralizedForm} consists of three generalized forces: a
2395 + system force (${\bf F}_{s}$), a frictional or dissipative force (${\bf
2396 + F}_{f}$) and a stochastic force (${\bf F}_{r}$). While the evolution
2397 + of the system in Newtonian mechanics is typically done in the lab
2398 + frame, it is convenient to handle the dynamics of rigid bodies in
2399 + body-fixed frames. Thus the friction and random forces on each
2400 + substructure are calculated in a body-fixed frame and may converted
2401 + back to the lab frame using that substructure's rotation matrix (${\bf
2402 + Q}$):
2403 + \begin{equation}
2404 + {\bf F}_{f,r} =
2405 + \left( \begin{array}{c}
2406 + {\bf f}_{f,r} \\
2407 + {\bf \tau}_{f,r}
2408 + \end{array} \right)
2409 + =
2410 + \left( \begin{array}{c}
2411 + {\bf Q}^{T} {\bf f}^{~b}_{f,r} \\
2412 + {\bf Q}^{T} {\bf \tau}^{~b}_{f,r}
2413 + \end{array} \right)
2414 + \end{equation}
2415 + The body-fixed friction force, ${\bf F}_{f}^{~b}$, is proportional to
2416 + the (body-fixed) velocity at the center of resistance
2417 + ${\bf v}_{R}^{~b}$ and the angular velocity ${\bf \omega}$
2418 + \begin{equation}
2419 + {\bf F}_{f}^{~b}(t) = \left( \begin{array}{l}
2420 + {\bf f}_{f}^{~b}(t) \\
2421 + {\bf \tau}_{f}^{~b}(t) \\
2422 + \end{array} \right) =  - \left( \begin{array}{*{20}c}
2423 +   \Xi_{R}^{tt} & \Xi_{R}^{rt} \\
2424 +   \Xi_{R}^{tr} & \Xi_{R}^{rr}    \\
2425 + \end{array} \right)\left( \begin{array}{l}
2426 + {\bf v}_{R}^{~b}(t) \\
2427 + {\bf \omega}(t) \\
2428 + \end{array} \right),
2429 + \end{equation}
2430 + while the random force, ${\bf F}_{r}$, is a Gaussian stochastic
2431 + variable with zero mean and variance,
2432 + \begin{equation}
2433 + \left\langle {{\bf F}_{r}(t) ({\bf F}_{r}(t'))^T } \right\rangle  =
2434 + \left\langle {{\bf F}_{r}^{~b} (t) ({\bf F}_{r}^{~b} (t'))^T } \right\rangle  =
2435 + 2 k_B T \Xi_R \delta(t - t'). \label{eq:randomForce}
2436 + \end{equation}
2437 + $\Xi_R$ is the $6\times6$ resistance tensor at the center of
2438 + resistance.  Once this tensor is known for a given rigid body (as
2439 + described in the previous section) obtaining a stochastic vector that
2440 + has the properties in Eq. (\ref{eq:randomForce}) can be done
2441 + efficiently by carrying out a one-time Cholesky decomposition to
2442 + obtain the square root matrix of the resistance tensor,
2443 + \begin{equation}
2444 + \Xi_R = {\bf S} {\bf S}^{T},
2445 + \label{eq:Cholesky}
2446 + \end{equation}
2447 + where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A
2448 + vector with the statistics required for the random force can then be
2449 + obtained by multiplying ${\bf S}$ onto a random 6-vector ${\bf Z}$ which
2450 + has elements chosen from a Gaussian distribution, such that:
2451 + \begin{equation}
2452 + \langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot
2453 + {\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij},
2454 + \end{equation}
2455 + where $\delta t$ is the timestep in use during the simulation. The
2456 + random force, ${\bf F}_{r}^{~b} = {\bf S} {\bf Z}$, can be shown to have the
2457 + correct properties required by Eq. (\ref{eq:randomForce}).
2458 +
2459 + The equation of motion for the translational velocity of the center of
2460 + mass (${\bf v}$) can be written as
2461 + \begin{equation}
2462 + m \dot{{\bf v}} (t) =  {\bf f}_{s}(t) + {\bf f}_{f}(t) +
2463 + {\bf f}_{r}(t)
2464 + \end{equation}
2465 + Since the frictional and random forces are applied at the center of
2466 + resistance, which generally does not coincide with the center of mass,
2467 + extra torques are exerted at the center of mass. Thus, the net
2468 + body-fixed torque at the center of mass, $\tau^{~b}(t)$,
2469 + is given by
2470 + \begin{equation}
2471 + \tau^{~b} \leftarrow \tau_{s}^{~b} + \tau_{f}^{~b} + \tau_{r}^{~b} + {\bf r}_{MR} \times \left( {\bf f}_{f}^{~b} + {\bf f}_{r}^{~b} \right)
2472 + \end{equation}
2473 + where ${\bf r}_{MR}$ is the vector from the center of mass to the center of
2474 + resistance. Instead of integrating the angular velocity in lab-fixed
2475 + frame, we consider the equation of motion for the angular momentum
2476 + (${\bf j}$) in the body-fixed frame
2477 + \begin{equation}
2478 + \frac{\partial}{\partial t}{\bf j}(t) = \tau^{~b}(t)
2479 + \end{equation}
2480 + Embedding the friction and random forces into the the total force and
2481 + torque, one can integrate the Langevin equations of motion for a rigid
2482 + body of arbitrary shape in a velocity-Verlet style 2-part algorithm,
2483 + where $h = \delta t$:
2484 +
2485 + {\tt move A:}
2486 + \begin{align*}
2487 + {\bf v}\left(t + h / 2\right)  &\leftarrow  {\bf v}(t)
2488 +    + \frac{h}{2} \left( {\bf f}(t) / m \right), \\
2489 + %
2490 + {\bf r}(t + h) &\leftarrow {\bf r}(t)
2491 +    + h  {\bf v}\left(t + h / 2 \right), \\
2492 + %
2493 + {\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t)
2494 +    + \frac{h}{2} {\bf \tau}^{~b}(t), \\
2495 + %
2496 + {\bf Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
2497 +    (t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right).
2498 + \end{align*}
2499 + In this context, $\overleftrightarrow{\mathsf{I}}$ is the diagonal
2500 + moment of inertia tensor, and the $\mathrm{rotate}$ function is the
2501 + reversible product of the three body-fixed rotations,
2502 + \begin{equation}
2503 + \mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
2504 + \mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y
2505 + / 2) \cdot \mathsf{G}_x(a_x /2),
2506 + \end{equation}
2507 + where each rotational propagator, $\mathsf{G}_\alpha(\theta)$,
2508 + rotates both the rotation matrix ($\mathbf{Q}$) and the body-fixed
2509 + angular momentum (${\bf j}$) by an angle $\theta$ around body-fixed
2510 + axis $\alpha$,
2511 + \begin{equation}
2512 + \mathsf{G}_\alpha( \theta ) = \left\{
2513 + \begin{array}{lcl}
2514 + \mathbf{Q}(t) & \leftarrow & \mathbf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
2515 + {\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf
2516 + j}(0).
2517 + \end{array}
2518 + \right.
2519 + \end{equation}
2520 + $\mathsf{R}_\alpha$ is a quadratic approximation to the single-axis
2521 + rotation matrix.  For example, in the small-angle limit, the
2522 + rotation matrix around the body-fixed x-axis can be approximated as
2523 + \begin{equation}
2524 + \mathsf{R}_x(\theta) \approx \left(
2525 + \begin{array}{ccc}
2526 + 1 & 0 & 0 \\
2527 + 0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}  & -\frac{\theta}{1+
2528 + \theta^2 / 4} \\
2529 + 0 & \frac{\theta}{1+ \theta^2 / 4} & \frac{1-\theta^2 / 4}{1 +
2530 + \theta^2 / 4}
2531 + \end{array}
2532 + \right).
2533 + \end{equation}
2534 + All other rotations follow in a straightforward manner. After the
2535 + first part of the propagation, the forces and body-fixed torques are
2536 + calculated at the new positions and orientations.  The system forces
2537 + and torques are derivatives of the total potential energy function
2538 + ($U$) with respect to the rigid body positions (${\bf r}$) and the
2539 + columns of the transposed rotation matrix ${\bf Q}^T = \left({\bf
2540 + u}_x, {\bf u}_y, {\bf u}_z \right)$:
2541  
2542 + {\tt Forces:}
2543 + \begin{align*}
2544 + {\bf f}_{s}(t + h) & \leftarrow
2545 +    - \left(\frac{\partial U}{\partial {\bf r}}\right)_{{\bf r}(t + h)} \\
2546 + %
2547 + {\bf \tau}_{s}(t + h) &\leftarrow {\bf u}(t + h)
2548 +    \times \frac{\partial U}{\partial {\bf u}} \\
2549 + %
2550 + {\bf v}^{b}_{R}(t+h) & \leftarrow \mathbf{Q}(t+h) \cdot \left({\bf v}(t+h) + {\bf \omega}(t+h) \times {\bf r}_{MR} \right) \\
2551 + %
2552 + {\bf f}_{R,f}^{b}(t+h) & \leftarrow - \Xi_{R}^{tt} \cdot
2553 + {\bf v}^{b}_{R}(t+h) - \Xi_{R}^{rt} \cdot {\bf \omega}(t+h) \\
2554 + %
2555 + {\bf \tau}_{R,f}^{b}(t+h) & \leftarrow - \Xi_{R}^{tr} \cdot
2556 + {\bf v}^{b}_{R}(t+h) - \Xi_{R}^{rr} \cdot {\bf \omega}(t+h) \\
2557 + %
2558 + Z & \leftarrow  {\tt GaussianNormal}(2 k_B T / h, 6) \\
2559 + {\bf F}_{R,r}^{b}(t+h) & \leftarrow {\bf S} \cdot Z \\
2560 + %
2561 + {\bf f}(t+h) & \leftarrow {\bf f}_{s}(t+h) + \mathbf{Q}^{T}(t+h)
2562 + \cdot \left({\bf f}_{R,f}^{~b} + {\bf f}_{R,r}^{~b} \right) \\
2563 + %
2564 + \tau(t+h) & \leftarrow \tau_{s}(t+h) + \mathbf{Q}^{T}(t+h) \cdot \left(\tau_{R,f}^{~b} + \tau_{R,r}^{~b} \right) + {\bf r}_{MR} \times \left({\bf f}_{f}(t+h) + {\bf f}_{r}(t+h) \right) \\
2565 + \tau^{~b}(t+h) & \leftarrow \mathbf{Q}(t+h) \cdot \tau(t+h) \\
2566 + \end{align*}
2567 + Frictional (and random) forces and torques must be computed at the
2568 + center of resistance, so there are additional steps required to find
2569 + the body-fixed velocity (${\bf v}_{R}^{~b}$) at this location.  Mapping
2570 + the frictional and random forces at the center of resistance back to
2571 + the center of mass also introduces an additional term in the torque
2572 + one obtains at the center of mass.
2573  
2574 + Once the forces and torques have been obtained at the new time step,
2575 + the velocities can be advanced to the same time value.
2576 +
2577 + {\tt move B:}
2578 + \begin{align*}
2579 + {\bf v}\left(t + h \right)  &\leftarrow  {\bf v}\left(t + h / 2
2580 + \right)
2581 +    + \frac{h}{2} \left( {\bf f}(t + h) / m \right), \\
2582 + %
2583 + {\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t + h / 2
2584 + \right)
2585 +    + \frac{h}{2} {\bf \tau}^{~b}(t + h) .
2586 + \end{align*}
2587 +
2588 + The viscosity of the implicit of solvents must be specified using {\tt
2589 + viscosity} keywords in the meta-data file to use langevin dynamics
2590 + integrator. For simple shaped particles (spheres and ellipsoids), no
2591 + further parameters are necessary. However, there are no analytical
2592 + solutions for composite shaped particles, the approximate methods have
2593 + to be applied to get the resistance tensor. The file which contains
2594 + the information about hydro properties is indicated by {\tt
2595 + HydroPropFile} keyword in meta-data file. The {\tt HydroPropFile} is
2596 + precalculated by {\tt Hydro}.
2597 +
2598 + \begin{longtable}[c]{ABCD}
2599 + \caption{Meta-data Keywords: Langevin Dynamics Parameters}
2600 + \\
2601 + {\bf keyword} & {\bf units} & {\bf use} & {\bf remarks}  \\ \hline
2602 + \endhead
2603 + \hline
2604 + \endfoot
2605 + {\tt viscosity} & centipoise & Sets the value of viscosity of implicit
2606 + solvents & \\ {\tt HydroPropFile} & string & specifies the resistance
2607 + tensor file & usually a {\tt .diff} file which is precalculated by
2608 + {\sc Hydro}. Not necessory for simple shaped particles (spheres and
2609 + ellipsoids) \\
2610 + {\tt beadSize} & $\mbox{\AA}$ & Sets the diameters of
2611 + beads when {\sc Rough Shell Model} is used to generate the resistance
2612 + tensor file. \\
2613 + \label{table:ldParameters}
2614 + \end{longtable}
2615 +
2616 +
2617   \section{\label{sec:constraints}Constraint Methods}
2618  
2619   \subsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond
# Line 3214 | Line 3457 | The options available for SimpleBuilder are as follows
3457      &  {\tt -{}-ny=INT}           &  number of unit cells in y\\
3458      &  {\tt -{}-nz=INT}            &  number of unit cells in z
3459   \end{longtable}
3460 +
3461 + \section{\label{oopseSec:Hydro}Hydro}
3462 + {\tt Hydro} generates {\tt .diff} file which is required when a Langevin
3463 + Dynamics simulation using approximate models (supports Bead Model and
3464 + Rough Shell Model) is performed. To generate the {\tt }.diff file, the
3465 + meta-data file is needed as the input file. The viscosity of the fluid
3466 + flow (solvent) and the temperature of the system have to be defined in
3467 + meta-data file. If the approximate model is {\tt Rough Shell Model},
3468 + the {\tt beadSize} which is the diameter of every beads must be
3469 + specified in meta-data file.
3470  
3471 + The options available for Hydro are as follows:
3472 + \begin{longtable}[c]{|EFG|}
3473 + \caption{Hydro Command-line Options}
3474 + \\ \hline
3475 + {\bf option} & {\bf verbose option} & {\bf behavior} \\ \hline
3476 + \endhead
3477 + \hline
3478 + \endfoot
3479 +  -h& {\tt -{}-help}               & Print help and exit\\
3480 +  -V& {\tt -{}-version}            & Print version and exit\\
3481 +  -i& {\tt -{}-input=filename}     & input MetaData (md) file\\
3482 +  -o& {\tt -{}-output=STRING}      & Output file name\\
3483 +   &  {\tt -{}-model=STRING}     & hydrodynamics model (supports
3484 + RoughShell and BeadModel)\\
3485 +  -b&  {\tt -{}-beads}            & generate the beads only,
3486 + hydrodynamics will be performed (default=off)\\
3487 + \end{longtable}
3488 +
3489 +
3490   \chapter{\label{oopseSec:parallelization} Parallel Simulation Implementation}
3491  
3492   Although processor power is continually improving, it is still

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines