| 312 |
|
\begin{figure} |
| 313 |
|
\centering |
| 314 |
|
\includegraphics[width=\linewidth]{lipidModel.eps} |
| 315 |
< |
\caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ % |
| 315 |
> |
\caption[A representation of a lipid model in {\sc duff}]{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ % |
| 316 |
|
is the bend angle, $\mu$ is the dipole moment of the head group, and n |
| 317 |
|
is the chain length.} |
| 318 |
|
\label{oopseFig:lipidModel} |
| 631 |
|
\phi_{ij}({\bf r}_{ij}) \\ |
| 632 |
|
\rho_{i} & = & \sum_{j \neq i} f_{j}({\bf r}_{ij}) |
| 633 |
|
\end{eqnarray} |
| 634 |
< |
where $F_{i} $ is the embedding function that equates the energy required to embed a |
| 635 |
< |
positively-charged core ion $i$ into a linear superposition of |
| 636 |
< |
spherically averaged atomic electron densities given by |
| 637 |
< |
$\rho_{i}$. $\phi_{ij}$ is a primarily repulsive pairwise interaction |
| 638 |
< |
between atoms $i$ and $j$. In the original formulation of |
| 639 |
< |
{\sc eam}\cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term, however |
| 640 |
< |
in later refinements to EAM have shown that non-uniqueness between $F$ |
| 641 |
< |
and $\phi$ allow for more general forms for $\phi$.\cite{Daw89} |
| 642 |
< |
There is a cutoff distance, $r_{cut}$, which limits the |
| 643 |
< |
summations in the {\sc eam} equation to the few dozen atoms |
| 634 |
> |
where $F_{i} $ is the embedding function that equates the energy |
| 635 |
> |
required to embed a positively-charged core ion $i$ into a linear |
| 636 |
> |
superposition of spherically averaged atomic electron densities given |
| 637 |
> |
by $\rho_{i}$. $\phi_{ij}$ is a primarily repulsive pairwise |
| 638 |
> |
interaction between atoms $i$ and $j$. In the original formulation of |
| 639 |
> |
{\sc eam}\cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term, |
| 640 |
> |
however in later refinements to {\sc eam} have shown that non-uniqueness |
| 641 |
> |
between $F$ and $\phi$ allow for more general forms for |
| 642 |
> |
$\phi$.\cite{Daw89} There is a cutoff distance, $r_{cut}$, which |
| 643 |
> |
limits the summations in the {\sc eam} equation to the few dozen atoms |
| 644 |
|
surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$ |
| 645 |
< |
interactions. Foiles et al. fit EAM potentials for fcc metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals\cite{FBD86}. These potential fits are in the DYNAMO 86 format and are included with {\sc oopse}. |
| 645 |
> |
interactions. Foiles \emph{et al}.~fit {\sc eam} potentials for the fcc |
| 646 |
> |
metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals.\cite{FBD86} |
| 647 |
> |
These fits, are included in {\sc oopse}. |
| 648 |
|
|
| 647 |
– |
|
| 649 |
|
\subsection{\label{oopseSec:pbc}Periodic Boundary Conditions} |
| 650 |
|
|
| 651 |
|
\newcommand{\roundme}{\operatorname{round}} |
| 659 |
|
simulation box is large enough to avoid ``feeling'' the symmetries of |
| 660 |
|
the periodic lattice, surface effects can be ignored. The available |
| 661 |
|
periodic cells in OOPSE are cubic, orthorhombic and parallelepiped. We |
| 662 |
< |
use a $3 \times 3$ matrix, $\mathbf{H}$, to describe the shape and |
| 663 |
< |
size of the simulation box. $\mathbf{H}$ is defined: |
| 662 |
> |
use a $3 \times 3$ matrix, $\mathsf{H}$, to describe the shape and |
| 663 |
> |
size of the simulation box. $\mathsf{H}$ is defined: |
| 664 |
|
\begin{equation} |
| 665 |
< |
\mathbf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z ) |
| 665 |
> |
\mathsf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z ) |
| 666 |
|
\end{equation} |
| 667 |
|
Where $\mathbf{h}_j$ is the column vector of the $j$th axis of the |
| 668 |
|
box. During the course of the simulation both the size and shape of |
| 669 |
< |
the box can be changed to allow volume fluctations when constraining |
| 669 |
> |
the box can be changed to allow volume fluctuations when constraining |
| 670 |
|
the pressure. |
| 671 |
|
|
| 672 |
|
A real space vector, $\mathbf{r}$ can be transformed in to a box space |
| 673 |
|
vector, $\mathbf{s}$, and back through the following transformations: |
| 674 |
|
\begin{align} |
| 675 |
< |
\mathbf{s} &= \mathbf{H}^{-1} \mathbf{r} \\ |
| 676 |
< |
\mathbf{r} &= \mathbf{H} \mathbf{s} |
| 675 |
> |
\mathbf{s} &= \mathsf{H}^{-1} \mathbf{r} \\ |
| 676 |
> |
\mathbf{r} &= \mathsf{H} \mathbf{s} |
| 677 |
|
\end{align} |
| 678 |
|
The vector $\mathbf{s}$ is now a vector expressed as the number of box |
| 679 |
|
lengths in the $\mathbf{h}_x$, $\mathbf{h}_y$, and $\mathbf{h}_z$ |
| 701 |
|
Finally, we obtain the minimum image coordinates $\mathbf{r}^{\prime}$ by |
| 702 |
|
transforming back to real space, |
| 703 |
|
\begin{equation} |
| 704 |
< |
\mathbf{r}^{\prime}=\mathbf{H}^{-1}\mathbf{s}^{\prime}% |
| 704 |
> |
\mathbf{r}^{\prime}=\mathsf{H}^{-1}\mathbf{s}^{\prime}% |
| 705 |
|
\end{equation} |
| 706 |
|
In this way, particles are allowed to diffuse freely in $\mathbf{r}$, |
| 707 |
|
but their minimum images, $\mathbf{r}^{\prime}$ are used to compute |
| 708 |
< |
the interatomic forces. |
| 708 |
> |
the inter-atomic forces. |
| 709 |
|
|
| 710 |
|
|
| 711 |
|
\section{\label{oopseSec:IOfiles}Input and Output Files} |
| 741 |
|
initialConfig = "./argon.init"; |
| 742 |
|
|
| 743 |
|
forceField = "LJ"; |
| 744 |
< |
ensemble = "NVE"; // specify the simulation enesemble |
| 744 |
> |
ensemble = "NVE"; // specify the simulation ensemble |
| 745 |
|
dt = 1.0; // the time step for integration |
| 746 |
|
runTime = 1e3; // the total simulation run time |
| 747 |
|
sampleTime = 100; // trajectory file frequency |
| 782 |
|
|
| 783 |
|
#include "argon.mdl" |
| 784 |
|
|
| 784 |
– |
molecule{ |
| 785 |
– |
name = "Ar"; |
| 786 |
– |
nAtoms = 1; |
| 787 |
– |
atom[0]{ |
| 788 |
– |
type="Ar"; |
| 789 |
– |
position( 0.0, 0.0, 0.0 ); |
| 790 |
– |
} |
| 791 |
– |
} |
| 792 |
– |
|
| 785 |
|
nComponents = 1; |
| 786 |
|
component{ |
| 787 |
|
type = "Ar"; |
| 811 |
|
entities are written out using quanternions, to save space in the |
| 812 |
|
output files. |
| 813 |
|
|
| 814 |
< |
\begin{lstlisting}[float,caption={[The format of the coordinate files]Shows the format of the coordinate files. The fist line is the number of atoms. The second line begins with the time stamp followed by the three $\mathbf{H}$ column vectors. The next lines are the atomic coordinates for all atoms in the system. First is the name followed by position, velocity, quanternions, and lastly angular momentum.},label=sch:dumpFormat] |
| 814 |
> |
\begin{lstlisting}[float,caption={[The format of the coordinate files]Shows the format of the coordinate files. The fist line is the number of atoms. The second line begins with the time stamp followed by the three $\mathsf{H}$ column vectors. It is important to note, that for extended system ensembles, additional information pertinent to the integrators may be stored on this line as well.. The next lines are the atomic coordinates for all atoms in the system. First is the name followed by position, velocity, quanternions, and lastly angular velocities.},label=sch:dumpFormat] |
| 815 |
|
|
| 816 |
|
nAtoms |
| 817 |
|
time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz; |
| 846 |
|
|
| 847 |
|
As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization |
| 848 |
|
file is needed to provide the starting coordinates for a |
| 849 |
< |
simulation. The {\sc oopse} package provides a program called |
| 850 |
< |
\texttt{sysBuilder} to aid in the creation of the \texttt{.init} |
| 851 |
< |
file. \texttt{sysBuilder} uses {\sc bass}, and will recognize |
| 849 |
> |
simulation. The {\sc oopse} package provides several system building |
| 850 |
> |
programs to aid in the creation of the \texttt{.init} |
| 851 |
> |
file. The programs use {\sc bass}, and will recognize |
| 852 |
|
arguments and parameters in the \texttt{.bass} file that would |
| 853 |
|
otherwise be ignored by the simulation. |
| 854 |
|
|
| 976 |
|
formulation of the {\sc shake} method\cite{ryckaert77} of iteratively |
| 977 |
|
solving the Lagrange multipliers of constraint. The system of lagrange |
| 978 |
|
multipliers allows one to reformulate the equations of motion with |
| 979 |
< |
explicit constraint forces on the equations of |
| 988 |
< |
motion.\cite{fowles99:lagrange} |
| 979 |
> |
explicit constraint forces.\cite{fowles99:lagrange} |
| 980 |
|
|
| 981 |
< |
Consider a system described by qoordinates $q_1$ and $q_2$ subject to an |
| 981 |
> |
Consider a system described by coordinates $q_1$ and $q_2$ subject to an |
| 982 |
|
equation of constraint: |
| 983 |
|
\begin{equation} |
| 984 |
|
\sigma(q_1, q_2,t) = 0 |
| 1050 |
|
In a simulation, this would involve the solution of a set of $(m + n)$ |
| 1051 |
|
number of equations. Where $m$ is the number of constraints, and $n$ |
| 1052 |
|
is the number of constrained coordinates. In practice, this is not |
| 1053 |
< |
done, as the matrix inversion neccassary to solve the system of |
| 1053 |
> |
done, as the matrix inversion necessary to solve the system of |
| 1054 |
|
equations would be very time consuming to solve. Additionally, the |
| 1055 |
|
numerical error in the solution of the set of $\lambda$'s would be |
| 1056 |
|
compounded by the error inherent in propagating by the Velocity Verlet |
| 1057 |
< |
algorithm ($\Delta t^4$). The verlet propagation error is negligible |
| 1058 |
< |
in an unconstrained system, as one is interested in the statisitics of |
| 1057 |
> |
algorithm ($\Delta t^4$). The Verlet propagation error is negligible |
| 1058 |
> |
in an unconstrained system, as one is interested in the statistics of |
| 1059 |
|
the run, and not that the run be numerically exact to the ``true'' |
| 1060 |
|
integration. This relates back to the ergodic hypothesis that a time |
| 1061 |
< |
integral of a valid trajectory will still give the correct enesemble |
| 1061 |
> |
integral of a valid trajectory will still give the correct ensemble |
| 1062 |
|
average. However, in the case of constraints, if the equations of |
| 1063 |
|
motion leave the ``true'' trajectory, they are departing from the |
| 1064 |
|
constrained surface. The method that is used, is to iteratively solve |
| 1077 |
|
Eq.~\ref{oopseEq:c1} is the set of bond constraints, where $d_{ij}$ is |
| 1078 |
|
the constrained distance between atom $i$ and |
| 1079 |
|
$j$. Eq.~\ref{oopseEq:c2} constrains the velocities of $i$ and $j$ to |
| 1080 |
< |
be perpindicular to the bond vector, so that the bond can neither grow |
| 1080 |
> |
be perpendicular to the bond vector, so that the bond can neither grow |
| 1081 |
|
nor shrink. The constrained dynamics equations become: |
| 1082 |
|
\begin{equation} |
| 1083 |
|
m_i \mathbf{\ddot{r}}_i = \mathbf{F}_i + \mathbf{\mathcal{G}}_i |
| 1084 |
|
\label{oopseEq:r1} |
| 1085 |
|
\end{equation} |
| 1086 |
< |
Where, |
| 1086 |
> |
Where,$\mathbf{\mathcal{G}}_i$ are the forces of constraint on $i$, |
| 1087 |
> |
and are defined: |
| 1088 |
|
\begin{equation} |
| 1089 |
|
\mathbf{\mathcal{G}}_i = - \sum_j \lambda_{ij}(t)\,\nabla \sigma_{ij} |
| 1090 |
|
\label{oopseEq:r2} |
| 1103 |
|
\mathbf{F}_i(t+h) + \mathbf{\mathcal{G}}_{Vi}(t+h) \Bigr] % |
| 1104 |
|
\label{oopseEq:vv2} |
| 1105 |
|
\end{align} |
| 1106 |
+ |
Where: |
| 1107 |
+ |
\begin{align} |
| 1108 |
+ |
\mathbf{\mathcal{G}}_{Ri}(t) &= |
| 1109 |
+ |
-2 \sum_j \lambda_{Rij}(t) \mathbf{r}_{ij}(t) \\ |
| 1110 |
+ |
% |
| 1111 |
+ |
\mathbf{\mathcal{G}}_{Vi}(t+h) &= |
| 1112 |
+ |
-2 \sum_j \lambda_{Vij}(t+h) \mathbf{r}(t+h) |
| 1113 |
+ |
\end{align} |
| 1114 |
+ |
Next, define: |
| 1115 |
+ |
\begin{align} |
| 1116 |
+ |
g_{ij} &= h \lambda_{Rij}(t) \\ |
| 1117 |
+ |
k_{ij} &= h \lambda_{Vij}(t+h) \\ |
| 1118 |
+ |
\mathbf{q}_i &= \mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} \mathbf{F}_i(t) |
| 1119 |
+ |
- \frac{1}{m_i}\sum_j g_{ij}\mathbf{r}_{ij}(t) |
| 1120 |
+ |
\end{align} |
| 1121 |
+ |
Using these definitions, Eq.~\ref{oopseEq:vv1} and \ref{oopseEq:vv2} |
| 1122 |
+ |
can be rewritten as, |
| 1123 |
+ |
\begin{align} |
| 1124 |
+ |
\mathbf{r}_i(t+h) &= \mathbf{r}_i(t) + h \mathbf{q}_i \\ |
| 1125 |
+ |
% |
| 1126 |
+ |
\mathbf{\dot{r}}(t+h) &= \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h) |
| 1127 |
+ |
-\frac{1}{m_i}\sum_j k_{ij} \mathbf{r}_{ij}(t+h) |
| 1128 |
+ |
\end{align} |
| 1129 |
|
|
| 1130 |
+ |
To integrate the equations of motion, the {\sc rattle} algorithm first |
| 1131 |
+ |
solves for $\mathbf{r}(t+h)$. Let, |
| 1132 |
+ |
\begin{equation} |
| 1133 |
+ |
\mathbf{q}_i = \mathbf{\dot{r}}(t) + \frac{h}{2m_i}\mathbf{F}_i(t) |
| 1134 |
+ |
\end{equation} |
| 1135 |
+ |
Here $\mathbf{q}_i$ corresponds to an initial unconstrained move. Next |
| 1136 |
+ |
pick a constraint $j$, and let, |
| 1137 |
+ |
\begin{equation} |
| 1138 |
+ |
\mathbf{s} = \mathbf{r}_i(t) + h\mathbf{q}_i(t) |
| 1139 |
+ |
- \mathbf{r}_j(t) + h\mathbf{q}_j(t) |
| 1140 |
+ |
\label{oopseEq:ra1} |
| 1141 |
+ |
\end{equation} |
| 1142 |
+ |
If |
| 1143 |
+ |
\begin{equation} |
| 1144 |
+ |
\Big| |\mathbf{s}|^2 - d_{ij}^2 \Big| > \text{tolerance}, |
| 1145 |
+ |
\end{equation} |
| 1146 |
+ |
then the constraint is unsatisfied, and corrections are made to the |
| 1147 |
+ |
positions. First we define a test corrected configuration as, |
| 1148 |
+ |
\begin{align} |
| 1149 |
+ |
\mathbf{r}_i^T(t+h) = \mathbf{r}_i(t) + h\biggl[\mathbf{q}_i - |
| 1150 |
+ |
g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_i} \biggr] \\ |
| 1151 |
+ |
% |
| 1152 |
+ |
\mathbf{r}_j^T(t+h) = \mathbf{r}_j(t) + h\biggl[\mathbf{q}_j + |
| 1153 |
+ |
g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_j} \biggr] |
| 1154 |
+ |
\end{align} |
| 1155 |
+ |
And we chose $g_{ij}$ such that, $|\mathbf{r}_i^T - \mathbf{r}_j^T|^2 |
| 1156 |
+ |
= d_{ij}^2$. Solving the quadratic for $g_{ij}$ we obtain the |
| 1157 |
+ |
approximation, |
| 1158 |
+ |
\begin{equation} |
| 1159 |
+ |
g_{ij} = \frac{(s^2 - d^2)}{2h[\mathbf{s}\cdot\mathbf{r}_{ij}(t)] |
| 1160 |
+ |
(\frac{1}{m_i} + \frac{1}{m_j})} |
| 1161 |
+ |
\end{equation} |
| 1162 |
+ |
Although not an exact solution for $g_{ij}$, as this is an iterative |
| 1163 |
+ |
scheme overall, the eventual solution will converge. With a trial |
| 1164 |
+ |
$g_{ij}$, the new $\mathbf{q}$'s become, |
| 1165 |
+ |
\begin{align} |
| 1166 |
+ |
\mathbf{q}_i &= \mathbf{q}^{\text{old}}_i - g_{ij}\, |
| 1167 |
+ |
\frac{\mathbf{r}_{ij}(t)}{m_i} \\ |
| 1168 |
+ |
% |
| 1169 |
+ |
\mathbf{q}_j &= \mathbf{q}^{\text{old}}_j + g_{ij}\, |
| 1170 |
+ |
\frac{\mathbf{r}_{ij}(t)}{m_j} |
| 1171 |
+ |
\end{align} |
| 1172 |
+ |
The whole algorithm is then repeated from Eq.~\ref{oopseEq:ra1} until |
| 1173 |
+ |
all constraints are satisfied. |
| 1174 |
|
|
| 1175 |
+ |
The second step of {\sc rattle}, is to then update the velocities. The |
| 1176 |
+ |
step starts with, |
| 1177 |
+ |
\begin{equation} |
| 1178 |
+ |
\mathbf{\dot{r}}_i(t+h) = \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h) |
| 1179 |
+ |
\end{equation} |
| 1180 |
+ |
Next we pick a constraint $j$, and calculate the dot product $\ell$. |
| 1181 |
+ |
\begin{equation} |
| 1182 |
+ |
\ell = \mathbf{r}_{ij}(t+h) \cdot \mathbf{\dot{r}}_{ij}(t+h) |
| 1183 |
+ |
\label{oopseEq:rv1} |
| 1184 |
+ |
\end{equation} |
| 1185 |
+ |
Here if constraint Eq.~\ref{oopseEq:c2} holds, $\ell$ should be |
| 1186 |
+ |
zero. Therefore if $\ell$ is greater than some tolerance, then |
| 1187 |
+ |
corrections are made to the $i$ and $j$ velocities. |
| 1188 |
+ |
\begin{align} |
| 1189 |
+ |
\mathbf{\dot{r}}_i^T &= \mathbf{\dot{r}}_i(t+h) - k_{ij} |
| 1190 |
+ |
\frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_i} \\ |
| 1191 |
+ |
% |
| 1192 |
+ |
\mathbf{\dot{r}}_j^T &= \mathbf{\dot{r}}_j(t+h) + k_{ij} |
| 1193 |
+ |
\frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_j} |
| 1194 |
+ |
\end{align} |
| 1195 |
+ |
Like in the previous step, we select a value for $k_{ij}$ such that |
| 1196 |
+ |
$\ell$ is zero. |
| 1197 |
+ |
\begin{equation} |
| 1198 |
+ |
k_{ij} = \frac{\ell}{d^2_{ij}(\frac{1}{m_i} + \frac{1}{m_j})} |
| 1199 |
+ |
\end{equation} |
| 1200 |
+ |
The test velocities, $\mathbf{\dot{r}}^T_i$ and |
| 1201 |
+ |
$\mathbf{\dot{r}}^T_j$, then replace their respective velocities, and |
| 1202 |
+ |
the algorithm is iterated from Eq.~\ref{oopseEq:rv1} until all |
| 1203 |
+ |
constraints are satisfied. |
| 1204 |
|
|
| 1205 |
+ |
|
| 1206 |
|
\subsection{\label{oopseSec:zcons}Z-Constraint Method} |
| 1207 |
|
|
| 1208 |
< |
Based on fluctuation-dissipation theorem, a force auto-correlation |
| 1209 |
< |
method was developed to investigate the dynamics of ions inside the ion |
| 1210 |
< |
channels.\cite{Roux91} Time-dependent friction coefficient can be calculated |
| 1211 |
< |
from the deviation of the instantaneous force from its mean force. |
| 1212 |
< |
|
| 1124 |
< |
% |
| 1125 |
< |
|
| 1208 |
> |
Based on the fluctuation-dissipation theorem, a force auto-correlation |
| 1209 |
> |
method was developed by Roux and Karplus to investigate the dynamics |
| 1210 |
> |
of ions inside ion channels.\cite{Roux91} The time-dependent friction |
| 1211 |
> |
coefficient can be calculated from the deviation of the instantaneous |
| 1212 |
> |
force from its mean force. |
| 1213 |
|
\begin{equation} |
| 1214 |
|
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T |
| 1215 |
|
\end{equation} |
| 1219 |
|
\end{equation} |
| 1220 |
|
|
| 1221 |
|
|
| 1222 |
< |
If the time-dependent friction decay rapidly, static friction coefficient can |
| 1223 |
< |
be approximated by% |
| 1137 |
< |
|
| 1222 |
> |
If the time-dependent friction decays rapidly, the static friction |
| 1223 |
> |
coefficient can be approximated by |
| 1224 |
|
\begin{equation} |
| 1225 |
|
\xi^{static}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt |
| 1226 |
|
\end{equation} |
| 1227 |
< |
|
| 1142 |
< |
|
| 1143 |
< |
Hence, diffusion constant can be estimated by |
| 1227 |
> |
Therefore, the diffusion constant can then be estimated by |
| 1228 |
|
\begin{equation} |
| 1229 |
|
D(z)=\frac{k_{B}T}{\xi^{static}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty |
| 1230 |
|
}\langle\delta F(z,t)\delta F(z,0)\rangle dt}% |
| 1231 |
|
\end{equation} |
| 1232 |
|
|
| 1233 |
< |
|
| 1234 |
< |
\bigskip Z-Constraint method, which fixed the z coordinates of the molecules |
| 1235 |
< |
with respect to the center of the mass of the system, was proposed to obtain |
| 1236 |
< |
the forces required in force auto-correlation method.\cite{Marrink94} However, |
| 1237 |
< |
simply resetting the coordinate will move the center of the mass of the whole |
| 1238 |
< |
system. To avoid this problem, a new method was used at {\sc oopse}. Instead of |
| 1239 |
< |
resetting the coordinate, we reset the forces of z-constraint molecules as |
| 1240 |
< |
well as subtract the total constraint forces from the rest of the system after |
| 1241 |
< |
force calculation at each time step. |
| 1233 |
> |
The Z-Constraint method, which fixes the z coordinates of the |
| 1234 |
> |
molecules with respect to the center of the mass of the system, has |
| 1235 |
> |
been a method suggested to obtain the forces required for the force |
| 1236 |
> |
auto-correlation calculation.\cite{Marrink94} However, simply resetting the |
| 1237 |
> |
coordinate will move the center of the mass of the whole system. To |
| 1238 |
> |
avoid this problem, a new method was used in {\sc oopse}. Instead of |
| 1239 |
> |
resetting the coordinate, we reset the forces of z-constraint |
| 1240 |
> |
molecules as well as subtract the total constraint forces from the |
| 1241 |
> |
rest of the system after force calculation at each time step. |
| 1242 |
|
\begin{align} |
| 1243 |
|
F_{\alpha i}&=0\\ |
| 1244 |
|
V_{\alpha i}&=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}\\ |
| 1246 |
|
V_{\alpha i}&=V_{\alpha i}-\frac{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}} |
| 1247 |
|
\end{align} |
| 1248 |
|
|
| 1249 |
< |
At the very beginning of the simulation, the molecules may not be at its |
| 1250 |
< |
constraint position. To move the z-constraint molecule to the specified |
| 1251 |
< |
position, a simple harmonic potential is used% |
| 1168 |
< |
|
| 1249 |
> |
At the very beginning of the simulation, the molecules may not be at their |
| 1250 |
> |
constrained positions. To move a z-constrained molecule to its specified |
| 1251 |
> |
position, a simple harmonic potential is used |
| 1252 |
|
\begin{equation} |
| 1253 |
< |
U(t)=\frac{1}{2}k_{Harmonic}(z(t)-z_{cons})^{2}% |
| 1253 |
> |
U(t)=\frac{1}{2}k_{\text{Harmonic}}(z(t)-z_{\text{cons}})^{2}% |
| 1254 |
|
\end{equation} |
| 1255 |
< |
where $k_{Harmonic}$\bigskip\ is the harmonic force constant, $z(t)$ is |
| 1256 |
< |
current z coordinate of the center of mass of the z-constraint molecule, and |
| 1257 |
< |
$z_{cons}$ is the restraint position. Therefore, the harmonic force operated |
| 1258 |
< |
on the z-constraint molecule at time $t$ can be calculated by% |
| 1255 |
> |
where $k_{\text{Harmonic}}$ is the harmonic force constant, $z(t)$ is the |
| 1256 |
> |
current $z$ coordinate of the center of mass of the constrained molecule, and |
| 1257 |
> |
$z_{\text{cons}}$ is the constrained position. The harmonic force operating |
| 1258 |
> |
on the z-constrained molecule at time $t$ can be calculated by |
| 1259 |
|
\begin{equation} |
| 1260 |
< |
F_{z_{Harmonic}}(t)=-\frac{\partial U(t)}{\partial z(t)}=-k_{Harmonic}% |
| 1261 |
< |
(z(t)-z_{cons}) |
| 1260 |
> |
F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}= |
| 1261 |
> |
-k_{\text{Harmonic}}(z(t)-z_{\text{cons}}) |
| 1262 |
|
\end{equation} |
| 1180 |
– |
Worthy of mention, other kinds of potential functions can also be used to |
| 1181 |
– |
drive the z-constraint molecule. |
| 1263 |
|
|
| 1264 |
|
\section{\label{oopseSec:props}Trajectory Analysis} |
| 1265 |
|
|
| 1266 |
|
\subsection{\label{oopseSec:staticProps}Static Property Analysis} |
| 1267 |
|
|
| 1268 |
|
The static properties of the trajectories are analyzed with the |
| 1269 |
< |
program \texttt{staticProps}. The code is capable of calculating the following |
| 1270 |
< |
pair correlations between species A and B: |
| 1271 |
< |
\begin{itemize} |
| 1272 |
< |
\item $g_{\text{AB}}(r)$: Eq.~\ref{eq:gofr} |
| 1192 |
< |
\item $g_{\text{AB}}(r, \cos \theta)$: Eq.~\ref{eq:gofrCosTheta} |
| 1193 |
< |
\item $g_{\text{AB}}(r, \cos \omega)$: Eq.~\ref{eq:gofrCosOmega} |
| 1194 |
< |
\item $g_{\text{AB}}(x, y, z)$: Eq.~\ref{eq:gofrXYZ} |
| 1195 |
< |
\item $\langle \cos \omega \rangle_{\text{AB}}(r)$: |
| 1196 |
< |
Eq.~\ref{eq:cosOmegaOfR} |
| 1197 |
< |
\end{itemize} |
| 1269 |
> |
program \texttt{staticProps}. The code is capable of calculating a |
| 1270 |
> |
number of pair correlations between species A and B. Some of which |
| 1271 |
> |
only apply to directional entities. The summary of pair correlations |
| 1272 |
> |
can be found in Table~\ref{oopseTb:gofrs} |
| 1273 |
|
|
| 1274 |
+ |
\begin{table} |
| 1275 |
+ |
\caption[The list of pair correlations in \texttt{staticProps}]{The different pair correlations in \texttt{staticProps} along with whether atom A or B must be directional.} |
| 1276 |
+ |
\label{oopseTb:gofrs} |
| 1277 |
+ |
\begin{center} |
| 1278 |
+ |
\begin{tabular}{|l|c|c|} |
| 1279 |
+ |
\hline |
| 1280 |
+ |
Name & Equation & Directional Atom \\ \hline |
| 1281 |
+ |
$g_{\text{AB}}(r)$ & Eq.~\ref{eq:gofr} & neither \\ \hline |
| 1282 |
+ |
$g_{\text{AB}}(r, \cos \theta)$ & Eq.~\ref{eq:gofrCosTheta} & A \\ \hline |
| 1283 |
+ |
$g_{\text{AB}}(r, \cos \omega)$ & Eq.~\ref{eq:gofrCosOmega} & both \\ \hline |
| 1284 |
+ |
$g_{\text{AB}}(x, y, z)$ & Eq.~\ref{eq:gofrXYZ} & neither \\ \hline |
| 1285 |
+ |
$\langle \cos \omega \rangle_{\text{AB}}(r)$ & Eq.~\ref{eq:cosOmegaOfR} &% |
| 1286 |
+ |
both \\ \hline |
| 1287 |
+ |
\end{tabular} |
| 1288 |
+ |
\end{center} |
| 1289 |
+ |
\end{table} |
| 1290 |
+ |
|
| 1291 |
|
The first pair correlation, $g_{\text{AB}}(r)$, is defined as follows: |
| 1292 |
|
\begin{equation} |
| 1293 |
|
g_{\text{AB}}(r) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle %% |
| 1365 |
|
correlation that gives the average correlation of two directional |
| 1366 |
|
entities as a function of their distance from each other. |
| 1367 |
|
|
| 1276 |
– |
All static properties are calculated on a frame by frame basis. The |
| 1277 |
– |
trajectory is read a single frame at a time, and the appropriate |
| 1278 |
– |
calculations are done on each frame. Once one frame is finished, the |
| 1279 |
– |
next frame is read in, and a running average of the property being |
| 1280 |
– |
calculated is accumulated in each frame. The program allows for the |
| 1281 |
– |
user to specify more than one property be calculated in single run, |
| 1282 |
– |
preventing the need to read a file multiple times. |
| 1283 |
– |
|
| 1368 |
|
\subsection{\label{dynamicProps}Dynamic Property Analysis} |
| 1369 |
|
|
| 1370 |
|
The dynamic properties of a trajectory are calculated with the program |
| 1371 |
< |
\texttt{dynamicProps}. The program will calculate the following properties: |
| 1371 |
> |
\texttt{dynamicProps}. The program calculates the following properties: |
| 1372 |
|
\begin{gather} |
| 1373 |
|
\langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle \label{eq:rms}\\ |
| 1374 |
|
\langle \mathbf{v}(t) \cdot \mathbf{v}(0) \rangle \label{eq:velCorr} \\ |
| 1375 |
|
\langle \mathbf{j}(t) \cdot \mathbf{j}(0) \rangle \label{eq:angularVelCorr} |
| 1376 |
|
\end{gather} |
| 1377 |
|
|
| 1378 |
< |
Eq.~\ref{eq:rms} is the root mean square displacement |
| 1379 |
< |
function. Eq.~\ref{eq:velCorr} and Eq.~\ref{eq:angularVelCorr} are the |
| 1380 |
< |
velocity and angular velocity correlation functions respectively. The |
| 1381 |
< |
latter is only applicable to directional species in the simulation. |
| 1378 |
> |
Eq.~\ref{eq:rms} is the root mean square displacement function. Which |
| 1379 |
> |
allows one to observe the average displacement of an atom as a |
| 1380 |
> |
function of time. The quantity is useful when calculating diffusion |
| 1381 |
> |
coefficients because of the Einstein Relation, which is valid at long |
| 1382 |
> |
times.\cite{allen87:csl} |
| 1383 |
> |
\begin{equation} |
| 1384 |
> |
2tD = \langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle |
| 1385 |
> |
\label{oopseEq:einstein} |
| 1386 |
> |
\end{equation} |
| 1387 |
|
|
| 1388 |
< |
The \texttt{dynamicProps} program handles he file in a manner different from |
| 1389 |
< |
\texttt{staticProps}. As the properties calculated by this program are time |
| 1390 |
< |
dependent, multiple frames must be read in simultaneously by the |
| 1391 |
< |
program. For small trajectories this is no problem, and the entire |
| 1392 |
< |
trajectory is read into memory. However, for long trajectories of |
| 1304 |
< |
large systems, the files can be quite large. In order to accommodate |
| 1305 |
< |
large files, \texttt{dynamicProps} adopts a scheme whereby two blocks of memory |
| 1306 |
< |
are allocated to read in several frames each. |
| 1307 |
< |
|
| 1308 |
< |
In this two block scheme, the correlation functions are first |
| 1309 |
< |
calculated within each memory block, then the cross correlations |
| 1310 |
< |
between the frames contained within the two blocks are |
| 1311 |
< |
calculated. Once completed, the memory blocks are incremented, and the |
| 1312 |
< |
process is repeated. A diagram illustrating the process is shown in |
| 1313 |
< |
Fig.~\ref{oopseFig:dynamicPropsMemory}. As was the case with |
| 1314 |
< |
\texttt{staticProps}, multiple properties may be calculated in a |
| 1315 |
< |
single run to avoid multiple reads on the same file. |
| 1316 |
< |
|
| 1317 |
< |
|
| 1388 |
> |
Eq.~\ref{eq:velCorr} and \ref{eq:angularVelCorr} are the translational |
| 1389 |
> |
velocity and angular velocity correlation functions respectively. The |
| 1390 |
> |
latter is only applicable to directional species in the |
| 1391 |
> |
simulation. The velocity autocorrelation functions are useful when |
| 1392 |
> |
determining vibrational information about the system of interest. |
| 1393 |
|
|
| 1394 |
|
\section{\label{oopseSec:design}Program Design} |
| 1395 |
|
|
| 1429 |
|
developed to utilize the routines provided by \texttt{libBASS} and |
| 1430 |
|
\texttt{libmdtools}. The main program of the package is \texttt{oopse} |
| 1431 |
|
and the corresponding parallel version \texttt{oopse\_MPI}. These two |
| 1432 |
< |
programs will take the \texttt{.bass} file, and create then integrate |
| 1432 |
> |
programs will take the \texttt{.bass} file, and create (and integrate) |
| 1433 |
|
the simulation specified in the script. The two analysis programs |
| 1434 |
|
\texttt{staticProps} and \texttt{dynamicProps} utilize the core |
| 1435 |
|
libraries to initialize and read in trajectories from previously |
| 1441 |
|
|
| 1442 |
|
\subsection{\label{oopseSec:parallelization} Parallelization of {\sc oopse}} |
| 1443 |
|
|
| 1444 |
< |
Although processor power is continually growing month by month, it is |
| 1445 |
< |
still unreasonable to simulate systems of more then a 1000 atoms on a |
| 1446 |
< |
single processor. To facilitate study of larger system sizes or |
| 1447 |
< |
smaller systems on long time scales in a reasonable period of time, |
| 1448 |
< |
parallel methods were developed allowing multiple CPU's to share the |
| 1449 |
< |
simulation workload. Three general categories of parallel |
| 1450 |
< |
decomposition method's have been developed including atomic, spatial |
| 1451 |
< |
and force decomposition methods. |
| 1444 |
> |
Although processor power is continually growing roughly following |
| 1445 |
> |
Moore's Law, it is still unreasonable to simulate systems of more then |
| 1446 |
> |
a 1000 atoms on a single processor. To facilitate study of larger |
| 1447 |
> |
system sizes or smaller systems on long time scales in a reasonable |
| 1448 |
> |
period of time, parallel methods were developed allowing multiple |
| 1449 |
> |
CPU's to share the simulation workload. Three general categories of |
| 1450 |
> |
parallel decomposition methods have been developed including atomic, |
| 1451 |
> |
spatial and force decomposition methods. |
| 1452 |
|
|
| 1453 |
< |
Algorithmically simplest of the three method's is atomic decomposition |
| 1453 |
> |
Algorithmically simplest of the three methods is atomic decomposition |
| 1454 |
|
where N particles in a simulation are split among P processors for the |
| 1455 |
|
duration of the simulation. Computational cost scales as an optimal |
| 1456 |
|
$O(N/P)$ for atomic decomposition. Unfortunately all processors must |
| 1457 |
< |
communicate positions and forces with all other processors leading |
| 1458 |
< |
communication to scale as an unfavorable $O(N)$ independent of the |
| 1459 |
< |
number of processors. This communication bottleneck led to the |
| 1460 |
< |
development of spatial and force decomposition methods in which |
| 1461 |
< |
communication among processors scales much more favorably. Spatial or |
| 1462 |
< |
domain decomposition divides the physical spatial domain into 3D boxes |
| 1463 |
< |
in which each processor is responsible for calculation of forces and |
| 1464 |
< |
positions of particles located in its box. Particles are reassigned to |
| 1465 |
< |
different processors as they move through simulation space. To |
| 1466 |
< |
calculate forces on a given particle, a processor must know the |
| 1467 |
< |
positions of particles within some cutoff radius located on nearby |
| 1468 |
< |
processors instead of the positions of particles on all |
| 1469 |
< |
processors. Both communication between processors and computation |
| 1470 |
< |
scale as $O(N/P)$ in the spatial method. However, spatial |
| 1457 |
> |
communicate positions and forces with all other processors at every |
| 1458 |
> |
force evaluation, leading communication costs to scale as an |
| 1459 |
> |
unfavorable $O(N)$, \emph{independent of the number of processors}. This |
| 1460 |
> |
communication bottleneck led to the development of spatial and force |
| 1461 |
> |
decomposition methods in which communication among processors scales |
| 1462 |
> |
much more favorably. Spatial or domain decomposition divides the |
| 1463 |
> |
physical spatial domain into 3D boxes in which each processor is |
| 1464 |
> |
responsible for calculation of forces and positions of particles |
| 1465 |
> |
located in its box. Particles are reassigned to different processors |
| 1466 |
> |
as they move through simulation space. To calculate forces on a given |
| 1467 |
> |
particle, a processor must know the positions of particles within some |
| 1468 |
> |
cutoff radius located on nearby processors instead of the positions of |
| 1469 |
> |
particles on all processors. Both communication between processors and |
| 1470 |
> |
computation scale as $O(N/P)$ in the spatial method. However, spatial |
| 1471 |
|
decomposition adds algorithmic complexity to the simulation code and |
| 1472 |
|
is not very efficient for small N since the overall communication |
| 1473 |
|
scales as the surface to volume ratio $(N/P)^{2/3}$ in three |
| 1474 |
|
dimensions. |
| 1475 |
|
|
| 1476 |
< |
Force decomposition assigns particles to processors based on a block |
| 1477 |
< |
decomposition of the force matrix. Processors are split into a |
| 1478 |
< |
optimally square grid forming row and column processor groups. Forces |
| 1479 |
< |
are calculated on particles in a given row by particles located in |
| 1480 |
< |
that processors column assignment. Force decomposition is less complex |
| 1481 |
< |
to implement then the spatial method but still scales computationally |
| 1482 |
< |
as $O(N/P)$ and scales as $(N/\sqrt{p})$ in communication |
| 1483 |
< |
cost. Plimpton also found that force decompositions scales more |
| 1484 |
< |
favorably then spatial decomposition up to 10,000 atoms and favorably |
| 1485 |
< |
competes with spatial methods for up to 100,000 atoms. |
| 1476 |
> |
The parallelization method used in {\sc oopse} is the force |
| 1477 |
> |
decomposition method. Force decomposition assigns particles to |
| 1478 |
> |
processors based on a block decomposition of the force |
| 1479 |
> |
matrix. Processors are split into an optimally square grid forming row |
| 1480 |
> |
and column processor groups. Forces are calculated on particles in a |
| 1481 |
> |
given row by particles located in that processors column |
| 1482 |
> |
assignment. Force decomposition is less complex to implement than the |
| 1483 |
> |
spatial method but still scales computationally as $O(N/P)$ and scales |
| 1484 |
> |
as $O(N/\sqrt{P})$ in communication cost. Plimpton has also found that |
| 1485 |
> |
force decompositions scale more favorably than spatial decompositions |
| 1486 |
> |
for systems up to 10,000 atoms and favorably compete with spatial |
| 1487 |
> |
methods up to 100,000 atoms.\cite{plimpton95} |
| 1488 |
|
|
| 1489 |
|
\subsection{\label{oopseSec:memAlloc}Memory Issues in Trajectory Analysis} |
| 1490 |
|
|
| 1491 |
|
For large simulations, the trajectory files can sometimes reach sizes |
| 1492 |
|
in excess of several gigabytes. In order to effectively analyze that |
| 1493 |
< |
amount of data+, two memory management schemes have been devised for |
| 1493 |
> |
amount of data, two memory management schemes have been devised for |
| 1494 |
|
\texttt{staticProps} and for \texttt{dynamicProps}. The first scheme, |
| 1495 |
|
developed for \texttt{staticProps}, is the simplest. As each frame's |
| 1496 |
|
statistics are calculated independent of each other, memory is |
| 1498 |
|
complete for the snapshot. To prevent multiple passes through a |
| 1499 |
|
potentially large file, \texttt{staticProps} is capable of calculating |
| 1500 |
|
all requested correlations per frame with only a single pair loop in |
| 1501 |
< |
each frame and a single read through of the file. |
| 1501 |
> |
each frame and a single read of the file. |
| 1502 |
|
|
| 1503 |
|
The second, more advanced memory scheme, is used by |
| 1504 |
|
\texttt{dynamicProps}. Here, the program must have multiple frames in |
| 1508 |
|
in blocks. The number of frames in each block is specified by the |
| 1509 |
|
user, and upon reading a block of the trajectory, |
| 1510 |
|
\texttt{dynamicProps} will calculate all of the time correlation frame |
| 1511 |
< |
pairs within the block. After in block correlations are complete, a |
| 1511 |
> |
pairs within the block. After in-block correlations are complete, a |
| 1512 |
|
second block of the trajectory is read, and the cross correlations are |
| 1513 |
|
calculated between the two blocks. this second block is then freed and |
| 1514 |
|
then incremented and the process repeated until the end of the |
| 1526 |
|
\label{oopseFig:dynamicPropsMemory} |
| 1527 |
|
\end{figure} |
| 1528 |
|
|
| 1452 |
– |
\subsection{\label{openSource}Open Source and Distribution License} |
| 1453 |
– |
|
| 1529 |
|
\section{\label{oopseSec:conclusion}Conclusion} |
| 1530 |
|
|
| 1531 |
|
We have presented the design and implementation of our open source |
| 1532 |
< |
simulation package {\sc oopse}. The package offers novel |
| 1533 |
< |
capabilities to the field of Molecular Dynamics simulation packages in |
| 1534 |
< |
the form of dipolar force fields, and symplectic integration of rigid |
| 1535 |
< |
body dynamics. It is capable of scaling across multiple processors |
| 1536 |
< |
through the use of MPI. It also implements several integration |
| 1537 |
< |
ensembles allowing the end user control over temperature and |
| 1538 |
< |
pressure. In addition, it is capable of integrating constrained |
| 1539 |
< |
dynamics through both the {\sc rattle} algorithm and the z-constraint |
| 1540 |
< |
method. |
| 1532 |
> |
simulation package {\sc oopse}. The package offers novel capabilities |
| 1533 |
> |
to the field of Molecular Dynamics simulation packages in the form of |
| 1534 |
> |
dipolar force fields, and symplectic integration of rigid body |
| 1535 |
> |
dynamics. It is capable of scaling across multiple processors through |
| 1536 |
> |
the use of force based decomposition using MPI. It also implements |
| 1537 |
> |
several advanced integrators allowing the end user control over |
| 1538 |
> |
temperature and pressure. In addition, it is capable of integrating |
| 1539 |
> |
constrained dynamics through both the {\sc rattle} algorithm and the |
| 1540 |
> |
z-constraint method. |
| 1541 |
|
|
| 1542 |
|
These features are all brought together in a single open-source |
| 1543 |
< |
development package. This allows researchers to not only benefit from |
| 1543 |
> |
program. Allowing researchers to not only benefit from |
| 1544 |
|
{\sc oopse}, but also contribute to {\sc oopse}'s development as |
| 1545 |
|
well.Documentation and source code for {\sc oopse} can be downloaded |
| 1546 |
|
from \texttt{http://www.openscience.org/oopse/}. |
