trunk/chuckDissertation/introduction.tex



%!TEX root = /Users/charles/Documents/chuckDissertation/dissertation.tex
\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND}

\section{\label{introSec:introintro}Introduction}

This dissertation presents research that I have performed and been involved with over the course of my studies at the University of Notre Dame. I have chosen to present this material in a chronological manner since the research naturally lends itself to this style. Introductory material relating to Molecular Dynamics techniques is presented in this opening chapter. The primary focus of this dissertation is the study of the structure and dynamics of phases exhibited by metallic systems, particularly metallic glass-formers and nanoparticles (NPs). Molecular Dynamics augments experimental information (both structural information and particle dynamics) about these systems and allows the probing of the fundamental nature of phase transitions in these materials.

Chapter \ref{chap:metalglass} explores transport dynamics in a known glass former (the eutectic mixture of silver and copper). $\mathrm{Ag}_6\mathrm{Cu}_4$ presents an interesting system to study because it closely resembles binary Lennard-Jones systems that have correlation functions that decay according to the Kohlrausch-Williams-Watts ({\sc kww}) law.\cite{Kohlrausch:1863zv,Williams:1970fk} This law is general such that it can be applied to not only relaxation in glasses, but also in NPs and will be explored in several chapters in this dissertation. In addition, the observation that dynamics in glassy materials can be described by a fractal waiting time distribution known as the Continuous Time Random Walk ({\sc CTRW}) will be investigated. 

Spontaneous alloying of bimetallic Au-Ag NPs after synthesis is the focus of study in Chapter \ref{chap:nanodiffusion}. Metallic nanoparticles exhibit many physical and chemical properties that differ significantly from their bulk counterparts. These unique properties may be partly due to the enormous surface area to volume ratio of NPs which leads to excess surface free energy. If the excess surface free energy is comparable in magnitude to the lattice or binding energies of the NPs, structural instability can result. 
Computational techniques will be used to explore whether the hypothesis that a small fraction of vacancies formed at the Au-Ag core-shell interface during synthesis can result in the alloying of the particle is consistent with experimental observations.

Chapter \ref{chap:bulkmod} reports on simulations of the transient
response of metallic nanoparticles to the nearly instantaneous heating
undergone when photons are absorbed during ultrafast laser excitation
experiments. The time scale for heating (determined by the electron-phonon
coupling constant) is faster than a single period of the breathing
mode for spherical nanoparticles. Hot-electron pressure and direct lattice heating contribute to the thermal excitation of the atomic degrees of freedom, and both
mechanisms are rapid enough to coherently excite the breathing mode of
the spherical particles. Molecular Dynamics is used to augment experimental data by replicating the laser-excitation event to probe the dynamics of the NPs after excitation.

Chapter \ref{chap:nanoglass} combines observations from the previous chapters to determine if it is plausible to construct a glassy nanobead from the $\mathrm{Ag}_6\mathrm{Cu}_4$ system. One of the side effects of absorption of laser radiation at these frequencies is the rapid (sub-picosecond) heating of the electronic degrees of freedom in the metal. Since these experiments are carried out in condensed phase surroundings, the large surface area to volume ratio makes the heat transfer to the surrounding solvent a relatively rapid process. In Chapter \ref{chap:bulkmod} we observe that the cooling rate for these particles (10$^{11}$-10$^{12}$ K/s) is in excess of the cooling rate required for glass formation in bulk metallic alloys.\cite{Greer:1995qy} Given this fact, it may be possible to use laser excitation to melt, alloy and quench metallic nanoparticles in order to form glassy nanobeads.

\section{\label{introSec:MolecularDynamics}Computer Simulation Methods} 

Computational techniques allow us to explore and validate hypotheses about experimentally observed behavior that may otherwise not be accessible through traditional experimental techniques. Additionally, computer simulations allow for theory to propose areas of interest to which physical experimental techniques may be applied. Two of the most common computer simulation methods, Monte Carlo ({\sc MC})and Molecular Dynamics ({\sc MD}) have been widely applied to a variety of model systems ranging from large proteins to metallic clusters. These methods typically (but not always) use classical approximations for the interactions between particles. This allows them to be applied to much larger and more complex systems than would be practical for more detailed {\it ab initio} quantum mechanical models. While these two techniques use similar descriptions for the interactions between particles, Monte Carlo and Molecular Dynamics differ in their fundamental technique for sampling configurations for the system of interest. Monte Carlo sampling is a stochastic approach based on randomly probing the configuration space commensurate with a given ensemble. Conversely, the Molecular Dynamics approach uses the potential energy surface to explicitly calculate the forces within a system and integrate those forces in time to produce a phase space trajectory. Both techniques allow for sampling of thermodynamic quantities within an ensemble, but only Molecular Dynamics allow us to calculate properties that are dynamical in nature. That is to say, properties that depend on knowing the trajectory of the system as a function of time through phase space. Such properties include diffusion constants, viscosities, thermal conductivities, and other time-dependent correlation functions. Since many of the interesting questions pertaining to systems presented in this dissertation involve dynamics, Molecular Dynamics was the method of choice for our computer simulations.\cite{Allen87,Frenkel02,Leach01} 

Many of the systems of interest in this dissertation present unique challenges for existing Molecular Dynamics computational codes. Because metallic potentials and parallel processing capabilities are found in few available Molecular Dynamics programs, we have developed an open source computational code that would address some of these deficiencies. The Object Oriented Parallel Simulation Engine ({\sc oopse}) was the result of the desire to create a Molecular Dynamics engine that incorporated modern object oriented syntax for building a molecular simulation with the ability to simulate large systems that require unique models for molecular interactions.\cite{Meineke:2004uq} We have released {\sc oopse} under a permissive open source license because of the belief that computational codes used to publish scientific material should be transparent and open to peer review. Computational codes used in the first two chapters of this dissertation were the precursor of what was eventually released as {\sc oopse}.

%TODO: Do I need to introduce the concepts of phase space and configuration space?
\subsection{\label{introSec:EmpiricalEfncs}Empirical Energy Functions} It is computationally advantageous to split the potential energy into intra-molecular and inter-molecular contributions to the total potential energy. Typically, intra-molecular potential energy functions include the short-ranged (bonded) portion of the potential and the long-range (non-bonded) interactions are inter-molecular. The general form of the potential is then, 
\begin{equation}
        V = V_{\mathrm{short-range}} + V_{\mathrm{long-range}}. 
\end{equation}
The short-ranged portion includes the explicit bonds, bends, torsions, and improper-torsions  for the molecules present in the simulation. The functional forms and parameters for these interactions are defined by the choice of force field for the simululation.

 For a system containing $N$ atoms, we can decompose the long-range potential into terms that contain the coordinates of individual atoms, pairs, triplets, etc. such that the total long-range potential is given by

\begin{equation}
        V_{\mathrm{long-range}} = \sum_{i}v_{1}(\mathbf{r}_{i}) + \sum_{i}\sum_{j>i}v_{2}(\mathbf{r}_{i},\mathbf{r}_{j}) + \sum_{i}\sum_{j>i}\sum_{k>j>i}v_{3}(\mathbf{r}_{i},\mathbf{r}_{j},\mathbf{r}_{k}) + \ldots .
\label{eq:genlongrange}
\end{equation}
The sum $\sum_{i}\sum_{j>i}$ indicates a summation over all pairs containing $i$ and $j$ which does not double count any pair ($ij$) and does not include self interactions. The first term in (\ref{eq:genlongrange}) is due to an external field acting on the particles in a system such as an external electric field or the interaction of particles with a container wall. The second term in (\ref{eq:genlongrange}) defines a pair potential which often depends only on the magnitude of the pair separation $r_{ij}=|\mathbf{r}_i-\mathbf{r}_j|$. The later terms presented in (\ref{eq:genlongrange}) correspond to three-body interactions, four body-interactions, and so on. For a system of three molecules, the contribution of the pair-pair interactions to the total potential energy would be $V_{\mathrm{AB}}+V_{\mathrm{BC}}+V_{\mathrm{AC}}$ where $V_{\mathrm{AB}}$ is evaluated as though molecule C were absent. This is often unphysical because the presence of molecule C {\em will} modify the interaction between molecules A and B. The correction due to the presence of molecule C can be included in our previous expression for the total potential energy such that there is a three-body correction $V_{\mathrm{long-range}} = V_{\mathrm{AB}}+V_{\mathrm{BC}}+V_{\mathrm{AC}}+V_{\mathrm{ABC}}$. If we have a system of four molecules, we would have to include a four-body correction term.  Because of the expense of computing terms higher than pair interactions and the fact that the contribution of higher order terms is generally small, pair interactions are the primary term for calculating the long-range potential. The pairwise potential can be modified to include average three-body effects forming an {\em effective} pair potential 


\begin{equation}
        V_{\mathrm{long-range}} = \sum_{i}v_{1}(\mathbf{r}_{i}) + \sum_{i}\sum_{j>i} v_2^{\mathrm{eff}}(r_{ij}).
        \label{eq:effpair}
\end{equation}

The functional form for pair-potential represented in computer simulations is typically an effective pair potential that approximates contributions from many-body interactions and includes them in the pair interaction.\cite{Allen87} These functional forms are typically fit to empirical data and quantum mechanical calculations. Effective pair potentials work best for describing the interactions in the noble gases increasing in effectiveness from solid to gas phases. As the system becomes more dense and polarizable, the effective pair approximation fails to adequately describe that system. For example, in a solid, the cohesive energy, $E_{\mathrm{coh}}$ would be calculated as a sum over the bond pairs as indicated by Equation  \ref{eq:genlongrange}. If we consider a set of crystal structures that differ only in coordination, $Z$, of the central atom , we would expect the cohesive energy to scale linearly with coordination. In reality, the cohesive energy scales more weakly due to multi-body screening effects and actually scales as $-Z^{1/2}$.\cite{Nieminen:1990hw,DAW:1993p1640} We will find in Section \ref{introSec:eam} that multi-body effects are very important in accurately describing metallic systems.

In general, the interaction energy between two molecules as a function of their inter-molecular distance takes the shape shown in Figure \ref{fig:ljform}. There is an attractive region at long-range where $-\partial V/\partial R < 0$ while at shorter distances, V, becomes steeply repulsive to account for the low compressibility of particles in condensed phase materials. We represent the distance at which the energy passes through zero and becomes deeply repulsive as $\sigma$. The minimum energy the pair configuration possesses, $-\epsilon$, occurs at the minimum energy interaction distance of $R_m$. 

\begin{figure}[htbp]
        %\centering
                \includegraphics[height=3in]{images/ljform.pdf}
        \caption{General shape of the pair interaction energy between molecules where $-\epsilon$ is the minimum energy, $\sigma$ is the close interaction distance and $R_m$ is the minimum energy separation.}
        \label{fig:ljform}
\end{figure}


As indicated in Equation \ref{eq:genlongrange}, calculation of the long-range (non-bonded) potential involves a sum over all pairs of atoms (as defined by the molecular model) (except for atom pairs which are involved in a short-ranged bonded interaction). It should be noted that atoms in this context refers to interaction sites as defined by the force-field model. These interaction sites may not necessarily correspond to the location of an atomic site in a molecule. The expense of directly calculating all of the long-range interactions for $N$ atoms would involve $N(N-1)/2$ evaluations of interatomic distances. To reduce the number of distance evaluations between pairs of atoms, it is common to use a switched cutoff distance beyond which interactions are not calculated.\cite{Allen87,Frenkel02,Leach01} Further, a Verlet neighbor list can be implemented which creates a look-up of interactions that are within a specific cutoff distance.\cite{Allen87} 

\section{\label{introSec:LJPot}The Lennard-Jones Force Field}

One of the most basic force fields is the Lennard-Jones force field, which mimics the van der Waals interaction at long distances and uses an empirical repulsion at short distances. The general functional form for the Lennard-Jones potential is presented in Figure \ref{fig:ljform} and is given by: 
\begin{equation}
        V_{\text{LJ}}(r_{ij}) = 4\epsilon_{ij} \biggl[ \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12} - \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6} \biggr], \label{eq:lennardJonesPot} 
\end{equation}
where $r_{ij}$ is the distance between particles $i$ and $j$, $\sigma_{ij}$ scales the length of the interaction, and $\epsilon_{ij}$ scales the well depth of the potential.

Interactions between dissimilar particles require parameters for $\sigma$ and $\epsilon$ to be generated. These parameters are usually determined using the Lorentz-Berthelot mixing rules:\cite{Allen87} 
\begin{equation}
        \sigma_{ij} = \frac{1}{2}[\sigma_{ii} + \sigma_{jj}], \label{eq:sigmaMix} 
\end{equation}
and 
\begin{equation}
        \epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}}. \label{eq:epsilonMix} 
\end{equation}

% TODO: General MD introduction, can grab from  OOPSE paper.
\section{\label{introSec:MetallicPotentials}Metallic Potentials} Bonding in metallic systems is significantly different from systems comprising elements that have a mix of covalent, coulombic, and van der Waals interactions. One of the first attempts to model bonding in metals was constructed by Drude in 1900 by extending the classical kinetic theory of gases to electrons in metallic systems.\cite{Ashcroft:1976zt,Drude:1900p1479,Drude:1900p1481} Just as in the kinetic theory of gases, the Drude model assumes that electrons behave as hard spheres that move in a straight line unless colliding with other particles. Electron collision events instantaneously alter the trajectory for an electron. Electronic thermal equilibrium is assumed to occur through electron collisions and obeys Maxwell-Boltzmann statistics. For metals, Drude assumed that electrons were mobile negative charge carriers and that there is a compensating positive charge center attached to much heavier immobile particles. A single atom of a metallic element consists of a nucleus surrounded by a core electron density that is tightly bound to the nucleus. The valence electrons are more weakly bound and form an electron gas when the isolated metallic atoms condense to form a metal. Despite using relatively simplistic assumptions, the Drude model provides a reasonable explanation for AC and DC electrical conductivity, the Hall effect and electron thermal conductivity.\cite{Ashcroft:1976zt}

Sommerfeld subsequently extended the Drude model by incorporating quantum effects.\cite{Ashcroft:1976zt,Kittel:1996fk} The Pauli exclusion principle requires that, since electrons are fermions, only one electron occupy any single electronic state at a time. Additionally, electrons are indistinguishable particles. Fermi-Dirac statistics account for both indistinguishability and anti-symmetric exchange between pairs of electrons by modifying the probability of finding an electron at a given energy 
\begin{equation}
        f(E)=\frac{1}{e^{(E-E_{F})/k_{B}T}+1} \label{eq:fermi-dirac} 
\end{equation}
where $E_{F}$ is the Fermi energy of the system at temperature $T$ and $k_B$ is Boltzmann's constant. $E_{F}$ is by definition the topmost filled level ($n_f$) in a N electron system that is in the ground state confined between $0$ and length $L$.
\begin{equation}
        E_{F} = \frac{\hbar^2}{2m}\left(\frac{n_f\pi}{L}\right)^{2}. 
\end{equation}
Sommerfeld's free electron model uses Fermi-Dirac statistics instead of the classical Maxwell-Boltzmann statistics in the Drude model producing a much more accurate model that successfully describes numerous phenomena in metallic systems. These include the Wiedemann-Franz law relating electrical and thermal conductivity, temperature dependence of the heat capacity, and the electronic density of states.\cite{Ashcroft:1976zt} However, the free electron model makes quantitative predictions that are directly contradicted by experimental observation. Ashcroft and Mermen have an extensive discussion in Chapter 3 of their book describing the failures of the free electron model due to the oversimplification of the role of ions and the fundamental nature of
interactions in metallic systems.\cite{Ashcroft:1976zt} Most significant of these failure is the inability to distinguish between metals, semi-metals, semiconductors and insulators and the relation of valence electrons in a free atom to the conduction band electrons in a metal. 

Many of these deficiencies can be corrected by accounting for ions contained in a periodic lattice with an electron density modulated by these lattice positions. A model pseudo-potential for the pair-pair interaction can be constructed that has the form

\begin{equation}
        v_2(r_{ij}) = v_d(r_{ij}) + v_l(r_{ij})
        \label{eq:metpseudo}
\end{equation}
where $v_d(r_{ij})$ represents the direct interaction between ion pairs and $v_l(r_{ij})$ is a indirect interaction due to ions being placed in a bath of electrons.\cite{Egelstaff:1992yb} Effective Medium Theory ({\sc EMT}) was one of the first pseudo-potential models that uses such a functional form for the potential.\cite{Nrskov:1980p1752,Nrskov:1982p1753,Stott:1980p1754} 

In {\sc EMT}, the real material is replaced by {\it jellium} which consists of a homogeneous electron gas. The metal ions are replaced by a constant positive background density. The first order approximation for the change in energy due to embedding an atom in jellium at position $\mathbf{r}$ is

\begin{equation}
        \Delta E(\mathbf{r}) = E_{\mathrm{atom+jellium}}-( E_{\mathrm{atom}}+E_{\mathrm{jellium}})
        \equiv \Delta E^{\mathrm{embed}}(\rho_0(\mathbf{r}))
\end{equation}
where $\Delta E^{\mathrm{embed}}(rho)$ is the {\it embedding energy} of an atom in a homogeneous electron gas with density $\rho$ and $\rho_0(\mathbf{r})$ is the electron density at $\mathbf{r}$. Puska et al. showed that the embedding energy is a universal function of the electron density.\cite{Puska:1981p1755} They found that for the noble gases, the embedding energy was linear for all values of $\rho$ because the closed electron shell only causes a repulsive interaction. Van der Waals interactions do provide attractive interactions, but these effects are not part of the embedded interaction. For other elements, there is a minimum in the potential due to the propensity of these elements to form bonds. Several different potential models for metals share the basic notion of {\sc EMT} but differ in the functional form for the embedding energy.

%TODO paragraph on Egelstaff metal potential, first minimum much deeper then 1/r^6, second minumum do to ion-ion pair interactions.
\subsection{\label{introSec:eam}Embedded Atom Method}

In section \ref{introSec:EmpiricalEfncs}, the pair potential approximation was introduced as a computationally expedient method for calculating the interaction energy for a system of molecules. However, it was pointed out that the effective pair approximation was inadequate for describing interactions in metallic solids and liquids. An example of the inadequacy of this approximation for metallic systems is presented by Daw {\it et al.} involving differences in elastic constants for cubic lattice solids composed of noble gases and fcc metals.\cite{DAW:1993p1640} It is known for a solid described by a purely pairwise interaction that the elastic stiffness constants $C_{12}$ and $C_{44}$ are equal. This is known as the Cauchy relation. We can also define a Cauchy pressure as the difference between the two elastic constants, $P_c=C_{12}-C_{44}$.\cite{Kittel:1996fk}  Both Ar and Kr have ratios very close to $1$, but fcc metals have ratios closer to 2. Pt and Au have ratios that are between 3 and 4. This can be seen in Table \ref{tab:manybod} (which is reproduced from Table 1 in Daw {\it et al.}). Also included in Table \ref{tab:manybod} is the ratio of vacancy formation energy ($E^f_v$) to cohesive energy($E_{\mathrm{coh}}$). In a solid with purely pairwise interactions, this ratio is once again $1.0$. The metals have values closer to $0.35$ with Au and Pt again as outliers significantly deviating from pairwise bonding. The last indicator presented in Table \ref{tab:manybod} is that of the ratio between cohesive energy and melting temperature. Once again, FCC metals differ significantly from a strictly pair-wise interaction model.

\begin{table}[tpb]
        \setlength{\capwidth}{0.7\textwidth}
        \begin{center}
        
\begin{minipage}{\textwidth}
        \renewcommand{\thefootnote}{\thempfootnote}
        \caption[DEVIATION FROM PAIRWISE BEHAVIOR IN FCC METALS COMPARED WITH NON-METALLIC SYSTEMS]{ DEVIATION FROM PAIRWISE BEHAVIOR IN FCC METALS COMPARED WITH NON-METALLIC SYSTEMS\footnote{Reproduced from Daw {\it et al.} \cite{DAW:1993p1640} Table 1.}}
        \centering
\begin{tabular}{cccc}
\toprule
\toprule
Solid & $\frac{C_{12}}{C_{44}}$ & $\frac{E^f_v}{E_{\mathrm{coh}}}$ & $\frac{E_{\mathrm{coh}}}{kT_m}$\\
\midrule
Pair Potential & & &\\
%\hline
LJ & 1.0 & 1.00 & 13\\
\hline
Rare Gases &&&\\
%\hline
Ar & 1.1 & 0.95 & 11\\
%\hline
Kr & 1.0 & 0.66 & 12\\
\hline
FCC Metals & & &\\
%\hline
Ni & 1.2 & 0.31 & 30\\
%\hline
Cu & 1.6 & 0.37 & 30\\
%\hline
Pd & 2.5 & 0.36 & 25\\
%\hline
Ag & 2.0 & 0.39 & 27\\
%\hline
Pt & 3.3 & 0.26 & 33\\
%\hline
Au & 3.7 & 0.23 & 34\\
\bottomrule
\end{tabular}
\renewcommand{\footnoterule}{}
\end{minipage}
\label{tab:manybod}

\end{center}
\end{table}

Accurate description of metallic bonding requires the consideration of coordination dependent bonding but at the same time must be computationally efficient to implement. The Embedded Atom Method ({\sc EAM}) of Daw and Baskes uses the idea in Equation \ref{eq:metpseudo} of separating metallic interactions into those due to ion-pair interactions and those due to embedding an ion in an electron gas.\cite{Daw84,DAW:1983ht} Density Functional Theory is a useful starting point to determine the embedding energy since one can prove using this theory that the ground state energy of a system of atoms can be written as a unique functional of the electronic density.\cite{Hohenberg:1964bs} Using density-functional theory as a starting point, they derived an approximate expression for the cohesive energy in a metallic system. The density-functional expression for the cohesive energy of a solid is given by

\begin{equation}
        E_{\mathrm{coh}} = F[\rho] + \frac{1}{2} \sum_{i,j,i\neq j} \frac{Z_i Z_j}{R_{ij}} - \sum_{i} \int \frac{Z_i \rho(\mathbf{r})}{|\mathbf{r}-\mathbf{R}_i|}\mathrm{d}\mathbf{r} + \frac{1}{2} \int\int \frac{\rho(\mathbf{r}_1)\rho(\mathbf{r}_2)}{r_{12}} \mathrm{d}\mathbf{r}_1 \mathrm{d}\mathbf{r}_2 - E_{\mathrm{atoms}}
        \label{eq:eamdensfunc}
\end{equation}
where i and j represent the nuclei of the solid and $Z_i$, $\mathbf{R}_i$ represent the charge and position of the $i$th nucleus and $\rho$ is the electron density. The collective energy of the isolated atoms is given by $E_{\mathrm{atoms}}$ and $F[\rho]$ is the kinetic, exchange, and correlation energy functional. 

If we assume that $F[\rho]$ can be described by 
\begin{equation}
        F[\rho] = \int F(\rho(\mathbf{r}),\nabla\rho(\mathbf{r}),\nabla^{2}\rho(\mathbf{r}),\ldots) \mathrm{d}\mathbf{r}
\end{equation}  
where $F$ is the density of the solid and is described by the local electron density. This assumption has been justified by studies of response functions in nearly free-electron gases. If we further assume that electron density in a solid can be described as a linear superposition of the electron density due to individual atoms such that $\rho_s(\mathbf{r}) = \sum_{i}\rho_{i}(\mathbf{r}-\mathbf{R}_i)$. Using the two assumptions, Daw and Baskes\cite{Daw84,DAW:1983ht} derived an expression for the cohesive energy in metals as

\begin{equation}
        V  =  \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i}
        \phi_{ij}( R_{ij}) + E_{\mathrm{err}}
        \label{introeq:eam}
\end{equation} 
where $\rho_{i}$ is given by

\begin{equation}
        \rho_{i} = \sum_{j\ne i} \rho_j^a(R_{ij}).
        \label{introeq:eamrho}
\end{equation}

The error associated with Equation \ref{introeq:eam} is represented by the term $E_{\mathrm{err}}$ and is a function of the background electron density $\bar\rho_i$. Daw derives an expression for the optimal background density by setting $E_{\mathrm{err}} = 0$.\cite{DAW:1989p1673} It is the embedding function that is the source of the many-body effects in {\sc EAM}. The embedding function is non-linear, reflecting saturation in the metallic bond by increasing background electron density. A useful way of viewing the many body effect of $F[\rho]$ is that as an atom makes more bonds, each new bond constructed is weaker than the previous bond. This is reasonable because constructing a new bond increases the total bonding energy but decreases the average energy per bond. The weakening of successive bonds will occur if $F(\rho)$ has positive curvature ($d^2F/d\rho^2 > 0$). Puska et al. demonstrated this from first principle calculations for embedding an atom in a homogenous electron gas.\cite{Puska:1981p1755} The functional forms for the various parameters are plotted in \ref{introfig:eam_all} for the metallic elements that are the focus of this dissertation. Positive curvature in $F[\rho]$ can be seen in the right panel in this Figure.

\begin{figure}[htbp]
        %\centering
                \includegraphics[width=\linewidth]{images/eam_all.pdf}
        \caption{Functional Forms for {\sc EAM} Parameters.}
        \label{introfig:eam_all}
\end{figure}

The pair potential in {\sc EAM} is defined to be 

\begin{equation}
        \phi_{ij}(R) = \left(\frac{1}{4\pi\epsilon_0}\right)\left(\frac{Z_i^a(R)Z_j^b(R)}{R}
        \right)
        \label{eq:eamphi}
\end{equation}
where the $Z^a_i(R)$ are effective screened charges of the nuclei ($i$) corresponding to atom of identity $a$. In the original formulation of {\sc EAM}, the pair-interaction was purely repulsive in nature.\cite{PhysRevB.33.7983} Later formulations of {\sc EAM} include an attractive component in the pair-interaction potential.\cite{Voter:95} $Z(R)$ was then fit to an empirical parameterized form

\begin{equation}
        Z(R) = Z_0(1+\beta R^\nu)e^{-\alpha R}
\end{equation}
where $Z_0$ is given by the number of outer electrons in the atom. The three remaining parameters, $\alpha$, $\beta$ and $\nu$ are determined empirically to give the desired screening to $Z(R)$. Likewise, $F[\rho]$ can be derived using a universal function that relates the sublimation energy of a metal $E(a)$, to the lattice constant, $a$.\cite{Rose:1984rw} The universal function

\begin{equation}
        E(a)=-E_{\mathrm{sub}}(1+a^*)e^{-a^*}
        \label{eq:rose_eqs}
\end{equation}
where $a^*$ is given by
\begin{equation}
a^*=(a/a_0-1)/(E_{\mathrm{sub}}/9B\Omega)^{1/2}.
\end{equation}
$E_{\mathrm{sub}}$ in the first expression corresponds to the absolute value of sublimation energy at zero temperature and pressure. The deviation from the equilibrium lattice constant, $a_0$, is determined by $a^*$ where $a$ is the fcc lattice constant. The bulk modulus for the system is given by $B$ and $\Omega$ is the equilibrium volume per atom. If the atomic electron densities as a function of $R$ and the pair interaction energy are both known, then the embedding energy can be determined by requiring the total energy computed using Equation \ref{introeq:eam} to agree with the equation of state given in Equation \ref{eq:rose_eqs}.


{\sc EAM} is appropriate for metallic systems where directional bonding is not important and thus fails to describe the behavior of elements that are semiconductors or are in the middle of the transition series. A directional form of {\sc EAM} has been developed to account for directional bonding effects.\cite{BASKES:1987p1743,BASKES:1989p1746,BASKES:1992p1735} However, the directional bonding potentials tend to be more computationally expensive because of their more complex functional forms for the electron density.


\subsection{\label{introSec:tightbind}Tight-Binding Formulation}

Other potentials have been developed that share a similar functional form to {\sc EMT} and therefore to {\sc EAM}.~\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} These are briefly reviewed by Voter and summarized in this section.\cite{Voter:95} Finnis and Sinclair developed a `N-body potential' based in part on the second moment approximation {\sc SMA} to tight-binding theory.\cite{Finnis84} In this theory, atomic orbitals interact via a one-electron Hamiltonian describing their interaction. Diagonalizing the system Hamiltonian matrix yields the set of orbital energies that can be populated with electrons. The distribution of energy levels can be represented by a density of states ({\sc DOS}) where the bond energy can be written as
\begin{equation}
        E_{\mathrm{bond}} = \int_{-\infty}^{E_f} \epsilon n(\epsilon) \mathrm{d}\epsilon
        \label{introeq:tb_dos}
\end{equation} 
where $n(\epsilon)$ is the density of states and $E_f$ is the Fermi level energy. Because it is relatively easy to compute Hamiltonian moments on an atom-by-atom basis, one can compute an approximate {\sc DOS} shape from lower moments of the distribution reducing the computational work-load. As an example, Voter assumes a single s orbital basis per atom $i$. The second moment of the {\sc DOS} is then

\begin{equation}
        \mu_{2i} = \left<\chi_i|H^2|\chi_i\right> = \sum_{j} \left<\chi_i|H|\chi_j\right> \left<\chi_j|H|\chi_i\right> = \sum_{j} h^2_{ij}.
        \label{introeq:tb_basis}
\end{equation}
In Equation \ref{introeq:tb_basis} $\chi_i$ is the basis function on atom $i$ and $h^2_{ij}=\left<\chi_i|H^2|\chi_i\right>$. Note that summation in Equation \ref{introeq:tb_basis} is due to the complete expansion in a given basis set. One can show that for any basic shape of the {\sc DOS} on atom $i$, the width of the {\sc DOS} is proportional to $(\mu_{2i})^{1/2}$. Evaluation of Equation \ref{introeq:tb_dos} given $\mu_{2i}$ and setting the first moment of the DOS, $\mu_{1i}=0$, gives
\begin{equation}
        E_{\mathrm{bond,i}} = A(\mu_{2i})^{1/2}
        \label{introeq:tb_form}
\end{equation}
where A is constant that is determined by the shape of the given {\sc DOS} and by the fractional electronic occupation of the orbital. 

In the spirit of {\sc EMT}, the tight-binding energy defined in Equation \ref{introeq:tb_form} should be augmented by a pairwise sum to account for ion-ion interactions. Using \ref{introeq:tb_basis} and \ref{introeq:tb_form} we can derive an expression for the metallic potential at atom $i$ due to all other atoms $j$ as

\begin{equation}
        E_{i} = \frac{1}{2}\sum_{j} \phi(r_{ij}) - A\left(\sum_{j} h^2_{ij}(r_{ij})\right)^{1/2}.
        \label{introeq:finnsincair}
\end{equation}
It should be obvious when comparing Equation \ref{introeq:finnsincair} to the functional form for {\sc EAM} that the electron density corresponds to 
\begin{equation}
        \rho(\mathbf{r}) \equiv h^2_{ij}(\mathbf{r})
        \label{introeq:eamrhoequiv}
\end{equation}
and
\begin{equation}
        F[\rho] \equiv -A(\rho)^{1/2}.
        \label{introeq:eamfrhoequiv}
\end{equation}

It should be noted that the negative square-root form for the embedding function has a positive second derivative for all $\rho$ and therefore is consistent with the trend in metallic bond energies. This particular embedding function also has a negative slope and therefore is a totally cohesive force and must be balanced by a repulsive pair potential. For alloy systems there is a slight difference between methods based on {\sc EAM} and {\sc SMA} due the physical origin of the electron density in their respective functional forms. Terms in the density summation for atom $i$ in the {\sc EAM} formulation depend on the atom type for neighbor $j$, but {\em not} on the identity of atom $i$ itself. The fixed electron density on atom $j$ is independent of the electron density associated with atom $i$. Equation \ref{introeq:tb_basis} presents the corresponding {\sc SMA} expression. Here the electron density is based on the square of matrix elements between atom $i$ and atom $j$ and therefore are dependent on the identity of both $i$ and $j$. {\sc SMA} therefore requires knowledge of $h_{ij}(r)$ to construct an alloy pair potential. 


The tight-binding model implemented in {\sc OOPSE} and used in this dissertation for the metallic glass simulations is due to Sutton and Chen ({\sc sc}).\cite{Chen90} The {\sc sc} potential takes a form similar to that of the Finnis-Sinclair potential in Equation \ref{introeq:finnsincair},
\begin{equation}
\label{introeq:SCP1}
U_{tot}=\sum _{i}\left[ \frac{1}{2}\sum _{j\neq
i}D_{ij}V^{pair}_{ij}(r_{ij})-c_{i}D_{ii}\sqrt{\rho_{i}}\right] ,
\end{equation}
 where the pair potential, $V^{pair}_{ij}$ is given by
\begin{equation}
\label{introeq:SCPPair}
V^{pair}_{ij}(r)=\left(
\frac{\alpha_{ij}}{r_{ij}}\right)^{n_{ij}} 
\end{equation}
once again accounts  for repulsion between ionic atomic cores. The local density term is then given by
\begin{equation}
        \rho_{i}=\sum_{j\neq i}\left(
        \frac{\alpha_{ij}}{r_{ij}}\right) ^{m_{ij}}.
\end{equation}
$\alpha$ is a length parameter scaling the particle spacing, c is a dimensionless parameter that scales the attractive portion of the potential. $\epsilon$ sets the entire energy scale. $n$, $m$ are integer parameters where $n>m$. The exponential parameters and c are fit to the equilibrium lattice parameter and $\epsilon$ is determined by the cohesive energy. $n$,$m$ can be determined from the Cauchy pressure, $P_c$, such that

\begin{equation}
        n = 3\sqrt{\left(\frac{\Omega B}{E_{\mathrm{coh}}}\right)\left(\frac{B}{2P_c}+1\right)}
\end{equation}

\begin{equation}
        m=6\sqrt{\frac{\left(\frac{\Omega B}{E_{\mathrm{coh}}}\right)}{\left(\frac{B}{2P_c}+1\right)}}
\end{equation}
where $\Omega$ once again is the volume per atom. 

The {\sc SC} potential was re-parameterized by Kimura et al. to include in addition to the experimental lattice parameter, the cohesive energy, bulk modulus, phonon frequencies, vacancy formation energy, and surface energies.\cite{kimura-quantum,PhysRevB.59.3527} The new parameterization was named Quantum Sutton-Chen ({\sc Q-SC}) because it takes into account zero-point energy. The functional forms for {\sc SC} and {\sc Q-SC} have computational advantages because they somewhat resemble the Lennard-Jones potential presented in Section \ref{introSec:LJPot}. In this dissertation the {\sc q-sc} formulation has been utilized for potential energies and forces. Combination rules for the alloy were taken to be the arithmetic average of the atomic parameters with the exception of $c_i$ since its
value is only dependent on the identity of the atom where the density
is evaluated.  For the {\sc q-sc} potential, cutoff distances are
traditionally taken to be $2\alpha_{ij}$ and include up to the sixth
coordination shell in FCC metals. Parameters for the various versions of Sutton-Chen  are presented in Table \ref{introtab:qsc}.

\begin{table}
        \begin{minipage}{\textwidth}
                \renewcommand{\thefootnote}{\thempfootnote}
                \caption[{\sc SC} AND {\sc Q-SC } PARAMETERS]{ {\sc SC} AND {\sc Q-SC } PARAMETERS\footnote{Reproduced from Kimera {\it et al.} \cite{kimura-quantum} Table 1.}}        
        \centering
\begin{tabular}{llccccc}
        \toprule
        \toprule
 & & n & m & $\epsilon$(eV) & c & $\alpha$ (\AA)\\
\midrule
Ni & ({\sc Q-SC}) & 10 & 5 & 7.3767E-3 & 84.745 & 3.5157\\

   & ({\sc SC})   & 9 & 6 & 1.5714E-2 & 39.756 & 3.5200\\

Cu & ({\sc Q-SC}) & 10 & 5 & 5.7921E-3 & 84.843 & 3.6030\\

   & ({\sc SC}) & 9 & 6 & 1.2351E-3 & 39.756 & 3.6100\\

Rh &({\sc Q-SC}) & 13 & 5 & 2.4612E-3 & 305.499 & 3.7984\\

   & ({\sc SC}) & 12 & 6 & 4.9371E-3 & 145.658 & 3.8000\\

Pd & ({\sc Q-SC}) & 12 & 6 & 3.2864E-3 & 148.205 & 3.8813\\

   & ({\sc SC}) & 12 & 7 & 4.1260E-3 & 108.526 & 3.8900\\

Ag & ({\sc Q-SC}) & 11 & 6 & 3.9450E-3 & 96.524 & 4.0691\\

   & ({\sc SC}) & 12 & 6 & 2.5330E-3 & 145.658 & 4.0900\\

Ir & ({\sc Q-SC}) & 13 & 6 & 3.7674E-3 & 224.815 & 3.8344\\

   & ({\sc SC}) & 14 & 6 & 2.4524E-3 & 337.831 & 3.8400\\

Pt & ({\sc Q-SC}) & 11 & 7 & 9,7894E-3 & 71.336 & 3.9136\\

   & ({\sc SC}) & 10 & 8 & 1.9768E-2 & 34.428 & 3.9200\\

Au & ({\sc Q-SC}) & 11 & 8 & 7.8052E-3 & 53.581 & 4.0651\\

   & ({\sc SC}) & 10 & 8 & 1.2896E-2 & 34.428 & 4.0800\\
\bottomrule
\end{tabular}
\renewcommand{\footnoterule}{}
\end{minipage}
\label{introtab:qsc}
\end{table}

\section{\label{introSec:eom}Integrating Equations of Motion}


Once a potential model governing the interaction between particles has been selected, the next step is to propagate the system forward in time according to the classical equations of motion. At any particular time, the system is defined by the positions ($\mathbf{q}$), velocities ($\dot{\mathbf{q}}$) (or momenta ($\mathbf{p}$))  and the forces on each particle comprising the system. We can formulate the equations of motion in the Lagrangian form 

\begin{equation}
L = T - V,
\label{introeq:lagrang}
\end{equation}
where $L$ is a function of the total kinetic ($T$) and potential energy ($V$). Using variational calculus and Hamilton's principle that the variation of the integral of the Lagrangian with respect to time is zero, we can derive a general set of equations of motion.\cite{Tolman:1938kl,Goldstein:2001uf,Allen87} Mathematically we represent this principle for a system of $f$ degrees of freedom by

\begin{equation}
\delta\int_{t_1}^{t_2}L(q_1\ldots q_f,\dot{q_1}\ldots \dot{q_f},t)dt = 0.
\label{introeq:variatLag}
\end{equation}
In terms of derivatives with respect to positions and velocities, we can rewrite \ref{introeq:variatLag} as
\begin{equation}
\int_{t_1}^{t_2}\sum_{i=1}^{f}\left(
        \frac{\partial L}{\partial \mathbf{q}_i}\delta \mathbf{q}_i
        + \frac{\partial L}{\partial \dot{\mathbf{q}}_i}\delta 
                \dot{\mathbf{q}}_i\right)dt = 0.
                \label{introeq:variatLag2}
\end{equation}
Integrating the term in $\dot{\mathbf{q}}$ (by parts) we can simplify \ref{introeq:variatLag2} to
\begin{equation}
\int_{t_1}^{t_2}\sum_{i=1}^{f}\left(
        \frac{d}{dt}\frac{\partial L}{\partial \dot{\mathbf{q}}_i}
        - \frac{\partial L}{\partial \mathbf{q}_i}\right)
                \delta {\mathbf{q}}_i dt = 0.
\end{equation}
Now, we have achieved a separation of variables in \ref{introeq:variatLag2}, we can further separate the contribution from each variable
\begin{equation}
\frac{d}{dt}\frac{\partial L}{\partial \dot{\mathbf{q}}_i}
        - \frac{\partial L}{\partial \mathbf{q}_i} = 0
                \quad\quad(i=1,2,\dots,f).
\label{introeq:eom1}
\end{equation}
Newton's equations of motion can be recovered from the above formulation by
substituting Equation \ref{introeq:lagrang} into \ref{introeq:eom1} and substituting the definition for kinetic energy,
$m\dot{\mathbf{q}}^2/2$, gives
\begin{equation}
\frac{d}{dt}\left(\frac{m\dot{\mathbf{q}}}{2}\right)+\frac{dV}{d\mathbf{q}}=0,
\end{equation}
or in the form of Newton's famous second law,
\begin{equation}
\mathbf{f} = m\mathbf{a}.
\end{equation}


\subsection{\label{introSec:verlet}Verlet Method of Intergration}
The evolution of time in Molecular Dynamics requires an integrator that is symplectic, meaning that it conserves the canonical circular structure of the Hamiltonian and guarantees that the error in the energy is bounded at each time step.\cite{McQuarrie:2000yt,Tolman:1938kl,Goldstein:2001uf} Given a Hamiltonian that is separable with respect to  $\mathbf{q}$ and ${\mathbf{p}}$
\begin{equation}
        H(\mathbf{q},\mathbf{p}) = T(\mathbf{p}) + V(\mathbf{q})
\end{equation}
and is described by the phase space $\mathbf{q}=\{q_\alpha\}$,$\mathbf{p}=\{p_\alpha\}$, $\alpha = 1,\ldots,f$. Our goal is to produce a transformation the preserves the flow in phase space generated by H such that
\begin{equation}
        (\mathbf{q}_0,\mathbf{p}_0)\text{ at time } t_0 \rightarrow (\mathbf{q},\mathbf{p})\text{at time } t
\end{equation}
and where we denote time step as $t-t_0 = \delta t$. Hamilton's equations in canonical form are
\begin{equation}
        \dot{\mathbf{p}} = -\frac{\partial H}{\partial \mathbf{q}}
        \label{eq:cannon1}
\end{equation}
and
\begin{equation}
        \dot{\mathbf{q}} = -\frac{\partial H}{\partial \mathbf{p}}.
        \label{eq:cannon2}
\end{equation}
Liouville's theorem describes how the phase space density, $\rho$ at a point of interest changes as a function of time
\begin{equation}
        \left(\frac{\partial \rho}{\partial t}\right)_{\mathbf{p},\mathbf{q}}=-\sum_i \left(\frac{\partial \rho}{\partial \mathbf{q}_i}\dot{\mathbf{q}_i}+\frac{\partial \rho}{\partial \mathbf{p}_i}\dot{\mathbf{p}_i}\right).
        \label{eq:liouville}
\end{equation}
Substituting \ref{eq:cannon1} and \ref{eq:cannon2} into \ref{eq:liouville}, we can rewrite Liouville's theorm as
\begin{equation}
        \left(\frac{\partial \rho}{\partial t}\right)_{\mathbf{p},\mathbf{q}}=-\sum_i \left(\frac{\partial \rho}{\partial \mathbf{q}_i}\frac{\partial H}{\partial \mathbf{p}_i}
        +\frac{\partial \rho}{\partial \mathbf{p}_i}\frac{\partial H}{\partial \mathbf{q}_i}\right).
        \label{eq:liouville2}
\end{equation}
If we define the Poisson bracket for two quantities $M$,$N$ as
\begin{equation}
        \{M,N\}= \sum_i \left(\frac{\partial M}{\partial \mathbf{q}_i}\frac{\partial N}{\partial \mathbf{p}_i}
        -\frac{\partial N}{\partial \mathbf{q}_i}\frac{\partial M}{\partial \mathbf{p}_i}\right).
        \label{eq:poisson}
\end{equation}
Equation \ref{eq:liouville2} is then
\begin{equation}
        \left(\frac{\partial \rho}{\partial t}\right)_{\mathbf{q},\mathbf{p}}=-\{\rho,H\}.
        \label{eq:liouville3}
\end{equation}
Using the Liouville operator, $L$, we can write the solution to Equation \ref{eq:liouville3} as 
\begin{equation}
        \rho(\mathbf{\dot{\mathbf{q}}},\mathbf{q},t) = U(t)\rho(\mathbf{\dot{\mathbf{q}}},\mathbf{q},0)
\end{equation}
where
\begin{equation}
        U(t) = e^{iLt}
        \label{eq:propagator}
\end{equation}
is the {\em classical propagator}. In relation to the Liouville equation, 
\begin{equation}
        iL = \{\cdots,H\}=\left[\frac{\vec{p}}{m}\frac{\partial}{\partial \vec{q}} + F(\vec{q})\frac{\partial}{\partial \vec{p}}\right].
\end{equation}
In numerical simulations $t$ is broken into discrete time steps $\Delta t$. To generate a discrete form for Equation \ref{eq:propagator}, we perform an operator expansion
\begin{equation}
        e^{iLt}=[e^{iL\Delta t}]^P,
        \label{eq:prop_expans1}
\end{equation}

\begin{equation}
        e^{iL\Delta t}=e^{iL_1(\Delta t/2)}e^{iL_2(\Delta t)}e^{iL_1(\Delta t/2)}+O(\Delta t^3)
        \label{eq:prop_expans3}
\end{equation}
where the operator $iL = iL_1+iL_2$ and $\Delta t=t/P$. The {\em Position} form of the Verlet integrator can be generated by substituting $iL_1=F(\mathbf{q})\partial/\partial \mathbf{p}$ and $iL_2=(\mathbf{p}/m)\partial/\partial \mathbf{q}$ and evaluating the evolution operator.\cite{swope:637,Verlet67,Allen87,Leach01,tuckerman:2278} The resultant integrator for the position is
\begin{equation}
        \mathbf{q}(\Delta t)=\mathbf{q}(0) + \frac{\Delta t}{m} \mathbf{p}(0) + \frac{\Delta t^2}{2m}F[\mathbf{q}(0)],
\end{equation}
and momentum
\begin{equation}
        \mathbf{p}(\Delta t) = \mathbf{p}(0)+\frac{\Delta t}{2}(F[\mathbf{q}(0)]+F[\mathbf{q}(\Delta t)]).
\end{equation}
The velocity Verlet algorithm was designed to improve the error associated with the calculation of velocities and therefore improve overall energy conservation. Velocity Verlet involves a half time step propagation of momenta
\begin{equation}
\mathbf{p}\left(t+\frac{1}{2}\Delta t\right) = \mathbf{p}(t) 
                                        + \frac{1}{2}\Delta t m \mathbf{a}(t).
\end{equation}
and a full time step propagation of the positions
\begin{equation}
\mathbf{q}(t+\Delta t) = \mathbf{q}(t) + \frac{\Delta t}{m}\mathbf{p}(t)
                                + \frac{1}{2}\Delta t^2\mathbf{a}(t),
\end{equation}
In many algorithms, this is referred to as {\tt moveA:}. Next, forces are calculated based on the new positions and then the momentum is propagated to a full time step in {\tt moveB:}
\begin{equation}
\mathbf{p}(t+\Delta t) = \mathbf{p}\left(t+\frac{1}{2}\Delta t\right) 
                                + \frac{m}{2}\Delta t\mathbf{a}(t+\Delta t).
\end{equation}
Velocity verlet reduces the error in calculation of the velocities to $\mathcal{O}(\Delta t^3)$, but at the expense of the error associated with the positions which increases to 
$\mathcal{O}(\Delta t^3)$.

\subsection{\label{introSec:LD}Langevin Dynamics}
Langevin dynamics is simplified model for the interaction of a particle with a surrounding solvent. This model uses a stochastic differential equation to approximate the effects of a surrounding solvent thus eliminating the need for an explicit solvent.\cite{Allen87,Frenkel02,Leach01,BROOKS:1985kx,BROOKS:1983uq,BRUNGER:1984fj,Schlick:2002hc} The Langevin equation of motion is given by 
\begin{equation}
        m\ddot{q} = -\mathbf{\nabla}_q V - \gamma\dot{q}(t) + \mathbf{R}(t)
        \label{introeq:langevin}
\end{equation}
where $-\mathbf{\nabla}_q V$ is the force due the the potential acting on a particle, $\gamma\dot{q}(t)$ is a frictional force and $\mathbf{R}(t)$ is a stochastic force. If the form of the potential, $V$ is quadratic or absent a potential, then Equation \ref{introeq:langevin} is known to be exact. For other potentials, the Langevin equation may not be exact. A more general form of the Langevin equation, is the {\em Generalized Langevin Equation}
\begin{equation}
        m\ddot{x} = -\frac{\partial U(x)}{\partial x} - \int_0^t \gamma(t)\dot{x}(t-\tau)\mathrm{d}\tau \mathbf+{R}(t)
        \label{introeq:GLE}
\end{equation} 
where $\gamma(t)$ is a friction kernel containing all solvent responses to a moving solute and $\mathbf{R}(t)$ is a random fluctuating force. $\mathbf{R}(t)$ often represents a Gaussian random process, so will have a time average of $0$
\begin{equation}
        \left<R(t)\right> = 0
\end{equation}
and the correlation between any two times must satisfy the fluctuation dissipation theorem: 
\begin{equation}
        \left<R(t)R(t')^T\right> = 2 m \gamma k_B T \delta(t-t')
\end{equation}
where $T$ is the temperature, $k_B$ is the Boltzmann constant and $\delta$ is the Dirac delta function. The strength of the inertial forces relative to the random force is determined by the magnitude of $\gamma$. $\gamma$ can be related to the viscosity of an implicit solvent, $\eta$, by Stokes' law
\begin{equation}
        \gamma = \frac{6\pi\eta a}{m}
\end{equation}
where $a$ is the hydrodynamic radius of the particle and $m$ is the particle's mass.

\section{Parallel Molecular Dynamics} 
Although processor power is continually improving, it is still
difficult to simulate systems of more than 10,000 atoms on a single
processor. To facilitate study of larger system sizes or smaller
systems for longer time scales, parallel methods were developed to
allow multiple CPUs to share the simulation workload. Three general
categories of parallel decomposition methods have been developed:
these are the atomic,\cite{Fox88} spatial~\cite{plimpton95} and
force~\cite{Paradyn} decomposition methods.

Algorithmically simplest of the three methods is atomic decomposition,
where $N$ particles in a simulation are split among $P$ processors for
the duration of the simulation. Computational cost scales as an
optimal $\mathcal{O}(N/P)$ for atomic decomposition. Unfortunately, all
processors must communicate positions and forces with all other
processors at every force evaluation, leading the communication costs
to scale as an unfavorable $\mathcal{O}(N)$, \emph{independent of the
number of processors}. This communication bottleneck led to the
development of spatial and force decomposition methods, in which
communication among processors scales much more favorably. Spatial or
domain decomposition divides the physical spatial domain into 3D boxes
in which each processor is responsible for calculation of forces and
positions of particles located in its box. Particles are reassigned to
different processors as they move through simulation space. To
calculate forces on a given particle, a processor must simply know the
positions of particles within some cutoff radius located on nearby
processors rather than the positions of particles on all
processors. Both communication between processors and computation
scale as $\mathcal{O}(N/P)$ in the spatial method. However, spatial
decomposition adds algorithmic complexity to the simulation code and
is not very efficient for small $N$, since the overall communication
scales as the surface to volume ratio $\mathcal{O}(N/P)^{2/3}$ in
three dimensions.

The parallelization method used in our codes is the force
decomposition method.\cite{hendrickson:95} Force decomposition assigns
particles to processors based on a block decomposition of the force
matrix. Processors are split into an optimally square grid forming row
and column processor groups. Forces are calculated on particles in a
given row by particles located in that processor's column
assignment. One deviation from the algorithm described by Hendrickson
{\it et al.} is the use of column ordering based on the row indexes
preventing the need for a transpose operation necessitating a second
communication step when gathering the final force components.  Force
decomposition is less complex to implement than the spatial method but
still scales computationally as $\mathcal{O}(N/P)$ and scales as
$\mathcal{O}(N/\sqrt{P})$ in communication cost. Plimpton has also
found that force decompositions scale more favorably than spatial
decompositions for systems up to 10,000 atoms and favorably compete
with spatial methods up to 100,000 atoms.\cite{plimpton95}
Revision:	3496
Committed:	Wed Apr 8 19:13:41 2009 UTC (16 years ago) by chuckv
Content type:	application/x-tex
File size:	51127 byte(s)
Log Message:	Final Version