ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/oopsePaper/oopsePaper.tex

Comparing trunk/oopsePaper/oopsePaper.tex (file contents):
Revision 1123 by mmeineke, Mon Apr 19 21:25:23 2004 UTC vs.
Revision 1425 by gezelter, Wed Jul 28 15:44:21 2004 UTC

#	Line 2 | Line 2
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}
#	Line 25 | Line 25 | Engine for Molecular Dynamics}
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}
#	Line 34 | Line 35 | We detail the capabilities of a new open-source parall
35		\maketitle
36
37		\begin{abstract}
38	<	We detail the capabilities of a new open-source parallel simulation
38	<	package ({\sc oopse}) that can perform molecular dynamics simulations
39	<	on atom types that are missing from other popular packages.  In
40	<	particular, {\sc oopse} is capable of performing orientational
41	<	dynamics on dipolar systems, and it can handle simulations of metallic
42	<	systems using the embedded atom method ({\sc eam}).
38	>	Need an abstract.
39		\end{abstract}
40
41		\section{\label{sec:intro}Introduction}
42
43	<	When choosing to simulate a chemical system with molecular dynamics,
48	<	there are a variety of options available. For simple systems, one
49	<	might consider writing one's own programming code. However, as systems
50	<	grow larger and more complex, building and maintaining code for the
51	<	simulations becomes a time consuming task. In such cases it is usually
52	<	more convenient for a researcher to turn to pre-existing simulation
53	<	packages. These packages, such as {\sc amber}\cite{pearlman:1995} and
54	<	{\sc charmm}\cite{Brooks83}, provide powerful tools for researchers to
55	<	conduct simulations of their systems without spending their time
56	<	developing a code base to conduct their research. This then frees them
57	<	to perhaps explore experimental analogues to their models.
43	>	UNDERWAY
44
59	–	Despite their utility, problems with these packages arise when
60	–	researchers try to develop techniques or energetic models that the
61	–	code was not originally designed to simulate. Examples of techniques
62	–	and energetics not commonly implemented include; dipole-dipole
63	–	interactions, rigid body dynamics, and metallic potentials. When faced
64	–	with these obstacles, a researcher must either develop their own code
65	–	or license and extend one of the commercial packages. What we have
66	–	elected to do is develop a body of simulation code capable of
67	–	implementing the types of models upon which our research is based.
45
46	<	In developing {\sc oopse}, we have adhered to the precepts of Open
47	<	Source development, and are releasing our source code with a
48	<	permissive license. It is our intent that by doing so, other
49	<	researchers might benefit from our work, and add their own
50	<	contributions to the package. The license under which {\sc oopse} is
51	<	distributed allows any researcher to download and modify the source
52	<	code for their own use. In this way further development of {\sc oopse}
53	<	is not limited to only the models of interest to ourselves, but also
54	<	those of the community of scientists who contribute back to the
78	<	project.
46	>	We have structured this paper to first discuss the underlying concepts
47	>	in this simulation package (Sec. \ref{oopseSec:IOfiles}).  The
48	>	empirical energy functions implemented are discussed in
49	>	Sec.~\ref{oopseSec:empiricalEnergy}.  Sec.~\ref{oopseSec:mechanics}
50	>	describes the various Molecular Dynamics algorithms {\sc oopse}
51	>	implements in the integration of the Newtonian equations of motion.
52	>	Program design considerations for parallel computing are presented in
53	>	Sec.~\ref{oopseSec:parallelization}. Concluding remarks are presented
54	>	in Sec.~\ref{oopseSec:conclusion}.
55
56	<	We have structured this chapter 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
89	<	considerations are presented in Sec.~\ref{oopseSec:design}. And
90	<	lastly, Sec.~\ref{oopseSec:conclusion} concludes the chapter.
56	>	\section{\label{oopseSec:IOfiles}Concepts \& Files}
57
58	<	\section{\label{oopseSec:empiricalEnergy}The Empirical Energy Functions}
58	>	A simulation in {\sc oopse} is built using a few fundamental
59	>	conceptual building blocks most of which are chemically intuitive.
60	>	The basic unit of a simulation is an {\tt atom}.  The parameters
61	>	describing an {\tt atom} have been generalized to make it as flexible
62	>	as possible; this means that in addition to translational degrees of
63	>	freedom, {\tt Atoms} may also have {\it orientational} degrees of freedom.
64
65	<	\subsection{\label{oopseSec:atomsMolecules}Atoms, Molecules and Rigid Bodies}
65	>	The fundamental (static) properties of {\tt atoms} are defined by the
66	>	{\tt forceField} chosen for the simulation.  The atomic properties
67	>	specified by a {\tt forceField} might include (but are not limited to)
68	>	charge, $\sigma$ and $\epsilon$ values for Lennard-Jones interactions,
69	>	the strength of the dipole moment ($\mu$), the mass, and the moments
70	>	of inertia.  Other more complicated properties of atoms might also be
71	>	specified by the {\tt forceField}.
72
73	<	The basic unit of an {\sc oopse} simulation is the atom. The
74	<	parameters describing the atom are generalized to make the atom as
75	<	flexible a representation as possible. They may represent specific
76	<	atoms of an element, or be used for collections of atoms such as
77	<	methyl and carbonyl groups. The atoms are also capable of having
78	<	directional components associated with them (\emph{e.g.}~permanent
79	<	dipoles). Charges, permanent dipoles, and Lennard-Jones parameters for
80	<	a given atom type are set in the force field parameter files.
73	>	{\tt Atoms} can be grouped together in many ways.  A {\tt rigidBody}
74	>	contains atoms that exert no forces on one another and which move as a
75	>	single rigid unit.  A {\tt cutoffGroup} may contain atoms which
76	>	function together as a (rigid {\it or} non-rigid) unit for potential
77	>	energy calculations,
78	>	\begin{equation}
79	>	V_{ab} = s(r_{ab}) \sum_{i \in a} \sum_{j \in b} V_{ij}(r_{ij})
80	>	\end{equation}
81	>	Here, $a$ and $b$ are two {\tt cutoffGroups} containing multiple atoms
82	>	($a = \left\{i\right\}$ and $b = \left\{j\right\}$).  $s(r_{ab})$ is a
83	>	generalized switching function which insures that the atoms in the two
84	>	{\tt cutoffGroups} are treated identically as the two groups enter or
85	>	leave an interaction region.
86
87	<	\begin{lstlisting}[float,caption={[Specifier for molecules and atoms] A sample specification of an Ar molecule},label=sch:AtmMole]
87	>	{\tt Atoms} may also be grouped in more traditional ways into {\tt
88	>	bonds}, {\tt bends}, and {\tt torsions}.  These groupings allow the
89	>	correct choice of interaction parameters for short-range interactions
90	>	to be chosen from the definitions in the {\tt forceField}.
91	>
92	>	All of these groups of {\tt atoms} are brought together in the {\tt
93	>	molecule}, which is the fundamental structure for setting up and {\sc
94	>	oopse} simulation.  {\tt Molecules} contain lists of {\tt atoms}
95	>	followed by listings of the other atomic groupings ({\tt bonds}, {\tt
96	>	bends}, {\tt torsions}, {\tt rigidBodies}, and {\tt cutoffGroups})
97	>	which relate the atoms to one another.
98	>
99	>	Simulations often involve heterogeneous collections of molecules.  To
100	>	specify a mixture of {\tt molecule} types, {\sc oopse} uses {\tt
101	>	components}.  Even simulations containing only one type of molecule
102	>	must specify a single {\tt component}.
103	>
104	>	Starting a simulation requires two types of information: {\it
105	>	meta-data}, which describes the types of objects present in the
106	>	simulation, and {\it configuration} information, which describes the
107	>	initial state of these objects.  The meta-data is given to {\sc oopse}
108	>	using a C-based syntax that is parsed at the beginning of the
109	>	simulation.  Configuration information is specified using an extended
110	>	XYZ file format.  Both the meta-data and configuration file formats
111	>	are described in the following sections.
112	>
113	>	\subsection{Meta-data Files}
114	>
115	>	{\sc oopse} uses a C-based script syntax to parse the meta-data files
116	>	at run time.  These files allow the user to completely describe the
117	>	system they wish to simulate, as well as tailor {\sc oopse}'s behavior
118	>	during the simulation.  Meta-data files are typically denoted with the
119	>	extension {\tt .md} (which can stand for Meta-Data or Molecular
120	>	Dynamics or Molecule Definition depending on the user's mood). An
121	>	example meta-data file is shown in Scheme~\ref{sch:mdExample}.
122	>
123	>	\begin{lstlisting}[float,caption={[An example of a complete meta-data
124	>	file] An example showing a complete meta-data
125	>	file.},label={sch:mdExample}]
126	>
127	>	molecule{
128	>	name = "Ar";
129	>	nAtoms = 1;
130	>	atom[0]{
131	>	type="Ar";
132	>	position( 0.0, 0.0, 0.0 );
133	>	}
134	>	}
135	>
136	>	nComponents = 1;
137	>	component{
138	>	type = "Ar";
139	>	nMol = 108;
140	>	}
141	>
142	>	initialConfig = "./argon.in";
143	>
144	>	forceField = "LJ";
145	>	ensemble = "NVE"; // specify the simulation ensemble
146	>	dt = 1.0;         // the time step for integration
147	>	runTime = 1e3;    // the total simulation run time
148	>	sampleTime = 100; // trajectory file frequency
149	>	statusTime = 50;  // statistics file frequency
150	>
151	>	\end{lstlisting}
152	>
153	>	Within the meta-data file it is necessary to provide a complete
154	>	description of the molecule before it is actually placed in the
155	>	simulation. {\sc oopse}'s meta-data syntax was originally developed
156	>	with this goal in mind, and allows for the use of {\it include files}
157	>	to specify all atoms in a molecular prototype, as well as any bonds,
158	>	bends, or torsions.  Include files allow the user to describe a
159	>	molecular prototype once, then simply include it into each simulation
160	>	containing that molecule. Returning to the example in
161	>	Scheme~\ref{sch:mdExample}, the include file's contents would be
162	>	Scheme~\ref{sch:mdIncludeExample}, and the new meta-data file would
163	>	become Scheme~\ref{sch:mdExPrime}.
164	>
165	>	\begin{lstlisting}[float,caption={An example molecule definition in an
166	>	include file.},label={sch:mdIncludeExample}]
167	>
168		molecule{
169		name = "Ar";
170		nAtoms = 1;
#	Line 111 | Line 173 | molecule{
173		position( 0.0, 0.0, 0.0 );
174		}
175		}
176	+
177		\end{lstlisting}
178
179	+	\begin{lstlisting}[float,caption={Revised meta-data example.},label={sch:mdExPrime}]
180
181	<	Atoms can be collected into secondary structures such as rigid bodies
118	<	or molecules. The molecule is a way for {\sc oopse} to keep track of
119	<	the atoms in a simulation in logical manner. Molecular units store the
120	<	identities of all the atoms and rigid bodies associated with
121	<	themselves, and are responsible for the evaluation of their own
122	<	internal interactions (\emph{i.e.}~bonds, bends, and torsions). Scheme
123	<	\ref{sch:AtmMole} shows how one creates a molecule in a ``model'' or
124	<	\texttt{.mdl} file. The position of the atoms given in the
125	<	declaration are relative to the origin of the molecule, and is used
126	<	when creating a system containing the molecule.
181	>	#include "argon.md"
182
183	<	As stated previously, one of the features that sets {\sc oopse} apart
184	<	from most of the current molecular simulation packages is the ability
185	<	to handle rigid body dynamics. Rigid bodies are non-spherical
186	<	particles or collections of particles that have a constant internal
183	>	nComponents = 1;
184	>	component{
185	>	type = "Ar";
186	>	nMol = 108;
187	>	}
188	>
189	>	initialConfig = "./argon.in";
190	>
191	>	forceField = "LJ";
192	>	ensemble = "NVE";
193	>	dt = 1.0;
194	>	runTime = 1e3;
195	>	sampleTime = 100;
196	>	statusTime = 50;
197	>
198	>	\end{lstlisting}
199	>
200	>	\subsection{\label{oopseSec:atomsMolecules}Atoms, Molecules, and other
201	>	ways of grouping atoms}
202	>
203	>	As mentioned above, the fundamental unit for an {\sc oopse} simulation
204	>	is the {\tt atom}.  Atoms can be collected into secondary structures
205	>	such as {\tt rigidBodies}, {\tt cutoffGroups}, or {\tt molecules}. The
206	>	{\tt molecule} is a way for {\sc oopse} to keep track of the atoms in
207	>	a simulation in logical manner. Molecular units store the identities
208	>	of all the atoms and rigid bodies associated with themselves, and they
209	>	are responsible for the evaluation of their own internal interactions
210	>	(\emph{i.e.}~bonds, bends, and torsions). Scheme
211	>	\ref{sch:mdIncludeExample} shows how one creates a molecule in an
212	>	included meta-data file. The positions of the atoms given in the
213	>	declaration are relative to the origin of the molecule, and the origin
214	>	is used when creating a system containing the molecule.
215	>
216	>	One of the features that sets {\sc oopse} apart from most of the
217	>	current molecular simulation packages is the ability to handle rigid
218	>	body dynamics. Rigid bodies are non-spherical particles or collections
219	>	of particles (e.g. $\mbox{C}_{60}$) that have a constant internal
220		potential and move collectively.\cite{Goldstein01} They are not
221		included in most simulation packages because of the algorithmic
222	<	complexity involved in propagating orientational degrees of
223	<	freedom. Until recently, integrators which propagate orientational
224	<	motion have been much worse than those available for translational
225	<	motion.
222	>	complexity involved in propagating orientational degrees of freedom.
223	>	Integrators which propagate orientational motion with an acceptable
224	>	level of energy conservation for molecular dynamics are relatively
225	>	new inventions.
226
227		Moving a rigid body involves determination of both the force and
228		torque applied by the surroundings, which directly affect the
#	Line 144 | Line 232 | the rigid body. The torque on rigid body $i$ is
232		the rigid body is simply the sum of these external forces.
233		Accumulation of the total torque on the rigid body is more complex
234		than the force because the torque is applied to the center of mass of
235	<	the rigid body. The torque on rigid body $i$ is
235	>	the rigid body. The space-fixed torque on rigid body $i$ is
236		\begin{equation}
237		\boldsymbol{\tau}_i=
238		\sum_{a}\biggl[(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia}
#	Line 171 | Line 259 | systems.\cite{Evans77}
259		performance enhancements, particularly for very small
260		systems.\cite{Evans77}
261
262	<	{\sc oopse} utilizes a relatively new scheme that propagates the
263	<	entire nine parameter rotation matrix. Further discussion
264	<	on this choice can be found in Sec.~\ref{oopseSec:integrate}. An example
265	<	definition of a rigid body can be seen in Scheme
266	<	\ref{sch:rigidBody}. The positions in the atom definitions are the
267	<	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.
262	>	Rather than use one of the previously stated methods, {\sc oopse}
263	>	utilizes a relatively new scheme that propagates the entire nine
264	>	parameter rotation matrix. Further discussion on this choice can be
265	>	found in Sec.~\ref{oopseSec:integrate}. An example definition of a
266	>	rigid body can be seen in Scheme
267	>	\ref{sch:rigidBody}.
268
269	<	\begin{lstlisting}[float,caption={[Defining rigid bodies]A sample definition of a rigid body},label={sch:rigidBody}]
269	>	\begin{lstlisting}[float,caption={[Defining rigid bodies]A sample
270	>	definition of a molecule containing a rigid body and a cutoff
271	>	group},label={sch:rigidBody}]
272		molecule{
273		name = "TIP3P";
274		nAtoms = 3;
#	Line 201 | Line 290 | molecule{
290		nMembers = 3;
291		members(0, 1, 2);
292		}
293	+
294	+	nCutoffGroups = 1;
295	+	cutoffGroup[0]{
296	+	nMembers = 3;
297	+	members(0, 1, 2);
298	+	}
299		}
300		\end{lstlisting}
301	+
302	+	\subsection{\label{sec:miscConcepts}Creating a Metadata File}
303	+
304	+	The actual creation of a metadata file requires several key
305	+	components. The first part of the file needs to be the declaration of
306	+	all of the molecule prototypes used in the simulation. This is
307	+	typically done through included meta-data files. Only the molecules
308	+	actually present in the simulation need to be declared; however, {\sc
309	+	oopse} allows for the declaration of more molecules than are
310	+	needed. This gives the user the ability to build up a library of
311	+	commonly used molecules into a single include file.
312	+
313	+	Once all prototypes are declared, the ordering of the rest of the
314	+	script is less stringent.  The molecular composition of the simulation
315	+	is specified with {\tt component} statements. Each different type of
316	+	molecule present in the simulation is considered a separate
317	+	component. The number of components must be declared before the first
318	+	component block statement (an example is shown in
319	+	Sch.~\ref{sch:mdExPrime}).  The component blocks tell {\sc oopse} the
320	+	number of molecules that will be in the simulation, and the order in
321	+	which the components blocks are declared sets the ordering of the real
322	+	atoms in the configuration file as well as in the output files. The
323	+	remainder of the script then sets the various simulation parameters
324	+	for the system of interest.
325	+
326	+	The required set of parameters that must be present in all simulations
327	+	is given in Table~\ref{table:reqParams}.  Since the user can use {\sc
328	+	oopse} to perform energy minimizations as well as molecular dynamics
329	+	simulations, one of the {\tt minimizer} or {\tt ensemble} keywords
330	+	must be present.  The {\tt ensemble} keyword is responsible for
331	+	selecting the integration method used for the calculation of the
332	+	equations of motion. An in depth discussion of the various methods
333	+	available in {\sc oopse} can be found in
334	+	Sec.~\ref{oopseSec:mechanics}.  The {\tt minimizer} keyword selects
335	+	which minimization method to use, and more details on the choices of
336	+	minimizer parameters can be found in
337	+	Sec.~\ref{oopseSec:minimizer}. The {\tt forceField} statement is
338	+	important for the selection of which forces will be used in the course
339	+	of the simulation. {\sc oopse} supports several force fields, as
340	+	outlined in Sec.~\ref{oopseSec:empericalEnergy}. The force fields are
341	+	interchangeable between simulations, with the only requirement being
342	+	that all atoms needed by the simulation are defined within the
343	+	selected force field.
344	+
345	+	For molecular dynamics simulations, the time step between force
346	+	evaluations is set with the {\tt dt} parameter, and {\tt runTime} will
347	+	set the time length of the simulation. Note, that {\tt runTime} is an
348	+	absolute time, meaning if the simulation is started at t = 10.0~ns
349	+	with a {\tt runTime} of 25.0~ns, the simulation will only run for an
350	+	additional 15.0~ns.
351	+
352	+	For energy minimizations, it is not necessary to specify {\tt dt} or
353	+	{\tt runTime}.
354	+
355	+	The final required parameter is the {\tt initialConfig}
356	+	statement. This will set the initial coordinates for the system, as
357	+	well as the initial time if the {\tt useInitalTime} flag is set to
358	+	{\tt true}. The format of the file specified in {\tt initialConfig},
359	+	is given in Sec.~\ref{oopseSec:coordFiles}. Additional parameters are
360	+	summarized in Table~\ref{table:genParams}.
361	+
362	+	It is important to note the fundamental units in all files which are
363	+	read and written by {\sc oopse}.  Energies are in $\mbox{kcal
364	+	mol}^{-1}$, distances are in $\mbox{\AA}$, times are in $\mbox{fs}$,
365	+	translational velocities are in $\mbox{\AA fs}^{-1}$, and masses are
366	+	in $\mbox{amu}$.  Orientational degrees of freedom are described using
367	+	quaternions (unitless, but $q_w^2 + q_x^2 + q_y^2 + q_z^2 = 1$),
368	+	body-fixed angular momenta ($\mbox{amu \AA}^{2} \mbox{radians
369	+	fs}^{-1}$), and body-fixed moments of inertia ($\mbox{amu \AA}^{2}$).
370	+
371	+	\begin{table}
372	+	\caption{Meta-data Keywords: Required Parameters}
373	+	\label{table:reqParams}
374	+	\begin{center}
375	+	% Note when adding or removing columns, the \hsize numbers must add up to the total number
376	+	% of columns.
377	+	\begin{tabularx}{\linewidth}%
378	+	{>{\setlength{\hsize}{1.00\hsize}}X%
379	+	>{\setlength{\hsize}{0.4\hsize}}X%
380	+	>{\setlength{\hsize}{1.2\hsize}}X%
381	+	>{\setlength{\hsize}{1.4\hsize}}X}
382	+
383	+	{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline
384	+
385	+	{\tt forceField} & string & Sets the force field. & Possible force fields are "DUFF", "LJ", and "EAM". \\
386	+	{\tt nComponents} & integer & Sets the number of components. & Needs to appear before the first {\tt Component} block. \\
387	+	{\tt initialConfig} & string & Sets the file containing the initial configuration. & Can point to any file containing the configuration in the correct order. \\
388	+	{\tt minimizer}& string & Chooses a minimizer & Possible minimizers
389	+	are "SD" and "CG". Either {\tt ensemble} or {\tt minimizer} must be specified. \\
390	+	{\tt ensemble} & string & Sets the ensemble. & Possible ensembles are
391	+	"NVE", "NVT", "NPTi", "NPTf", and "NPTxyz".  Either {\tt ensemble}
392	+	or {\tt minimizer} must be specified. \\
393	+	{\tt dt} & fs & Sets the time step. & Selection of {\tt dt} should be
394	+	small enough to sample the fastest motion of the simulation. (required
395	+	for molecular dynamics simulations)\\
396	+	{\tt runTime} & fs & Sets the time at which the simulation should
397	+	end. & This is an absolute time, and will end the simulation when the
398	+	current time meets or exceeds the {\tt runTime}. (required for
399	+	molecular dynamics simulations)\\
400	+
401	+	\end{tabularx}
402	+	\end{center}
403	+	\end{table}
404	+
405	+	\begin{table}
406	+	\caption{Meta-data Keywords: General Parameters}
407	+	\label{table:genParams}
408	+	\begin{center}
409	+	% Note when adding or removing columns, the \hsize numbers must add up to the total number
410	+	% of columns.
411	+	\begin{tabularx}{\linewidth}%
412	+	{>{\setlength{\hsize}{1.00\hsize}}X%
413	+	>{\setlength{\hsize}{0.4\hsize}}X%
414	+	>{\setlength{\hsize}{1.2\hsize}}X%
415	+	>{\setlength{\hsize}{1.4\hsize}}X}
416	+
417	+	{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline
418	+
419	+	{\tt finalConfig} & string & Sets the name of the final
420	+	output file. & Useful when stringing simulations together. Defaults to
421	+	the root name of the initial meta-data file but with an {\tt .eor}
422	+	extension. \\
423	+	{\tt useInitialTime} & logical & Sets whether the initial time is taken from the {\tt .in} file. & Useful when recovering a simulation from a crashed processor. Default is false. \\
424	+	{\tt sampleTime} & fs & Sets the frequency at which the {\tt .dump} file is written. & Default sets the frequency to the {\tt runTime}. \\
425	+	{\tt statusTime} & fs & Sets the frequency at which the {\tt .stat} file is written. & Defaults set the frequency to the {\tt sampleTime}. \\
426	+	{\tt cutoffRadius} & $\mbox{\AA}$ & Manually sets the cutoffRadius & Defaults to
427	+	$15\mbox{\AA}$ for systems containing charges or dipoles or to $2.5
428	+	\sigma_{L}$, where $\sigma_{L}$ is the largest LJ $\sigma$ in the
429	+	simulation. \\
430	+	{\tt switchingRadius} & $\mbox{\AA}$  & Manually sets the inner radius for the switching function. & Defaults to 95~\% of the {\tt cutoffRadius}. \\
431	+	{\tt useReactionField} & logical & Turns the reaction field correction on/off. & Default is "false". \\
432	+	{\tt dielectric} & unitless & Sets the dielectric constant for reaction field. & If {\tt useReactionField} is true, then {\tt dielectric} must be set. \\
433	+	{\tt usePeriodicBoundaryConditions} & & & \\
434	+	& logical & Turns periodic boundary conditions on/off. & Default is "true". \\
435	+	{\tt seed } & integer & Sets the seed value for the random number
436	+	generator. & The seed needs to be at least 9 digits long. The default
437	+	is to take the seed from the CPU clock. \\
438	+	{\tt forceFieldVariant} & string & Sets the name of the variant of the
439	+	force field.  ({\sc eam} has three variants: {\tt u3}, {\tt u6}, and
440	+	{\tt VC}.
441	+
442	+	\end{tabularx}
443	+	\end{center}
444	+	\end{table}
445	+
446	+
447	+	\subsection{\label{oopseSec:coordFiles}Coordinate Files}
448	+
449	+	The standard format for storage of a systems coordinates is a modified
450	+	xyz-file syntax, the exact details of which can be seen in
451	+	Scheme~\ref{sch:dumpFormat}. As all bonding and molecular information
452	+	is stored in the meta-data files, the coordinate files contain only
453	+	the coordinates of the objects which move independently during the
454	+	simulation.  It is important to note that {\it not all atoms} are
455	+	capable of independent motion.  Atoms which are part of rigid bodies
456	+	are not ``integrable objects'' in the equations of motion; the rigid
457	+	bodies themselves are the integrable objects.  Therefore, the
458	+	coordinate file contains coordinates of all the {\tt
459	+	integrableObjects} in the system.  For systems without rigid bodies,
460	+	this is simply the coordinates of all the atoms.
461	+
462	+	It is important to note that although the simulation propagates the
463	+	complete rotation matrix, directional entities are written out using
464	+	quaternions to save space in the output files.  All objects (atoms,
465	+	orientational atoms, and rigid bodies) are given quaternions and
466	+	angular momenta in coordinate files which are output by OOPSE, but it
467	+	is not necessary for the user to specify the quaternions or angular
468	+	momenta for atoms without orientational degrees of freedom.
469	+
470	+	\begin{lstlisting}[float,caption={[The format of the coordinate
471	+	files] An example of the format of the coordinate files. The fist line
472	+	is the number of {\tt integrableObjects} (freely-moving atoms and
473	+	rigid bodies). The second line begins with the time stamp followed by
474	+	the three $\mathsf{H}$ column vectors. It is important to note that
475	+	for extended system ensembles, additional information pertinent to the
476	+	integrators may be stored on this line as well. The next lines are the
477	+	coordinates for all integrable objects in the system.  The name of the
478	+	integrable object is followed by position, velocity, quaternions, and
479	+	lastly, body fixed angular momentum.},label=sch:dumpFormat]
480	+
481	+	nIntegrable
482	+	time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz;
483	+	Name1 x y z vx vy vz qw qx qy qz jx jy jz
484	+	Name2 x y z vx vy vz qw qx qy qz jx jy jz
485	+	etc...
486	+
487	+	\end{lstlisting}
488	+
489	+	The {\tt name} field for atoms is simply the atom type as specified in
490	+	the meta-data file.  The {\tt name} field for a rigid body is
491	+	specified as {\tt MOLTYPE\_RB\_N}, to specify that this is {\tt
492	+	rigidBody} N in a molecule of type MOLTYPE.  In simulations with rigid
493	+	body models of water, a sample coordinate line might be:
494	+
495	+	\begin{tt}
496	+	TIP3P\_RB\_0  x y z vx vy vz qw qx qy qz jx jy jz
497	+	\end{tt}
498	+
499	+	which tells the program that the rigid body representing a TIP3P
500	+	molecule (rigid body \# 0) is listed on that line.
501	+
502	+	There are three files used by {\sc oopse} which are written in the
503	+	coordinate format.  They are: the initial coordinate file
504	+	(\texttt{.in}), the simulation trajectory file (\texttt{.dump}), and
505	+	the final coordinates or ``end-of-run'' for the simulation
506	+	(\texttt{.eor}). The initial coordinate file is necessary for {\sc
507	+	oopse} to start the simulation with the proper coordinates, and this
508	+	file must be generated by the user before the simulation run. The
509	+	trajectory (or ``dump'') file is updated during simulation and is used
510	+	to store snapshots of the coordinates at regular intervals. The first
511	+	frame is a duplication of the
512	+	\texttt{.in} file, and each subsequent frame is appended to the file
513	+	at an interval specified in the meta-data file with the
514	+	\texttt{sampleTime} flag. The final coordinate file is the
515	+	``end-of-run'' file.  The \texttt{.eor} file stores the final
516	+	configuration of the system for a given simulation. The file is
517	+	updated at the same time as the \texttt{.dump} file, but it only
518	+	contains the most recent frame. In this way, an \texttt{.eor} file may
519	+	be used to initialize a second simulation should it be necessary to
520	+	recover from a crash or power outage.
521	+
522	+	\subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates}
523	+
524	+	As was stated in Sec.~\ref{oopseSec:coordFiles}, an initial coordinate
525	+	file is needed to provide the starting coordinates for a simulation.
526	+	Since each simulation is different, system creation is left to the end
527	+	user; however, we have included a few sample programs which make some
528	+	specialized structures.  The {\tt .in} file must list the integrable
529	+	objects in the correct order.  The ordering of the integrable objects
530	+	relies on the ordering of molecules within the meta-data file. {\sc
531	+	oopse} expects the order to comply with the following guidelines:
532	+	\begin{enumerate}
533	+	\item All of the molecules of the first declared component are given
534	+	before proceeding to the molecules of the second component, and so on
535	+	for all subsequently declared components.
536	+	\item The ordering of the atoms for each molecule follows the order
537	+	declared in the molecule's declaration within the model file.
538	+	\item Only atoms which are not members of a {\tt rigidBody} are
539	+	included
540	+	\item Rigid Body coordinates for a molecule are listed immediately
541	+	after the the other atoms in a molecule.  Some molecules may be
542	+	entirely rigid, in which case, only the rigid body coordinates are
543	+	given.
544	+	\end{enumerate}
545	+	An example is given in the meta-data file in Scheme~\ref{sch:initEx1}
546	+	which results in the {\tt .in} file shown in Scheme~\ref{sch:initEx2}.
547	+
548	+	\begin{lstlisting}[float,caption={Example declaration of the
549	+	$\text{I}_2$ molecule and the HCl molecule. The two molecules are then
550	+	included into a simulation.}, label=sch:initEx1]
551	+
552	+	molecule{
553	+	name = "I2";
554	+	nAtoms = 2;
555	+	atom[0]{
556	+	type = "I";
557	+	}
558	+	atom[1]{
559	+	type = "I";
560	+	}
561	+	nBonds = 1;
562	+	bond[0]{
563	+	members( 0, 1);
564	+	}
565	+	}
566	+
567	+	molecule{
568	+	name = "HCl"
569	+	nAtoms = 2;
570	+	atom[0]{
571	+	type = "H";
572	+	}
573	+	atom[1]{
574	+	type = "Cl";
575	+	}
576	+	nBonds = 1;
577	+	bond[0]{
578	+	members( 0, 1);
579	+	}
580	+	}
581	+
582	+	nComponents = 2;
583	+	component{
584	+	type = "HCl";
585	+	nMol = 4;
586	+	}
587	+	component{
588	+	type = "I2";
589	+	nMol = 1;
590	+	}
591	+
592	+	initialConfig = "mixture.in";
593	+
594	+	\end{lstlisting}
595	+
596	+	\begin{lstlisting}[float,caption={The contents of the {\tt
597	+	mixture.in} file matching the declarations in
598	+	Scheme~\ref{sch:initEx1}. Note that even though $\text{I}_2$ is
599	+	declared before HCl, the {\tt .in} file follows the order {\it in
600	+	which the components were included}.},label=sch:initEx2]
601	+
602	+	10
603	+	0.0;  10.0  0.0  0.0;  0.0  10.0  0.0;  0.0  0.0  10.0;
604	+	H  ...
605	+	Cl ...
606	+	H  ...
607	+	Cl ...
608	+	H  ...
609	+	Cl ...
610	+	H  ...
611	+	Cl ...
612	+	I  ...
613	+	I  ...
614	+
615	+	\end{lstlisting}
616	+
617	+
618	+	\subsection{The Statistics File}
619	+
620	+	The last output file generated by {\sc oopse} is the statistics
621	+	file. This file records such statistical quantities as the
622	+	instantaneous temperature (in $K$), volume (in $\mbox{\AA}^{3}$),
623	+	pressure (in $\mbox{atm}$), etc. It is written out with the frequency
624	+	specified in the meta-data file with the
625	+	\texttt{statusTime} keyword. The file allows the user to observe the
626	+	system variables as a function of simulation time while the simulation
627	+	is in progress. One useful function the statistics file serves is to
628	+	monitor the conserved quantity of a given simulation ensemble,
629	+	allowing the user to gauge the stability of the integrator. The
630	+	statistics file is denoted with the \texttt{.stat} file extension.
631	+
632	+	\section{\label{oopseSec:empiricalEnergy}The Empirical Energy
633	+	Functions}
634	+
635	+	Like many simulation packages, {\sc oopse} splits the potential energy
636	+	into the short-ranged (bonded) portion and a long-range (non-bonded)
637	+	potential,
638	+	\begin{equation}
639	+	V = V_{\mathrm{short-range}} + V_{\mathrm{long-range}}.
640	+	\end{equation}
641	+	The short-ranged portion includes explicit bonds, bends and torsions,
642	+	which have been defined in the meta-data file for the molecules which
643	+	present in the simulation.  The functional forms and parameters for
644	+	these interactions are defined by the force field which is chosen.
645
646	+	Calculating long-range (non-bonded) potential involves a sum over all
647	+	pairs of atoms (except for those atoms which are involved in a bond,
648	+	bend, or torsion with each other).  If done poorly, calculating the
649	+	the long-range interactions for $N$ atoms would involve $N^2$
650	+	evaluations of atomic distance.  To reduce the number of distance
651	+	evaluations between pairs of atoms, {\sc oopse} uses a switched cutoff
652	+	with Verlet neighbor lists.\cite{allen87:csl} It is well known that
653	+	neutral groups which contain charges will exhibit pathological forces
654	+	unless the cutoff is applied to the neutral groups evenly instead of
655	+	to the individual atoms.\cite{leach01:mm} {\sc oopse} allows users to
656	+	specify cutoff groups which may contain an arbitrary number of atoms
657	+	in the molecule.  Atoms in a cutoff group are treated as a single unit
658	+	for the evaluation of the switching function:
659	+	\begin{equation}
660	+	V_{\mathrm{long-range}} = \sum_{a} \sum_{b>a} s(r_{ab}) \sum_{i \in a} \sum_{j \in b} V_{ij}(r_{ij}),
661	+	\end{equation}
662	+	where $r_{ab}$ is the distance between the centers of mass of the two
663	+	cutoff groups ($a$ and $b$).
664	+
665	+	The sums over $a$ and $b$ are over the cutoffGroups that are present
666	+	in the simulation.  Atoms which are not explicitly defined as members
667	+	of a {\tt cutoffGroup} are treated as a group consisting of only one
668	+	atom.  The switching function, $s(r)$ is the standard cubic switching
669	+	function,
670	+	\begin{equation}
671	+	S(r) =
672	+	\begin{cases}
673	+	1 & \text{if $r \le r_{\text{sw}}$},\\
674	+	\frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2}
675	+	{(r_{\text{cut}} - r_{\text{sw}})^2}
676	+	& \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\
677	+	0 & \text{if $r > r_{\text{cut}}$.}
678	+	\end{cases}
679	+	\label{eq:dipoleSwitching}
680	+	\end{equation}
681	+	Here, $r_{\text{sw}}$ is the {\tt switchingRadius}, or the distance
682	+	beyond which interactions are reduced, and $r_{\text{cut}}$ is the
683	+	{\tt cutoffRadius}, or the distance at which interactions are
684	+	truncated.
685	+
686	+	Users of {\sc oopse} do not need to specify the {\tt cutoffRadius} or
687	+	{\tt switchingRadius}.  In simulations containing only Lennard-Jones
688	+	atoms, the cutoff radius has a default value of $2.5\sigma_{ii}$,
689	+	where $\sigma_{ii}$ is the largest Lennard-Jones length parameter
690	+	present in the simulation.  In simulations containing charged or
691	+	dipolar atoms, the default cutoff Radius is $15 \mbox{\AA}$.
692	+
693	+	The {\tt switchingRadius} is set to a default value of 95\% of the
694	+	{\tt cutoffRadius}.  In the special case of a simulation containing
695	+	{\it only} Lennard-Jones atoms, the default switching radius takes the
696	+	same value as the cutoff radius, and {\sc oopse} will use a shifted
697	+	potential to remove discontinuities in the potential at the cutoff.
698	+	Both radii may be specified in the meta-data file.
699	+
700	+	Force fields can easily be added to {\sc oopse}, although it comes
701	+	with a few simple examples (Lennard-Jones, {\sc duff}, {\sc water},
702	+	and {\sc eam}) which are explained in the following sections.
703	+
704		\subsection{\label{sec:LJPot}The Lennard Jones Force Field}
705
706		The most basic force field implemented in {\sc oopse} is the
707	<	Lennard-Jones force field, which mimics the van der Waals interaction at
708	<	long distances, and uses an empirical repulsion at short
707	>	Lennard-Jones force field, which mimics the van der Waals interaction
708	>	at long distances and uses an empirical repulsion at short
709		distances. The Lennard-Jones potential is given by:
710		\begin{equation}
711		V_{\text{LJ}}(r_{ij}) =
#	Line 221 | Line 718 | $\epsilon_{ij}$ scales the well depth of the potential
718		where $r_{ij}$ is the distance between particles $i$ and $j$,
719		$\sigma_{ij}$ scales the length of the interaction, and
720		$\epsilon_{ij}$ scales the well depth of the potential. Scheme
721	<	\ref{sch:LJFF} gives an example \texttt{.bass} file that
721	>	\ref{sch:LJFF} gives an example meta-data file that
722		sets up a system of 108 Ar particles to be simulated using the
723		Lennard-Jones force field.
724
725	<	\begin{lstlisting}[float,caption={[Invocation of the Lennard-Jones force field] A sample system using the Lennard-Jones force field.},label={sch:LJFF}]
725	>	\begin{lstlisting}[float,caption={[Invocation of the Lennard-Jones
726	>	force field] A sample meta-data file for a small Lennard-Jones
727	>	simulation.},label={sch:LJFF}]
728
729	<	#include "argon.mdl"
729	>	#include "argon.md"
730
731		nComponents = 1;
732		component{
#	Line 235 | Line 734 | initialConfig = "./argon.init";
734		nMol = 108;
735		}
736
737	<	initialConfig = "./argon.init";
737	>	initialConfig = "./argon.in";
738
739		forceField = "LJ";
740		\end{lstlisting}
741
243	–	Because this potential is calculated between all pairs, the force
244	–	evaluation can become computationally expensive for large systems. To
245	–	keep the pair evaluations to a manageable number, {\sc oopse} employs
246	–	a cut-off radius.\cite{allen87:csl} The cutoff radius can either be
247	–	specified in the \texttt{.bass} file, or left as its default value of
248	–	$2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest Lennard-Jones
249	–	length parameter present in the simulation. Truncating the calculation
250	–	at $r_{\text{cut}}$ introduces a discontinuity into the potential
251	–	energy and the force. To offset this discontinuity in the potential,
252	–	the energy value at $r_{\text{cut}}$ is subtracted from the
253	–	potential. This causes the potential to go to zero smoothly at the
254	–	cut-off radius, and preserves conservation of energy in integrating
255	–	the equations of motion. There still remains a discontinuity in the derivative (the forces), however, this does not significantly affect the dynamics.
256	–
742		Interactions between dissimilar particles requires the generation of
743	<	cross term parameters for $\sigma$ and $\epsilon$. These are
744	<	calculated through the Lorentz-Berthelot mixing
743	>	cross term parameters for $\sigma$ and $\epsilon$. These parameters
744	>	are determined using the Lorentz-Berthelot mixing
745		rules:\cite{allen87:csl}
746		\begin{equation}
747		\sigma_{ij} = \frac{1}{2}[\sigma_{ii} + \sigma_{jj}],
#	Line 271 | Line 756 | simulate lipid bilayers. The simulations require a mod
756		\subsection{\label{oopseSec:DUFF}Dipolar Unified-Atom Force Field}
757
758		The dipolar unified-atom force field ({\sc duff}) was developed to
759	<	simulate lipid bilayers. The simulations require a model capable of
760	<	forming bilayers, while still being sufficiently computationally
761	<	efficient to allow large systems ($\sim$100's of phospholipids,
762	<	$\sim$1000's of waters) to be simulated for long times
763	<	($\sim$10's of nanoseconds).
759	>	simulate lipid bilayers. These types of simulations require a model
760	>	capable of forming bilayers, while still being sufficiently
761	>	computationally efficient to allow large systems ($\sim$100's of
762	>	phospholipids, $\sim$1000's of waters) to be simulated for long times
763	>	($\sim$10's of nanoseconds). With this goal in mind, {\sc duff} has no
764	>	point charges. Charge-neutral distributions are replaced with dipoles,
765	>	while most atoms and groups of atoms are reduced to Lennard-Jones
766	>	interaction sites. This simplification reduces the length scale of
767	>	long range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$,
768	>	removing the need for the computationally expensive Ewald
769	>	sum. Instead, Verlet neighbor-lists and cutoff radii are used for the
770	>	dipolar interactions, and, if desired, a reaction field may be added
771	>	to mimic longer range interactions.
772
280	–	With this goal in mind, {\sc duff} has no point
281	–	charges. Charge-neutral distributions were replaced with dipoles,
282	–	while most atoms and groups of atoms were reduced to Lennard-Jones
283	–	interaction sites. This simplification cuts the length scale of long
284	–	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.
288	–
773		As an example, lipid head-groups in {\sc duff} are represented as
774	<	point dipole interaction sites. By placing a dipole at the head
775	<	group's center of mass, our model mimics the charge separation found
776	<	in common phospholipid head groups such as
777	<	phosphatidylcholine.\cite{Cevc87} Additionally, a large Lennard-Jones
778	<	site is located at the pseudoatom's center of mass. The model is
779	<	illustrated by the red atom in Fig.~\ref{oopseFig:lipidModel}. The
780	<	water model we use to complement the dipoles of the lipids is our
781	<	reparameterization of the soft sticky dipole (SSD) model of Ichiye
774	>	point dipole interaction sites.  Placing a dipole at the head group's
775	>	center of mass mimics the charge separation found in common
776	>	phospholipid head groups such as phosphatidylcholine.\cite{Cevc87}
777	>	Additionally, a large Lennard-Jones site is located at the
778	>	pseudoatom's center of mass. The model is illustrated by the red atom
779	>	in Fig.~\ref{oopseFig:lipidModel}. The water model we use to
780	>	complement the dipoles of the lipids is a
781	>	reparameterization\cite{fennell04} of the soft sticky dipole (SSD)
782	>	model of Ichiye
783		\emph{et al.}\cite{liu96:new_model}
784
785		\begin{figure}
786		\centering
787	<	\includegraphics[width=\linewidth]{twoChainFig.pdf}
788	<	\caption[A representation of a lipid model in {\sc duff}]{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ %
789	<	is the bend angle, and $\mu$ is the dipole moment of the head group.}
787	>	\includegraphics[width=\linewidth]{lipidModel.eps}
788	>	\caption[A representation of a lipid model in {\sc duff}]{A
789	>	representation of the lipid model. $\phi$ is the torsion angle,
790	>	$\theta$ is the bend angle, and $\mu$ is the dipole moment of the head
791	>	group.}
792		\label{oopseFig:lipidModel}
793		\end{figure}
794
795	<	We have used a set of scalable parameters to model the alkyl groups
796	<	with Lennard-Jones sites. For this, we have borrowed parameters from
797	<	the TraPPE force field of Siepmann
798	<	\emph{et al}.\cite{Siepmann1998} TraPPE is a unified-atom
799	<	representation of n-alkanes, which is parametrized against phase
800	<	equilibria using Gibbs ensemble Monte Carlo simulation
801	<	techniques.\cite{Siepmann1998} One of the advantages of TraPPE is that
802	<	it generalizes the types of atoms in an alkyl chain to keep the number
803	<	of pseudoatoms to a minimum; the parameters for a unified atom such as
804	<	$\text{CH}_2$ do not change depending on what species are bonded to
318	<	it.
795	>	A set of scalable parameters has been used to model the alkyl groups
796	>	with Lennard-Jones sites. For this, parameters from the TraPPE force
797	>	field of Siepmann \emph{et al.}\cite{Siepmann1998} have been
798	>	utilized. TraPPE is a unified-atom representation of n-alkanes which
799	>	is parametrized against phase equilibria using Gibbs ensemble Monte
800	>	Carlo simulation techniques.\cite{Siepmann1998} One of the advantages
801	>	of TraPPE is that it generalizes the types of atoms in an alkyl chain
802	>	to keep the number of pseudoatoms to a minimum; thus, the parameters
803	>	for a unified atom such as $\text{CH}_2$ do not change depending on
804	>	what species are bonded to it.
805
806	<	TraPPE also constrains all bonds to be of fixed length. Typically,
807	<	bond vibrations are the fastest motions in a molecular dynamic
808	<	simulation. Small time steps between force evaluations must be used to
809	<	ensure adequate energy conservation in the bond degrees of freedom. By
810	<	constraining the bond lengths, larger time steps may be used when
811	<	integrating the equations of motion. A simulation using {\sc duff} is
812	<	illustrated in Scheme \ref{sch:DUFF}.
806	>	As is required by TraPPE, {\sc duff} also constrains all bonds to be
807	>	of fixed length. Typically, bond vibrations are the fastest motions in
808	>	a molecular dynamic simulation.  With these vibrations present, small
809	>	time steps between force evaluations must be used to ensure adequate
810	>	energy conservation in the bond degrees of freedom. By constraining
811	>	the bond lengths, larger time steps may be used when integrating the
812	>	equations of motion. A simulation using {\sc duff} is illustrated in
813	>	Scheme \ref{sch:DUFF}.
814
815	<	\begin{lstlisting}[float,caption={[Invocation of {\sc duff}]A portion of a \texttt{.bass} file showing a simulation utilizing {\sc duff}},label={sch:DUFF}]
815	>	\begin{lstlisting}[float,caption={[Invocation of {\sc duff}]A portion
816	>	of a meta-data file showing a simulation utilizing {\sc
817	>	duff}},label={sch:DUFF}]
818
819	<	#include "water.mdl"
820	<	#include "lipid.mdl"
819	>	#include "water.md"
820	>	#include "lipid.md"
821
822		nComponents = 2;
823		component{
#	Line 341 | Line 830 | initialConfig = "bilayer.init";
830		nMol = 1936;
831		}
832
833	<	initialConfig = "bilayer.init";
833	>	initialConfig = "bilayer.in";
834
835		forceField = "DUFF";
836
837		\end{lstlisting}
838
839	<	\subsection{\label{oopseSec:energyFunctions}{\sc duff} Energy Functions}
839	>	\subsubsection{\label{oopseSec:energyFunctions}{\sc duff} Energy Functions}
840
841		The total potential energy function in {\sc duff} is
842		\begin{equation}
#	Line 367 | Line 856 | within the molecule $I$, and $V_{\text{torsion}}$ is t
856		\label{eq:internalPotential}
857		\end{equation}
858		Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs
859	<	within the molecule $I$, and $V_{\text{torsion}}$ is the torsion potential
860	<	for all 1, 4 bonded pairs. The pairwise portions of the internal
861	<	potential are excluded for atom pairs that are involved in the same bond, bend, or torsion. All other atom pairs within the molecule are subject to the LJ pair potential.
859	>	within the molecule $I$, and $V_{\text{torsion}}$ is the torsion
860	>	potential for all 1, 4 bonded pairs.  The pairwise portions of the
861	>	non-bonded interactions are excluded for atom pairs that are involved
862	>	in the smae bond, bend, or torsion. All other atom pairs within a
863	>	molecule are subject to the LJ pair potential.
864
374	–
865		The bend potential of a molecule is represented by the following function:
866		\begin{equation}
867	<	V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 )^2, \label{eq:bendPot}
867	>	V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0
868	>	)^2, \label{eq:bendPot}
869		\end{equation}
870		where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$
871		(see Fig.~\ref{oopseFig:lipidModel}), $\theta_0$ is the equilibrium
#	Line 414 | Line 905 | evaluations are avoided during the calculation of the
905		k_3 &= 4c_3.
906		\end{align*}
907		By recasting the potential as a power series, repeated trigonometric
908	<	evaluations are avoided during the calculation of the potential energy.
908	>	evaluations are avoided during the calculation of the potential
909	>	energy.
910
911
912	<	The cross potential between molecules $I$ and $J$, $V^{IJ}_{\text{Cross}}$, is
913	<	as follows:
912	>	The cross potential between molecules $I$ and $J$,
913	>	$V^{IJ}_{\text{Cross}}$, is as follows:
914		\begin{equation}
915		V^{IJ}_{\text{Cross}} =
916		\sum_{i \in I} \sum_{j \in J}
#	Line 448 | Line 940 | respectively. $|\mu_i|$ is the magnitude of the dipole
940		Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing
941		towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$
942		are the orientational degrees of freedom for atoms $i$ and $j$
943	<	respectively. $|\mu_i|$ is the magnitude of the dipole moment of atom
944	<	$i$, $\boldsymbol{\hat{u}}_i$ is the standard unit orientation vector
945	<	of $\boldsymbol{\Omega}_i$, and $\boldsymbol{\hat{r}}_{ij}$ is the
946	<	unit vector pointing along $\mathbf{r}_{ij}$
943	>	respectively. The magnitude of the dipole moment of atom $i$ is
944	>	$|\mu_i|$, $\boldsymbol{\hat{u}}_i$ is the standard unit orientation
945	>	vector of $\boldsymbol{\Omega}_i$, and $\boldsymbol{\hat{r}}_{ij}$ is
946	>	the unit vector pointing along $\mathbf{r}_{ij}$
947		($\boldsymbol{\hat{r}}_{ij}=\mathbf{r}_{ij}/|\mathbf{r}_{ij}|$).
948
949	<	To improve computational efficiency of the dipole-dipole interactions,
950	<	{\sc oopse} employs an electrostatic cutoff radius. This parameter can
459	<	be set in the \texttt{.bass} file, and controls the length scale over
460	<	which dipole interactions are felt. To compensate for the
461	<	discontinuity in the potential and the forces at the cutoff radius, we
462	<	have implemented a switching function to smoothly scale the
463	<	dipole-dipole interaction at the cutoff.
464	<	\begin{equation}
465	<	S(r_{ij}) =
466	<	\begin{cases}
467	<	1 & \text{if $r_{ij} \le r_t$},\\
468	<	\frac{(r_{\text{cut}} + 2r_{ij} - 3r_t)(r_{\text{cut}} - r_{ij})^2}
469	<	{(r_{\text{cut}} - r_t)^2}
470	<	& \text{if $r_t < r_{ij} \le r_{\text{cut}}$}, \\
471	<	0 & \text{if $r_{ij} > r_{\text{cut}}$.}
472	<	\end{cases}
473	<	\label{eq:dipoleSwitching}
474	<	\end{equation}
475	<	Here $S(r_{ij})$ scales the potential at a given $r_{ij}$, and $r_t$
476	<	is the taper radius some given thickness less than the electrostatic
477	<	cutoff. The switching thickness can be set in the \texttt{.bass} file.
949	>	\subsubsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E
950	>	and SSD/RF}
951
479	–	\subsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF}
480	–
952		In the interest of computational efficiency, the default solvent used
953		by {\sc oopse} is the extended Soft Sticky Dipole (SSD/E) water
954		model.\cite{fennell04} The original SSD was developed by Ichiye
#	Line 536 | Line 1007 | Since SSD/E is a single-point {\it dipolar} model, the
1007		can be found in the original SSD
1008		articles.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md,Ichiye03}
1009
1010	+	\begin{figure}
1011	+	\centering
1012	+	\includegraphics[width=\linewidth]{waterAngle.eps}
1013	+	\caption[Coordinate definition for the SSD/E water model]{Coordinates
1014	+	for the interaction between two SSD/E water molecules.  $\theta_{ij}$
1015	+	is the angle that $r_{ij}$ makes with the $\hat{z}$ vector in the
1016	+	body-fixed frame for molecule $i$.  The $\hat{z}$ vector bisects the
1017	+	HOH angle in each water molecule. }
1018	+	\label{oopseFig:ssd}
1019	+	\end{figure}
1020	+
1021	+
1022		Since SSD/E is a single-point {\it dipolar} model, the force
1023		calculations are simplified significantly relative to the standard
1024		{\it charged} multi-point models. In the original Monte Carlo
1025		simulations using this model, Ichiye {\it et al.} reported that using
1026		SSD decreased computer time by a factor of 6-7 compared to other
1027	<	models.\cite{liu96:new_model} What is most impressive is that these savings
1028	<	did not come at the expense of accurate depiction of the liquid state
1029	<	properties.  Indeed, SSD/E maintains reasonable agreement with the Head-Gordon
1030	<	diffraction data for the structural features of liquid
1031	<	water.\cite{hura00,liu96:new_model} Additionally, the dynamical properties
1032	<	exhibited by SSD/E agree with experiment better than those of more
1033	<	computationally expensive models (like TIP3P and
1034	<	SPC/E).\cite{chandra99:ssd_md} The combination of speed and accurate depiction
1035	<	of solvent properties makes SSD/E a very attractive model for the
1036	<	simulation of large scale biochemical simulations.
1027	>	models.\cite{liu96:new_model} What is most impressive is that these
1028	>	savings did not come at the expense of accurate depiction of the
1029	>	liquid state properties.  Indeed, SSD/E maintains reasonable agreement
1030	>	with the Head-Gordon diffraction data for the structural features of
1031	>	liquid water.\cite{hura00,liu96:new_model} Additionally, the dynamical
1032	>	properties exhibited by SSD/E agree with experiment better than those
1033	>	of more computationally expensive models (like TIP3P and
1034	>	SPC/E).\cite{chandra99:ssd_md} The combination of speed and accurate
1035	>	depiction of solvent properties makes SSD/E a very attractive model
1036	>	for the simulation of large scale biochemical simulations.
1037
1038		Recent constant pressure simulations revealed issues in the original
1039		SSD model that led to lower than expected densities at all target
#	Line 559 | Line 1042 | model (an SSD variant  parameterized for reaction fiel
1042		exhibits improved liquid structure and transport behavior. If the use
1043		of a reaction field long-range interaction correction is desired, it
1044		is recommended that the parameters be modified to those of the SSD/RF
1045	<	model (an SSD variant  parameterized for reaction field). Solvent parameters can be easily modified in an accompanying
1046	<	\texttt{.bass} file as illustrated in the scheme below. A table of the
1047	<	parameter values and the drawbacks and benefits of the different
1048	<	density corrected SSD models can be found in
1049	<	reference~\cite{fennell04}.
567	<
568	<	\begin{lstlisting}[float,caption={[A simulation of {\sc ssd} water]A portion of a \texttt{.bass} file showing a simulation including {\sc ssd} water.},label={sch:ssd}]
569	<
570	<	#include "water.mdl"
571	<
572	<	nComponents = 1;
573	<	component{
574	<	type = "SSD_water";
575	<	nMol = 864;
576	<	}
577	<
578	<	initialConfig = "liquidWater.init";
1045	>	model (an SSD variant parameterized for reaction field). These solvent
1046	>	parameters are listed and can be easily modified in the {\sc duff}
1047	>	force field file ({\tt DUFF.frc}).  A table of the parameter values
1048	>	and the drawbacks and benefits of the different density corrected SSD
1049	>	models can be found in reference~\citen{fennell04}.
1050
580	–	forceField = "DUFF";
581	–
582	–	/*
583	–	* The following two flags set the cutoff
584	–	* radius for the electrostatic forces
585	–	* as well as the skin thickness of the switching
586	–	* function.
587	–	*/
588	–
589	–	electrostaticCutoffRadius  = 9.2;
590	–	electrostaticSkinThickness = 1.38;
591	–
592	–	\end{lstlisting}
593	–
594	–
1051		\subsection{\label{oopseSec:eam}Embedded Atom Method}
1052
1053	<	There are Molecular Dynamics packages which have the
1054	<	capacity to simulate metallic systems, including some that have
1055	<	parallel computational abilities\cite{plimpton93}. Potentials that
1056	<	describe bonding transition metal
1057	<	systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have an
602	<	attractive interaction which models  ``Embedding''
603	<	a positively charged metal ion in the electron density due to the
1053	>	{\sc oopse} implements a potential that describes bonding in
1054	>	transition metal
1055	>	systems.~\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} This
1056	>	potential has an attractive interaction which models ``Embedding'' a
1057	>	positively charged pseudo-atom core in the electron density due to the
1058		free valance ``sea'' of electrons created by the surrounding atoms in
1059	<	the system. A mostly-repulsive pairwise part of the potential
1060	<	describes the interaction of the positively charged metal core ions
1061	<	with one another. A particular potential description called the
1062	<	Embedded Atom Method\cite{Daw84,FBD86,johnson89,Lu97}({\sc eam}) that has
1063	<	particularly wide adoption has been selected for inclusion in {\sc oopse}. A
1064	<	good review of {\sc eam} and other metallic potential formulations was written
1065	<	by Voter.\cite{voter}
1059	>	the system.  A pairwise part of the potential (which is primarily
1060	>	repulsive) describes the interaction of the positively charged metal
1061	>	core ions with one another.  The Embedded Atom Method ({\sc
1062	>	eam})~\cite{Daw84,FBD86,johnson89,Lu97} has been widely adopted in the
1063	>	materials science community and has been included in {\sc oopse}. A
1064	>	good review of {\sc eam} and other formulations of metallic potentials
1065	>	was given by Voter.\cite{Voter:95}
1066
1067		The {\sc eam} potential has the form:
1068	<	\begin{eqnarray}
1069	<	V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i}
1070	<	\phi_{ij}({\bf r}_{ij}),  \\
1071	<	\rho_{i}  & = & \sum_{j \neq i} f_{j}({\bf r}_{ij}),
1072	<	\end{eqnarray}
619	<	where $F_{i} $ is the embedding function that equates the energy
1068	>	\begin{equation}
1069	>	V  =  \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i}
1070	>	\phi_{ij}({\bf r}_{ij})
1071	>	\end{equation}
1072	>	where $F_{i} $ is an embedding functional that approximates the energy
1073		required to embed a positively-charged core ion $i$ into a linear
1074		superposition of spherically averaged atomic electron densities given
1075	<	by $\rho_{i}$.  $\phi_{ij}$ is a primarily repulsive pairwise
623	<	interaction between atoms $i$ and $j$. In the original formulation of
624	<	{\sc eam}\cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term,
625	<	however in later refinements to {\sc eam} have shown that non-uniqueness
626	<	between $F$ and $\phi$ allow for more general forms for
627	<	$\phi$.\cite{Daw89} There is a cutoff distance, $r_{cut}$, which
628	<	limits the summations in the {\sc eam} equation to the few dozen atoms
629	<	surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$
630	<	interactions. Foiles \emph{et al}.~fit {\sc eam} potentials for the fcc
631	<	metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals.\cite{FBD86}
632	<	These fits are included in {\sc oopse}.
633	<
634	<	\subsection{\label{oopseSec:pbc}Periodic Boundary Conditions}
635	<
636	<	\newcommand{\roundme}{\operatorname{round}}
637	<
638	<	\textit{Periodic boundary conditions} are widely used to simulate bulk properties with a relatively small number of particles. The
639	<	simulation box is replicated throughout space to form an infinite
640	<	lattice.  During the simulation, when a particle moves in the primary
641	<	cell, its image in other cells move in exactly the same direction with
642	<	exactly the same orientation. Thus, as a particle leaves the primary
643	<	cell, one of its images will enter through the opposite face. If the
644	<	simulation box is large enough to avoid ``feeling'' the symmetries of
645	<	the periodic lattice, surface effects can be ignored. The available
646	<	periodic cells in OOPSE are cubic, orthorhombic and parallelepiped. We
647	<	use a $3 \times 3$ matrix, $\mathsf{H}$, to describe the shape and
648	<	size of the simulation box. $\mathsf{H}$ is defined:
1075	>	by $\rho_{i}$,
1076		\begin{equation}
1077	<	\mathsf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z ),
1077	>	\rho_{i}   =  \sum_{j \neq i} f_{j}({\bf r}_{ij}),
1078		\end{equation}
1079	<	where $\mathbf{h}_{\alpha}$ is the column vector of the $\alpha$ axis of the
1080	<	box.  During the course of the simulation both the size and shape of
1081	<	the box can be changed to allow volume fluctuations when constraining
1082	<	the pressure.
1079	>	Since the density at site $i$ ($\rho_i$) must be computed before the
1080	>	embedding functional can be evaluated, {\sc eam} and the related
1081	>	transition metal potentials require two loops through the atom pairs
1082	>	to compute the inter-atomic forces.
1083
1084	<	A real space vector, $\mathbf{r}$ can be transformed in to a box space
1085	<	vector, $\mathbf{s}$, and back through the following transformations:
1086	<	\begin{align}
1087	<	\mathbf{s} &= \mathsf{H}^{-1} \mathbf{r}, \\
1088	<	\mathbf{r} &= \mathsf{H} \mathbf{s}.
1089	<	\end{align}
1090	<	The vector $\mathbf{s}$ is now a vector expressed as the number of box
1091	<	lengths in the $\mathbf{h}_x$, $\mathbf{h}_y$, and $\mathbf{h}_z$
1092	<	directions. To find the minimum image of a vector $\mathbf{r}$, we
1093	<	first convert it to its corresponding vector in box space, and then,
1094	<	cast each element to lie in the range $[-0.5,0.5]$:
668	<	\begin{equation}
669	<	s_{i}^{\prime}=s_{i}-\roundme(s_{i}),
670	<	\end{equation}
671	<	where $s_i$ is the $i$th element of $\mathbf{s}$, and
672	<	$\roundme(s_i)$ is given by
673	<	\begin{equation}
674	<	\roundme(x) =
675	<	\begin{cases}
676	<	\lfloor x+0.5 \rfloor & \text{if $x \ge 0$,} \\
677	<	\lceil x-0.5 \rceil & \text{if $x < 0$.}
678	<	\end{cases}
679	<	\end{equation}
680	<	Here $\lfloor x \rfloor$ is the floor operator, and gives the largest
681	<	integer value that is not greater than $x$, and $\lceil x \rceil$ is
682	<	the ceiling operator, and gives the smallest integer that is not less
683	<	than $x$.  For example, $\roundme(3.6)=4$, $\roundme(3.1)=3$,
684	<	$\roundme(-3.6)=-4$, $\roundme(-3.1)=-3$.
1084	>	The pairwise portion of the potential, $\phi_{ij}$, is a primarily
1085	>	repulsive interaction between atoms $i$ and $j$. In the original
1086	>	formulation of {\sc eam}\cite{Daw84}, $\phi_{ij}$ was an entirely
1087	>	repulsive term; however later refinements to {\sc eam} allowed for
1088	>	more general forms for $\phi$.\cite{Daw89} The effective cutoff
1089	>	distance, $r_{{\text cut}}$ is the distance at which the values of
1090	>	$f(r)$ and $\phi(r)$ drop to zero for all atoms present in the
1091	>	simulation.  In practice, this distance is fairly small, limiting the
1092	>	summations in the {\sc eam} equation to the few dozen atoms
1093	>	surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$
1094	>	interactions.
1095
1096	<	Finally, we obtain the minimum image coordinates $\mathbf{r}^{\prime}$ by
1097	<	transforming back to real space,
1098	<	\begin{equation}
1099	<	\mathbf{r}^{\prime}=\mathsf{H}^{-1}\mathbf{s}^{\prime}.%
1100	<	\end{equation}
1101	<	In this way, particles are allowed to diffuse freely in $\mathbf{r}$,
1102	<	but their minimum images, $\mathbf{r}^{\prime}$ are used to compute
1103	<	the inter-atomic forces.
1096	>	In computing forces for alloys, mixing rules as outlined by
1097	>	Johnson~\cite{johnson89} are used to compute the heterogenous pair
1098	>	potential,
1099	>	\begin{eqnarray}
1100	>	\label{eq:johnson}
1101	>	\phi_{ab}(r)=\frac{1}{2}\left(
1102	>	\frac{f_{b}(r)}{f_{a}(r)}\phi_{aa}(r)+
1103	>	\frac{f_{a}(r)}{f_{b}(r)}\phi_{bb}(r)
1104	>	\right).
1105	>	\end{eqnarray}
1106	>	No mixing rule is needed for the densities, since the density at site
1107	>	$i$ is simply the linear sum of density contributions of all the other
1108	>	atoms.
1109
1110	+	The {\sc eam} force field illustrates an additional feature of {\sc
1111	+	oopse}.  Foiles, Baskes and Daw fit {\sc eam} potentials for Cu, Ag,
1112	+	Au, Ni, Pd, Pt and alloys of these metals.\cite{FBD86} These fits are
1113	+	included in {\sc oopse} as the {\tt u3} variant of the {\sc eam} force
1114	+	field.  Voter and Chen reparamaterized a set of {\sc eam} functions
1115	+	which do a better job of predicting melting points.\cite{Voter:87}
1116	+	These functions are included in {\sc oopse} as the {\tt VC} variant of
1117	+	the {\sc eam} force field.  An additional set of functions (the
1118	+	``Universal 6'' functions) are included in {\sc oopse} as the {\tt u6}
1119	+	variant of {\sc eam}.  For example, to specify the Voter-Chen variant
1120	+	of the {\sc eam} force field, the user would add the {\tt
1121	+	forceFieldVariant = "VC";} line to the meta-data file.
1122
1123	<	\section{\label{oopseSec:IOfiles}Input and Output Files}
1123	>	The potential files used by the {\sc eam} force field are in the
1124	>	standard {\tt funcfl} format, which is the format utilized by a number
1125	>	of other codes (e.g. ParaDyn~\cite{Paradyn}, {\sc dynamo 86}).  It
1126	>	should be noted that the energy units in these files are in eV, not
1127	>	$\mbox{kcal mol}^{-1}$ as in the rest of the {\sc oopse} force field
1128	>	files.
1129
1130	<	\subsection{{\sc bass} and Model Files}
1130	>	\subsection{\label{oopseSec:pbc}Periodic Boundary Conditions}
1131
1132	<	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}.
1132	>	\newcommand{\roundme}{\operatorname{round}}
1133
1134	<	\begin{lstlisting}[float,caption={[An example of a complete {\sc bass} file] An example showing a complete {\sc bass} file.},label={sch:bassExample}]
1135	<
1136	<	molecule{
1137	<	name = "Ar";
1138	<	nAtoms = 1;
1139	<	atom[0]{
1140	<	type="Ar";
1141	<	position( 0.0, 0.0, 0.0 );
1142	<	}
1143	<	}
1144	<
1145	<	nComponents = 1;
1146	<	component{
1147	<	type = "Ar";
1148	<	nMol = 108;
1149	<	}
1150	<
1151	<	initialConfig = "./argon.init";
1152	<
1153	<	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}
1134	>	\textit{Periodic boundary conditions} are widely used to simulate bulk
1135	>	properties with a relatively small number of particles. In this method
1136	>	the simulation box is replicated throughout space to form an infinite
1137	>	lattice.  During the simulation, when a particle moves in the primary
1138	>	cell, its image in other cells move in exactly the same direction with
1139	>	exactly the same orientation. Thus, as a particle leaves the primary
1140	>	cell, one of its images will enter through the opposite face. If the
1141	>	simulation box is large enough to avoid ``feeling'' the symmetries of
1142	>	the periodic lattice, surface effects can be ignored. The available
1143	>	periodic cells in {\sc oopse} are cubic, orthorhombic and
1144	>	parallelepiped.  {\sc oopse} use a $3 \times 3$ matrix, $\mathsf{H}$,
1145	>	to describe the shape and size of the simulation box. $\mathsf{H}$ is
1146	>	defined:
1147	>	\begin{equation}
1148	>	\mathsf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z ),
1149	>	\end{equation}
1150	>	where $\mathbf{h}_{\alpha}$ is the column vector of the $\alpha$ axis of the
1151	>	box.  During the course of the simulation both the size and shape of
1152	>	the box can be changed to allow volume fluctuations when constraining
1153	>	the pressure.
1154
1155	<
1156	<	There are three major files used by {\sc oopse} written in the
1157	<	coordinate format, they are as follows: the initialization file
1158	<	(\texttt{.init}), the simulation trajectory file (\texttt{.dump}), and
1159	<	the final coordinates of the simulation. The initialization file is
1160	<	necessary for {\sc oopse} to start the simulation with the proper
1161	<	coordinates, and is generated before the simulation run. The
1162	<	trajectory file is created at the beginning of the simulation, and is
1163	<	used to store snapshots of the simulation at regular intervals. The
1164	<	first frame is a duplication of the
1165	<	\texttt{.init} file, and each subsequent frame is appended to the file
1166	<	at an interval specified in the \texttt{.bass} file with the
1167	<	\texttt{sampleTime} flag. The final coordinate file is the end of run file. The
1168	<	\texttt{.eor} file stores the final configuration of the system for a
1169	<	given simulation. The file is updated at the same time as the
1170	<	\texttt{.dump} file, however, it only contains the most recent
1171	<	frame. In this way, an \texttt{.eor} file may be used as the
1172	<	initialization file to a second simulation in order to continue a
1173	<	simulation or recover one from a processor that has crashed during the
1174	<	course of the run.
1175	<
1176	<	\subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates}
1155	>	A real space vector, $\mathbf{r}$ can be transformed in to a box space
1156	>	vector, $\mathbf{s}$, and back through the following transformations:
1157	>	\begin{align}
1158	>	\mathbf{s} &= \mathsf{H}^{-1} \mathbf{r}, \\
1159	>	\mathbf{r} &= \mathsf{H} \mathbf{s}.
1160	>	\end{align}
1161	>	The vector $\mathbf{s}$ is now a vector expressed as the number of box
1162	>	lengths in the $\mathbf{h}_x$, $\mathbf{h}_y$, and $\mathbf{h}_z$
1163	>	directions. To find the minimum image of a vector $\mathbf{r}$, {\sc
1164	>	oopse} first converts it to its corresponding vector in box space, and
1165	>	then casts each element to lie in the range $[-0.5,0.5]$:
1166	>	\begin{equation}
1167	>	s_{i}^{\prime}=s_{i}-\roundme(s_{i}),
1168	>	\end{equation}
1169	>	where $s_i$ is the $i$th element of $\mathbf{s}$, and
1170	>	$\roundme(s_i)$ is given by
1171	>	\begin{equation}
1172	>	\roundme(x) =
1173	>	\begin{cases}
1174	>	\lfloor x+0.5 \rfloor & \text{if $x \ge 0$,} \\
1175	>	\lceil x-0.5 \rceil & \text{if $x < 0$.}
1176	>	\end{cases}
1177	>	\end{equation}
1178	>	Here $\lfloor x \rfloor$ is the floor operator, and gives the largest
1179	>	integer value that is not greater than $x$, and $\lceil x \rceil$ is
1180	>	the ceiling operator, and gives the smallest integer that is not less
1181	>	than $x$.
1182
1183	<	As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization
1184	<	file is needed to provide the starting coordinates for a
1185	<	simulation. The {\sc oopse} package provides several system building
1186	<	programs to aid in the creation of the \texttt{.init}
1187	<	file. The programs use {\sc bass}, and will recognize
1188	<	arguments and parameters in the \texttt{.bass} file that would
1189	<	otherwise be ignored by the simulation.
1183	>	Finally, the minimum image coordinates $\mathbf{r}^{\prime}$ are
1184	>	obtained by transforming back to real space,
1185	>	\begin{equation}
1186	>	\mathbf{r}^{\prime}=\mathsf{H}^{-1}\mathbf{s}^{\prime}.%
1187	>	\end{equation}
1188	>	In this way, particles are allowed to diffuse freely in $\mathbf{r}$,
1189	>	but their minimum images, or $\mathbf{r}^{\prime}$, are used to compute
1190	>	the inter-atomic forces.
1191
840	–	\subsection{The Statistics File}
1192
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.
1193
1194		\section{\label{oopseSec:mechanics}Mechanics}
1195
#	Line 866 | Line 1207 | motion for $\phi$ and $\psi$,\cite{allen87:csl} leadin
1207		Previous integration methods for orientational motion have problems
1208		that are avoided in the DLM method.  Direct propagation of the Euler
1209		angles has a known $1/\sin\theta$ divergence in the equations of
1210	<	motion for $\phi$ and $\psi$,\cite{allen87:csl} leading to
1211	<	numerical instabilities any time one of the directional atoms or rigid
1212	<	bodies has an orientation near $\theta=0$ or $\theta=\pi$.  More
1213	<	modern quaternion-based integration methods have relatively poor
1214	<	energy conservation.  While quaternions work well for orientational
1215	<	motion in other ensembles, the microcanonical ensemble has a
1216	<	constant energy requirement that is quite sensitive to errors in the
1217	<	equations of motion.  An earlier implementation of {\sc oopse}
1218	<	utilized quaternions for propagation of rotational motion; however, a
1219	<	detailed investigation showed that they resulted in a steady drift in
879	<	the total energy, something that has been observed by
880	<	Laird {\it et al.}\cite{Laird97}
1210	>	motion for $\phi$ and $\psi$,\cite{allen87:csl} leading to numerical
1211	>	instabilities any time one of the directional atoms or rigid bodies
1212	>	has an orientation near $\theta=0$ or $\theta=\pi$.  Quaternion-based
1213	>	integration methods work well for propagating orientational motion;
1214	>	however, energy conservation concerns arise when using the
1215	>	microcanonical (NVE) ensemble.  An earlier implementation of {\sc
1216	>	oopse} utilized quaternions for propagation of rotational motion;
1217	>	however, a detailed investigation showed that they resulted in a
1218	>	steady drift in the total energy, something that has been observed by
1219	>	Laird {\it et al.}\cite{Laird97}
1220
1221		The key difference in the integration method proposed by Dullweber
1222		\emph{et al.} is that the entire $3 \times 3$ rotation matrix is
#	Line 951 | Line 1290 | is equivalent to the more familiar body-fixed forms,
1290		represented by ${\bf j}$.  This equation of motion for angular momenta
1291		is equivalent to the more familiar body-fixed forms,
1292		\begin{eqnarray}
1293	<	\dot{j_{x}} & = & \tau^b_x(t)  +
1294	<	\left(\overleftrightarrow{\mathsf{I}}_{yy} - \overleftrightarrow{\mathsf{I}}_{zz} \right) j_y j_z, \\
1295	<	\dot{j_{y}} & = & \tau^b_y(t) +
1296	<	\left(\overleftrightarrow{\mathsf{I}}_{zz} - \overleftrightarrow{\mathsf{I}}_{xx} \right) j_z j_x,\\
1297	<	\dot{j_{z}} & = & \tau^b_z(t) +
1298	<	\left(\overleftrightarrow{\mathsf{I}}_{xx} - \overleftrightarrow{\mathsf{I}}_{yy} \right) j_x j_y,
1293	>	\dot{j_{x}} & = & \tau^b_x(t)  -
1294	>	\left(\overleftrightarrow{\mathsf{I}}_{yy}^{-1} - \overleftrightarrow{\mathsf{I}}_{zz}^{-1} \right) j_y j_z, \\
1295	>	\dot{j_{y}} & = & \tau^b_y(t) -
1296	>	\left(\overleftrightarrow{\mathsf{I}}_{zz}^{-1} - \overleftrightarrow{\mathsf{I}}_{xx}^{-1} \right) j_z j_x,\\
1297	>	\dot{j_{z}} & = & \tau^b_z(t) -
1298	>	\left(\overleftrightarrow{\mathsf{I}}_{xx}^{-1} - \overleftrightarrow{\mathsf{I}}_{yy}^{-1} \right) j_x j_y,
1299		\end{eqnarray}
1300		which utilize the body-fixed torques, ${\bf \tau}^b$. Torques are
1301		most easily derived in the space-fixed frame,
#	Line 1078 | Line 1417 | Fig.~\ref{timestep}.
1417
1418		\begin{figure}
1419		\centering
1420	<	\includegraphics[width=\linewidth]{timeStep.pdf}
1420	>	\includegraphics[width=\linewidth]{timeStep.eps}
1421		\caption[Energy conservation for quaternion versus DLM dynamics]{Energy conservation using quaternion based integration versus
1422		the method proposed by Dullweber \emph{et al.} with increasing time
1423		step. For each time step, the dotted line is total energy using the
#	Line 1106 | Line 1445 | and that is the time step for integrating the equation
1445
1446		\begin{center}
1447		\begin{tabular}{llll}
1448	<	{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf
1448	>	{\bf variable} & {\bf Meta-data keyword} & {\bf units} & {\bf
1449		default value} \\
1450		$h$ & {\tt dt = 2.0;} & fs & none
1451		\end{tabular}
#	Line 1116 | Line 1455 | integrator can selected with the {\tt ensemble} keywor
1455
1456		{\sc oopse} implements a number of extended system integrators for
1457		sampling from other ensembles relevant to chemical physics.  The
1458	<	integrator can selected with the {\tt ensemble} keyword in the
1459	<	{\tt .bass} file:
1458	>	integrator can be selected with the {\tt ensemble} keyword in the
1459	>	meta-data file:
1460
1461		\begin{center}
1462		\begin{tabular}{lll}
1463	<	{\bf Integrator} & {\bf Ensemble} & {\bf {\tt .bass} line} \\
1463	>	{\bf Integrator} & {\bf Ensemble} & {\bf Meta-data instruction} \\
1464		NVE & microcanonical & {\tt ensemble = NVE; } \\
1465		NVT & canonical & {\tt ensemble = NVT; } \\
1466		NPTi & isobaric-isothermal & {\tt ensemble = NPTi;} \\
#	Line 1136 | Line 1475 | system, and has been shown to sample the canonical dis
1475		The relatively well-known Nos\'e-Hoover thermostat\cite{Hoover85} is
1476		implemented in {\sc oopse}'s NVT integrator.  This method couples an
1477		extra degree of freedom (the thermostat) to the kinetic energy of the
1478	<	system, and has been shown to sample the canonical distribution in the
1479	<	system degrees of freedom while conserving a quantity that is, to
1478	>	system and it has been shown to sample the canonical distribution in
1479	>	the system degrees of freedom while conserving a quantity that is, to
1480		within a constant, the Helmholtz free energy.\cite{melchionna93}
1481
1482		NPT algorithms attempt to maintain constant pressure in the system by
#	Line 1161 | Line 1500 | variables.
1500
1501		\begin{center}
1502		\begin{tabular}{llll}
1503	<	{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf
1503	>	{\bf variable} & {\bf Meta-data instruction} & {\bf units} & {\bf
1504		default value} \\
1505		$T_{\mathrm{target}}$ & {\tt targetTemperature = 300;} &  K & none \\
1506		$P_{\mathrm{target}}$ & {\tt targetPressure = 1;} & atm & none \\
#	Line 1219 | Line 1558 | $\tau_T$ or $T_{\mathrm{target}}$ in a simulation, one
1558		In eq.(\ref{eq:nosehooverext}), $\tau_T$ is the time constant for
1559		relaxation of the temperature to the target value.  To set values for
1560		$\tau_T$ or $T_{\mathrm{target}}$ in a simulation, one would use the
1561	<	{\tt tauThermostat} and {\tt targetTemperature} keywords in the {\tt
1562	<	.bass} file.  The units for {\tt tauThermostat} are fs, and the units
1563	<	for the {\tt targetTemperature} are degrees K.   The integration of
1564	<	the equations of motion is carried out in a velocity-Verlet style 2
1561	>	{\tt tauThermostat} and {\tt targetTemperature} keywords in the
1562	>	meta-data file.  The units for {\tt tauThermostat} are fs, and the
1563	>	units for the {\tt targetTemperature} are degrees K.   The integration
1564	>	of the equations of motion is carried out in a velocity-Verlet style 2
1565		part algorithm:
1566
1567		{\tt moveA:}
#	Line 1284 | Line 1623 | Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are requir
1623		\chi(t + h) \right) .
1624		\end{align*}
1625
1626	<	Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to caclculate
1626	>	Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to calculate
1627		$T(t + h)$ as well as $\chi(t + h)$, they indirectly depend on their
1628		own values at time $t + h$.  {\tt moveB} is therefore done in an
1629		iterative fashion until $\chi(t + h)$ becomes self-consistent.  The
#	Line 1312 | Line 1651 | To carry out isobaric-isothermal ensemble calculations
1651		\subsection{\label{sec:NPTi}Constant-pressure integration with
1652		isotropic box deformations (NPTi)}
1653
1654	<	To carry out isobaric-isothermal ensemble calculations {\sc oopse}
1654	>	To carry out isobaric-isothermal ensemble calculations, {\sc oopse}
1655		implements the Melchionna modifications to the Nos\'e-Hoover-Andersen
1656	<	equations of motion,\cite{melchionna93}
1656	>	equations of motion.\cite{melchionna93} The equations of motion are
1657	>	the same as NVT with the following exceptions:
1658
1659		\begin{eqnarray}
1660		\dot{{\bf r}} & = & {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right), \\
1661		\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) ,\\
1662		\dot{\eta} & = & \frac{1}{\tau_{B}^2 f k_B T_{\mathrm{target}}} V \left( P -
1663		P_{\mathrm{target}} \right), \\
1664		\dot{\mathcal{V}} & = & 3 \mathcal{V} \eta .
#	Line 1352 | Line 1685 | outer} product of the velocities denoted by the $\otim
1685		\overleftrightarrow{\mathsf{W}}(t).
1686		\end{equation}
1687		The kinetic contribution to the pressure tensor utilizes the {\it
1688	<	outer} product of the velocities denoted by the $\otimes$ symbol.  The
1688	>	outer} product of the velocities, denoted by the $\otimes$ symbol.  The
1689		stress tensor is calculated from another outer product of the
1690		inter-atomic separation vectors (${\bf r}_{ij} = {\bf r}_j - {\bf
1691		r}_i$) with the forces between the same two atoms,
#	Line 1360 | Line 1693 | The instantaneous pressure is then simply obtained fro
1693		\overleftrightarrow{\mathsf{W}}(t) = \sum_{i} \sum_{j>i} {\bf r}_{ij}(t)
1694		\otimes {\bf f}_{ij}(t).
1695		\end{equation}
1696	+	In systems containing cutoff groups, the stress tensor is computed
1697	+	between the centers-of-mass of the cutoff groups:
1698	+	\begin{equation}
1699	+	\overleftrightarrow{\mathsf{W}}(t) = \sum_{a} \sum_{b} {\bf r}_{ab}(t)
1700	+	\otimes {\bf f}_{ab}(t).
1701	+	\end{equation}
1702	+	where ${\bf r}_{ab}$ is the distance between the centers of mass, and
1703	+	\begin{equation}
1704	+	{\bf f}_{ab} = s(r_{ab}) \sum_{i \in a} \sum_{j \in b} {\bf f}_{ij} +
1705	+	s\prime(r_{ab}) \frac{{\bf r}_{ab}}{|r_{ab}|} \sum_{i \in a} \sum_{j
1706	+	\in b} V_{ij}({\bf r}_{ij}).
1707	+	\end{equation}
1708	+
1709		The instantaneous pressure is then simply obtained from the trace of
1710	<	the Pressure tensor,
1710	>	the pressure tensor,
1711		\begin{equation}
1712		P(t) = \frac{1}{3} \mathrm{Tr} \left( \overleftrightarrow{\mathsf{P}}(t).
1713		\right)
#	Line 1370 | Line 1716 | $\tau_B$ or $P_{\mathrm{target}}$ in a simulation, one
1716		In eq.(\ref{eq:melchionna1}), $\tau_B$ is the time constant for
1717		relaxation of the pressure to the target value.  To set values for
1718		$\tau_B$ or $P_{\mathrm{target}}$ in a simulation, one would use the
1719	<	{\tt tauBarostat} and {\tt targetPressure} keywords in the {\tt .bass}
1719	>	{\tt tauBarostat} and {\tt targetPressure} keywords in the meta-data
1720		file.  The units for {\tt tauBarostat} are fs, and the units for the
1721		{\tt targetPressure} are atmospheres.  Like in the NVT integrator, the
1722		integration of the equations of motion is carried out in a
1723	<	velocity-Verlet style 2 part algorithm:
1723	>	velocity-Verlet style 2 part algorithm with only the following differences:
1724
1725		{\tt moveA:}
1726		\begin{align*}
1381	–	T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\
1382	–	%
1727		P(t) &\leftarrow \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\} ,\\
1728		%
1729		{\bf v}\left(t + h / 2\right)  &\leftarrow {\bf v}(t)
1730		+ \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t)
1731		\left(\chi(t) + \eta(t) \right) \right), \\
1732		%
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	–	%
1733		\eta(t + h / 2) &\leftarrow \eta(t) + \frac{h
1734		\mathcal{V}(t)}{2 N k_B T(t) \tau_B^2} \left( P(t)
1735		- P_{\mathrm{target}} \right), \\
#	Line 1411 | Line 1743 | Most of these equations are identical to their counter
1743		\mathsf{H}(t).
1744		\end{align*}
1745
1746	<	Most of these equations are identical to their counterparts in the NVT
1415	<	integrator, but the propagation of positions to time $t + h$
1746	>	The propagation of positions to time $t + h$
1747		depends on the positions at the same time.  {\sc oopse} carries out
1748		this step iteratively (with a limit of 5 passes through the iterative
1749		loop).  Also, the simulation box $\mathsf{H}$ is scaled uniformly for
#	Line 1421 | Line 1752 | the box by
1752		h / 2$.  Reshaping the box uniformly also scales the volume of
1753		the box by
1754		\begin{equation}
1755	<	\mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)}.
1755	>	\mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)} \times
1756		\mathcal{V}(t)
1757		\end{equation}
1758
#	Line 1432 | Line 1763 | T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
1763
1764		{\tt moveB:}
1765		\begin{align*}
1435	–	T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
1436	–	\left\{{\bf j}(t + h)\right\} ,\\
1437	–	%
1766		P(t + h) &\leftarrow  \left\{{\bf r}(t + h)\right\},
1767		\left\{{\bf v}(t + h)\right\}, \\
1768		%
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	–	%
1769		\eta(t + h) &\leftarrow \eta(t + h / 2) +
1770		\frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h)
1771		\tau_B^2} \left( P(t + h) - P_{\mathrm{target}} \right), \\
#	Line 1458 | Line 1782 | to caclculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, a
1782		\end{align*}
1783
1784		Once again, since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required
1785	<	to caclculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, and $\eta(t +
1785	>	to calculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, and $\eta(t +
1786		h)$, they indirectly depend on their own values at time $t + h$.  {\tt
1787		moveB} is therefore done in an iterative fashion until $\chi(t + h)$
1788		and $\eta(t + h)$ become self-consistent.  The relative tolerance for
#	Line 1498 | Line 1822 | the box shape.  The equations of motion for this metho
1822		{\it shape} as well as in the volume of the box.  This method utilizes
1823		the full $3 \times 3$ pressure tensor and introduces a tensor of
1824		extended variables ($\overleftrightarrow{\eta}$) to control changes to
1825	<	the box shape.  The equations of motion for this method are
1825	>	the box shape.  The equations of motion for this method differ from
1826	>	those of NPTi as follows:
1827		\begin{eqnarray}
1828		\dot{{\bf r}} & = & {\bf v} + \overleftrightarrow{\eta} \cdot \left( {\bf r} - {\bf R}_0 \right), \\
1829		\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\overleftrightarrow{\eta} +
1830		\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) ,\\
1831		\dot{\overleftrightarrow{\eta}} & = & \frac{1}{\tau_{B}^2 f k_B
1832		T_{\mathrm{target}}} V \left( \overleftrightarrow{\mathsf{P}} - P_{\mathrm{target}}\mathsf{1} \right) ,\\
1833		\dot{\mathsf{H}} & = &  \overleftrightarrow{\eta} \cdot \mathsf{H} .
#	Line 1525 | Line 1843 | T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\b
1843
1844		{\tt moveA:}
1845		\begin{align*}
1528	–	T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\
1529	–	%
1846		\overleftrightarrow{\mathsf{P}}(t) &\leftarrow \left\{{\bf r}(t)\right\},
1847		\left\{{\bf v}(t)\right\} ,\\
1848		%
#	Line 1535 | Line 1851 | T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\b
1851		\left(\chi(t)\mathsf{1} + \overleftrightarrow{\eta}(t) \right) \cdot
1852		{\bf v}(t) \right), \\
1853		%
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	–	%
1854		\overleftrightarrow{\eta}(t + h / 2) &\leftarrow
1855		\overleftrightarrow{\eta}(t) + \frac{h \mathcal{V}(t)}{2 N k_B
1856		T(t) \tau_B^2} \left( \overleftrightarrow{\mathsf{P}}(t)
#	Line 1568 | Line 1872 | T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
1872
1873		{\tt moveB:}
1874		\begin{align*}
1571	–	T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
1572	–	\left\{{\bf j}(t + h)\right\}, \\
1573	–	%
1875		\overleftrightarrow{\mathsf{P}}(t + h) &\leftarrow \left\{{\bf r}
1876		(t + h)\right\}, \left\{{\bf v}(t
1877		+ h)\right\}, \left\{{\bf f}(t + h)\right\} ,\\
1878		%
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	–	%
1879		\overleftrightarrow{\eta}(t + h) &\leftarrow
1880		\overleftrightarrow{\eta}(t + h / 2) +
1881		\frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h)
#	Line 1590 | Line 1887 | T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
1887		\frac{{\bf f}(t + h)}{m} -
1888		(\chi(t + h)\mathsf{1} + \overleftrightarrow{\eta}(t
1889		+ 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) .
1890		\end{align*}
1891
1892		The iterative schemes for both {\tt moveA} and {\tt moveB} are
#	Line 1611 | Line 1904 | elongated and sheared geometries which become smaller
1904		This integrator must be used with care, particularly in liquid
1905		simulations.  Liquids have very small restoring forces in the
1906		off-diagonal directions, and the simulation box can very quickly form
1907	<	elongated and sheared geometries which become smaller than the
1908	<	electrostatic or Lennard-Jones cutoff radii.  The NPTf integrator
1909	<	finds most use in simulating crystals or liquid crystals which assume
1617	<	non-orthorhombic geometries.
1907	>	elongated and sheared geometries which become smaller than the cutoff
1908	>	radius.  The NPTf integrator finds most use in simulating crystals or
1909	>	liquid crystals which assume non-orthorhombic geometries.
1910
1911		\subsection{\label{nptxyz}Constant pressure in 3 axes (NPTxyz)}
1912
#	Line 1637 | Line 1929 | simulations).
1929		orientational anisotropy in the system (i.e. in lipid bilayer
1930		simulations).
1931
1932	<	\subsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond
1932	>	\subsection{\label{sec:constraints}Constraint Methods}
1933	>
1934	>	\subsubsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond
1935		Constraints}
1936
1937		In order to satisfy the constraints of fixed bond lengths within {\sc
1938		oopse}, we have implemented the {\sc rattle} algorithm of
1939	<	Andersen.\cite{andersen83} The algorithm is a velocity verlet
1940	<	formulation of the {\sc shake} method\cite{ryckaert77} of iteratively
1941	<	solving the Lagrange multipliers of constraint. The system of Lagrange
1942	<	multipliers allows one to reformulate the equations of motion with
1943	<	explicit constraint forces.\cite{fowles99:lagrange}
1939	>	Andersen.\cite{andersen83} {\sc rattle} is a velocity-Verlet
1940	>	formulation of the {\sc shake} method\cite{ryckaert77} for iteratively
1941	>	solving the Lagrange multipliers which maintain the holonomic
1942	>	constraints.  Both methods are covered in depth in the
1943	>	literature,\cite{leach01:mm,allen87:csl} and a detailed description of
1944	>	this method would be redundant.
1945
1946	<	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}
1946	>	\subsubsection{\label{oopseSec:zcons}The Z-Constraint Method}
1947
1948	<	In a simulation, this would involve the solution of a set of $(m + n)$
1949	<	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}
1762	<
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	<
1878	<	Based on the fluctuation-dissipation theorem, a force auto-correlation
1879	<	method was developed by Roux and Karplus to investigate the dynamics
1948	>	A force auto-correlation method based on the fluctuation-dissipation
1949	>	theorem was developed by Roux and Karplus to investigate the dynamics
1950		of ions inside ion channels.\cite{Roux91} The time-dependent friction
1951		coefficient can be calculated from the deviation of the instantaneous
1952	<	force from its mean force.
1952	>	force from its mean value:
1953		\begin{equation}
1954		\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T,
1955		\end{equation}
#	Line 1888 | Line 1958 | where%
1958		\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle.
1959		\end{equation}
1960
1891	–
1961		If the time-dependent friction decays rapidly, the static friction
1962		coefficient can be approximated by
1963		\begin{equation}
1964		\xi_{\text{static}}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt.
1965		\end{equation}
1966	<	Allowing diffusion constant to then be calculated through the
1966	>
1967	>	This allows the diffusion constant to then be calculated through the
1968		Einstein relation:\cite{Marrink94}
1969		\begin{equation}
1970		D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty
1971		}\langle\delta F(z,t)\delta F(z,0)\rangle dt}.%
1972		\end{equation}
1973
1974	<	The Z-Constraint method, which fixes the z coordinates of the
1975	<	molecules with respect to the center of the mass of the system, has
1976	<	been a method suggested to obtain the forces required for the force
1977	<	auto-correlation calculation.\cite{Marrink94} However, simply resetting the
1978	<	coordinate will move the center of the mass of the whole system. To
1979	<	avoid this problem, a new method was used in {\sc oopse}. Instead of
1980	<	resetting the coordinate, we reset the forces of z-constrained
1981	<	molecules as well as subtract the total constraint forces from the
1982	<	rest of the system after the force calculation at each time step.
1974	>	The Z-Constraint method, which fixes the $z$ coordinates of a few
1975	>	``tagged'' molecules with respect to the center of the mass of the
1976	>	system is a technique that was proposed to obtain the forces required
1977	>	for the force auto-correlation calculation.\cite{Marrink94} However,
1978	>	simply resetting the coordinate will move the center of the mass of
1979	>	the whole system. To avoid this problem, we have developed a new
1980	>	method that is utilized in {\sc oopse}. Instead of resetting the
1981	>	coordinates, we reset the forces of $z$-constrained molecules and
1982	>	subtract the total constraint forces from the rest of the system after
1983	>	the force calculation at each time step.
1984
1985	<	After the force calculation, define $G_\alpha$ as
1985	>	After the force calculation, the total force on molecule $\alpha$,
1986		\begin{equation}
1987		G_{\alpha} = \sum_i F_{\alpha i},
1988		\label{oopseEq:zc1}
1989		\end{equation}
1990	<	where $F_{\alpha i}$ is the force in the z direction of atom $i$ in
1991	<	z-constrained molecule $\alpha$. The forces of the z constrained
1992	<	molecule are then set to:
1990	>	where $F_{\alpha i}$ is the force in the $z$ direction on atom $i$ in
1991	>	$z$-constrained molecule $\alpha$. The forces on the atoms in the
1992	>	$z$-constrained molecule are then adjusted to remove the total force
1993	>	on molecule $\alpha$:
1994		\begin{equation}
1995		F_{\alpha i} = F_{\alpha i} -
1996		\frac{m_{\alpha i} G_{\alpha}}{\sum_i m_{\alpha i}}.
1997		\end{equation}
1998	<	Here, $m_{\alpha i}$ is the mass of atom $i$ in the z-constrained
1999	<	molecule. Having rescaled the forces, the velocities must also be
2000	<	rescaled to subtract out any center of mass velocity in the z
2001	<	direction.
1998	>	Here, $m_{\alpha i}$ is the mass of atom $i$ in the $z$-constrained
1999	>	molecule.  After the forces have been adjusted, the velocities must
2000	>	also be modified to subtract out molecule $\alpha$'s center-of-mass
2001	>	velocity in the $z$ direction.
2002		\begin{equation}
2003		v_{\alpha i} = v_{\alpha i} -
2004		\frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}},
2005		\end{equation}
2006		where $v_{\alpha i}$ is the velocity of atom $i$ in the z direction.
2007	<	Lastly, all of the accumulated z constrained forces must be subtracted
2008	<	from the system to keep the system center of mass from drifting.
2007	>	Lastly, all of the accumulated constraint forces must be subtracted
2008	>	from the rest of the unconstrained system to keep the system center of
2009	>	mass of the entire system from drifting.
2010		\begin{equation}
2011		F_{\beta i} = F_{\beta i} - \frac{m_{\beta i} \sum_{\alpha} G_{\alpha}}
2012		{\sum_{\beta}\sum_i m_{\beta i}},
2013		\end{equation}
2014	<	where $\beta$ are all of the unconstrained molecules in the
2014	>	where $\beta$ denotes all {\it unconstrained} molecules in the
2015		system. Similarly, the velocities of the unconstrained molecules must
2016	<	also be scaled.
2016	>	also be scaled:
2017		\begin{equation}
2018	<	v_{\beta i} = v_{\beta i} + \sum_{\alpha}
2019	<	\frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}.
2018	>	v_{\beta i} = v_{\beta i} + \sum_{\alpha} \frac{\sum_i m_{\alpha i}
2019	>	v_{\alpha i}}{\sum_i m_{\alpha i}}.
2020		\end{equation}
2021
2022	<	At the very beginning of the simulation, the molecules may not be at their
2023	<	constrained positions. To move a z-constrained molecule to its specified
2024	<	position, a simple harmonic potential is used
2022	>	This method will pin down the centers-of-mass of all of the
2023	>	$z$-constrained molecules, and will also keep the entire system fixed
2024	>	at the original system center-of-mass location.
2025	>
2026	>	At the very beginning of the simulation, the molecules may not be at
2027	>	their desired positions. To steer a $z$-constrained molecule to its
2028	>	specified position, a simple harmonic potential is used:
2029		\begin{equation}
2030		U(t)=\frac{1}{2}k_{\text{Harmonic}}(z(t)-z_{\text{cons}})^{2},%
2031		\end{equation}
2032	<	where $k_{\text{Harmonic}}$ is the harmonic force constant, $z(t)$ is the
2033	<	current $z$ coordinate of the center of mass of the constrained molecule, and
2034	<	$z_{\text{cons}}$ is the constrained position. The harmonic force operating
2035	<	on the z-constrained molecule at time $t$ can be calculated by
2032	>	where $k_{\text{Harmonic}}$ is an harmonic force constant, $z(t)$ is
2033	>	the current $z$ coordinate of the center of mass of the constrained
2034	>	molecule, and $z_{\text{cons}}$ is the desired constraint
2035	>	position. The harmonic force operating on the $z$-constrained molecule
2036	>	at time $t$ can be calculated by
2037		\begin{equation}
2038		F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}=
2039		-k_{\text{Harmonic}}(z(t)-z_{\text{cons}}).
2040		\end{equation}
2041
2042	<	\section{\label{oopseSec:props}Trajectory Analysis}
2042	>	The user may also specify the use of a constant velocity method
2043	>	(steered molecular dynamics) to move the molecules to their desired
2044	>	initial positions.
2045
2046	<	\subsection{\label{oopseSec:staticProps}Static Property Analysis}
2046	>	To use of the $z$-constraint method in an {\sc oopse} simulation, the
2047	>	molecules must be specified using the {\tt nZconstraints} keyword in
2048	>	the meta-data file.  The other parameters for modifying the behavior
2049	>	of the $z$-constraint method are listed in table~\ref{table:zconParams}.
2050
1968	–	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}
1973	–
2051		\begin{table}
2052	<	\caption{THE DIFFERENT PAIR CORRELATIONS IN \texttt{staticProps}}
2053	<	\label{oopseTb:gofrs}
2052	>	\caption{The Global Keywords: Z-Constraint Parameters}
2053	>	\label{table:zconParams}
2054		\begin{center}
2055	<	\begin{tabular}{|l|c|c|}
2056	<	\hline
2057	<	Name      & Equation & Directional Atom \\ \hline
2058	<	$g_{\text{AB}}(r)$              & Eq.~\ref{eq:gofr}         & neither \\ \hline
2059	<	$g_{\text{AB}}(r, \cos \theta)$ & Eq.~\ref{eq:gofrCosTheta} & A \\ \hline
2060	<	$g_{\text{AB}}(r, \cos \omega)$ & Eq.~\ref{eq:gofrCosOmega} & both \\ \hline
2061	<	$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}
1993	<	\end{center}
1994	<	\end{table}
2055	>	% Note when adding or removing columns, the \hsize numbers must add up to the total number
2056	>	% of columns.
2057	>	\begin{tabularx}{\linewidth}%
2058	>	{>{\setlength{\hsize}{1.00\hsize}}X%
2059	>	>{\setlength{\hsize}{0.4\hsize}}X%
2060	>	>{\setlength{\hsize}{1.2\hsize}}X%
2061	>	>{\setlength{\hsize}{1.4\hsize}}X}
2062
2063	<	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$.
2063	>	{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline
2064
2065	<	The next two pair correlations, $g_{\text{AB}}(r, \cos \theta)$ and
2066	<	$g_{\text{AB}}(r, \cos \omega)$, are similar in that they are both two
2067	<	dimensional histograms. Both use $r$ for the primary axis then a
2068	<	$\cos$ for the secondary axis ($\cos \theta$ for
2069	<	Eq.~\ref{eq:gofrCosTheta} and $\cos \omega$ for
2070	<	Eq.~\ref{eq:gofrCosOmega}). This allows for the investigator to
2071	<	correlate alignment on directional entities. $g_{\text{AB}}(r, \cos
2072	<	\theta)$ is defined as follows:
2073	<	\begin{equation}
2074	<	g_{\text{AB}}(r, \cos \theta) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle
2075	<	\sum_{i \in \text{A}} \sum_{j \in \text{B}}
2076	<	\delta( \cos \theta - \cos \theta_{ij})
2077	<	\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}$.
2065	>	{\tt nZconstraints} & integer &  The number of zconstraint molecules& If using zconstraint method, {\tt nZconstraints} must be set \\
2066	>	{\tt zconsTime} & fs & Sets the frequency at which the {\tt .fz} file is written &  \\
2067	>	{\tt zconsForcePolicy} & string & The strategy of subtracting
2068	>	zconstraint force from the unconstrained molecules & Possible
2069	>	strategies are {\tt BYMASS} and {\tt BYNUMBER}. Default
2070	>	strategy is set to {\tt BYMASS}\\
2071	>	{\tt zconsGap} & $\mbox{\AA}$ & Set the distance between two adjacent
2072	>	constraint positions& Used mainly in moving molecules through a simulation \\
2073	>	{\tt zconsFixtime} & fs & Sets how long the zconstraint molecule is
2074	>	fixed & {\tt zconsFixtime} must be set if {\tt zconsGap} is set\\
2075	>	{\tt zconsUsingSMD} &logical & Flag for using Steered Molecular
2076	>	Dynamics or Harmonic Forces to move the molecule  & Harmonic Forces are
2077	>	used by default\\
2078
2079	<	The second two dimensional histogram is of the form:
2080	<	\begin{equation}
2081	<	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}
2039	<	\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$.
2079	>	\end{tabularx}
2080	>	\end{center}
2081	>	\end{table}
2082
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}$.
2083
2084	<	The final pair correlation is similar to
2061	<	Eq.~\ref{eq:gofrCosOmega}. $\langle \cos \omega
2062	<	\rangle_{\text{AB}}(r)$ is calculated in the following way:
2063	<	\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.
2084	>	\section{\label{sec:minimize}Energy Minimization}
2085
2086	<	\subsection{\label{dynamicProps}Dynamic Property Analysis}
2086	>	As one of the basic procedures of molecular modeling, energy
2087	>	minimization is used to identify local configurations that are stable
2088	>	points on the potential energy surface. There is a vast literature on
2089	>	energy minimization algorithms have been developed to search for the
2090	>	global energy minimum as well as to find local structures which are
2091	>	stable fixed points on the surface.  We have included two simple
2092	>	minimization algorithms: steepest descent, ({\sc sd}) and conjugate
2093	>	gradient ({\sc cg}) to help users find reasonable local minima from
2094	>	their initial configurations.
2095
2096	<	The dynamic properties of a trajectory are calculated with the program
2097	<	\texttt{dynamicProps}. The program calculates the following properties:
2098	<	\begin{gather}
2099	<	\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}
2096	>	Since {\sc oopse} handles atoms and rigid bodies which have
2097	>	orientational coordinates as well as translational coordinates, there
2098	>	is some subtlety to the choice of parameters for minimization
2099	>	algorithms.
2100
2101	<	Eq.~\ref{eq:rms} is the root mean square displacement function. Which
2102	<	allows one to observe the average displacement of an atom as a
2103	<	function of time. The quantity is useful when calculating diffusion
2104	<	coefficients because of the Einstein Relation, which is valid at long
2105	<	times.\cite{allen87:csl}
2101	>	Given a coordinate set $x_{k}$ and a search direction $d_{k}$, a line
2102	>	search algorithm is performed along $d_{k}$ to produce
2103	>	$x_{k+1}=x_{k}+$ $\lambda _{k}d_{k}$.
2104	>
2105	>	In the steepest descent ({\sc sd}) algorithm,%
2106		\begin{equation}
2107	<	2tD = \langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle.
2090	<	\label{oopseEq:einstein}
2107	>	d_{k}=-\nabla V(x_{k})
2108		\end{equation}
2109	<
2110	<	Eq.~\ref{eq:velCorr} and \ref{eq:angularVelCorr} are the translational
2111	<	velocity and angular velocity correlation functions respectively. The
2112	<	latter is only applicable to directional species in the
2113	<	simulation. The velocity autocorrelation functions are useful when
2114	<	determining vibrational information about the system of interest.
2109	>	The gradient and the direction of next step are always orthogonal.
2110	>	This may cause oscillatory behavior in narrow valleys.  To overcome
2111	>	this problem, the Fletcher-Reeves variant~\cite{FletcherReeves} of the
2112	>	conjugate gradient ({\sc cg}) algorithm is used to generate $d_{k+1}$
2113	>	via simple recursion:
2114	>	\begin{align}
2115	>	d_{k+1}  &  =-\nabla V(x_{k+1})+\gamma_{k}d_{k}\\
2116	>	\gamma_{k}  &  =\frac{\nabla V(x_{k+1})^{T}\nabla V(x_{k+1})}{\nabla
2117	>	V(x_{k})^{T}\nabla V(x_{k})}%
2118	>	\end{align}
2119
2120	<	\section{\label{oopseSec:design}Program Design}
2120	>	The Polak-Ribiere variant~\cite{PolakRibiere} of the conjugate
2121	>	gradient ($\gamma_{k}$) is defined as%
2122	>	\begin{equation}
2123	>	\gamma_{k}=\frac{[\nabla V(x_{k+1})-\nabla V(x)]^{T}\nabla V(x_{k+1})}{\nabla
2124	>	V(x_{k})^{T}\nabla V(x_{k})}%
2125	>	\end{equation}
2126
2127	<	\subsection{\label{sec:architecture} {\sc oopse} Architecture}
2127	>	The conjugate gradient method assumes that the conformation is close
2128	>	enough to a local minimum that the potential energy surface is very
2129	>	nearly quadratic.  When the initial structure is far from the minimum,
2130	>	the steepest descent method can be superior to the conjugate gradient
2131	>	method. Hence, the steepest descent method is often used for the first
2132	>	10-100 steps of minimization. Another useful feature of minimization
2133	>	methods in {\sc oopse} is that a modified {\sc shake} algorithm can be
2134	>	applied during the minimization to constraint the bond lengths if this
2135	>	is required by the force field. Meta-data parameters concerning the
2136	>	minimizer are given in Table~\ref{table:minimizeParams}
2137
2138	<	The core of OOPSE is divided into two main object libraries:
2139	<	\texttt{libBASS} and \texttt{libmdtools}. \texttt{libBASS} is the
2140	<	library developed around the parsing engine and \texttt{libmdtools}
2141	<	is the software library developed around the simulation engine. These
2142	<	two libraries are designed to encompass all the basic functions and
2143	<	tools that {\sc oopse} provides. Utility programs, such as the
2144	<	property analyzers, need only link against the software libraries to
2145	<	gain access to parsing, force evaluation, and input / output
2146	<	routines.
2138	>	\begin{table}
2139	>	\caption{The Global Keywords: Energy Minimizer Parameters}
2140	>	\label{table:minimizeParams}
2141	>	\begin{center}
2142	>	% Note when adding or removing columns, the \hsize numbers must add up to the total number
2143	>	% of columns.
2144	>	\begin{tabularx}{\linewidth}%
2145	>	{>{\setlength{\hsize}{1.2\hsize}}X%
2146	>	>{\setlength{\hsize}{0.6\hsize}}X%
2147	>	>{\setlength{\hsize}{1.1\hsize}}X%
2148	>	>{\setlength{\hsize}{1.1\hsize}}X}
2149
2150	<	Contained in \texttt{libBASS} are all the routines associated with
2114	<	reading and parsing the \texttt{.bass} input files. Given a
2115	<	\texttt{.bass} file, \texttt{libBASS} will open it and any associated
2116	<	\texttt{.mdl} files; then create structures in memory that are
2117	<	templates of all the molecules specified in the input files. In
2118	<	addition, any simulation parameters set in the \texttt{.bass} file
2119	<	will be placed in a structure for later query by the controlling
2120	<	program.
2150	>	{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline
2151
2152	<	Located in \texttt{libmdtools} are all other routines necessary to a
2153	<	Molecular Dynamics simulation. The library uses the main data
2154	<	structures returned by \texttt{libBASS} to initialize the various
2155	<	parts of the simulation: the atom structures and positions, the force
2156	<	field, the integrator, \emph{et cetera}. After initialization, the
2157	<	library can be used to perform a variety of tasks: integrate a
2158	<	Molecular Dynamics trajectory, query phase space information from a
2159	<	specific frame of a completed trajectory, or even recalculate force or
2160	<	energetic information about specific frames from a completed
2161	<	trajectory.
2152	>	{\tt minimizer} & string &  selects the minimization method to be used
2153	>	& either {\tt CG} (conjugate gradient) or {\tt SD} (steepest
2154	>	descent) \\
2155	>	{\tt minimizerMaxIter} & steps & Sets the maximum iteration number in the energy minimization & Default value is 200\\
2156	>	{\tt minimizerWriteFreq} & steps & Sets the frequency at which the {\tt .dump} and {\tt .stat} files are writtern during energy minimization & \\
2157	>	{\tt minimizerStepSize} & $\mbox{\AA}$ &  Set the step size of line search & Default value is 0.01\\
2158	>	{\tt minimizerFTol} & $\mbox{kcal mol}^{-1}$  & Sets energy tolerance  & Default value is $10^{-8}$\\
2159	>	{\tt minimizerGTol} & $\mbox{kcal mol}^{-1}\mbox{\AA}^{-1}$ & Sets gradient tolerance & Default value is $10^{-8}$\\
2160	>	{\tt minimizerLSTol} &  $\mbox{kcal mol}^{-1}$ & Sets line search tolerance & Default value is $10^{-8}$\\
2161	>	{\tt minimizerLSMaxIter} & steps &  Sets the maximum iteration of line searching & Default value is 50\\
2162
2163	<	With these core libraries in place, several programs have been
2164	<	developed to utilize the routines provided by \texttt{libBASS} and
2165	<	\texttt{libmdtools}. The main program of the package is \texttt{oopse}
2136	<	and the corresponding parallel version \texttt{oopse\_MPI}. These two
2137	<	programs will take the \texttt{.bass} file, and create (and integrate)
2138	<	the simulation specified in the script. The two analysis programs
2139	<	\texttt{staticProps} and \texttt{dynamicProps} utilize the core
2140	<	libraries to initialize and read in trajectories from previously
2141	<	completed simulations, in addition to the ability to use functionality
2142	<	from \texttt{libmdtools} to recalculate forces and energies at key
2143	<	frames in the trajectories. Lastly, the family of system building
2144	<	programs (Sec.~\ref{oopseSec:initCoords}) also use the libraries to
2145	<	store and output the system configurations they create.
2163	>	\end{tabularx}
2164	>	\end{center}
2165	>	\end{table}
2166
2167	<	\subsection{\label{oopseSec:parallelization} Parallelization of {\sc oopse}}
2167	>	\section{\label{oopseSec:parallelization} Parallel Simulation Implementation}
2168
2169	<	Although processor power is continually growing roughly following
2170	<	Moore's Law, it is still unreasonable to simulate systems of more then
2171	<	a 1000 atoms on a single processor. To facilitate study of larger
2172	<	system sizes or smaller systems on long time scales in a reasonable
2173	<	period of time, parallel methods were developed allowing multiple
2174	<	CPU's to share the simulation workload. Three general categories of
2175	<	parallel decomposition methods have been developed including atomic,
2176	<	spatial and force decomposition methods.
2169	>	Although processor power is continually improving, it is still
2170	>	unreasonable to simulate systems of more then a 1000 atoms on a single
2171	>	processor. To facilitate study of larger system sizes or smaller
2172	>	systems for longer time scales, parallel methods were developed to
2173	>	allow multiple CPU's to share the simulation workload. Three general
2174	>	categories of parallel decomposition methods have been developed:
2175	>	these are the atomic,\cite{Fox88} spatial~\cite{plimpton95} and
2176	>	force~\cite{Paradyn} decomposition methods.
2177
2178	<	Algorithmically simplest of the three methods is atomic decomposition
2179	<	where N particles in a simulation are split among P processors for the
2180	<	duration of the simulation. Computational cost scales as an optimal
2181	<	$\mathcal{O}(N/P)$ for atomic decomposition. Unfortunately all
2178	>	Algorithmically simplest of the three methods is atomic decomposition,
2179	>	where $N$ particles in a simulation are split among $P$ processors for
2180	>	the duration of the simulation. Computational cost scales as an
2181	>	optimal $\mathcal{O}(N/P)$ for atomic decomposition. Unfortunately all
2182		processors must communicate positions and forces with all other
2183	<	processors at every force evaluation, leading communication costs to
2184	<	scale as an unfavorable $\mathcal{O}(N)$, \emph{independent of the
2183	>	processors at every force evaluation, leading the communication costs
2184	>	to scale as an unfavorable $\mathcal{O}(N)$, \emph{independent of the
2185		number of processors}. This communication bottleneck led to the
2186	<	development of spatial and force decomposition methods in which
2186	>	development of spatial and force decomposition methods, in which
2187		communication among processors scales much more favorably. Spatial or
2188		domain decomposition divides the physical spatial domain into 3D boxes
2189		in which each processor is responsible for calculation of forces and
2190		positions of particles located in its box. Particles are reassigned to
2191		different processors as they move through simulation space. To
2192	<	calculate forces on a given particle, a processor must know the
2192	>	calculate forces on a given particle, a processor must simply know the
2193		positions of particles within some cutoff radius located on nearby
2194	<	processors instead of the positions of particles on all
2194	>	processors rather than the positions of particles on all
2195		processors. Both communication between processors and computation
2196		scale as $\mathcal{O}(N/P)$ in the spatial method. However, spatial
2197		decomposition adds algorithmic complexity to the simulation code and
2198	<	is not very efficient for small N since the overall communication
2198	>	is not very efficient for small $N$, since the overall communication
2199		scales as the surface to volume ratio $\mathcal{O}(N/P)^{2/3}$ in
2200		three dimensions.
2201
#	Line 2184 | Line 2204 | given row by particles located in that processors colu
2204		processors based on a block decomposition of the force
2205		matrix. Processors are split into an optimally square grid forming row
2206		and column processor groups. Forces are calculated on particles in a
2207	<	given row by particles located in that processors column
2207	>	given row by particles located in that processor's column
2208		assignment. Force decomposition is less complex to implement than the
2209		spatial method but still scales computationally as $\mathcal{O}(N/P)$
2210		and scales as $\mathcal{O}(N/\sqrt{P})$ in communication
#	Line 2193 | Line 2213 | atoms.\cite{plimpton95}
2213		and favorably compete with spatial methods up to 100,000
2214		atoms.\cite{plimpton95}
2215
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	–
2216		\section{\label{oopseSec:conclusion}Conclusion}
2217
2218	<	We have presented the design and implementation of our open source
2219	<	simulation package {\sc oopse}. The package offers novel capabilities
2220	<	to the field of Molecular Dynamics simulation packages in the form of
2221	<	dipolar force fields, and symplectic integration of rigid body
2222	<	dynamics. It is capable of scaling across multiple processors through
2223	<	the use of force based decomposition using MPI. It also implements
2224	<	several advanced integrators allowing the end user control over
2225	<	temperature and pressure. In addition, it is capable of integrating
2226	<	constrained dynamics through both the {\sc rattle} algorithm and the
2227	<	z-constraint method.
2218	>	We have presented a new open source parallel simulation program {\sc
2219	>	oopse}. This program offers some novel capabilities, but mostly makes
2220	>	available a library of modern object-oriented code for the scientific
2221	>	community to use freely.  Notably, {\sc oopse} can handle symplectic
2222	>	integration of objects (atoms and rigid bodies) which have
2223	>	orientational degrees of freedom.  It can also work with transition
2224	>	metal force fields and point-dipoles. It is capable of scaling across
2225	>	multiple processors through the use of force based decomposition. It
2226	>	also implements several advanced integrators allowing the end user
2227	>	control over temperature and pressure. In addition, it is capable of
2228	>	integrating constrained dynamics through both the {\sc rattle}
2229	>	algorithm and the $z$-constraint method.
2230
2231	<	These features are all brought together in a single open-source
2232	<	program. This allows researchers to not only benefit from
2233	<	{\sc oopse}, but also contribute to {\sc oopse}'s development as
2234	<	well.
2231	>	We encourage other researchers to download and apply this program to
2232	>	their own research problems.  By making the code available, we hope to
2233	>	encourage other researchers to contribute their own code and make it a
2234	>	more powerful package for everyone in the molecular dynamics community
2235	>	to use.  All source code for {\sc oopse} is available for download at
2236	>	{\tt http://oopse.org}.
2237
2254	–
2238		\newpage
2239		\section{Acknowledgments}
2257	–	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.
2240
2241	+	Development of {\sc oopse} was funded by a New Faculty Award from the
2242	+	Camille and Henry Dreyfus Foundation and by the National Science
2243	+	Foundation under grant CHE-0134881. Computation time was provided by
2244	+	the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant
2245	+	DMR-0079647.
2246	+
2247		\bibliographystyle{achemso}
2248		\bibliography{oopsePaper}
2249

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines