--- trunk/oopsePaper/EmpericalEnergy.tex 2004/01/12 18:43:56 925 +++ trunk/oopsePaper/EmpericalEnergy.tex 2004/01/13 20:03:21 933 @@ -9,7 +9,7 @@ are not currently suporrted by {\sc oopse}. element, or be used for collections of atoms such as a methyl group. The atoms are also capable of having a directional component associated with them, often in the form of a dipole. Charges on atoms -are not currently suporrted by {\sc oopse}. +are not currently suported by {\sc oopse}. The second most basic building block of a simulation is the molecule. The molecule is a way for {\sc oopse} to keep track of the @@ -68,12 +68,9 @@ potential. The Lennard-Jones potential mimics the attr \subsection{\label{sec:LJPot}The Lennard Jones Potential} The most basic force field implemented in OOPSE is the Lennard-Jones -potential. The Lennard-Jones potential mimics the attractive forces -two charge neutral particles experience when spontaneous dipoles are -induced on each other. This is the predominant intermolecular force in -systems of such as noble gases and simple alkanes. - -The Lennard-Jones potential is given by: +potential. The Lennard-Jones potential. Which mimics the Van der Waals +interaction at long distances, and uses an emperical repulsion at +short distances. The Lennard-Jones potential is given by: \begin{equation} V_{\text{LJ}}(r_{ij}) = 4\epsilon_{ij} \biggl[ @@ -82,20 +79,20 @@ Where $r_ij$ is the distance between particle $i$ and \biggr] \label{eq:lennardJonesPot} \end{equation} -Where $r_ij$ is the distance between particle $i$ and $j$, $\sigma_{ij}$ -scales the length of the interaction, and $\epsilon_{ij}$ scales the -energy well depth of the potential. +Where $r_{ij}$ is the distance between particle $i$ and $j$, +$\sigma_{ij}$ scales the length of the interaction, and +$\epsilon_{ij}$ scales the well depth of the potential. -Because the potential is calculated between all pairs, the force +Because this potential is calculated between all pairs, the force evaluation can become computationally expensive for large systems. To -keep the pair evaluation to a manegable number, OOPSE employs the use -of a cut-off radius.\cite{allen87:csl} The cutoff radius is set to be -$2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest self self length +keep the pair evaluation to a manegable number, OOPSE employs a +cut-off radius.\cite{allen87:csl} The cutoff radius is set to be +$2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest Lennard-Jones length parameter in the system. Truncating the calculation at $r_{\text{cut}}$ introduces a discontinuity into the potential energy. To offset this discontinuity, the energy value at $r_{\text{cut}}$ is subtracted from the entire potential. This causes -the equation to go to zero at the cut-off radius. +the potential to go to zero at the cut-off radius. Interactions between dissimilar particles requires the generation of cross term parameters for $\sigma$ and $\epsilon$. These are @@ -114,51 +111,40 @@ The \underline{D}ipolar \underline{U}nified-Atom \subsection{\label{sec:DUFF}Dipolar Unified-Atom Force Field} -The \underline{D}ipolar \underline{U}nified-Atom -\underline{F}orce \underline{F}ield ({\sc duff}) was developed to -simulate lipid bilayers. We needed a model capable of forming +The Dipolar Unified-atom Force Field ({\sc duff}) was developed to +simulate lipid bilayers. The systems require a model capable of forming bilayers, while still being sufficiently computationally efficient to allow simulations of large systems ($\approx$100's of phospholipids, $\approx$1000's of waters) for long times ($\approx$10's of nanoseconds). -With this goal in mind, we have eliminated all point charges. Charge -distributions were replaced with dipoles, and charge-neutral -distributions were reduced to Lennard-Jones interaction sites. This +With this goal in mind, {\sc duff} has no point charges. Charge +neutral distributions were replaced with dipoles, while most atoms and +groups of atoms were reduced to Lennard-Jones interaction sites. This simplification cuts the length scale of long range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, allowing us to avoid the -computationally expensive Ewald-Sum. Instead, we can use -neighbor-lists and cutoff radii for the dipolar interactions. +computationally expensive Ewald sum. Instead, we can use +neighbor-lists, reaction field, and cutoff radii for the dipolar +interactions. -As an example, lipid head groups in {\sc duff} are represented as point -dipole interaction sites.PC and PE Lipid head groups are typically -zwitterionic in nature, with charges separated by as much as -6~$\mbox{\AA}$. By placing a dipole of 20.6~Debye at the head group -center of mass, our model mimics the head group of PC.\cite{Cevc87} -Additionally, a Lennard-Jones site is located at the -pseudoatom's center of mass. The model is illustrated by the dark grey -atom in Fig.~\ref{fig:lipidModel}. +As an example, lipid head-groups in {\sc duff} are represented as +point dipole interaction sites. By placing a dipole of 20.6~Debye at +the head group center of mass, our model mimics the head group of +phosphatidylcholine.\cite{Cevc87} Additionally, a Lennard-Jones site +is located at the pseudoatom's center of mass. The model is +illustrated by the dark grey atom in Fig.~\ref{fig:lipidModel}. The water model we use to complement the dipoles of the lipids is out +repaarameterization of the soft sticky dipole (SSD) model of Ichiye +\emph{et al.}\cite{liu96:new_model} \begin{figure} +\epsfxsize=\linewidth \epsfbox{lipidModel.eps} \caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ % -is the bend angle, $\mu$ is the dipole moment of the head group, and n is the chain length.} +is the bend angle, $\mu$ is the dipole moment of the head group, and n +is the chain length.} \label{fig:lipidModel} \end{figure} -The water model we use to complement the dipoles of the lipids is -the soft sticky dipole (SSD) model of Ichiye \emph{et -al.}\cite{liu96:new_model} This model is discussed in greater detail -in Sec.~\ref{sec:SSD}. In all cases we reduce water to a single -Lennard-Jones interaction site. The site also contains a dipole to -mimic the partial charges on the hydrogens and the oxygen. However, -what makes the SSD model unique is the inclusion of a tetrahedral -short range potential to recover the hydrogen bonding of water, an -important factor when modeling bilayers, as it has been shown that -hydrogen bond network formation is a leading contribution to the -entropic driving force towards lipid bilayer formation.\cite{Cevc87} - - Turning to the tails of the phospholipids, we have used a set of scalable parameters to model the alkyl groups with Lennard-Jones sites. For this, we have used the TraPPE force field of Siepmann @@ -184,14 +170,14 @@ V_{\text{Total}} = \sum^{N}_{I=1} V^{I}_{\text{Interna The total energy of function in {\sc duff} is given by the following: \begin{equation} V_{\text{Total}} = \sum^{N}_{I=1} V^{I}_{\text{Internal}} - + \sum^{N}_{I=1} \sum^{N}_{J=1} V^{IJ}_{\text{Cross}} + + \sum^{N}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}} \label{eq:totalPotential} \end{equation} Where $V^{I}_{\text{Internal}}$ is the internal potential of a molecule: \begin{equation} V^{I}_{\text{Internal}} = \sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk}) - + \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\theta_{ijkl}) + + \sum_{\phi_{ijkl} \in I} V_{\text{torsion}}(\phi_{ijkl}) + \sum_{i \in I} \sum_{(j>i+4) \in I} \biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}} (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) @@ -199,26 +185,12 @@ within in the molecule. $V_{\text{torsion}}$ is the to \label{eq:internalPotential} \end{equation} Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs -within in the molecule. $V_{\text{torsion}}$ is the torsion The -pairwise portions of the internal potential are excluded for pairs -that ar closer than three bonds, i.e.~atom pairs farther away than a -torsion are included in the pair-wise loop. +within the molecule, and $V_{\text{torsion}}$ is the torsion potential +for all 1, 4 bonded pairs. The pairwise portions of the internal +potential are excluded for pairs that are closer than three bonds, +i.e.~atom pairs farther away than a torsion are included in the +pair-wise loop. -The cross portion of the total potential, $V^{IJ}_{\text{Cross}}$, is -as follows: -\begin{equation} -V^{IJ}_{\text{Cross}} = - \sum_{i \in I} \sum_{j \in J} - \biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}} - (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) - + V_{\text{sticky}} - (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) - \biggr] -\label{eq:crossPotentail} -\end{equation} -Where $V_{\text{LJ}}$ is the Lennard Jones potential, -$V_{\text{dipole}}$ is the dipole dipole potential, and -$V_{\text{sticky}}$ is the sticky potential defined by the SSD model. The bend potential of a molecule is represented by the following function: \begin{equation} @@ -233,16 +205,16 @@ V_{\text{torsion}}(\phi_{ijkl}) = c_1[1 + \cos \phi] The torsion potential and parameters are also taken from TraPPE. It is of the form: \begin{equation} -V_{\text{torsion}}(\phi_{ijkl}) = c_1[1 + \cos \phi] +V_{\text{torsion}}(\phi) = c_1[1 + \cos \phi] + c_2[1 + \cos(2\phi)] + c_3[1 + \cos(3\phi)] \label{eq:origTorsionPot} \end{equation} Here $\phi_{ijkl}$ is the angle defined by four bonded neighbors $i$, -$j$, $k$, and $l$ (again, see Fig.~\ref{fig:lipidModel}). However, -for computaional efficency, the torsion potentail has been recast -after the method of CHARMM\cite{charmm1983} whereby the angle series -is converted to a power series of the form: +$j$, $k$, and $l$ (again, see Fig.~\ref{fig:lipidModel}). For +computaional efficency, the torsion potential has been recast after +the method of CHARMM\cite{charmm1983} whereby the angle series is +converted to a power series of the form: \begin{equation} V_{\text{torsion}}(\phi_{ijkl}) = k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0 @@ -259,23 +231,40 @@ evaluations are avoided during the calculation of the evaluations are avoided during the calculation of the potential. +The cross portion of the total potential, $V^{IJ}_{\text{Cross}}$, is +as follows: +\begin{equation} +V^{IJ}_{\text{Cross}} = + \sum_{i \in I} \sum_{j \in J} + \biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}} + (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) + + V_{\text{sticky}} + (\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) + \biggr] +\label{eq:crossPotentail} +\end{equation} +Where $V_{\text{LJ}}$ is the Lennard Jones potential, +$V_{\text{dipole}}$ is the dipole dipole potential, and +$V_{\text{sticky}}$ is the sticky potential defined by the SSD +model. Note that not all atom types include all interactions. The dipole-dipole potential has the following form: \begin{equation} V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i}, - \boldsymbol{\Omega}_{j}) = \frac{1}{4\pi\epsilon_{0}} \biggl[ - \frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}} + \boldsymbol{\Omega}_{j}) = \frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[ + \boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j} - - \frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) % - (\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) } - {r^{5}_{ij}} \biggr] + \frac{3(\boldsymbol{\hat{u}}_i \cdot \mathbf{r}_{ij}) % + (\boldsymbol{\hat{u}}_j \cdot \mathbf{r}_{ij}) } + {r^{2}_{ij}} \biggr] \label{eq:dipolePot} \end{equation} Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ -are the Euler angles of atom $i$ and $j$ -respectively. $\boldsymbol{\mu}_i$ is the dipole vector of atom -$i$ it takes its orientation from $\boldsymbol{\Omega}_i$. +are the Euler angles of atom $i$ and $j$ respectively. $|\mu_i|$ is +the magnitude of the dipole moment of atom $i$ and +$\boldsymbol{\hat{u}}_i$ is the standard unit orientation vector of +$\boldsymbol{\Omega}_i$. \subsection{\label{sec:SSD}The {\sc DUFF} Water Models: SSD/E and SSD/RF} @@ -367,22 +356,22 @@ models can be found in reference \ref{Gezelter04}. !!!Place a {\sc BASS} scheme showing SSD parameters around here!!! -\subsection{\label{sec:eam}Embedded Atom Model} +\subsection{\label{sec:eam}Embedded Atom Method} -Several molecular dynamics codes\cite{dynamo86} exist which have the +Several other molecular dynamics packages\cite{dynamo86} exist which have the capacity to simulate metallic systems, including some that have parallel computational abilities\cite{plimpton93}. Potentials that describe bonding transition metal systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have a -attractive interaction which models the stabilization of ``Embedding'' -a positively charged metal ion in an electron density created by the +attractive interaction which models ``Embedding'' +a positively charged metal ion in the electron density due to the free valance ``sea'' of electrons created by the surrounding atoms in the system. A mostly repulsive pairwise part of the potential describes the interaction of the positively charged metal core ions with one another. A particular potential description called the -Embedded Atom Method\cite{Daw84,FBD86,johnson89,Lu97}(EAM) that has +Embedded Atom Method\cite{Daw84,FBD86,johnson89,Lu97}({\sc eam}) that has particularly wide adoption has been selected for inclusion in OOPSE. A -good review of EAM and other metallic potential formulations was done +good review of {\sc eam} and other metallic potential formulations was done by Voter.\cite{voter} The {\sc eam} potential has the form: @@ -390,73 +379,78 @@ V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i} \phi_{ij}({\bf r}_{ij}) \\ \rho_{i} & = & \sum_{j \neq i} f_{j}({\bf r}_{ij}) -\end{eqnarray} +\end{eqnarray}S -where $\phi_{ij}$ is a primarily repulsive pairwise interaction -between atoms $i$ and $j$.In the origional formulation of -EAM\cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term, however -in later refinements to EAM have shown that nonuniqueness between $F$ -and $\phi$ allow for more general forms for $\phi$.\cite{Daw89} The -embedding function $F_{i}$ is the energy required to embedded an -positively-charged core ion $i$ into a linear supeposition of +where $F_{i} $ is the embedding function that equates the energy required to embed a +positively-charged core ion $i$ into a linear superposition of sperically averaged atomic electron densities given by -$\rho_{i}$. There is a cutoff distance, $r_{cut}$, which limits the +$\rho_{i}$. $\phi_{ij}$ is a primarily repulsive pairwise interaction +between atoms $i$ and $j$. In the original formulation of +{\sc eam} cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term, however +in later refinements to EAM have shown that non-uniqueness between $F$ +and $\phi$ allow for more general forms for $\phi$.\cite{Daw89} + There is a cutoff distance, $r_{cut}$, which limits the summations in the {\sc eam} equation to the few dozen atoms surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$ -interactions. +interactions. Foiles et al. fit EAM potentials for fcc metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals\cite{FDB86}. These potential fits are in the DYNAMO 86 format and are included with {\sc oopse}. + \subsection{\label{Sec:pbc}Periodic Boundary Conditions} -\textit{Periodic boundary conditions} are widely used to simulate truly -macroscopic systems with a relatively small number of particles. Simulation -box is replicated throughout space to form an infinite lattice. During the -simulation, when a particle moves in the primary cell, its periodic image -particles in other boxes move in exactly the same direction with exactly the -same orientation.Thus, as a particle leaves the primary cell, one of its -images will enter through the opposite face.If the simulation box is large -enough to avoid "feeling" the symmetric of the periodic lattice,the surface -effect could be ignored. Cubic and parallelepiped are the available periodic -cells. \bigskip In OOPSE, we use the matrix instead of the vector to describe -the property of the simulation box. Therefore, not only the size of the -simulation box could be changed during the simulation, but also the shape of -it. The transformation from box space vector $\overrightarrow{s}$ to its -corresponding real space vector $\overrightarrow{r}$ is defined by -\begin{equation} -\overrightarrow{r}=H\overrightarrow{s}% -\end{equation} - - -where $H=(h_{x},h_{y},h_{z})$ is a transformation matrix made up of the three -box axis vectors. $h_{x},h_{y}$ and $h_{z}$ represent the sides of the -simulation box respectively. - -To find the minimum image, we need to convert the real vector to its -corresponding vector in box space first, \bigskip% -\begin{equation} -\overrightarrow{s}=H^{-1}\overrightarrow{r}% -\end{equation} -And then, each element of $\overrightarrow{s}$ is casted to lie between -0.5 -to 0.5, -\begin{equation} -s_{i}^{\prime}=s_{i}-round(s_{i}) -\end{equation} -where% - -\begin{equation} -round(x)=\lfloor{x+0.5}\rfloor\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }if\text{ -}x\geqslant0 -\end{equation} -% - -\begin{equation} -round(x)=\lceil{x-0.5}\rceil\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ }if\text{ }x<0 -\end{equation} - - -For example, $round(3.6)=4$,$round(3.1)=3$, $round(-3.6)=-4$, $round(-3.1)=-3$. - -Finally, we could get the minimum image by transforming back to real space,% - -\begin{equation} -\overrightarrow{r^{\prime}}=H^{-1}\overrightarrow{s^{\prime}}% -\end{equation} +\textit{Periodic boundary conditions} are widely used to simulate truly +macroscopic systems with a relatively small number of particles. The +simulation box is replicated throughout space to form an infinite +lattice. During the simulation, when a particle moves in the primary +cell, its image in other boxes move in exactly the same direction with +exactly the same orientation.Thus, as a particle leaves the primary +cell, one of its images will enter through the opposite face.If the +simulation box is large enough to avoid "feeling" the symmetries of +the periodic lattice, surface effects can be ignored. Cubic, +orthorhombic and parallelepiped are the available periodic cells In +OOPSE. We use a matrix to describe the property of the simulation +box. Therefore, both the size and shape of the simulation box can be +changed during the simulation. The transformation from box space +vector $\mathbf{s}$ to its corresponding real space vector +$\mathbf{r}$ is defined by +\begin{equation} +\mathbf{r}=\underline{\underline{H}}\cdot\mathbf{s}% +\end{equation} + + +where $H=(h_{x},h_{y},h_{z})$ is a transformation matrix made up of +the three box axis vectors. $h_{x},h_{y}$ and $h_{z}$ represent the +three sides of the simulation box respectively. + +To find the minimum image, we convert the real vector to its +corresponding vector in box space first, \bigskip% +\begin{equation} +\mathbf{s}=\underline{\underline{H}}^{-1}\cdot\mathbf{r}% +\end{equation} +And then, each element of $\mathbf{s}$ is wrapped to lie between -0.5 to 0.5, +\begin{equation} +s_{i}^{\prime}=s_{i}-round(s_{i}) +\end{equation} +where + +% + +\begin{equation} +round(x)=\left\{ +\begin{array}[c]{c}% +\lfloor{x+0.5}\rfloor & \text{if \ }x\geqslant0\\ +\lceil{x-0.5}\rceil & \text{otherwise}% +\end{array} +\right. +\end{equation} + + +For example, $round(3.6)=4$,$round(3.1)=3$, $round(-3.6)=-4$, +$round(-3.1)=-3$. + +Finally, we obtain the minimum image coordinates by transforming back +to real space,% + +\begin{equation} +\mathbf{r}^{\prime}=\underline{\underline{H}}^{-1}\cdot\mathbf{s}^{\prime}% +\end{equation} +