| 190 |
|
dipoles). Charges, permanent dipoles, and Lennard-Jones parameters for |
| 191 |
|
a given atom type are set in the force field parameter files. |
| 192 |
|
|
| 193 |
– |
\begin{lstlisting}[float,caption={[Specifier for molecules and atoms] A sample specification of an Ar molecule},label=sch:AtmMole] |
| 194 |
– |
molecule{ |
| 195 |
– |
name = "Ar"; |
| 196 |
– |
nAtoms = 1; |
| 197 |
– |
atom[0]{ |
| 198 |
– |
type="Ar"; |
| 199 |
– |
position( 0.0, 0.0, 0.0 ); |
| 200 |
– |
} |
| 201 |
– |
} |
| 202 |
– |
\end{lstlisting} |
| 203 |
– |
|
| 204 |
– |
|
| 193 |
|
Atoms can be collected into secondary structures such as rigid bodies |
| 194 |
|
or molecules. The molecule is a way for {\sc oopse} to keep track of |
| 195 |
|
the atoms in a simulation in logical manner. Molecular units store the |
| 196 |
|
identities of all the atoms and rigid bodies associated with |
| 197 |
|
themselves, and are responsible for the evaluation of their own |
| 198 |
|
internal interactions (\emph{i.e.}~bonds, bends, and torsions). Scheme |
| 199 |
< |
\ref{sch:AtmMole} shows how one creates a molecule in a ``model'' or |
| 199 |
> |
\ref{sch:mdlExample} shows how one creates a molecule in a ``model'' or |
| 200 |
|
\texttt{.mdl} file. The position of the atoms given in the |
| 201 |
|
declaration are relative to the origin of the molecule, and is used |
| 202 |
|
when creating a system containing the molecule. |
| 209 |
|
included in most simulation packages because of the algorithmic |
| 210 |
|
complexity involved in propagating orientational degrees of |
| 211 |
|
freedom. Until recently, integrators which propagate orientational |
| 212 |
< |
motion have been much worse than those available for translational |
| 212 |
> |
motion have had energy conservation problems when compared to those available for translational |
| 213 |
|
motion. |
| 214 |
|
|
| 215 |
|
Moving a rigid body involves determination of both the force and |
| 251 |
|
entire nine parameter rotation matrix. Further discussion |
| 252 |
|
on this choice can be found in Sec.~\ref{oopseSec:integrate}. An example |
| 253 |
|
definition of a rigid body can be seen in Scheme |
| 254 |
< |
\ref{sch:rigidBody}. The positions in the atom definitions are the |
| 267 |
< |
placements of the atoms relative to the origin of the rigid body, |
| 268 |
< |
which itself has a position relative to the origin of the molecule. |
| 254 |
> |
\ref{sch:rigidBody}. |
| 255 |
|
|
| 256 |
< |
\begin{lstlisting}[float,caption={[Defining rigid bodies]A sample definition of a rigid body},label={sch:rigidBody}] |
| 256 |
> |
\begin{lstlisting}[float,caption={[Defining rigid bodies]A sample definition of a molecule containing a rigid body},label={sch:rigidBody}] |
| 257 |
|
molecule{ |
| 258 |
|
name = "TIP3P"; |
| 259 |
|
nAtoms = 3; |
| 282 |
|
|
| 283 |
|
The actual creation of a {\sc bass} script requires several key components. The first part of the script needs to be the declaration of all of the molecule prototypes used in the simulation. This is typically done through the inclusion of {\tt .mdl} files. Only the molecules actually present in the simulation need to be declared, however {\sc bass} allows for the declaration of more molecules than are needed. This gives the user the ability to build up a library of commonly used molecules into a single {\tt .mdl} file. |
| 284 |
|
|
| 285 |
< |
Once all prototypes are declared, the ordering of the rest of the script is less stringent. Typically, the next to follow the molecular prototypes are the component statements. These statements specify which molecules are present within the simulation. The number of components must first be declared before the first component block statement (an example is seen in Sch.~\ref{sch:bassExPrime}). |
| 285 |
> |
Once all prototypes are declared, the ordering of the rest of the script is less stringent. Typically, the next to follow the molecular prototypes are the component statements. These statements specify which molecules are present within the simulation. The number of components must first be declared before the first component block statement (an example is seen in Sch.~\ref{sch:bassExPrime}). The component blocks tell {\sc oopse} the number of molecules that will be in the simulation, and the order in which the components blocks are declared sets the ordering of the real atoms within the simulation as well as in the output files. |
| 286 |
> |
|
| 287 |
> |
The remainder of the script then sets the various simulation parameters for the system of interest. The required set of parameters that must be present in all simulations is given in Table~\ref{table:reqParams}. The {\tt ensemble} statement is responsible for selecting the integration method used for the calculation of the equations of motion. An in depth discussion of the various methods available in {\sc oopse} can be found in Sec.~\ref{oopseSec:mechanics}. The {\tt forceField} statement is important for the selection of which forces will be used in the course of the simulation. {\sc oopse} supports several force fields, as outlined in Sec.~\ref{oopseSec:empericalEnergy}. The force fields are interchangeable between simulations, with the only requirement being that all atoms needed by the simulation are defined within the selected force field. The time step between force evaluations is set with the {\tt dt} parameter, and {\tt runTime} will set the time length of the simulation. Note, that {\tt runTime} is an absolute time, meaning if the simulation is started at t = 10.0~ns with a {\tt runTime} of 25.0~ns, the simulation will only run for an additional 15.0~ns. The final required parameter, is the {\tt initialConfig} statement. This will set the initial coordinates for the system, as well as the initial time if the {\tt useInitalTime = true;} flag is given. The format of the file specified in {\tt initialConfig}, is given in Sec.~\ref{oopseSec:coordFiles}. Additional parameters are summarized in Table~\ref{table:genParams}. |
| 288 |
> |
|
| 289 |
> |
\begin{table} |
| 290 |
> |
\caption{The Global Keywords: Required Parameters} |
| 291 |
> |
\label{table:reqParams} |
| 292 |
> |
\begin{center} |
| 293 |
> |
% Note when adding or removing columns, the \hsize numbers must add up to the total number |
| 294 |
> |
% of columns. |
| 295 |
> |
\begin{tabularx}{\linewidth}% |
| 296 |
> |
{>{\setlength{\hsize}{1.00\hsize}}X% |
| 297 |
> |
>{\setlength{\hsize}{0.4\hsize}}X% |
| 298 |
> |
>{\setlength{\hsize}{1.2\hsize}}X% |
| 299 |
> |
>{\setlength{\hsize}{1.4\hsize}}X} |
| 300 |
|
|
| 301 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 302 |
+ |
|
| 303 |
+ |
{\tt forceField} & string & Sets the force field. & Possible force fields are "DUFF", "LJ", and "EAM". \\ |
| 304 |
+ |
{\tt ensemble} & string & Sets the ensemble. & Possible ensembles are "NVE", "NVT", "NPTi", "NPTf", and "NPTxyz".\\ |
| 305 |
+ |
{\tt dt} & fs & Sets the time step. & Selection of {\tt dt} should be small enough to sample the fastest motion of the simulation. \\ |
| 306 |
+ |
{\tt nComponents} & integer & Sets the number of components. & Needs to appear before the first {\tt Component} block. \\ |
| 307 |
+ |
{\tt initialConfig} & string & Sets the file containing the initial configuration. & Can point to any file containing the configuration in the correct order. \\ |
| 308 |
+ |
{\tt runTime} & fs & Sets the time at which the simulation should end. & This is an absolute time, and will end the simulation when the current time meets or exceeds the {\tt runTime}. \\ |
| 309 |
+ |
|
| 310 |
+ |
|
| 311 |
+ |
\end{tabularx} |
| 312 |
+ |
\end{center} |
| 313 |
+ |
\end{table} |
| 314 |
+ |
|
| 315 |
+ |
\begin{table} |
| 316 |
+ |
\caption{The Global Keywords: General Parameters} |
| 317 |
+ |
\label{table:genParams} |
| 318 |
+ |
\begin{center} |
| 319 |
+ |
% Note when adding or removing columns, the \hsize numbers must add up to the total number |
| 320 |
+ |
% of columns. |
| 321 |
+ |
\begin{tabularx}{\linewidth}% |
| 322 |
+ |
{>{\setlength{\hsize}{1.00\hsize}}X% |
| 323 |
+ |
>{\setlength{\hsize}{0.4\hsize}}X% |
| 324 |
+ |
>{\setlength{\hsize}{1.2\hsize}}X% |
| 325 |
+ |
>{\setlength{\hsize}{1.4\hsize}}X} |
| 326 |
+ |
|
| 327 |
+ |
{\bf keyword} & {\bf units} & {\bf use} & {\bf remarks} \\ \hline |
| 328 |
+ |
|
| 329 |
+ |
{\tt finalConfig} & string & Option to set the name of the final output file. & Useful when stringing simulations together. Defaults to the {\tt .bass} file with an {\tt .eor} extension. \\ |
| 330 |
+ |
{\tt useInitialTime} & logical & Sets whether the initial time is taken from the {\tt .init} file. & Useful when recovering a simulation from a crashed processor. Default is false. \\ |
| 331 |
+ |
{\tt sampleTime} & fs & Sets the frequency at which the {\tt .dump} file is written. & Default sets the frequency to the {\tt runTime}. \\ |
| 332 |
+ |
{\tt statusTime} & fs & Sets the frequency at which the {\tt .stat} file is written. & Defaults sets the frequency to the {\tt sampleTime}. \\ |
| 333 |
+ |
{\tt LJrcut} & $\mbox{\AA}$ & Manually sets the Lennard-Jones cutoff. & Defaults to $2.5\sigma_L$, where $\sigma_L$ is the largest LJ $\sigma$ in the simulation. \\ |
| 334 |
+ |
{\tt electrostaticCutoffRadius}& & & \\ |
| 335 |
+ |
& $\mbox{\AA}$ & Manually sets the cutoff used by the electrostatic potentials. & Defaults to $15\mbox{\AA}$ \\ |
| 336 |
+ |
{\tt electrostaticSkinThickness} & & & \\ |
| 337 |
+ |
& $\mbox{\AA}$ & Manually sets the skin thickness for the electrostatic switching function. & Defaults to 5~\% of the {\tt electrostaticSkinThickness}. \\ |
| 338 |
+ |
{\tt useReactionField} & logical & Turns the reaction field correction on/off. & Default is "false". \\ |
| 339 |
+ |
{\tt dielectric} & unitless & Sets the dielectric constant for reaction field. & If {\tt useReactionField} is true, then {\tt dielectric} must be set. \\ |
| 340 |
+ |
{\tt usePeriodicBoundaryConditions} & & & \\ |
| 341 |
+ |
& logical & Turns periodic boundary conditions on/off. & Default is "true". \\ |
| 342 |
+ |
{\tt seed } & integer & Sets the seed value for the random number generator. & The seed needs to be at least 9 digits long. The default is to take the seed from the CPU clock. |
| 343 |
+ |
|
| 344 |
+ |
\end{tabularx} |
| 345 |
+ |
\end{center} |
| 346 |
+ |
\end{table} |
| 347 |
+ |
|
| 348 |
+ |
|
| 349 |
+ |
|
| 350 |
|
\subsection{\label{oopseSec:coordFiles}Coordinate Files} |
| 351 |
|
|
| 352 |
|
The standard format for storage of a systems coordinates is a modified |
| 561 |
|
charges. Charge-neutral distributions were replaced with dipoles, |
| 562 |
|
while most atoms and groups of atoms were reduced to Lennard-Jones |
| 563 |
|
interaction sites. This simplification cuts the length scale of long |
| 564 |
< |
range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, and allows |
| 516 |
< |
us to avoid the computationally expensive Ewald sum. Instead, we can |
| 517 |
< |
use neighbor-lists and cutoff radii for the dipolar interactions, or |
| 518 |
< |
include a reaction field to mimic larger range interactions. |
| 564 |
> |
range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, removing the need for the computationally expensive Ewald sum. Instead, we Verlet neighbor-lists and cutoff radii are used for the dipolar interactions, or a reaction field is added to mimic longer range interactions. |
| 565 |
|
|
| 566 |
|
As an example, lipid head-groups in {\sc duff} are represented as |
| 567 |
|
point dipole interaction sites. By placing a dipole at the head |
| 571 |
|
site is located at the pseudoatom's center of mass. The model is |
| 572 |
|
illustrated by the red atom in Fig.~\ref{oopseFig:lipidModel}. The |
| 573 |
|
water model we use to complement the dipoles of the lipids is our |
| 574 |
< |
reparameterization of the soft sticky dipole (SSD) model of Ichiye |
| 574 |
> |
reparameterization\cite{fennell04} of the soft sticky dipole (SSD) model of Ichiye |
| 575 |
|
\emph{et al.}\cite{liu96:new_model} |
| 576 |
|
|
| 577 |
|
\begin{figure} |
| 594 |
|
$\text{CH}_2$ do not change depending on what species are bonded to |
| 595 |
|
it. |
| 596 |
|
|
| 597 |
< |
TraPPE also constrains all bonds to be of fixed length. Typically, |
| 597 |
> |
TraPPE and {\sc duff} also constrain all bonds to be of fixed length. Typically, |
| 598 |
|
bond vibrations are the fastest motions in a molecular dynamic |
| 599 |
|
simulation. Small time steps between force evaluations must be used to |
| 600 |
|
ensure adequate energy conservation in the bond degrees of freedom. By |
| 624 |
|
|
| 625 |
|
\end{lstlisting} |
| 626 |
|
|
| 627 |
< |
\subsection{\label{oopseSec:energyFunctions}{\sc duff} Energy Functions} |
| 627 |
> |
\subsubsection{\label{oopseSec:energyFunctions}{\sc duff} Energy Functions} |
| 628 |
|
|
| 629 |
|
The total potential energy function in {\sc duff} is |
| 630 |
|
\begin{equation} |
| 753 |
|
is the taper radius some given thickness less than the electrostatic |
| 754 |
|
cutoff. The switching thickness can be set in the \texttt{.bass} file. |
| 755 |
|
|
| 756 |
< |
\subsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF} |
| 756 |
> |
\subsubsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF} |
| 757 |
|
|
| 758 |
|
In the interest of computational efficiency, the default solvent used |
| 759 |
|
by {\sc oopse} is the extended Soft Sticky Dipole (SSD/E) water |
| 871 |
|
|
| 872 |
|
\subsection{\label{oopseSec:eam}Embedded Atom Method} |
| 873 |
|
|
| 874 |
< |
There are Molecular Dynamics packages which have the |
| 875 |
< |
capacity to simulate metallic systems, including some that have |
| 876 |
< |
parallel computational abilities\cite{plimpton93}. Potentials that |
| 831 |
< |
describe bonding transition metal |
| 832 |
< |
systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have an |
| 833 |
< |
attractive interaction which models ``Embedding'' |
| 874 |
> |
{\sc oopse} implements a potential that |
| 875 |
> |
describes bonding transition metal |
| 876 |
> |
systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} and has attractive interaction which models ``Embedding'' |
| 877 |
|
a positively charged metal ion in the electron density due to the |
| 878 |
|
free valance ``sea'' of electrons created by the surrounding atoms in |
| 879 |
|
the system. A mostly-repulsive pairwise part of the potential |
