| 2 |
|
\usepackage{amsmath} |
| 3 |
|
\usepackage{amssymb} |
| 4 |
|
\usepackage{endfloat} |
| 5 |
– |
%\usepackage{berkeley} |
| 5 |
|
\usepackage{listings} |
| 6 |
+ |
\usepackage{palatino} |
| 7 |
|
\usepackage{graphicx} |
| 8 |
|
\usepackage[ref]{overcite} |
| 9 |
|
\usepackage{setspace} |
| 25 |
|
\title{{\sc oopse}: An Open Source Object-Oriented Parallel Simulation |
| 26 |
|
Engine for Molecular Dynamics} |
| 27 |
|
|
| 28 |
< |
\author{Matthew A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher J. Fennell and J. Daniel Gezelter\\ |
| 28 |
> |
\author{Matthew A. Meineke, Charles F. Vardeman II, Teng Lin,\\ |
| 29 |
> |
Christopher J. Fennell and J. Daniel Gezelter\\ |
| 30 |
|
Department of Chemistry and Biochemistry\\ |
| 31 |
|
University of Notre Dame\\ |
| 32 |
|
Notre Dame, Indiana 46556} |
| 36 |
|
|
| 37 |
|
\begin{abstract} |
| 38 |
|
We detail the capabilities of a new open-source parallel simulation |
| 39 |
< |
package ({\sc oopse}) that can perform molecular dynamics simulations |
| 40 |
< |
on atom types that are missing from other popular packages. In |
| 41 |
< |
particular, {\sc oopse} is capable of performing orientational |
| 41 |
< |
dynamics on dipolar systems, and it can handle simulations of metallic |
| 39 |
> |
progrm for MD ({\sc oopse}) that can work with atom types that are missing from other popular packages. In |
| 40 |
> |
particular, {\sc oopse} is capable of performing efficient orientational |
| 41 |
> |
dynamics on dipolar or rigid body systems, and it can handle simulations of metallic |
| 42 |
|
systems using the embedded atom method ({\sc eam}). |
| 43 |
|
\end{abstract} |
| 44 |
|
|
| 77 |
|
those of the community of scientists who contribute back to the |
| 78 |
|
project. |
| 79 |
|
|
| 80 |
< |
We have structured this chapter to first discuss the empirical energy |
| 80 |
> |
We have structured this paper to first discuss the empirical energy |
| 81 |
|
functions that {\sc oopse } implements in |
| 82 |
|
Sec.~\ref{oopseSec:empiricalEnergy}. Following that is a discussion of |
| 83 |
|
the various input and output files associated with the package |
| 84 |
|
(Sec.~\ref{oopseSec:IOfiles}). Sec.~\ref{oopseSec:mechanics} |
| 85 |
|
elucidates the various Molecular Dynamics algorithms {\sc oopse} |
| 86 |
|
implements in the integration of the Newtonian equations of |
| 87 |
< |
motion. Basic analysis of the trajectories obtained from the |
| 88 |
< |
simulation is discussed in Sec.~\ref{oopseSec:props}. Program design |
| 87 |
> |
motion. Program design |
| 88 |
|
considerations are presented in Sec.~\ref{oopseSec:design}. And |
| 89 |
|
lastly, Sec.~\ref{oopseSec:conclusion} concludes the chapter. |
| 90 |
|
|
| 91 |
< |
\section{\label{oopseSec:empiricalEnergy}The Empirical Energy Functions} |
| 91 |
> |
\section{\label{oopseSec:IOfiles}Concepts \& Files} |
| 92 |
|
|
| 93 |
< |
\subsection{\label{oopseSec:atomsMolecules}Atoms, Molecules and Rigid Bodies} |
| 93 |
> |
\subsection{{\sc bass} and Model Files} |
| 94 |
|
|
| 95 |
< |
The basic unit of an {\sc oopse} simulation is the atom. The |
| 96 |
< |
parameters describing the atom are generalized to make the atom as |
| 97 |
< |
flexible a representation as possible. They may represent specific |
| 98 |
< |
atoms of an element, or be used for collections of atoms such as |
| 99 |
< |
methyl and carbonyl groups. The atoms are also capable of having |
| 100 |
< |
directional components associated with them (\emph{e.g.}~permanent |
| 101 |
< |
dipoles). Charges, permanent dipoles, and Lennard-Jones parameters for |
| 102 |
< |
a given atom type are set in the force field parameter files. |
| 95 |
> |
Every {\sc oopse} simulation begins with a Bizarre Atom Simulation |
| 96 |
> |
Syntax ({\sc bass}) file. {\sc bass} is a script syntax that is parsed |
| 97 |
> |
by {\sc oopse} at runtime. The {\sc bass} file allows for the user to |
| 98 |
> |
completely describe the system they wish to simulate, as well as tailor |
| 99 |
> |
{\sc oopse}'s behavior during the simulation. {\sc bass} files are |
| 100 |
> |
denoted with the extension |
| 101 |
> |
\texttt{.bass}, an example file is shown in |
| 102 |
> |
Scheme~\ref{sch:bassExample}. |
| 103 |
|
|
| 104 |
< |
\begin{lstlisting}[float,caption={[Specifier for molecules and atoms] A sample specification of an Ar molecule},label=sch:AtmMole] |
| 104 |
> |
\begin{lstlisting}[float,caption={[An example of a complete {\sc bass} file] An example showing a complete {\sc bass} file.},label={sch:bassExample}] |
| 105 |
> |
|
| 106 |
|
molecule{ |
| 107 |
|
name = "Ar"; |
| 108 |
|
nAtoms = 1; |
| 111 |
|
position( 0.0, 0.0, 0.0 ); |
| 112 |
|
} |
| 113 |
|
} |
| 114 |
+ |
|
| 115 |
+ |
nComponents = 1; |
| 116 |
+ |
component{ |
| 117 |
+ |
type = "Ar"; |
| 118 |
+ |
nMol = 108; |
| 119 |
+ |
} |
| 120 |
+ |
|
| 121 |
+ |
initialConfig = "./argon.init"; |
| 122 |
+ |
|
| 123 |
+ |
forceField = "LJ"; |
| 124 |
+ |
ensemble = "NVE"; // specify the simulation ensemble |
| 125 |
+ |
dt = 1.0; // the time step for integration |
| 126 |
+ |
runTime = 1e3; // the total simulation run time |
| 127 |
+ |
sampleTime = 100; // trajectory file frequency |
| 128 |
+ |
statusTime = 50; // statistics file frequency |
| 129 |
+ |
|
| 130 |
|
\end{lstlisting} |
| 131 |
|
|
| 132 |
+ |
Within the \texttt{.bass} file it is necessary to provide a complete |
| 133 |
+ |
description of the molecule before it is actually placed in the |
| 134 |
+ |
simulation. The {\sc bass} syntax was originally developed with this |
| 135 |
+ |
goal in mind, and allows for the specification of all the atoms in a |
| 136 |
+ |
molecular prototype, as well as any bonds, bends, or torsions. These |
| 137 |
+ |
descriptions can become lengthy for complex molecules, and it would be |
| 138 |
+ |
inconvenient to duplicate the simulation at the beginning of each {\sc |
| 139 |
+ |
bass} script. Addressing this issue {\sc bass} allows for the |
| 140 |
+ |
inclusion of model files at the top of a \texttt{.bass} file. These |
| 141 |
+ |
model files, denoted with the \texttt{.mdl} extension, allow the user |
| 142 |
+ |
to describe a molecular prototype once, then simply include it into |
| 143 |
+ |
each simulation containing that molecule. Returning to the example in |
| 144 |
+ |
Scheme~\ref{sch:bassExample}, the \texttt{.mdl} file's contents would |
| 145 |
+ |
be Scheme~\ref{sch:mdlExample}, and the new \texttt{.bass} file would |
| 146 |
+ |
become Scheme~\ref{sch:bassExPrime}. |
| 147 |
|
|
| 148 |
+ |
\begin{lstlisting}[float,caption={An example \texttt{.mdl} file.},label={sch:mdlExample}] |
| 149 |
+ |
|
| 150 |
+ |
molecule{ |
| 151 |
+ |
name = "Ar"; |
| 152 |
+ |
nAtoms = 1; |
| 153 |
+ |
atom[0]{ |
| 154 |
+ |
type="Ar"; |
| 155 |
+ |
position( 0.0, 0.0, 0.0 ); |
| 156 |
+ |
} |
| 157 |
+ |
} |
| 158 |
+ |
|
| 159 |
+ |
\end{lstlisting} |
| 160 |
+ |
|
| 161 |
+ |
\begin{lstlisting}[float,caption={Revised {\sc bass} example.},label={sch:bassExPrime}] |
| 162 |
+ |
|
| 163 |
+ |
#include "argon.mdl" |
| 164 |
+ |
|
| 165 |
+ |
nComponents = 1; |
| 166 |
+ |
component{ |
| 167 |
+ |
type = "Ar"; |
| 168 |
+ |
nMol = 108; |
| 169 |
+ |
} |
| 170 |
+ |
|
| 171 |
+ |
initialConfig = "./argon.init"; |
| 172 |
+ |
|
| 173 |
+ |
forceField = "LJ"; |
| 174 |
+ |
ensemble = "NVE"; |
| 175 |
+ |
dt = 1.0; |
| 176 |
+ |
runTime = 1e3; |
| 177 |
+ |
sampleTime = 100; |
| 178 |
+ |
statusTime = 50; |
| 179 |
+ |
|
| 180 |
+ |
\end{lstlisting} |
| 181 |
+ |
|
| 182 |
+ |
\subsection{\label{oopseSec:atomsMolecules}Atoms, Molecules and Rigid Bodies} |
| 183 |
+ |
|
| 184 |
+ |
The basic unit of an {\sc oopse} simulation is the atom. The |
| 185 |
+ |
parameters describing the atom are generalized to make the atom as |
| 186 |
+ |
flexible a representation as possible. They may represent specific |
| 187 |
+ |
atoms of an element, or be used for collections of atoms such as |
| 188 |
+ |
methyl and carbonyl groups. The atoms are also capable of having |
| 189 |
+ |
directional components associated with them (\emph{e.g.}~permanent |
| 190 |
+ |
dipoles). Charges, permanent dipoles, and Lennard-Jones parameters for |
| 191 |
+ |
a given atom type are set in the force field parameter files. |
| 192 |
+ |
|
| 193 |
|
Atoms can be collected into secondary structures such as rigid bodies |
| 194 |
|
or molecules. The molecule is a way for {\sc oopse} to keep track of |
| 195 |
|
the atoms in a simulation in logical manner. Molecular units store the |
| 196 |
|
identities of all the atoms and rigid bodies associated with |
| 197 |
|
themselves, and are responsible for the evaluation of their own |
| 198 |
|
internal interactions (\emph{i.e.}~bonds, bends, and torsions). Scheme |
| 199 |
< |
\ref{sch:AtmMole} shows how one creates a molecule in a ``model'' or |
| 199 |
> |
\ref{sch:mdlExample} shows how one creates a molecule in a ``model'' or |
| 200 |
|
\texttt{.mdl} file. The position of the atoms given in the |
| 201 |
|
declaration are relative to the origin of the molecule, and is used |
| 202 |
|
when creating a system containing the molecule. |
| 209 |
|
included in most simulation packages because of the algorithmic |
| 210 |
|
complexity involved in propagating orientational degrees of |
| 211 |
|
freedom. Until recently, integrators which propagate orientational |
| 212 |
< |
motion have been much worse than those available for translational |
| 212 |
> |
motion have had energy conservation problems when compared to those available for translational |
| 213 |
|
motion. |
| 214 |
|
|
| 215 |
|
Moving a rigid body involves determination of both the force and |
| 251 |
|
entire nine parameter rotation matrix. Further discussion |
| 252 |
|
on this choice can be found in Sec.~\ref{oopseSec:integrate}. An example |
| 253 |
|
definition of a rigid body can be seen in Scheme |
| 254 |
< |
\ref{sch:rigidBody}. The positions in the atom definitions are the |
| 179 |
< |
placements of the atoms relative to the origin of the rigid body, |
| 180 |
< |
which itself has a position relative to the origin of the molecule. |
| 254 |
> |
\ref{sch:rigidBody}. |
| 255 |
|
|
| 256 |
< |
\begin{lstlisting}[float,caption={[Defining rigid bodies]A sample definition of a rigid body},label={sch:rigidBody}] |
| 256 |
> |
\begin{lstlisting}[float,caption={[Defining rigid bodies]A sample definition of a molecule containing a rigid body},label={sch:rigidBody}] |
| 257 |
|
molecule{ |
| 258 |
|
name = "TIP3P"; |
| 259 |
|
nAtoms = 3; |
| 276 |
|
members(0, 1, 2); |
| 277 |
|
} |
| 278 |
|
} |
| 279 |
+ |
\end{lstlisting} |
| 280 |
+ |
|
| 281 |
+ |
\subsection{\label{sec:miscConcepts}Putting a Script Together} |
| 282 |
+ |
|
| 283 |
+ |
The actual creation of a {\sc bass} script requires several key components. The first part of the script needs to be the declaration of all of the molecule prototypes used in the simulation. This is typically done through the inclusion of {\tt .mdl} files. Only the molecules actually present in the simulation need to be declared, however {\sc bass} allows for the declaration of more molecules than are needed. This gives the user the ability to build up a library of commonly used molecules into a single {\tt .mdl} file. |
| 284 |
+ |
|
| 285 |
+ |
Once all prototypes are declared, the ordering of the rest of the script is less stringent. Typically, the next to follow the molecular prototypes are the component statements. These statements specify which molecules are present within the simulation. The number of components must first be declared before the first component block statement (an example is seen in Sch.~\ref{sch:bassExPrime}). The component blocks tell {\sc oopse} the number of molecules that will be in the simulation, and the order in which the components blocks are declared sets the ordering of the real atoms within the simulation as well as in the output files. |
| 286 |
+ |
|
| 287 |
+ |
The remainder of the script then sets the various simulation parameters for the system of interest. The required set of parameters that must be present in all simulations is given in Table~\ref{table:reqParams}. The {\tt ensemble} statement is responsible for selecting the integration method used for the calculation of the equations of motion. An in depth discussion of the various methods available in {\sc oopse} can be found in Sec.~\ref{oopseSec:mechanics}. The {\tt forceField} statement is important for the selection of which forces will be used in the course of the simulation. {\sc oopse} supports several force fields, as outlined in Sec.~\ref{oopseSec:empericalEnergy}. The force fields are interchangeable between simulations, with the only requirement being that all atoms needed by the simulation are defined within the selected force field. The time step between force evaluations is set with the {\tt dt} parameter, and {\tt runTime} will set the time length of the simulation. Note, that {\tt runTime} is an absolute time, meaning if the simulation is started at t = 10.0~ns with a {\tt runTime} of 25.0~ns, the simulation will only run for an additional 15.0~ns. The final required parameter, is the {\tt initialConfig} statement. This will set the initial coordinates for the system, as well as the initial time if the {\tt useInitalTime = true;} flag is given. The format of the file specified in {\tt initialConfig}, is given in Sec.~\ref{oopseSec:coordFiles}. Additional parameters are summarized in Table~\ref{table:genParams}. |
| 288 |
+ |
|
| 289 |
+ |
\begin{table} |
| 290 |
+ |
\caption{The Global Keywords: Required Parameters} |
| 291 |
+ |
\label{table:reqParams} |
| 292 |
+ |
\begin{center} |
| 293 |
+ |
% Note when adding or removing columns, the \hsize numbers must add up to the total number |
| 294 |
+ |
% of columns. |
| 295 |
+ |
\begin{tabularx}{\linewidth}% |
| 296 |
+ |
{>{\setlength{\hsize}{1.00\hsize}}X% |
| 297 |
+ |
>{\setlength{\hsize}{0.4\hsize}}X% |
| 298 |
+ |
>{\setlength{\hsize}{1.2\hsize}}X% |
| 299 |
+ |
>{\setlength{\hsize}{1.4\hsize}}X} |
| 300 |
+ |
|
| 301 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 302 |
+ |
|
| 303 |
+ |
{\tt forceField} & string & Sets the force field. & Possible force fields are "DUFF", "LJ", and "EAM". \\ |
| 304 |
+ |
{\tt ensemble} & string & Sets the ensemble. & Possible ensembles are "NVE", "NVT", "NPTi", "NPTf", and "NPTxyz".\\ |
| 305 |
+ |
{\tt dt} & fs & Sets the time step. & Selection of {\tt dt} should be small enough to sample the fastest motion of the simulation. \\ |
| 306 |
+ |
{\tt nComponents} & integer & Sets the number of components. & Needs to appear before the first {\tt Component} block. \\ |
| 307 |
+ |
{\tt initialConfig} & string & Sets the file containing the initial configuration. & Can point to any file containing the configuration in the correct order. \\ |
| 308 |
+ |
{\tt runTime} & fs & Sets the time at which the simulation should end. & This is an absolute time, and will end the simulation when the current time meets or exceeds the {\tt runTime}. \\ |
| 309 |
+ |
|
| 310 |
+ |
|
| 311 |
+ |
\end{tabularx} |
| 312 |
+ |
\end{center} |
| 313 |
+ |
\end{table} |
| 314 |
+ |
|
| 315 |
+ |
\begin{table} |
| 316 |
+ |
\caption{The Global Keywords: General Parameters} |
| 317 |
+ |
\label{table:genParams} |
| 318 |
+ |
\begin{center} |
| 319 |
+ |
% Note when adding or removing columns, the \hsize numbers must add up to the total number |
| 320 |
+ |
% of columns. |
| 321 |
+ |
\begin{tabularx}{\linewidth}% |
| 322 |
+ |
{>{\setlength{\hsize}{1.00\hsize}}X% |
| 323 |
+ |
>{\setlength{\hsize}{0.4\hsize}}X% |
| 324 |
+ |
>{\setlength{\hsize}{1.2\hsize}}X% |
| 325 |
+ |
>{\setlength{\hsize}{1.4\hsize}}X} |
| 326 |
+ |
|
| 327 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 328 |
+ |
|
| 329 |
+ |
{\tt finalConfig} & string & Option to set the name of the final output file. & Useful when stringing simulations together. Defaults to the {\tt .bass} file with an {\tt .eor} extension. \\ |
| 330 |
+ |
{\tt useInitialTime} & logical & Sets whether the initial time is taken from the {\tt .init} file. & Useful when recovering a simulation from a crashed processor. Default is false. \\ |
| 331 |
+ |
{\tt sampleTime} & fs & Sets the frequency at which the {\tt .dump} file is written. & Default sets the frequency to the {\tt runTime}. \\ |
| 332 |
+ |
{\tt statusTime} & fs & Sets the frequency at which the {\tt .stat} file is written. & Defaults sets the frequency to the {\tt sampleTime}. \\ |
| 333 |
+ |
{\tt LJrcut} & $\mbox{\AA}$ & Manually sets the Lennard-Jones cutoff. & Defaults to $2.5\sigma_L$, where $\sigma_L$ is the largest LJ $\sigma$ in the simulation. \\ |
| 334 |
+ |
{\tt electrostaticCutoffRadius}& & & \\ |
| 335 |
+ |
& $\mbox{\AA}$ & Manually sets the cutoff used by the electrostatic potentials. & Defaults to $15\mbox{\AA}$ \\ |
| 336 |
+ |
{\tt electrostaticSkinThickness} & & & \\ |
| 337 |
+ |
& $\mbox{\AA}$ & Manually sets the skin thickness for the electrostatic switching function. & Defaults to 5~\% of the {\tt electrostaticSkinThickness}. \\ |
| 338 |
+ |
{\tt useReactionField} & logical & Turns the reaction field correction on/off. & Default is "false". \\ |
| 339 |
+ |
{\tt dielectric} & unitless & Sets the dielectric constant for reaction field. & If {\tt useReactionField} is true, then {\tt dielectric} must be set. \\ |
| 340 |
+ |
{\tt usePeriodicBoundaryConditions} & & & \\ |
| 341 |
+ |
& logical & Turns periodic boundary conditions on/off. & Default is "true". \\ |
| 342 |
+ |
{\tt seed } & integer & Sets the seed value for the random number generator. & The seed needs to be at least 9 digits long. The default is to take the seed from the CPU clock. |
| 343 |
+ |
|
| 344 |
+ |
\end{tabularx} |
| 345 |
+ |
\end{center} |
| 346 |
+ |
\end{table} |
| 347 |
+ |
|
| 348 |
+ |
|
| 349 |
+ |
|
| 350 |
+ |
\subsection{\label{oopseSec:coordFiles}Coordinate Files} |
| 351 |
+ |
|
| 352 |
+ |
The standard format for storage of a systems coordinates is a modified |
| 353 |
+ |
xyz-file syntax, the exact details of which can be seen in |
| 354 |
+ |
Scheme~\ref{sch:dumpFormat}. As all bonding and molecular information |
| 355 |
+ |
is stored in the \texttt{.bass} and \texttt{.mdl} files, the |
| 356 |
+ |
coordinate files are simply the complete set of coordinates for each |
| 357 |
+ |
atom at a given simulation time. One important note, although the |
| 358 |
+ |
simulation propagates the complete rotation matrix, directional |
| 359 |
+ |
entities are written out using quanternions, to save space in the |
| 360 |
+ |
output files. |
| 361 |
+ |
|
| 362 |
+ |
\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, body fixed angular momentum.},label=sch:dumpFormat] |
| 363 |
+ |
|
| 364 |
+ |
nAtoms |
| 365 |
+ |
time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz; |
| 366 |
+ |
Name1 x y z vx vy vz q0 q1 q2 q3 jx jy jz |
| 367 |
+ |
Name2 x y z vx vy vz q0 q1 q2 q3 jx jy jz |
| 368 |
+ |
etc... |
| 369 |
+ |
|
| 370 |
|
\end{lstlisting} |
| 371 |
|
|
| 372 |
+ |
|
| 373 |
+ |
There are three major files used by {\sc oopse} written in the |
| 374 |
+ |
coordinate format, they are as follows: the initialization file |
| 375 |
+ |
(\texttt{.init}), the simulation trajectory file (\texttt{.dump}), and |
| 376 |
+ |
the final coordinates of the simulation (\texttt{.eor}). The initialization file is |
| 377 |
+ |
necessary for {\sc oopse} to start the simulation with the proper |
| 378 |
+ |
coordinates, and is generated before the simulation run. The |
| 379 |
+ |
trajectory file is created at the beginning of the simulation, and is |
| 380 |
+ |
used to store snapshots of the simulation at regular intervals. The |
| 381 |
+ |
first frame is a duplication of the |
| 382 |
+ |
\texttt{.init} file, and each subsequent frame is appended to the file |
| 383 |
+ |
at an interval specified in the \texttt{.bass} file with the |
| 384 |
+ |
\texttt{sampleTime} flag. The final coordinate file is the end of run file. The |
| 385 |
+ |
\texttt{.eor} file stores the final configuration of the system for a |
| 386 |
+ |
given simulation. The file is updated at the same time as the |
| 387 |
+ |
\texttt{.dump} file, however, it only contains the most recent |
| 388 |
+ |
frame. In this way, an \texttt{.eor} file may be used as the |
| 389 |
+ |
initialization file to a second simulation in order to continue a |
| 390 |
+ |
simulation or recover one from a processor that has crashed during the |
| 391 |
+ |
course of the run. |
| 392 |
+ |
|
| 393 |
+ |
\subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates} |
| 394 |
+ |
|
| 395 |
+ |
As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization |
| 396 |
+ |
file is needed to provide the starting coordinates for a |
| 397 |
+ |
simulation. Several helper programs are provided with {\sc oopse} to illustrate possible build routes. However, as each simulation is different, system creation is left to the end user. The {\tt .init} file must list the atoms in the correct order or {\sc oopse} will give an atom mismatch error. |
| 398 |
+ |
|
| 399 |
+ |
The correct ordering of the atoms relies on the ordering of atoms and molecules within the model and {\sc bass} scripts. {\sc oopse} expects the order to comply with the following guidelines: |
| 400 |
+ |
\begin{enumerate} |
| 401 |
+ |
\item All of the molecules of the first declared component are given before proceeding to the molecules of the second component, and so on for all declared components. |
| 402 |
+ |
\item The ordering of the atoms for each molecule follows the order declared in the molecule's declaration within the model file. |
| 403 |
+ |
\end{enumerate} |
| 404 |
+ |
An example is given in Scheme~\ref{sch:initEx1} resulting in the {\tt .init} file shown in Scheme~\ref{sch:initEx2}. |
| 405 |
+ |
|
| 406 |
+ |
\begin{lstlisting}[float,caption={This scheme illustrates the declaration of the $\text{I}_2$ molecule and the HCl molecule. The two molecules are then included into a simulation.}, label=sch:initEx1] |
| 407 |
+ |
|
| 408 |
+ |
molecule{ |
| 409 |
+ |
name = "I2"; |
| 410 |
+ |
nAtoms = 2; |
| 411 |
+ |
atom[0]{ |
| 412 |
+ |
type = "I"; |
| 413 |
+ |
} |
| 414 |
+ |
atom[1]{ |
| 415 |
+ |
type = "I"; |
| 416 |
+ |
} |
| 417 |
+ |
nBonds = 1; |
| 418 |
+ |
bond[0]{ |
| 419 |
+ |
members( 0, 1); |
| 420 |
+ |
} |
| 421 |
+ |
} |
| 422 |
+ |
|
| 423 |
+ |
molecule{ |
| 424 |
+ |
name = "HCl" |
| 425 |
+ |
nAtoms = 2; |
| 426 |
+ |
atom[0]{ |
| 427 |
+ |
type = "H"; |
| 428 |
+ |
} |
| 429 |
+ |
atom[1]{ |
| 430 |
+ |
type = "Cl"; |
| 431 |
+ |
} |
| 432 |
+ |
nBonds = 1; |
| 433 |
+ |
bond[0]{ |
| 434 |
+ |
members( 0, 1); |
| 435 |
+ |
} |
| 436 |
+ |
} |
| 437 |
+ |
|
| 438 |
+ |
nComponents = 2; |
| 439 |
+ |
component{ |
| 440 |
+ |
type = "HCl"; |
| 441 |
+ |
nMol = 4; |
| 442 |
+ |
} |
| 443 |
+ |
component{ |
| 444 |
+ |
type = "I2"; |
| 445 |
+ |
nMol = 1; |
| 446 |
+ |
} |
| 447 |
+ |
|
| 448 |
+ |
initialConfig = "mixture.init"; |
| 449 |
+ |
|
| 450 |
+ |
\end{lstlisting} |
| 451 |
+ |
|
| 452 |
+ |
\begin{lstlisting}[float,caption={This is the contents of the {\tt mixture.init} file matching the declarations in Scheme~\ref{sch:initEx1}. Note that even though $\text{I}_2$ is declared before HCl, the {\tt .init} file follows the order in which the components were included.},label=sch:initEx2] |
| 453 |
+ |
|
| 454 |
+ |
10 |
| 455 |
+ |
0.0; 10.0 0.0 0.0; 0.0 10.0 0.0; 0.0 0.0 10.0; |
| 456 |
+ |
H ... |
| 457 |
+ |
Cl ... |
| 458 |
+ |
H ... |
| 459 |
+ |
Cl ... |
| 460 |
+ |
H ... |
| 461 |
+ |
Cl ... |
| 462 |
+ |
H ... |
| 463 |
+ |
Cl ... |
| 464 |
+ |
I ... |
| 465 |
+ |
I ... |
| 466 |
+ |
|
| 467 |
+ |
\end{lstlisting} |
| 468 |
+ |
|
| 469 |
+ |
|
| 470 |
+ |
\subsection{The Statistics File} |
| 471 |
+ |
|
| 472 |
+ |
The last output file generated by {\sc oopse} is the statistics |
| 473 |
+ |
file. This file records such statistical quantities as the |
| 474 |
+ |
instantaneous temperature, volume, pressure, etc. It is written out |
| 475 |
+ |
with the frequency specified in the \texttt{.bass} file with the |
| 476 |
+ |
\texttt{statusTime} keyword. The file allows the user to observe the |
| 477 |
+ |
system variables as a function of simulation time while the simulation |
| 478 |
+ |
is in progress. One useful function the statistics file serves is to |
| 479 |
+ |
monitor the conserved quantity of a given simulation ensemble, this |
| 480 |
+ |
allows the user to observe the stability of the integrator. The |
| 481 |
+ |
statistics file is denoted with the \texttt{.stat} file extension. |
| 482 |
+ |
|
| 483 |
+ |
|
| 484 |
+ |
\section{\label{oopseSec:empiricalEnergy}The Empirical Energy Functions} |
| 485 |
+ |
|
| 486 |
+ |
\ |
| 487 |
|
\subsection{\label{sec:LJPot}The Lennard Jones Force Field} |
| 488 |
|
|
| 489 |
|
The most basic force field implemented in {\sc oopse} is the |
| 561 |
|
charges. Charge-neutral distributions were replaced with dipoles, |
| 562 |
|
while most atoms and groups of atoms were reduced to Lennard-Jones |
| 563 |
|
interaction sites. This simplification cuts the length scale of long |
| 564 |
< |
range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, and allows |
| 285 |
< |
us to avoid the computationally expensive Ewald sum. Instead, we can |
| 286 |
< |
use neighbor-lists and cutoff radii for the dipolar interactions, or |
| 287 |
< |
include a reaction field to mimic larger range interactions. |
| 564 |
> |
range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, removing the need for the computationally expensive Ewald sum. Instead, we Verlet neighbor-lists and cutoff radii are used for the dipolar interactions, or a reaction field is added to mimic longer range interactions. |
| 565 |
|
|
| 566 |
|
As an example, lipid head-groups in {\sc duff} are represented as |
| 567 |
|
point dipole interaction sites. By placing a dipole at the head |
| 571 |
|
site is located at the pseudoatom's center of mass. The model is |
| 572 |
|
illustrated by the red atom in Fig.~\ref{oopseFig:lipidModel}. The |
| 573 |
|
water model we use to complement the dipoles of the lipids is our |
| 574 |
< |
reparameterization of the soft sticky dipole (SSD) model of Ichiye |
| 574 |
> |
reparameterization\cite{fennell04} of the soft sticky dipole (SSD) model of Ichiye |
| 575 |
|
\emph{et al.}\cite{liu96:new_model} |
| 576 |
|
|
| 577 |
|
\begin{figure} |
| 594 |
|
$\text{CH}_2$ do not change depending on what species are bonded to |
| 595 |
|
it. |
| 596 |
|
|
| 597 |
< |
TraPPE also constrains all bonds to be of fixed length. Typically, |
| 597 |
> |
TraPPE and {\sc duff} also constrain all bonds to be of fixed length. Typically, |
| 598 |
|
bond vibrations are the fastest motions in a molecular dynamic |
| 599 |
|
simulation. Small time steps between force evaluations must be used to |
| 600 |
|
ensure adequate energy conservation in the bond degrees of freedom. By |
| 624 |
|
|
| 625 |
|
\end{lstlisting} |
| 626 |
|
|
| 627 |
< |
\subsection{\label{oopseSec:energyFunctions}{\sc duff} Energy Functions} |
| 627 |
> |
\subsubsection{\label{oopseSec:energyFunctions}{\sc duff} Energy Functions} |
| 628 |
|
|
| 629 |
|
The total potential energy function in {\sc duff} is |
| 630 |
|
\begin{equation} |
| 753 |
|
is the taper radius some given thickness less than the electrostatic |
| 754 |
|
cutoff. The switching thickness can be set in the \texttt{.bass} file. |
| 755 |
|
|
| 756 |
< |
\subsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF} |
| 756 |
> |
\subsubsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF} |
| 757 |
|
|
| 758 |
|
In the interest of computational efficiency, the default solvent used |
| 759 |
|
by {\sc oopse} is the extended Soft Sticky Dipole (SSD/E) water |
| 871 |
|
|
| 872 |
|
\subsection{\label{oopseSec:eam}Embedded Atom Method} |
| 873 |
|
|
| 874 |
< |
There are Molecular Dynamics packages which have the |
| 875 |
< |
capacity to simulate metallic systems, including some that have |
| 876 |
< |
parallel computational abilities\cite{plimpton93}. Potentials that |
| 600 |
< |
describe bonding transition metal |
| 601 |
< |
systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have an |
| 602 |
< |
attractive interaction which models ``Embedding'' |
| 874 |
> |
{\sc oopse} implements a potential that |
| 875 |
> |
describes bonding transition metal |
| 876 |
> |
systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} and has attractive interaction which models ``Embedding'' |
| 877 |
|
a positively charged metal ion in the electron density due to the |
| 878 |
|
free valance ``sea'' of electrons created by the surrounding atoms in |
| 879 |
|
the system. A mostly-repulsive pairwise part of the potential |
| 965 |
|
In this way, particles are allowed to diffuse freely in $\mathbf{r}$, |
| 966 |
|
but their minimum images, $\mathbf{r}^{\prime}$ are used to compute |
| 967 |
|
the inter-atomic forces. |
| 694 |
– |
|
| 695 |
– |
|
| 696 |
– |
\section{\label{oopseSec:IOfiles}Input and Output Files} |
| 697 |
– |
|
| 698 |
– |
\subsection{{\sc bass} and Model Files} |
| 699 |
– |
|
| 700 |
– |
Every {\sc oopse} simulation begins with a Bizarre Atom Simulation |
| 701 |
– |
Syntax ({\sc bass}) file. {\sc bass} is a script syntax that is parsed |
| 702 |
– |
by {\sc oopse} at runtime. The {\sc bass} file allows for the user to |
| 703 |
– |
completely describe the system they wish to simulate, as well as tailor |
| 704 |
– |
{\sc oopse}'s behavior during the simulation. {\sc bass} files are |
| 705 |
– |
denoted with the extension |
| 706 |
– |
\texttt{.bass}, an example file is shown in |
| 707 |
– |
Scheme~\ref{sch:bassExample}. |
| 708 |
– |
|
| 709 |
– |
\begin{lstlisting}[float,caption={[An example of a complete {\sc bass} file] An example showing a complete {\sc bass} file.},label={sch:bassExample}] |
| 710 |
– |
|
| 711 |
– |
molecule{ |
| 712 |
– |
name = "Ar"; |
| 713 |
– |
nAtoms = 1; |
| 714 |
– |
atom[0]{ |
| 715 |
– |
type="Ar"; |
| 716 |
– |
position( 0.0, 0.0, 0.0 ); |
| 717 |
– |
} |
| 718 |
– |
} |
| 719 |
– |
|
| 720 |
– |
nComponents = 1; |
| 721 |
– |
component{ |
| 722 |
– |
type = "Ar"; |
| 723 |
– |
nMol = 108; |
| 724 |
– |
} |
| 725 |
– |
|
| 726 |
– |
initialConfig = "./argon.init"; |
| 727 |
– |
|
| 728 |
– |
forceField = "LJ"; |
| 729 |
– |
ensemble = "NVE"; // specify the simulation ensemble |
| 730 |
– |
dt = 1.0; // the time step for integration |
| 731 |
– |
runTime = 1e3; // the total simulation run time |
| 732 |
– |
sampleTime = 100; // trajectory file frequency |
| 733 |
– |
statusTime = 50; // statistics file frequency |
| 734 |
– |
|
| 735 |
– |
\end{lstlisting} |
| 736 |
– |
|
| 737 |
– |
Within the \texttt{.bass} file it is necessary to provide a complete |
| 738 |
– |
description of the molecule before it is actually placed in the |
| 739 |
– |
simulation. The {\sc bass} syntax was originally developed with this |
| 740 |
– |
goal in mind, and allows for the specification of all the atoms in a |
| 741 |
– |
molecular prototype, as well as any bonds, bends, or torsions. These |
| 742 |
– |
descriptions can become lengthy for complex molecules, and it would be |
| 743 |
– |
inconvenient to duplicate the simulation at the beginning of each {\sc |
| 744 |
– |
bass} script. Addressing this issue {\sc bass} allows for the |
| 745 |
– |
inclusion of model files at the top of a \texttt{.bass} file. These |
| 746 |
– |
model files, denoted with the \texttt{.mdl} extension, allow the user |
| 747 |
– |
to describe a molecular prototype once, then simply include it into |
| 748 |
– |
each simulation containing that molecule. Returning to the example in |
| 749 |
– |
Scheme~\ref{sch:bassExample}, the \texttt{.mdl} file's contents would |
| 750 |
– |
be Scheme~\ref{sch:mdlExample}, and the new \texttt{.bass} file would |
| 751 |
– |
become Scheme~\ref{sch:bassExPrime}. |
| 752 |
– |
|
| 753 |
– |
\begin{lstlisting}[float,caption={An example \texttt{.mdl} file.},label={sch:mdlExample}] |
| 754 |
– |
|
| 755 |
– |
molecule{ |
| 756 |
– |
name = "Ar"; |
| 757 |
– |
nAtoms = 1; |
| 758 |
– |
atom[0]{ |
| 759 |
– |
type="Ar"; |
| 760 |
– |
position( 0.0, 0.0, 0.0 ); |
| 761 |
– |
} |
| 762 |
– |
} |
| 763 |
– |
|
| 764 |
– |
\end{lstlisting} |
| 765 |
– |
|
| 766 |
– |
\begin{lstlisting}[float,caption={Revised {\sc bass} example.},label={sch:bassExPrime}] |
| 767 |
– |
|
| 768 |
– |
#include "argon.mdl" |
| 769 |
– |
|
| 770 |
– |
nComponents = 1; |
| 771 |
– |
component{ |
| 772 |
– |
type = "Ar"; |
| 773 |
– |
nMol = 108; |
| 774 |
– |
} |
| 775 |
– |
|
| 776 |
– |
initialConfig = "./argon.init"; |
| 777 |
– |
|
| 778 |
– |
forceField = "LJ"; |
| 779 |
– |
ensemble = "NVE"; |
| 780 |
– |
dt = 1.0; |
| 781 |
– |
runTime = 1e3; |
| 782 |
– |
sampleTime = 100; |
| 783 |
– |
statusTime = 50; |
| 784 |
– |
|
| 785 |
– |
\end{lstlisting} |
| 786 |
– |
|
| 787 |
– |
\subsection{\label{oopseSec:coordFiles}Coordinate Files} |
| 788 |
– |
|
| 789 |
– |
The standard format for storage of a systems coordinates is a modified |
| 790 |
– |
xyz-file syntax, the exact details of which can be seen in |
| 791 |
– |
Scheme~\ref{sch:dumpFormat}. As all bonding and molecular information |
| 792 |
– |
is stored in the \texttt{.bass} and \texttt{.mdl} files, the |
| 793 |
– |
coordinate files are simply the complete set of coordinates for each |
| 794 |
– |
atom at a given simulation time. One important note, although the |
| 795 |
– |
simulation propagates the complete rotation matrix, directional |
| 796 |
– |
entities are written out using quanternions, to save space in the |
| 797 |
– |
output files. |
| 798 |
– |
|
| 799 |
– |
\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] |
| 800 |
– |
|
| 801 |
– |
nAtoms |
| 802 |
– |
time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz; |
| 803 |
– |
Name1 x y z vx vy vz q0 q1 q2 q3 jx jy jz |
| 804 |
– |
Name2 x y z vx vy vz q0 q1 q2 q3 jx jy jz |
| 805 |
– |
etc... |
| 806 |
– |
|
| 807 |
– |
\end{lstlisting} |
| 808 |
– |
|
| 809 |
– |
|
| 810 |
– |
There are three major files used by {\sc oopse} written in the |
| 811 |
– |
coordinate format, they are as follows: the initialization file |
| 812 |
– |
(\texttt{.init}), the simulation trajectory file (\texttt{.dump}), and |
| 813 |
– |
the final coordinates of the simulation. The initialization file is |
| 814 |
– |
necessary for {\sc oopse} to start the simulation with the proper |
| 815 |
– |
coordinates, and is generated before the simulation run. The |
| 816 |
– |
trajectory file is created at the beginning of the simulation, and is |
| 817 |
– |
used to store snapshots of the simulation at regular intervals. The |
| 818 |
– |
first frame is a duplication of the |
| 819 |
– |
\texttt{.init} file, and each subsequent frame is appended to the file |
| 820 |
– |
at an interval specified in the \texttt{.bass} file with the |
| 821 |
– |
\texttt{sampleTime} flag. The final coordinate file is the end of run file. The |
| 822 |
– |
\texttt{.eor} file stores the final configuration of the system for a |
| 823 |
– |
given simulation. The file is updated at the same time as the |
| 824 |
– |
\texttt{.dump} file, however, it only contains the most recent |
| 825 |
– |
frame. In this way, an \texttt{.eor} file may be used as the |
| 826 |
– |
initialization file to a second simulation in order to continue a |
| 827 |
– |
simulation or recover one from a processor that has crashed during the |
| 828 |
– |
course of the run. |
| 968 |
|
|
| 830 |
– |
\subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates} |
| 969 |
|
|
| 832 |
– |
As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization |
| 833 |
– |
file is needed to provide the starting coordinates for a |
| 834 |
– |
simulation. The {\sc oopse} package provides several system building |
| 835 |
– |
programs to aid in the creation of the \texttt{.init} |
| 836 |
– |
file. The programs use {\sc bass}, and will recognize |
| 837 |
– |
arguments and parameters in the \texttt{.bass} file that would |
| 838 |
– |
otherwise be ignored by the simulation. |
| 970 |
|
|
| 840 |
– |
\subsection{The Statistics File} |
| 841 |
– |
|
| 842 |
– |
The last output file generated by {\sc oopse} is the statistics |
| 843 |
– |
file. This file records such statistical quantities as the |
| 844 |
– |
instantaneous temperature, volume, pressure, etc. It is written out |
| 845 |
– |
with the frequency specified in the \texttt{.bass} file with the |
| 846 |
– |
\texttt{statusTime} keyword. The file allows the user to observe the |
| 847 |
– |
system variables as a function of simulation time while the simulation |
| 848 |
– |
is in progress. One useful function the statistics file serves is to |
| 849 |
– |
monitor the conserved quantity of a given simulation ensemble, this |
| 850 |
– |
allows the user to observe the stability of the integrator. The |
| 851 |
– |
statistics file is denoted with the \texttt{.stat} file extension. |
| 852 |
– |
|
| 971 |
|
\section{\label{oopseSec:mechanics}Mechanics} |
| 972 |
|
|
| 973 |
|
\subsection{\label{oopseSec:integrate}Integrating the Equations of Motion: the |
| 1432 |
|
|
| 1433 |
|
To carry out isobaric-isothermal ensemble calculations {\sc oopse} |
| 1434 |
|
implements the Melchionna modifications to the Nos\'e-Hoover-Andersen |
| 1435 |
< |
equations of motion,\cite{melchionna93} |
| 1435 |
> |
equations of motion.\cite{melchionna93} The equations of motion are the same as NVT with the following exceptions: |
| 1436 |
|
|
| 1437 |
|
\begin{eqnarray} |
| 1438 |
|
\dot{{\bf r}} & = & {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right), \\ |
| 1439 |
|
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\eta + \chi) {\bf v}, \\ |
| 1322 |
– |
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 1323 |
– |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right),\\ |
| 1324 |
– |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 1325 |
– |
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 1326 |
– |
V}{\partial \mathsf{A}} \right) - \chi {\bf j}, \\ |
| 1327 |
– |
\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left( |
| 1328 |
– |
\frac{T}{T_{\mathrm{target}}} - 1 \right) ,\\ |
| 1440 |
|
\dot{\eta} & = & \frac{1}{\tau_{B}^2 f k_B T_{\mathrm{target}}} V \left( P - |
| 1441 |
|
P_{\mathrm{target}} \right), \\ |
| 1442 |
|
\dot{\mathcal{V}} & = & 3 \mathcal{V} \eta . |
| 1485 |
|
file. The units for {\tt tauBarostat} are fs, and the units for the |
| 1486 |
|
{\tt targetPressure} are atmospheres. Like in the NVT integrator, the |
| 1487 |
|
integration of the equations of motion is carried out in a |
| 1488 |
< |
velocity-Verlet style 2 part algorithm: |
| 1488 |
> |
velocity-Verlet style 2 part algorithm with only the following differences: |
| 1489 |
|
|
| 1490 |
|
{\tt moveA:} |
| 1491 |
|
\begin{align*} |
| 1381 |
– |
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\ |
| 1382 |
– |
% |
| 1492 |
|
P(t) &\leftarrow \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\} ,\\ |
| 1493 |
|
% |
| 1494 |
|
{\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t) |
| 1495 |
|
+ \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) |
| 1496 |
|
\left(\chi(t) + \eta(t) \right) \right), \\ |
| 1497 |
|
% |
| 1389 |
– |
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
| 1390 |
– |
+ \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 1391 |
– |
\chi(t) \right), \\ |
| 1392 |
– |
% |
| 1393 |
– |
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h * |
| 1394 |
– |
{\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1} |
| 1395 |
– |
\right) ,\\ |
| 1396 |
– |
% |
| 1397 |
– |
\chi\left(t + h / 2 \right) &\leftarrow \chi(t) + |
| 1398 |
– |
\frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1 |
| 1399 |
– |
\right) ,\\ |
| 1400 |
– |
% |
| 1498 |
|
\eta(t + h / 2) &\leftarrow \eta(t) + \frac{h |
| 1499 |
|
\mathcal{V}(t)}{2 N k_B T(t) \tau_B^2} \left( P(t) |
| 1500 |
|
- P_{\mathrm{target}} \right), \\ |
| 1508 |
|
\mathsf{H}(t). |
| 1509 |
|
\end{align*} |
| 1510 |
|
|
| 1511 |
< |
Most of these equations are identical to their counterparts in the NVT |
| 1415 |
< |
integrator, but the propagation of positions to time $t + h$ |
| 1511 |
> |
The propagation of positions to time $t + h$ |
| 1512 |
|
depends on the positions at the same time. {\sc oopse} carries out |
| 1513 |
|
this step iteratively (with a limit of 5 passes through the iterative |
| 1514 |
|
loop). Also, the simulation box $\mathsf{H}$ is scaled uniformly for |
| 1528 |
|
|
| 1529 |
|
{\tt moveB:} |
| 1530 |
|
\begin{align*} |
| 1435 |
– |
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, |
| 1436 |
– |
\left\{{\bf j}(t + h)\right\} ,\\ |
| 1437 |
– |
% |
| 1531 |
|
P(t + h) &\leftarrow \left\{{\bf r}(t + h)\right\}, |
| 1532 |
|
\left\{{\bf v}(t + h)\right\}, \\ |
| 1533 |
|
% |
| 1441 |
– |
\chi\left(t + h \right) &\leftarrow \chi\left(t + h / |
| 1442 |
– |
2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+h)} |
| 1443 |
– |
{T_{\mathrm{target}}} - 1 \right), \\ |
| 1444 |
– |
% |
| 1534 |
|
\eta(t + h) &\leftarrow \eta(t + h / 2) + |
| 1535 |
|
\frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h) |
| 1536 |
|
\tau_B^2} \left( P(t + h) - P_{\mathrm{target}} \right), \\ |
| 1587 |
|
{\it shape} as well as in the volume of the box. This method utilizes |
| 1588 |
|
the full $3 \times 3$ pressure tensor and introduces a tensor of |
| 1589 |
|
extended variables ($\overleftrightarrow{\eta}$) to control changes to |
| 1590 |
< |
the box shape. The equations of motion for this method are |
| 1590 |
> |
the box shape. The equations of motion for this method differ from those of NPTi as follows: |
| 1591 |
|
\begin{eqnarray} |
| 1592 |
|
\dot{{\bf r}} & = & {\bf v} + \overleftrightarrow{\eta} \cdot \left( {\bf r} - {\bf R}_0 \right), \\ |
| 1593 |
|
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\overleftrightarrow{\eta} + |
| 1594 |
|
\chi \cdot \mathsf{1}) {\bf v}, \\ |
| 1506 |
– |
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 1507 |
– |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) ,\\ |
| 1508 |
– |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 1509 |
– |
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 1510 |
– |
V}{\partial \mathsf{A}} \right) - \chi {\bf j} ,\\ |
| 1511 |
– |
\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left( |
| 1512 |
– |
\frac{T}{T_{\mathrm{target}}} - 1 \right) ,\\ |
| 1595 |
|
\dot{\overleftrightarrow{\eta}} & = & \frac{1}{\tau_{B}^2 f k_B |
| 1596 |
|
T_{\mathrm{target}}} V \left( \overleftrightarrow{\mathsf{P}} - P_{\mathrm{target}}\mathsf{1} \right) ,\\ |
| 1597 |
|
\dot{\mathsf{H}} & = & \overleftrightarrow{\eta} \cdot \mathsf{H} . |
| 1607 |
|
|
| 1608 |
|
{\tt moveA:} |
| 1609 |
|
\begin{align*} |
| 1528 |
– |
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\ |
| 1529 |
– |
% |
| 1610 |
|
\overleftrightarrow{\mathsf{P}}(t) &\leftarrow \left\{{\bf r}(t)\right\}, |
| 1611 |
|
\left\{{\bf v}(t)\right\} ,\\ |
| 1612 |
|
% |
| 1615 |
|
\left(\chi(t)\mathsf{1} + \overleftrightarrow{\eta}(t) \right) \cdot |
| 1616 |
|
{\bf v}(t) \right), \\ |
| 1617 |
|
% |
| 1538 |
– |
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
| 1539 |
– |
+ \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 1540 |
– |
\chi(t) \right), \\ |
| 1541 |
– |
% |
| 1542 |
– |
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h * |
| 1543 |
– |
{\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1} |
| 1544 |
– |
\right), \\ |
| 1545 |
– |
% |
| 1546 |
– |
\chi\left(t + h / 2 \right) &\leftarrow \chi(t) + |
| 1547 |
– |
\frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} |
| 1548 |
– |
- 1 \right), \\ |
| 1549 |
– |
% |
| 1618 |
|
\overleftrightarrow{\eta}(t + h / 2) &\leftarrow |
| 1619 |
|
\overleftrightarrow{\eta}(t) + \frac{h \mathcal{V}(t)}{2 N k_B |
| 1620 |
|
T(t) \tau_B^2} \left( \overleftrightarrow{\mathsf{P}}(t) |
| 1636 |
|
|
| 1637 |
|
{\tt moveB:} |
| 1638 |
|
\begin{align*} |
| 1571 |
– |
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, |
| 1572 |
– |
\left\{{\bf j}(t + h)\right\}, \\ |
| 1573 |
– |
% |
| 1639 |
|
\overleftrightarrow{\mathsf{P}}(t + h) &\leftarrow \left\{{\bf r} |
| 1640 |
|
(t + h)\right\}, \left\{{\bf v}(t |
| 1641 |
|
+ h)\right\}, \left\{{\bf f}(t + h)\right\} ,\\ |
| 1642 |
|
% |
| 1578 |
– |
\chi\left(t + h \right) &\leftarrow \chi\left(t + h / |
| 1579 |
– |
2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+ |
| 1580 |
– |
h)}{T_{\mathrm{target}}} - 1 \right), \\ |
| 1581 |
– |
% |
| 1643 |
|
\overleftrightarrow{\eta}(t + h) &\leftarrow |
| 1644 |
|
\overleftrightarrow{\eta}(t + h / 2) + |
| 1645 |
|
\frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h) |
| 1651 |
|
\frac{{\bf f}(t + h)}{m} - |
| 1652 |
|
(\chi(t + h)\mathsf{1} + \overleftrightarrow{\eta}(t |
| 1653 |
|
+ h)) \right) \cdot {\bf v}(t + h), \\ |
| 1593 |
– |
% |
| 1594 |
– |
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t |
| 1595 |
– |
+ h / 2 \right) + \frac{h}{2} \left( {\bf \tau}^b(t |
| 1596 |
– |
+ h) - {\bf j}(t + h) \chi(t + h) \right) . |
| 1654 |
|
\end{align*} |
| 1655 |
|
|
| 1656 |
|
The iterative schemes for both {\tt moveA} and {\tt moveB} are |
| 1694 |
|
orientational anisotropy in the system (i.e. in lipid bilayer |
| 1695 |
|
simulations). |
| 1696 |
|
|
| 1697 |
< |
\subsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond |
| 1697 |
> |
\subsection{\label{sec:constraints}Constraint Methods} |
| 1698 |
> |
|
| 1699 |
> |
\subsubsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond |
| 1700 |
|
Constraints} |
| 1701 |
|
|
| 1702 |
|
In order to satisfy the constraints of fixed bond lengths within {\sc |
| 1703 |
|
oopse}, we have implemented the {\sc rattle} algorithm of |
| 1704 |
|
Andersen.\cite{andersen83} The algorithm is a velocity verlet |
| 1705 |
|
formulation of the {\sc shake} method\cite{ryckaert77} of iteratively |
| 1706 |
< |
solving the Lagrange multipliers of constraint. The system of Lagrange |
| 1648 |
< |
multipliers allows one to reformulate the equations of motion with |
| 1649 |
< |
explicit constraint forces.\cite{fowles99:lagrange} |
| 1706 |
> |
solving the Lagrange multipliers of constraint. |
| 1707 |
|
|
| 1708 |
< |
Consider a system described by coordinates $q_1$ and $q_2$ subject to an |
| 1652 |
< |
equation of constraint: |
| 1653 |
< |
\begin{equation} |
| 1654 |
< |
\sigma(q_1, q_2,t) = 0 |
| 1655 |
< |
\label{oopseEq:lm1} |
| 1656 |
< |
\end{equation} |
| 1657 |
< |
The Lagrange formulation of the equations of motion can be written: |
| 1658 |
< |
\begin{equation} |
| 1659 |
< |
\delta\int_{t_1}^{t_2}L\, dt = |
| 1660 |
< |
\int_{t_1}^{t_2} \sum_i \biggl [ \frac{\partial L}{\partial q_i} |
| 1661 |
< |
- \frac{d}{dt}\biggl(\frac{\partial L}{\partial \dot{q}_i} |
| 1662 |
< |
\biggr ) \biggr] \delta q_i \, dt = 0. |
| 1663 |
< |
\label{oopseEq:lm2} |
| 1664 |
< |
\end{equation} |
| 1665 |
< |
Here, $\delta q_i$ is not independent for each $q$, as $q_1$ and $q_2$ |
| 1666 |
< |
are linked by $\sigma$. However, $\sigma$ is fixed at any given |
| 1667 |
< |
instant of time, giving: |
| 1668 |
< |
\begin{align} |
| 1669 |
< |
\delta\sigma &= \biggl( \frac{\partial\sigma}{\partial q_1} \delta q_1 % |
| 1670 |
< |
+ \frac{\partial\sigma}{\partial q_2} \delta q_2 \biggr) = 0 ,\\ |
| 1671 |
< |
% |
| 1672 |
< |
\frac{\partial\sigma}{\partial q_1} \delta q_1 &= % |
| 1673 |
< |
- \frac{\partial\sigma}{\partial q_2} \delta q_2, \\ |
| 1674 |
< |
% |
| 1675 |
< |
\delta q_2 &= - \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
| 1676 |
< |
\frac{\partial\sigma}{\partial q_2} \biggr) \delta q_1. |
| 1677 |
< |
\end{align} |
| 1678 |
< |
Substituted back into Eq.~\ref{oopseEq:lm2}, |
| 1679 |
< |
\begin{equation} |
| 1680 |
< |
\int_{t_1}^{t_2}\biggl [ \biggl(\frac{\partial L}{\partial q_1} |
| 1681 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1682 |
< |
\biggr) |
| 1683 |
< |
- \biggl( \frac{\partial L}{\partial q_1} |
| 1684 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1685 |
< |
\biggr) \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
| 1686 |
< |
\frac{\partial\sigma}{\partial q_2} \biggr)\biggr] \delta q_1 \, dt = 0. |
| 1687 |
< |
\label{oopseEq:lm3} |
| 1688 |
< |
\end{equation} |
| 1689 |
< |
Leading to, |
| 1690 |
< |
\begin{equation} |
| 1691 |
< |
\frac{\biggl(\frac{\partial L}{\partial q_1} |
| 1692 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1693 |
< |
\biggr)}{\frac{\partial\sigma}{\partial q_1}} = |
| 1694 |
< |
\frac{\biggl(\frac{\partial L}{\partial q_2} |
| 1695 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_2} |
| 1696 |
< |
\biggr)}{\frac{\partial\sigma}{\partial q_2}}. |
| 1697 |
< |
\label{oopseEq:lm4} |
| 1698 |
< |
\end{equation} |
| 1699 |
< |
This relation can only be statisfied, if both are equal to a single |
| 1700 |
< |
function $-\lambda(t)$, |
| 1701 |
< |
\begin{align} |
| 1702 |
< |
\frac{\biggl(\frac{\partial L}{\partial q_1} |
| 1703 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1704 |
< |
\biggr)}{\frac{\partial\sigma}{\partial q_1}} &= -\lambda(t), \\ |
| 1705 |
< |
% |
| 1706 |
< |
\frac{\partial L}{\partial q_1} |
| 1707 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} &= |
| 1708 |
< |
-\lambda(t)\,\frac{\partial\sigma}{\partial q_1} ,\\ |
| 1709 |
< |
% |
| 1710 |
< |
\frac{\partial L}{\partial q_1} |
| 1711 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1712 |
< |
+ \mathcal{G}_i &= 0, |
| 1713 |
< |
\end{align} |
| 1714 |
< |
where $\mathcal{G}_i$, the force of constraint on $i$, is: |
| 1715 |
< |
\begin{equation} |
| 1716 |
< |
\mathcal{G}_i = \lambda(t)\,\frac{\partial\sigma}{\partial q_1}. |
| 1717 |
< |
\label{oopseEq:lm5} |
| 1718 |
< |
\end{equation} |
| 1719 |
< |
|
| 1720 |
< |
In a simulation, this would involve the solution of a set of $(m + n)$ |
| 1721 |
< |
number of equations. Where $m$ is the number of constraints, and $n$ |
| 1722 |
< |
is the number of constrained coordinates. In practice, this is not |
| 1723 |
< |
done, as the matrix inversion necessary to solve the system of |
| 1724 |
< |
equations would be very time consuming to solve. Additionally, the |
| 1725 |
< |
numerical error in the solution of the set of $\lambda$'s would be |
| 1726 |
< |
compounded by the error inherent in propagating by the Velocity Verlet |
| 1727 |
< |
algorithm ($\Delta t^4$). The Verlet propagation error is negligible |
| 1728 |
< |
in an unconstrained system, as one is interested in the statistics of |
| 1729 |
< |
the run, and not that the run be numerically exact to the ``true'' |
| 1730 |
< |
integration. This relates back to the ergodic hypothesis that a time |
| 1731 |
< |
integral of a valid trajectory will still give the correct ensemble |
| 1732 |
< |
average. However, in the case of constraints, if the equations of |
| 1733 |
< |
motion leave the ``true'' trajectory, they are departing from the |
| 1734 |
< |
constrained surface. The method that is used, is to iteratively solve |
| 1735 |
< |
for $\lambda(t)$ at each time step. |
| 1736 |
< |
|
| 1737 |
< |
In {\sc rattle} the equations of motion are modified subject to the |
| 1738 |
< |
following two constraints: |
| 1739 |
< |
\begin{align} |
| 1740 |
< |
\sigma_{ij}[\mathbf{r}(t)] \equiv |
| 1741 |
< |
[ \mathbf{r}_i(t) - \mathbf{r}_j(t)]^2 - d_{ij}^2 &= 0 % |
| 1742 |
< |
\label{oopseEq:c1}, \\ |
| 1743 |
< |
% |
| 1744 |
< |
[\mathbf{\dot{r}}_i(t) - \mathbf{\dot{r}}_j(t)] \cdot |
| 1745 |
< |
[\mathbf{r}_i(t) - \mathbf{r}_j(t)] &= 0 .\label{oopseEq:c2} |
| 1746 |
< |
\end{align} |
| 1747 |
< |
Eq.~\ref{oopseEq:c1} is the set of bond constraints, where $d_{ij}$ is |
| 1748 |
< |
the constrained distance between atom $i$ and |
| 1749 |
< |
$j$. Eq.~\ref{oopseEq:c2} constrains the velocities of $i$ and $j$ to |
| 1750 |
< |
be perpendicular to the bond vector, so that the bond can neither grow |
| 1751 |
< |
nor shrink. The constrained dynamics equations become: |
| 1752 |
< |
\begin{equation} |
| 1753 |
< |
m_i \mathbf{\ddot{r}}_i = \mathbf{F}_i + \mathbf{\mathcal{G}}_i, |
| 1754 |
< |
\label{oopseEq:r1} |
| 1755 |
< |
\end{equation} |
| 1756 |
< |
where,$\mathbf{\mathcal{G}}_i$ are the forces of constraint on $i$, |
| 1757 |
< |
and are defined: |
| 1758 |
< |
\begin{equation} |
| 1759 |
< |
\mathbf{\mathcal{G}}_i = - \sum_j \lambda_{ij}(t)\,\nabla \sigma_{ij}. |
| 1760 |
< |
\label{oopseEq:r2} |
| 1761 |
< |
\end{equation} |
| 1708 |
> |
\subsubsection{\label{oopseSec:zcons}Z-Constraint Method} |
| 1709 |
|
|
| 1763 |
– |
In Velocity Verlet, if $\Delta t = h$, the propagation can be written: |
| 1764 |
– |
\begin{align} |
| 1765 |
– |
\mathbf{r}_i(t+h) &= |
| 1766 |
– |
\mathbf{r}_i(t) + h\mathbf{\dot{r}}(t) + |
| 1767 |
– |
\frac{h^2}{2m_i}\,\Bigl[ \mathbf{F}_i(t) + |
| 1768 |
– |
\mathbf{\mathcal{G}}_{Ri}(t) \Bigr] \label{oopseEq:vv1}, \\ |
| 1769 |
– |
% |
| 1770 |
– |
\mathbf{\dot{r}}_i(t+h) &= |
| 1771 |
– |
\mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} |
| 1772 |
– |
\Bigl[ \mathbf{F}_i(t) + \mathbf{\mathcal{G}}_{Ri}(t) + |
| 1773 |
– |
\mathbf{F}_i(t+h) + \mathbf{\mathcal{G}}_{Vi}(t+h) \Bigr] ,% |
| 1774 |
– |
\label{oopseEq:vv2} |
| 1775 |
– |
\end{align} |
| 1776 |
– |
where: |
| 1777 |
– |
\begin{align} |
| 1778 |
– |
\mathbf{\mathcal{G}}_{Ri}(t) &= |
| 1779 |
– |
-2 \sum_j \lambda_{Rij}(t) \mathbf{r}_{ij}(t) ,\\ |
| 1780 |
– |
% |
| 1781 |
– |
\mathbf{\mathcal{G}}_{Vi}(t+h) &= |
| 1782 |
– |
-2 \sum_j \lambda_{Vij}(t+h) \mathbf{r}(t+h). |
| 1783 |
– |
\end{align} |
| 1784 |
– |
Next, define: |
| 1785 |
– |
\begin{align} |
| 1786 |
– |
g_{ij} &= h \lambda_{Rij}(t) ,\\ |
| 1787 |
– |
k_{ij} &= h \lambda_{Vij}(t+h), \\ |
| 1788 |
– |
\mathbf{q}_i &= \mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} \mathbf{F}_i(t) |
| 1789 |
– |
- \frac{1}{m_i}\sum_j g_{ij}\mathbf{r}_{ij}(t). |
| 1790 |
– |
\end{align} |
| 1791 |
– |
Using these definitions, Eq.~\ref{oopseEq:vv1} and \ref{oopseEq:vv2} |
| 1792 |
– |
can be rewritten as, |
| 1793 |
– |
\begin{align} |
| 1794 |
– |
\mathbf{r}_i(t+h) &= \mathbf{r}_i(t) + h \mathbf{q}_i ,\\ |
| 1795 |
– |
% |
| 1796 |
– |
\mathbf{\dot{r}}(t+h) &= \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h) |
| 1797 |
– |
-\frac{1}{m_i}\sum_j k_{ij} \mathbf{r}_{ij}(t+h). |
| 1798 |
– |
\end{align} |
| 1799 |
– |
|
| 1800 |
– |
To integrate the equations of motion, the {\sc rattle} algorithm first |
| 1801 |
– |
solves for $\mathbf{r}(t+h)$. Let, |
| 1802 |
– |
\begin{equation} |
| 1803 |
– |
\mathbf{q}_i = \mathbf{\dot{r}}(t) + \frac{h}{2m_i}\mathbf{F}_i(t). |
| 1804 |
– |
\end{equation} |
| 1805 |
– |
Here $\mathbf{q}_i$ corresponds to an initial unconstrained move. Next |
| 1806 |
– |
pick a constraint $j$, and let, |
| 1807 |
– |
\begin{equation} |
| 1808 |
– |
\mathbf{s} = \mathbf{r}_i(t) + h\mathbf{q}_i(t) |
| 1809 |
– |
- \mathbf{r}_j(t) + h\mathbf{q}_j(t). |
| 1810 |
– |
\label{oopseEq:ra1} |
| 1811 |
– |
\end{equation} |
| 1812 |
– |
If |
| 1813 |
– |
\begin{equation} |
| 1814 |
– |
\Big| |\mathbf{s}|^2 - d_{ij}^2 \Big| > \text{tolerance}, |
| 1815 |
– |
\end{equation} |
| 1816 |
– |
then the constraint is unsatisfied, and corrections are made to the |
| 1817 |
– |
positions. First we define a test corrected configuration as, |
| 1818 |
– |
\begin{align} |
| 1819 |
– |
\mathbf{r}_i^T(t+h) = \mathbf{r}_i(t) + h\biggl[\mathbf{q}_i - |
| 1820 |
– |
g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_i} \biggr] ,\\ |
| 1821 |
– |
% |
| 1822 |
– |
\mathbf{r}_j^T(t+h) = \mathbf{r}_j(t) + h\biggl[\mathbf{q}_j + |
| 1823 |
– |
g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_j} \biggr]. |
| 1824 |
– |
\end{align} |
| 1825 |
– |
And we chose $g_{ij}$ such that, $|\mathbf{r}_i^T - \mathbf{r}_j^T|^2 |
| 1826 |
– |
= d_{ij}^2$. Solving the quadratic for $g_{ij}$ we obtain the |
| 1827 |
– |
approximation, |
| 1828 |
– |
\begin{equation} |
| 1829 |
– |
g_{ij} = \frac{(s^2 - d^2)}{2h[\mathbf{s}\cdot\mathbf{r}_{ij}(t)] |
| 1830 |
– |
(\frac{1}{m_i} + \frac{1}{m_j})}. |
| 1831 |
– |
\end{equation} |
| 1832 |
– |
Although not an exact solution for $g_{ij}$, as this is an iterative |
| 1833 |
– |
scheme overall, the eventual solution will converge. With a trial |
| 1834 |
– |
$g_{ij}$, the new $\mathbf{q}$'s become, |
| 1835 |
– |
\begin{align} |
| 1836 |
– |
\mathbf{q}_i &= \mathbf{q}^{\text{old}}_i - g_{ij}\, |
| 1837 |
– |
\frac{\mathbf{r}_{ij}(t)}{m_i} ,\\ |
| 1838 |
– |
% |
| 1839 |
– |
\mathbf{q}_j &= \mathbf{q}^{\text{old}}_j + g_{ij}\, |
| 1840 |
– |
\frac{\mathbf{r}_{ij}(t)}{m_j} . |
| 1841 |
– |
\end{align} |
| 1842 |
– |
The whole algorithm is then repeated from Eq.~\ref{oopseEq:ra1} until |
| 1843 |
– |
all constraints are satisfied. |
| 1844 |
– |
|
| 1845 |
– |
The second step of {\sc rattle}, is to then update the velocities. The |
| 1846 |
– |
step starts with, |
| 1847 |
– |
\begin{equation} |
| 1848 |
– |
\mathbf{\dot{r}}_i(t+h) = \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h). |
| 1849 |
– |
\end{equation} |
| 1850 |
– |
Next we pick a constraint $j$, and calculate the dot product $\ell$. |
| 1851 |
– |
\begin{equation} |
| 1852 |
– |
\ell = \mathbf{r}_{ij}(t+h) \cdot \mathbf{\dot{r}}_{ij}(t+h). |
| 1853 |
– |
\label{oopseEq:rv1} |
| 1854 |
– |
\end{equation} |
| 1855 |
– |
Here if constraint Eq.~\ref{oopseEq:c2} holds, $\ell$ should be |
| 1856 |
– |
zero. Therefore if $\ell$ is greater than some tolerance, then |
| 1857 |
– |
corrections are made to the $i$ and $j$ velocities. |
| 1858 |
– |
\begin{align} |
| 1859 |
– |
\mathbf{\dot{r}}_i^T &= \mathbf{\dot{r}}_i(t+h) - k_{ij} |
| 1860 |
– |
\frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_i}, \\ |
| 1861 |
– |
% |
| 1862 |
– |
\mathbf{\dot{r}}_j^T &= \mathbf{\dot{r}}_j(t+h) + k_{ij} |
| 1863 |
– |
\frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_j}. |
| 1864 |
– |
\end{align} |
| 1865 |
– |
Like in the previous step, we select a value for $k_{ij}$ such that |
| 1866 |
– |
$\ell$ is zero. |
| 1867 |
– |
\begin{equation} |
| 1868 |
– |
k_{ij} = \frac{\ell}{d^2_{ij}(\frac{1}{m_i} + \frac{1}{m_j})}. |
| 1869 |
– |
\end{equation} |
| 1870 |
– |
The test velocities, $\mathbf{\dot{r}}^T_i$ and |
| 1871 |
– |
$\mathbf{\dot{r}}^T_j$, then replace their respective velocities, and |
| 1872 |
– |
the algorithm is iterated from Eq.~\ref{oopseEq:rv1} until all |
| 1873 |
– |
constraints are satisfied. |
| 1874 |
– |
|
| 1875 |
– |
|
| 1876 |
– |
\subsection{\label{oopseSec:zcons}Z-Constraint Method} |
| 1877 |
– |
|
| 1710 |
|
Based on the fluctuation-dissipation theorem, a force auto-correlation |
| 1711 |
|
method was developed by Roux and Karplus to investigate the dynamics |
| 1712 |
|
of ions inside ion channels.\cite{Roux91} The time-dependent friction |
| 1792 |
|
F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}= |
| 1793 |
|
-k_{\text{Harmonic}}(z(t)-z_{\text{cons}}). |
| 1794 |
|
\end{equation} |
| 1795 |
+ |
Parameters concerning the z-constraint method are summarized in Table~\ref{table:zconParams}. |
| 1796 |
|
|
| 1797 |
< |
\section{\label{oopseSec:props}Trajectory Analysis} |
| 1798 |
< |
|
| 1799 |
< |
\subsection{\label{oopseSec:staticProps}Static Property Analysis} |
| 1797 |
> |
\begin{table} |
| 1798 |
> |
\caption{The Global Keywords: Z-Constraint Parameters} |
| 1799 |
> |
\label{table:zconParams} |
| 1800 |
> |
\begin{center} |
| 1801 |
> |
% Note when adding or removing columns, the \hsize numbers must add up to the total number |
| 1802 |
> |
% of columns. |
| 1803 |
> |
\begin{tabularx}{\linewidth}% |
| 1804 |
> |
{>{\setlength{\hsize}{1.00\hsize}}X% |
| 1805 |
> |
>{\setlength{\hsize}{0.4\hsize}}X% |
| 1806 |
> |
>{\setlength{\hsize}{1.2\hsize}}X% |
| 1807 |
> |
>{\setlength{\hsize}{1.4\hsize}}X} |
| 1808 |
|
|
| 1809 |
< |
The static properties of the trajectories are analyzed with the |
| 1969 |
< |
program \texttt{staticProps}. The code is capable of calculating a |
| 1970 |
< |
number of pair correlations between species A and B. Some of which |
| 1971 |
< |
only apply to directional entities. The summary of pair correlations |
| 1972 |
< |
can be found in Table~\ref{oopseTb:gofrs} |
| 1809 |
> |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 1810 |
|
|
| 1811 |
< |
\begin{table} |
| 1812 |
< |
\caption{THE DIFFERENT PAIR CORRELATIONS IN \texttt{staticProps}} |
| 1813 |
< |
\label{oopseTb:gofrs} |
| 1814 |
< |
\begin{center} |
| 1815 |
< |
\begin{tabular}{|l|c|c|} |
| 1816 |
< |
\hline |
| 1817 |
< |
Name & Equation & Directional Atom \\ \hline |
| 1818 |
< |
$g_{\text{AB}}(r)$ & Eq.~\ref{eq:gofr} & neither \\ \hline |
| 1982 |
< |
$g_{\text{AB}}(r, \cos \theta)$ & Eq.~\ref{eq:gofrCosTheta} & A \\ \hline |
| 1983 |
< |
$g_{\text{AB}}(r, \cos \omega)$ & Eq.~\ref{eq:gofrCosOmega} & both \\ \hline |
| 1984 |
< |
$g_{\text{AB}}(x, y, z)$ & Eq.~\ref{eq:gofrXYZ} & neither \\ \hline |
| 1985 |
< |
$\langle \cos \omega \rangle_{\text{AB}}(r)$ & Eq.~\ref{eq:cosOmegaOfR} &% |
| 1986 |
< |
both \\ \hline |
| 1987 |
< |
\end{tabular} |
| 1988 |
< |
\begin{minipage}{\linewidth} |
| 1989 |
< |
\centering |
| 1990 |
< |
\vspace{2mm} |
| 1991 |
< |
The third column specifies which atom, if any, need be a directional entity. |
| 1992 |
< |
\end{minipage} |
| 1811 |
> |
{\tt zconsTime} & fs & Sets the frequency at which the {\tt .fz} file is written. & Default sets the frequency to the {\tt runTime} \\ |
| 1812 |
> |
{\tt nZconstraints} & integer & The number of zconstraint molecules& If using zconstraint method, {\tt nZconstraints} must be set \\ |
| 1813 |
> |
{\tt zconsForcePolicy} & string& The strategy of subtracting zconstraint force from unconstraint molecules & Possible strategies are BYMASS and BYNUMBER. Default strategy is set to BYMASS\\ |
| 1814 |
> |
{\tt zconsGap} & \r(A) & Set the distance between two adjacent constraint positions& Used mainly in constraining molecules sequentially \\ |
| 1815 |
> |
{\tt zconsFixtime} & fs & Sets how long the zconstraint molecule is fixed & {\tt zconsGap} must be set if {\tt zconsGap} is already set.\\ |
| 1816 |
> |
{\tt zconsUsingSMD} &logical & Flag of using Steered Molecular Dynamics or Harmonic Force to move the molecule & Using harmonic force by default\\ |
| 1817 |
> |
|
| 1818 |
> |
\end{tabularx} |
| 1819 |
|
\end{center} |
| 1820 |
|
\end{table} |
| 1821 |
|
|
| 1996 |
– |
The first pair correlation, $g_{\text{AB}}(r)$, is defined as follows: |
| 1997 |
– |
\begin{equation} |
| 1998 |
– |
g_{\text{AB}}(r) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle %% |
| 1999 |
– |
\sum_{i \in \text{A}} \sum_{j \in \text{B}} %% |
| 2000 |
– |
\delta( r - |\mathbf{r}_{ij}|) \rangle, \label{eq:gofr} |
| 2001 |
– |
\end{equation} |
| 2002 |
– |
where $\mathbf{r}_{ij}$ is the vector |
| 2003 |
– |
\begin{equation*} |
| 2004 |
– |
\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i, \notag |
| 2005 |
– |
\end{equation*} |
| 2006 |
– |
and $\frac{V}{N_{\text{A}}N_{\text{B}}}$ normalizes the average over |
| 2007 |
– |
the expected pair density at a given $r$. |
| 1822 |
|
|
| 2009 |
– |
The next two pair correlations, $g_{\text{AB}}(r, \cos \theta)$ and |
| 2010 |
– |
$g_{\text{AB}}(r, \cos \omega)$, are similar in that they are both two |
| 2011 |
– |
dimensional histograms. Both use $r$ for the primary axis then a |
| 2012 |
– |
$\cos$ for the secondary axis ($\cos \theta$ for |
| 2013 |
– |
Eq.~\ref{eq:gofrCosTheta} and $\cos \omega$ for |
| 2014 |
– |
Eq.~\ref{eq:gofrCosOmega}). This allows for the investigator to |
| 2015 |
– |
correlate alignment on directional entities. $g_{\text{AB}}(r, \cos |
| 2016 |
– |
\theta)$ is defined as follows: |
| 2017 |
– |
\begin{equation} |
| 2018 |
– |
g_{\text{AB}}(r, \cos \theta) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
| 2019 |
– |
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2020 |
– |
\delta( \cos \theta - \cos \theta_{ij}) |
| 2021 |
– |
\delta( r - |\mathbf{r}_{ij}|) \rangle. |
| 2022 |
– |
\label{eq:gofrCosTheta} |
| 2023 |
– |
\end{equation} |
| 2024 |
– |
Here |
| 2025 |
– |
\begin{equation*} |
| 2026 |
– |
\cos \theta_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{r}}_{ij}, |
| 2027 |
– |
\end{equation*} |
| 2028 |
– |
where $\mathbf{\hat{i}}$ is the unit directional vector of species $i$ |
| 2029 |
– |
and $\mathbf{\hat{r}}_{ij}$ is the unit vector associated with vector |
| 2030 |
– |
$\mathbf{r}_{ij}$. |
| 1823 |
|
|
| 1824 |
< |
The second two dimensional histogram is of the form: |
| 1824 |
> |
\section{\label{sec:minimize}Energy Minimization} |
| 1825 |
> |
|
| 1826 |
> |
|
| 1827 |
> |
As one of the basic procedures of molecular modeling, energy minimization |
| 1828 |
> |
method is used to identify configurations that are stable points on the energy |
| 1829 |
> |
surface by adjusting the atomic coordinates. Given a potential energy function |
| 1830 |
> |
$V$ which depends on a set of coordinates, energy minimization algorithm is |
| 1831 |
> |
developed to find its minimun value. Different from other packages, the |
| 1832 |
> |
coordinates in OOPSE not only include cartesian coordinates but also euler |
| 1833 |
> |
angle if directional atom or rigidbody is involved. Unfortunately, due to the |
| 1834 |
> |
number of local minima and the cost of computation, in most cases, it is |
| 1835 |
> |
always impossible to identify the global minimum. OOPSE provides two |
| 1836 |
> |
frequently used first-derivative algorithms, steepest descents and conjugate |
| 1837 |
> |
gradient, to find a reasonable local minima. |
| 1838 |
> |
|
| 1839 |
> |
Given a coordinate set $x_{k}$ and a search direction $d_{k}$, a line search |
| 1840 |
> |
algorithm is performed along $d_{k}$ to produce $x_{k+1}=x_{k}+$ $\lambda |
| 1841 |
> |
_{k}d_{k}$. |
| 1842 |
> |
|
| 1843 |
> |
In steepest descent algorithm,% |
| 1844 |
> |
|
| 1845 |
|
\begin{equation} |
| 1846 |
< |
g_{\text{AB}}(r, \cos \omega) = |
| 2035 |
< |
\frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
| 2036 |
< |
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2037 |
< |
\delta( \cos \omega - \cos \omega_{ij}) |
| 2038 |
< |
\delta( r - |\mathbf{r}_{ij}|) \rangle. \label{eq:gofrCosOmega} |
| 1846 |
> |
d_{k}=-\nabla V(x_{k}) |
| 1847 |
|
\end{equation} |
| 2040 |
– |
Here |
| 2041 |
– |
\begin{equation*} |
| 2042 |
– |
\cos \omega_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{j}}. |
| 2043 |
– |
\end{equation*} |
| 2044 |
– |
Again, $\mathbf{\hat{i}}$ and $\mathbf{\hat{j}}$ are the unit |
| 2045 |
– |
directional vectors of species $i$ and $j$. |
| 1848 |
|
|
| 2047 |
– |
The static analysis code is also cable of calculating a three |
| 2048 |
– |
dimensional pair correlation of the form: |
| 2049 |
– |
\begin{equation}\label{eq:gofrXYZ} |
| 2050 |
– |
g_{\text{AB}}(x, y, z) = |
| 2051 |
– |
\frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
| 2052 |
– |
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2053 |
– |
\delta( x - x_{ij}) |
| 2054 |
– |
\delta( y - y_{ij}) |
| 2055 |
– |
\delta( z - z_{ij}) \rangle, |
| 2056 |
– |
\end{equation} |
| 2057 |
– |
where $x_{ij}$, $y_{ij}$, and $z_{ij}$ are the $x$, $y$, and $z$ |
| 2058 |
– |
components respectively of vector $\mathbf{r}_{ij}$. |
| 1849 |
|
|
| 1850 |
< |
The final pair correlation is similar to |
| 1851 |
< |
Eq.~\ref{eq:gofrCosOmega}. $\langle \cos \omega |
| 1852 |
< |
\rangle_{\text{AB}}(r)$ is calculated in the following way: |
| 1853 |
< |
\begin{equation}\label{eq:cosOmegaOfR} |
| 2064 |
< |
\langle \cos \omega \rangle_{\text{AB}}(r) = |
| 2065 |
< |
\langle \sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2066 |
< |
(\cos \omega_{ij}) \delta( r - |\mathbf{r}_{ij}|) \rangle. |
| 2067 |
< |
\end{equation} |
| 2068 |
< |
Here $\cos \omega_{ij}$ is defined in the same way as in |
| 2069 |
< |
Eq.~\ref{eq:gofrCosOmega}. This equation is a single dimensional pair |
| 2070 |
< |
correlation that gives the average correlation of two directional |
| 2071 |
< |
entities as a function of their distance from each other. |
| 1850 |
> |
Therefore, the gradient and the direction of next step are always orthogonal |
| 1851 |
> |
which may causes oscillatory behavior in narrow valleys. To overcome this |
| 1852 |
> |
problem, the Fletcher-Reeves variant of the conjugate algorithm generates |
| 1853 |
> |
$d_{k+1}$ from the simple recursion% |
| 1854 |
|
|
| 1855 |
< |
\subsection{\label{dynamicProps}Dynamic Property Analysis} |
| 1855 |
> |
\begin{align} |
| 1856 |
> |
d_{k+1} & =-\nabla V(x_{k+1})+\gamma_{k}d_{k}\\ |
| 1857 |
> |
\gamma_{k} & =\frac{\nabla V(x_{k+1})^{T}\nabla V(x_{k+1})}{\nabla |
| 1858 |
> |
V(x_{k})^{T}\nabla V(x_{k})}% |
| 1859 |
> |
\end{align} |
| 1860 |
|
|
| 2075 |
– |
The dynamic properties of a trajectory are calculated with the program |
| 2076 |
– |
\texttt{dynamicProps}. The program calculates the following properties: |
| 2077 |
– |
\begin{gather} |
| 2078 |
– |
\langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle, \label{eq:rms}\\ |
| 2079 |
– |
\langle \mathbf{v}(t) \cdot \mathbf{v}(0) \rangle, \label{eq:velCorr} \\ |
| 2080 |
– |
\langle \mathbf{j}(t) \cdot \mathbf{j}(0) \rangle. \label{eq:angularVelCorr} |
| 2081 |
– |
\end{gather} |
| 1861 |
|
|
| 1862 |
< |
Eq.~\ref{eq:rms} is the root mean square displacement function. Which |
| 1863 |
< |
allows one to observe the average displacement of an atom as a |
| 2085 |
< |
function of time. The quantity is useful when calculating diffusion |
| 2086 |
< |
coefficients because of the Einstein Relation, which is valid at long |
| 2087 |
< |
times.\cite{allen87:csl} |
| 1862 |
> |
The Polak-Ribiere variant of conjugate gradient defines as% |
| 1863 |
> |
|
| 1864 |
|
\begin{equation} |
| 1865 |
< |
2tD = \langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle. |
| 1866 |
< |
\label{oopseEq:einstein} |
| 1865 |
> |
\gamma_{k}=\frac{[\nabla V(x_{k+1})-\nabla V(x)]^{T}\nabla V(x_{k+1})}{\nabla |
| 1866 |
> |
V(x_{k})^{T}\nabla V(x_{k})}% |
| 1867 |
|
\end{equation} |
| 1868 |
|
|
| 2093 |
– |
Eq.~\ref{eq:velCorr} and \ref{eq:angularVelCorr} are the translational |
| 2094 |
– |
velocity and angular velocity correlation functions respectively. The |
| 2095 |
– |
latter is only applicable to directional species in the |
| 2096 |
– |
simulation. The velocity autocorrelation functions are useful when |
| 2097 |
– |
determining vibrational information about the system of interest. |
| 1869 |
|
|
| 1870 |
+ |
The conjugate gradient method assumes that the conformation is close enough to |
| 1871 |
+ |
a local minimum that the potential energy surface is very nearly quadratic. |
| 1872 |
+ |
When initial structure is far from the minimimum, the steepest descent method |
| 1873 |
+ |
can be superiror to conjugate gradient. Hence, steepest descents may generally |
| 1874 |
+ |
be used for the first 10-100 steps of minimization. Another useful feature of |
| 1875 |
+ |
minimization methods in OOPSE is that a modified SHAKE algorithm can be |
| 1876 |
+ |
applied duing the minimization to constraint the bond length. {\tt bass} parameters concerning the minimizer are given in Table~\ref{table:minimizeParams} |
| 1877 |
+ |
|
| 1878 |
+ |
\begin{table} |
| 1879 |
+ |
\caption{The Global Keywords: Energy Minimizer Parameters} |
| 1880 |
+ |
\label{table:minimizeParams} |
| 1881 |
+ |
\begin{center} |
| 1882 |
+ |
% Note when adding or removing columns, the \hsize numbers must add up to the total number |
| 1883 |
+ |
% of columns. |
| 1884 |
+ |
\begin{tabularx}{\linewidth}% |
| 1885 |
+ |
{>{\setlength{\hsize}{1.00\hsize}}X% |
| 1886 |
+ |
>{\setlength{\hsize}{0.4\hsize}}X% |
| 1887 |
+ |
>{\setlength{\hsize}{1.2\hsize}}X% |
| 1888 |
+ |
>{\setlength{\hsize}{1.4\hsize}}X} |
| 1889 |
+ |
|
| 1890 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 1891 |
+ |
|
| 1892 |
+ |
{\tt minimizer} & & & \\ |
| 1893 |
+ |
{\tt minMaxIter} & integer & Sets the maximum iteration in energy minimization & Default value is 200\\ |
| 1894 |
+ |
{\tt minWriteFreq} & interger & Sets the frequency at which the {\tt .dump} and {\tt .stat} files are writtern in energy minimization & \\ |
| 1895 |
+ |
{\tt minStepSize} & double & Set the step size of line search & Default value is 0.01\\ |
| 1896 |
+ |
{\tt minFTol} & double & Sets energy tolerance & Default value is $10^(-8)$\\ |
| 1897 |
+ |
{\tt minGTol} & double & Sets gradient tolerance & Default value is $10^(-8)$\\ |
| 1898 |
+ |
{\tt minLSTol} & double & Sets line search tolerance & Default value is $10^(-8)$\\ |
| 1899 |
+ |
{\tt minLSMaxIter} & integer & Sets the maximum iteration of line searching & Default value is 50\\ |
| 1900 |
+ |
|
| 1901 |
+ |
\end{tabularx} |
| 1902 |
+ |
\end{center} |
| 1903 |
+ |
\end{table} |
| 1904 |
+ |
|
| 1905 |
+ |
|
| 1906 |
|
\section{\label{oopseSec:design}Program Design} |
| 1907 |
|
|
| 1908 |
|
\subsection{\label{sec:architecture} {\sc oopse} Architecture} |
| 2000 |
|
and favorably compete with spatial methods up to 100,000 |
| 2001 |
|
atoms.\cite{plimpton95} |
| 2002 |
|
|
| 2196 |
– |
\subsection{\label{oopseSec:memAlloc}Memory Issues in Trajectory Analysis} |
| 2197 |
– |
|
| 2198 |
– |
For large simulations, the trajectory files can sometimes reach sizes |
| 2199 |
– |
in excess of several gigabytes. In order to effectively analyze that |
| 2200 |
– |
amount of data, two memory management schemes have been devised for |
| 2201 |
– |
\texttt{staticProps} and for \texttt{dynamicProps}. The first scheme, |
| 2202 |
– |
developed for \texttt{staticProps}, is the simplest. As each frame's |
| 2203 |
– |
statistics are calculated independent of each other, memory is |
| 2204 |
– |
allocated for each frame, then freed once correlation calculations are |
| 2205 |
– |
complete for the snapshot. To prevent multiple passes through a |
| 2206 |
– |
potentially large file, \texttt{staticProps} is capable of calculating |
| 2207 |
– |
all requested correlations per frame with only a single pair loop in |
| 2208 |
– |
each frame and a single read of the file. |
| 2209 |
– |
|
| 2210 |
– |
The second, more advanced memory scheme, is used by |
| 2211 |
– |
\texttt{dynamicProps}. Here, the program must have multiple frames in |
| 2212 |
– |
memory to calculate time dependent correlations. In order to prevent a |
| 2213 |
– |
situation where the program runs out of memory due to large |
| 2214 |
– |
trajectories, the user is able to specify that the trajectory be read |
| 2215 |
– |
in blocks. The number of frames in each block is specified by the |
| 2216 |
– |
user, and upon reading a block of the trajectory, |
| 2217 |
– |
\texttt{dynamicProps} will calculate all of the time correlation frame |
| 2218 |
– |
pairs within the block. After in-block correlations are complete, a |
| 2219 |
– |
second block of the trajectory is read, and the cross correlations are |
| 2220 |
– |
calculated between the two blocks. This second block is then freed and |
| 2221 |
– |
then incremented and the process repeated until the end of the |
| 2222 |
– |
trajectory. Once the end is reached, the first block is freed then |
| 2223 |
– |
incremented, and the again the internal time correlations are |
| 2224 |
– |
calculated. The algorithm with the second block is then repeated with |
| 2225 |
– |
the new origin block, until all frame pairs have been correlated in |
| 2226 |
– |
time. This process is illustrated in |
| 2227 |
– |
Fig.~\ref{oopseFig:dynamicPropsMemory}. |
| 2228 |
– |
|
| 2229 |
– |
%\begin{figure} |
| 2230 |
– |
%\centering |
| 2231 |
– |
%\includegraphics[width=\linewidth]{dynamicPropsMem.eps} |
| 2232 |
– |
%\caption[A representation of the block correlations in \texttt{dynamicProps}]{This diagram illustrates the memory management used by \texttt{dynamicProps}, which follows the scheme: $\sum^{N_{\text{memory blocks}}}_{i=1}[ \operatorname{self}(i) + \sum^{N_{\text{memory blocks}}}_{j>i} \operatorname{cross}(i,j)]$. The shaded region represents the self correlation of the memory block, and the open blocks are read one at a time and the cross correlations between blocks are calculated.} |
| 2233 |
– |
%\label{oopseFig:dynamicPropsMemory} |
| 2234 |
– |
%\end{figure} |
| 2235 |
– |
|
| 2003 |
|
\section{\label{oopseSec:conclusion}Conclusion} |
| 2004 |
|
|
| 2005 |
|
We have presented the design and implementation of our open source |
| 2021 |
|
|
| 2022 |
|
\newpage |
| 2023 |
|
\section{Acknowledgments} |
| 2024 |
< |
The authors would like to thank espresso for fueling this work, and |
| 2258 |
< |
would also like to send a special acknowledgement to single malt |
| 2259 |
< |
scotch for its wonderful calming effects and its ability to make the |
| 2260 |
< |
troubles of the world float away. |
| 2024 |
> |
The authors would like to thank the Notre Dame BoB computer cluster where much of this project was tested. Additionally, the authors would like to acknowledge their funding from {\LARGE FIX ME}. |
| 2025 |
|
|
| 2026 |
|
\bibliographystyle{achemso} |
| 2027 |
|
\bibliography{oopsePaper} |
