| 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 |
– |
|
| 2374 |
– |
\section{\label{sec:constraints}Constraint Methods} |
| 2372 |
|
|
| 2373 |
< |
\subsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond |
| 2377 |
< |
Constraints} |
| 2378 |
< |
|
| 2379 |
< |
In order to satisfy the constraints of fixed bond lengths within {\sc |
| 2380 |
< |
oopse}, we have implemented the {\sc rattle} algorithm of |
| 2381 |
< |
Andersen.\cite{andersen83} {\sc rattle} is a velocity-Verlet |
| 2382 |
< |
formulation of the {\sc shake} method\cite{ryckaert77} for iteratively |
| 2383 |
< |
solving the Lagrange multipliers which maintain the holonomic |
| 2384 |
< |
constraints. Both methods are covered in depth in the |
| 2385 |
< |
literature,\cite{leach01:mm,Allen87} and a detailed description of |
| 2386 |
< |
this method would be redundant. |
| 2387 |
< |
|
| 2388 |
< |
\subsection{\label{oopseSec:zcons}The Z-Constraint Method} |
| 2389 |
< |
|
| 2390 |
< |
A force auto-correlation method based on the fluctuation-dissipation |
| 2391 |
< |
theorem was developed by Roux and Karplus to investigate the dynamics |
| 2392 |
< |
of ions inside ion channels.\cite{Roux91} The time-dependent friction |
| 2393 |
< |
coefficient can be calculated from the deviation of the instantaneous |
| 2394 |
< |
force from its mean value: |
| 2395 |
< |
\begin{equation} |
| 2396 |
< |
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T, |
| 2397 |
< |
\end{equation} |
| 2398 |
< |
where% |
| 2399 |
< |
\begin{equation} |
| 2400 |
< |
\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle. |
| 2401 |
< |
\end{equation} |
| 2402 |
< |
|
| 2403 |
< |
If the time-dependent friction decays rapidly, the static friction |
| 2404 |
< |
coefficient can be approximated by |
| 2405 |
< |
\begin{equation} |
| 2406 |
< |
\xi_{\text{static}}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt. |
| 2407 |
< |
\end{equation} |
| 2408 |
< |
|
| 2409 |
< |
This allows the diffusion constant to then be calculated through the |
| 2410 |
< |
Einstein relation:\cite{Marrink94} |
| 2411 |
< |
\begin{equation} |
| 2412 |
< |
D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty |
| 2413 |
< |
}\langle\delta F(z,t)\delta F(z,0)\rangle dt}.% |
| 2414 |
< |
\end{equation} |
| 2373 |
> |
\section{Langevin Dynamics (LANGEVINDYNAMICS)\label{LDRB}} |
| 2374 |
|
|
| 2375 |
< |
The Z-Constraint method, which fixes the $z$ coordinates of a few |
| 2376 |
< |
``tagged'' molecules with respect to the center of the mass of the |
| 2377 |
< |
system is a technique that was proposed to obtain the forces required |
| 2378 |
< |
for the force auto-correlation calculation.\cite{Marrink94} However, |
| 2379 |
< |
simply resetting the coordinate will move the center of the mass of |
| 2380 |
< |
the whole system. To avoid this problem, we have developed a new |
| 2381 |
< |
method that is utilized in {\sc oopse}. Instead of resetting the |
| 2423 |
< |
coordinates, we reset the forces of $z$-constrained molecules and |
| 2424 |
< |
subtract the total constraint forces from the rest of the system after |
| 2425 |
< |
the force calculation at each time step. |
| 2426 |
< |
|
| 2427 |
< |
After the force calculation, the total force on molecule $\alpha$ is: |
| 2428 |
< |
\begin{equation} |
| 2429 |
< |
G_{\alpha} = \sum_i F_{\alpha i}, |
| 2430 |
< |
\label{oopseEq:zc1} |
| 2431 |
< |
\end{equation} |
| 2432 |
< |
where $F_{\alpha i}$ is the force in the $z$ direction on atom $i$ in |
| 2433 |
< |
$z$-constrained molecule $\alpha$. The forces on the atoms in the |
| 2434 |
< |
$z$-constrained molecule are then adjusted to remove the total force |
| 2435 |
< |
on molecule $\alpha$: |
| 2436 |
< |
\begin{equation} |
| 2437 |
< |
F_{\alpha i} = F_{\alpha i} - |
| 2438 |
< |
\frac{m_{\alpha i} G_{\alpha}}{\sum_i m_{\alpha i}}. |
| 2439 |
< |
\end{equation} |
| 2440 |
< |
Here, $m_{\alpha i}$ is the mass of atom $i$ in the $z$-constrained |
| 2441 |
< |
molecule. After the forces have been adjusted, the velocities must |
| 2442 |
< |
also be modified to subtract out molecule $\alpha$'s center-of-mass |
| 2443 |
< |
velocity in the $z$ direction. |
| 2444 |
< |
\begin{equation} |
| 2445 |
< |
v_{\alpha i} = v_{\alpha i} - |
| 2446 |
< |
\frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}, |
| 2447 |
< |
\end{equation} |
| 2448 |
< |
where $v_{\alpha i}$ is the velocity of atom $i$ in the $z$ direction. |
| 2449 |
< |
Lastly, all of the accumulated constraint forces must be subtracted |
| 2450 |
< |
from the rest of the unconstrained system to keep the system center of |
| 2451 |
< |
mass of the entire system from drifting. |
| 2452 |
< |
\begin{equation} |
| 2453 |
< |
F_{\beta i} = F_{\beta i} - \frac{m_{\beta i} \sum_{\alpha} G_{\alpha}} |
| 2454 |
< |
{\sum_{\beta}\sum_i m_{\beta i}}, |
| 2455 |
< |
\end{equation} |
| 2456 |
< |
where $\beta$ denotes all {\it unconstrained} molecules in the |
| 2457 |
< |
system. Similarly, the velocities of the unconstrained molecules must |
| 2458 |
< |
also be scaled: |
| 2459 |
< |
\begin{equation} |
| 2460 |
< |
v_{\beta i} = v_{\beta i} + \sum_{\alpha} \frac{\sum_i m_{\alpha i} |
| 2461 |
< |
v_{\alpha i}}{\sum_i m_{\alpha i}}. |
| 2462 |
< |
\end{equation} |
| 2463 |
< |
|
| 2464 |
< |
This method will pin down the centers-of-mass of all of the |
| 2465 |
< |
$z$-constrained molecules, and will also keep the entire system fixed |
| 2466 |
< |
at the original system center-of-mass location. |
| 2467 |
< |
|
| 2468 |
< |
At the very beginning of the simulation, the molecules may not be at |
| 2469 |
< |
their desired positions. To steer a $z$-constrained molecule to its |
| 2470 |
< |
specified position, a simple harmonic potential is used: |
| 2471 |
< |
\begin{equation} |
| 2472 |
< |
U(t)=\frac{1}{2}k_{\text{Harmonic}}(z(t)-z_{\text{cons}})^{2},% |
| 2473 |
< |
\end{equation} |
| 2474 |
< |
where $k_{\text{Harmonic}}$ is an harmonic force constant, $z(t)$ is |
| 2475 |
< |
the current $z$ coordinate of the center of mass of the constrained |
| 2476 |
< |
molecule, and $z_{\text{cons}}$ is the desired constraint |
| 2477 |
< |
position. The harmonic force operating on the $z$-constrained molecule |
| 2478 |
< |
at time $t$ can be calculated by |
| 2479 |
< |
\begin{equation} |
| 2480 |
< |
F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}= |
| 2481 |
< |
-k_{\text{Harmonic}}(z(t)-z_{\text{cons}}). |
| 2482 |
< |
\end{equation} |
| 2483 |
< |
|
| 2484 |
< |
The user may also specify the use of a constant velocity method |
| 2485 |
< |
(steered molecular dynamics) to move the molecules to their desired |
| 2486 |
< |
initial positions. Based on concepts from atomic force microscopy, |
| 2487 |
< |
{\sc smd} has been used to study many processes which occur via rare |
| 2488 |
< |
events on the time scale of a few hundreds of picoseconds. For |
| 2489 |
< |
example,{\sc smd} has been used to observe the dissociation of |
| 2490 |
< |
Streptavidin-biotin Complex.\cite{smd} |
| 2491 |
< |
|
| 2492 |
< |
To use of the $z$-constraint method in an {\sc oopse} simulation, the |
| 2493 |
< |
molecules must be specified using the {\tt nZconstraints} keyword in |
| 2494 |
< |
the meta-data file. The other parameters for modifying the behavior |
| 2495 |
< |
of the $z$-constraint method are listed in table~\ref{table:zconParams}. |
| 2496 |
< |
|
| 2497 |
< |
\begin{longtable}[c]{ABCD} |
| 2498 |
< |
\caption{Meta-data Keywords: Z-Constraint Parameters} |
| 2499 |
< |
\\ |
| 2500 |
< |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 2501 |
< |
\endhead |
| 2502 |
< |
\hline |
| 2503 |
< |
\endfoot |
| 2504 |
< |
{\tt zconsTime} & fs & Sets the frequency at which the {\tt .fz} file |
| 2505 |
< |
is written & \\ |
| 2506 |
< |
{\tt zconsForcePolicy} & string & The strategy for subtracting |
| 2507 |
< |
the $z$-constraint force from the {\it unconstrained} molecules & Possible |
| 2508 |
< |
strategies are {\tt BYMASS} and {\tt BYNUMBER}. The default |
| 2509 |
< |
strategy is {\tt BYMASS}\\ |
| 2510 |
< |
{\tt zconsGap} & $\mbox{\AA}$ & Sets the distance between two adjacent |
| 2511 |
< |
constraint positions&Used mainly to move molecules through a |
| 2512 |
< |
simulation to estimate potentials of mean force. \\ |
| 2513 |
< |
{\tt zconsFixtime} & fs & Sets the length of time the $z$-constraint |
| 2514 |
< |
molecule is held fixed & {\tt zconsFixtime} must be set if {\tt |
| 2515 |
< |
zconsGap} is set\\ |
| 2516 |
< |
{\tt zconsUsingSMD} & logical & Flag for using Steered Molecular |
| 2517 |
< |
Dynamics to move the molecules to the correct constrained positions & |
| 2518 |
< |
Harmonic Forces are used by default |
| 2519 |
< |
\label{table:zconParams} |
| 2520 |
< |
\end{longtable} |
| 2521 |
< |
|
| 2522 |
< |
\section{Langevin Dynamics for Rigid Particles of Arbitrary Shape (LANGEVINDYNAMICS)\label{LDRB}} |
| 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 |
< |
{\sc oopse} implements langevin dynamics integrator to perform the |
| 2525 |
< |
molecular dynamics simulations at implicit solvents environment to |
| 2526 |
< |
speed up the simulation when the properties of solvents are not |
| 2527 |
< |
important. Consider the Langevin equations of motion in generalized |
| 2528 |
< |
coordinates |
| 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) |
| 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 |
| 2620 |
+ |
Constraints} |
| 2621 |
+ |
|
| 2622 |
+ |
In order to satisfy the constraints of fixed bond lengths within {\sc |
| 2623 |
+ |
oopse}, we have implemented the {\sc rattle} algorithm of |
| 2624 |
+ |
Andersen.\cite{andersen83} {\sc rattle} is a velocity-Verlet |
| 2625 |
+ |
formulation of the {\sc shake} method\cite{ryckaert77} for iteratively |
| 2626 |
+ |
solving the Lagrange multipliers which maintain the holonomic |
| 2627 |
+ |
constraints. Both methods are covered in depth in the |
| 2628 |
+ |
literature,\cite{leach01:mm,Allen87} and a detailed description of |
| 2629 |
+ |
this method would be redundant. |
| 2630 |
+ |
|
| 2631 |
+ |
\subsection{\label{oopseSec:zcons}The Z-Constraint Method} |
| 2632 |
+ |
|
| 2633 |
+ |
A force auto-correlation method based on the fluctuation-dissipation |
| 2634 |
+ |
theorem was developed by Roux and Karplus to investigate the dynamics |
| 2635 |
+ |
of ions inside ion channels.\cite{Roux91} The time-dependent friction |
| 2636 |
+ |
coefficient can be calculated from the deviation of the instantaneous |
| 2637 |
+ |
force from its mean value: |
| 2638 |
+ |
\begin{equation} |
| 2639 |
+ |
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T, |
| 2640 |
+ |
\end{equation} |
| 2641 |
+ |
where% |
| 2642 |
+ |
\begin{equation} |
| 2643 |
+ |
\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle. |
| 2644 |
+ |
\end{equation} |
| 2645 |
+ |
|
| 2646 |
+ |
If the time-dependent friction decays rapidly, the static friction |
| 2647 |
+ |
coefficient can be approximated by |
| 2648 |
+ |
\begin{equation} |
| 2649 |
+ |
\xi_{\text{static}}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt. |
| 2650 |
+ |
\end{equation} |
| 2651 |
+ |
|
| 2652 |
+ |
This allows the diffusion constant to then be calculated through the |
| 2653 |
+ |
Einstein relation:\cite{Marrink94} |
| 2654 |
+ |
\begin{equation} |
| 2655 |
+ |
D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty |
| 2656 |
+ |
}\langle\delta F(z,t)\delta F(z,0)\rangle dt}.% |
| 2657 |
+ |
\end{equation} |
| 2658 |
+ |
|
| 2659 |
+ |
The Z-Constraint method, which fixes the $z$ coordinates of a few |
| 2660 |
+ |
``tagged'' molecules with respect to the center of the mass of the |
| 2661 |
+ |
system is a technique that was proposed to obtain the forces required |
| 2662 |
+ |
for the force auto-correlation calculation.\cite{Marrink94} However, |
| 2663 |
+ |
simply resetting the coordinate will move the center of the mass of |
| 2664 |
+ |
the whole system. To avoid this problem, we have developed a new |
| 2665 |
+ |
method that is utilized in {\sc oopse}. Instead of resetting the |
| 2666 |
+ |
coordinates, we reset the forces of $z$-constrained molecules and |
| 2667 |
+ |
subtract the total constraint forces from the rest of the system after |
| 2668 |
+ |
the force calculation at each time step. |
| 2669 |
+ |
|
| 2670 |
+ |
After the force calculation, the total force on molecule $\alpha$ is: |
| 2671 |
+ |
\begin{equation} |
| 2672 |
+ |
G_{\alpha} = \sum_i F_{\alpha i}, |
| 2673 |
+ |
\label{oopseEq:zc1} |
| 2674 |
+ |
\end{equation} |
| 2675 |
+ |
where $F_{\alpha i}$ is the force in the $z$ direction on atom $i$ in |
| 2676 |
+ |
$z$-constrained molecule $\alpha$. The forces on the atoms in the |
| 2677 |
+ |
$z$-constrained molecule are then adjusted to remove the total force |
| 2678 |
+ |
on molecule $\alpha$: |
| 2679 |
+ |
\begin{equation} |
| 2680 |
+ |
F_{\alpha i} = F_{\alpha i} - |
| 2681 |
+ |
\frac{m_{\alpha i} G_{\alpha}}{\sum_i m_{\alpha i}}. |
| 2682 |
+ |
\end{equation} |
| 2683 |
+ |
Here, $m_{\alpha i}$ is the mass of atom $i$ in the $z$-constrained |
| 2684 |
+ |
molecule. After the forces have been adjusted, the velocities must |
| 2685 |
+ |
also be modified to subtract out molecule $\alpha$'s center-of-mass |
| 2686 |
+ |
velocity in the $z$ direction. |
| 2687 |
+ |
\begin{equation} |
| 2688 |
+ |
v_{\alpha i} = v_{\alpha i} - |
| 2689 |
+ |
\frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}, |
| 2690 |
+ |
\end{equation} |
| 2691 |
+ |
where $v_{\alpha i}$ is the velocity of atom $i$ in the $z$ direction. |
| 2692 |
+ |
Lastly, all of the accumulated constraint forces must be subtracted |
| 2693 |
+ |
from the rest of the unconstrained system to keep the system center of |
| 2694 |
+ |
mass of the entire system from drifting. |
| 2695 |
+ |
\begin{equation} |
| 2696 |
+ |
F_{\beta i} = F_{\beta i} - \frac{m_{\beta i} \sum_{\alpha} G_{\alpha}} |
| 2697 |
+ |
{\sum_{\beta}\sum_i m_{\beta i}}, |
| 2698 |
+ |
\end{equation} |
| 2699 |
+ |
where $\beta$ denotes all {\it unconstrained} molecules in the |
| 2700 |
+ |
system. Similarly, the velocities of the unconstrained molecules must |
| 2701 |
+ |
also be scaled: |
| 2702 |
+ |
\begin{equation} |
| 2703 |
+ |
v_{\beta i} = v_{\beta i} + \sum_{\alpha} \frac{\sum_i m_{\alpha i} |
| 2704 |
+ |
v_{\alpha i}}{\sum_i m_{\alpha i}}. |
| 2705 |
+ |
\end{equation} |
| 2706 |
+ |
|
| 2707 |
+ |
This method will pin down the centers-of-mass of all of the |
| 2708 |
+ |
$z$-constrained molecules, and will also keep the entire system fixed |
| 2709 |
+ |
at the original system center-of-mass location. |
| 2710 |
+ |
|
| 2711 |
+ |
At the very beginning of the simulation, the molecules may not be at |
| 2712 |
+ |
their desired positions. To steer a $z$-constrained molecule to its |
| 2713 |
+ |
specified position, a simple harmonic potential is used: |
| 2714 |
+ |
\begin{equation} |
| 2715 |
+ |
U(t)=\frac{1}{2}k_{\text{Harmonic}}(z(t)-z_{\text{cons}})^{2},% |
| 2716 |
+ |
\end{equation} |
| 2717 |
+ |
where $k_{\text{Harmonic}}$ is an harmonic force constant, $z(t)$ is |
| 2718 |
+ |
the current $z$ coordinate of the center of mass of the constrained |
| 2719 |
+ |
molecule, and $z_{\text{cons}}$ is the desired constraint |
| 2720 |
+ |
position. The harmonic force operating on the $z$-constrained molecule |
| 2721 |
+ |
at time $t$ can be calculated by |
| 2722 |
+ |
\begin{equation} |
| 2723 |
+ |
F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}= |
| 2724 |
+ |
-k_{\text{Harmonic}}(z(t)-z_{\text{cons}}). |
| 2725 |
+ |
\end{equation} |
| 2726 |
+ |
|
| 2727 |
+ |
The user may also specify the use of a constant velocity method |
| 2728 |
+ |
(steered molecular dynamics) to move the molecules to their desired |
| 2729 |
+ |
initial positions. Based on concepts from atomic force microscopy, |
| 2730 |
+ |
{\sc smd} has been used to study many processes which occur via rare |
| 2731 |
+ |
events on the time scale of a few hundreds of picoseconds. For |
| 2732 |
+ |
example,{\sc smd} has been used to observe the dissociation of |
| 2733 |
+ |
Streptavidin-biotin Complex.\cite{smd} |
| 2734 |
+ |
|
| 2735 |
+ |
To use of the $z$-constraint method in an {\sc oopse} simulation, the |
| 2736 |
+ |
molecules must be specified using the {\tt nZconstraints} keyword in |
| 2737 |
+ |
the meta-data file. The other parameters for modifying the behavior |
| 2738 |
+ |
of the $z$-constraint method are listed in table~\ref{table:zconParams}. |
| 2739 |
+ |
|
| 2740 |
+ |
\begin{longtable}[c]{ABCD} |
| 2741 |
+ |
\caption{Meta-data Keywords: Z-Constraint Parameters} |
| 2742 |
+ |
\\ |
| 2743 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 2744 |
+ |
\endhead |
| 2745 |
+ |
\hline |
| 2746 |
+ |
\endfoot |
| 2747 |
+ |
{\tt zconsTime} & fs & Sets the frequency at which the {\tt .fz} file |
| 2748 |
+ |
is written & \\ |
| 2749 |
+ |
{\tt zconsForcePolicy} & string & The strategy for subtracting |
| 2750 |
+ |
the $z$-constraint force from the {\it unconstrained} molecules & Possible |
| 2751 |
+ |
strategies are {\tt BYMASS} and {\tt BYNUMBER}. The default |
| 2752 |
+ |
strategy is {\tt BYMASS}\\ |
| 2753 |
+ |
{\tt zconsGap} & $\mbox{\AA}$ & Sets the distance between two adjacent |
| 2754 |
+ |
constraint positions&Used mainly to move molecules through a |
| 2755 |
+ |
simulation to estimate potentials of mean force. \\ |
| 2756 |
+ |
{\tt zconsFixtime} & fs & Sets the length of time the $z$-constraint |
| 2757 |
+ |
molecule is held fixed & {\tt zconsFixtime} must be set if {\tt |
| 2758 |
+ |
zconsGap} is set\\ |
| 2759 |
+ |
{\tt zconsUsingSMD} & logical & Flag for using Steered Molecular |
| 2760 |
+ |
Dynamics to move the molecules to the correct constrained positions & |
| 2761 |
+ |
Harmonic Forces are used by default |
| 2762 |
+ |
\label{table:zconParams} |
| 2763 |
+ |
\end{longtable} |
| 2764 |
+ |
|
| 2765 |
|
\chapter{\label{oopseSec:thermInt}Thermodynamic Integration} |
| 2766 |
|
|
| 2767 |
|
Thermodynamic integration is an established technique that has been |
