start/chuckDissertation/introduction.tex



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

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

% TODO: Write a general introduction tying everything together.
This dissertation presents research that I have performed and been involved with over my course of 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 on metallic system is presented in the opening chapter.

\section{\label{introSec:MolecularDynamics}COMPUTER SIMULATION METHODS} 

Computational techniques allow us to explore and validate hypotheses for 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 use classical approximations for interactions between particles that allow them to be applied to much larger and more complex systems than would be practical for more detailed ab-initio quantum based 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 phase space configurations for the system of interest. Monte Carlo is a stochastic approach based upon exploring a systems potential energy surface by randomly probing the systems configuration space within a given ensemble. Conversely, the Molecular Dynamic 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, viscosity, thermal conductivity, 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 model 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-torsons which have been defined in the {\sc oopse} meta-data file 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}_{j}) + \ldots .
\label{eq:genlongrange}
\end{equation}
The sum $\sum_{i}\sum_{j>i}$ indicates a summation over all pairs containing $i$ and $j$ while not double counting any pair $ij$ twice and not including a self interaction. The first term in equation \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 the pair potential depending 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 clearly unphysical because the presence of molecule C {\em will} modify the interaction between molecules A and B. The correction do to the presence of molecule C can be included in out 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 then pair-pair interactions and that the contribution of higher order terms is generally small, pair-pair interactions are included as the main term for calculating the long-range potential. The pairwise potential can 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-pair interaction.\cite{Allen87} These functional forms are typically fit from empirical data and quantum mechanical calculations. Effective pair potentials work well for describing the interactions in the noble gases progressively well in sold, liquid and gas phases. As the system becomes more dense and polarizable, the effective pair approximation begins to fail 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 weekly do 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 form shown in Figure \ref{fig:ljform}. There is an attractive region at long-range where $-\partial V/\partial R < 0$ and becomes steeply repulsive at small $R$ to account for compressibility in condensed phase material. 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 form of the pair-pair interaction energy between molecule where $-\epsilon$ is the minimum energy, $\sigma$ is the close interaction distance and $R_m$ is the minimum energy seperation.}
        \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 (except for atom pairs which are involved in an intra-molecular interaction). It should be noted that atoms in this context refer to long-range 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 atomic distances. To reduce the number of distance evaluations between pairs of atoms, {\sc oopse} uses a switched cutoff distance beyond which interactions are not calculated. Further, a Verlet neighbor list is implemented which creates a look-up of interactions that are within a specific cutoff distance.\cite{Allen87} A complication of using cutoff distances with neutral groups which contain charges is that the system will exhibit pathological forces unless the cutoff is applied to the neutral groups evenly instead of to the individual long-range interaction sites.\cite{Leach01} {\sc oopse} allows users to specify cutoff groups which may contain an arbitrary number of interaction sites with-in a molecule. Thus, members in a cutoff group are treated as a single unit for the purpose of evaluating the switching function: 
\begin{equation}
        V_{\mathrm{long-range}} = \sum_{a} \sum_{b>a} s(r_{ab}) \sum_{i \in a} \sum_{j \in b} V_{ij}(r_{ij}), 
\end{equation}
where $r_{ab}$ is the distance between the centers of mass of the two cutoff groups ($a$ and $b$). The sums over $a$ and $b$ are over the cutoff groups that are present in the simulation. Atoms which are not explicitly defined as members of a {\tt cutoffGroup} are treated as a group consisting of only one atom. In simulations containing only Lennard-Jones atoms, the cutoff radius has a default value of $2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest Lennard-Jones length parameter present in the simulation. In simulations containing charged or dipolar atoms, the default cutoff radius is $15 \mbox{\AA}$. 


\section{\label{introSec:LJPot}THE LENNARD-JONES FORCE FIELD}

The most basic force field implemented in {\sc oopse} 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} 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 cross term parameters for $\sigma$ and $\epsilon$ to be generated. These parameters are 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 comprised elements that have a covalently dominated binding mechanism. One of the first attempts to construct a model for 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 gasses, the Drude model assumes that electrons behave as a 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 valance electrons are more weekly bound and form an electron gas when the isolated metallic atoms condense to form a metal. This electron gas is typically referred to as the conduction band in a condensed system. Despite using relatively simplistic assumptions, the Drude model provides a reasonable explanation for AC and DC electrical conductivity, Hall effect and electron thermal conductivity.\cite{Ashcroft:1976zt}

Sommerfeld subsequently extended the Drude model by incorporating quantum effects into the model.\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 both in intrinsic physical properties and, do the Heisenberg uncertainty principle, with respect to the individual position and momentum an electron possesses. 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. $E_{F}$ can be written as a function of the number density of electrons, 
\begin{equation}
        E_{F} = \frac{\hbar^2}{2m}\left(3\pi^2 Z\rho\right)^{3/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 describe numerous phenomenon in metallic systems. These include Wiedemann-Franz law relating electrical and thermal conductivity, temperature dependence of the heat capacity, and the electronic density of states. 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 do to the oversimplification of the role of ions and the fundamental nature 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 valance 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_0(\mathbf{r}))$ 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 us a universal function of the electron density.\cite{Puska:1981p1755} They found that for the noble gasses, 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 do 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 et al. involving differences in elastic constants for cubic lattice solids composed of noble gasses and fcc metals.\cite{DAW:1993p1640} It is known for a solid described by a purely pairwise interaction that the elastic constants $C_{12}$ and $C_{44}$ are equal and is known as the Chauchy relation. We can also define a Chauchy pressure as the difference between the two elastic constants, $P_c=C_{12}-C_{44}$.\cite{Kittel:1996fk} Thus, the ideal ratio of $C_{12}:C_{44}$ is $1$. Both Ar, and Kr have ratios very close to $1$, but fcc metals have ratios closer to 2. Both Pt, and Au have ratios that are between 3 and 4. This can be seen in Table \ref{tab:manybod} reproduced from Table 1 in Daw et al. Also included in this 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}
        \caption{ DEVIATION FROM PAIRWISE BEHAVIOR IN FCC METALS COMPARED WITH NON-METALLIC SYSTEMS. REPRODUCED FROM REFERENCE \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 Gasses &&&\\
%\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}
\label{tab:manybod}
\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 as in Equation \ref{eq:metpseudo} of separating metallic interactions into those do to ion-pair interactions and those do 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 energy of a system of atoms can be written as a 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}} = G[\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 $\mathbf{R}_i$ represents the charge and position of the $i$th nucleus. The collective energy of the isolated atoms is given by $E_{\mathrm{atoms}}$ and $G[\rho]$ is the kinetic, exchange, and correlation energy functional. 

If we assume that $G[\rho]$ can be described by 
\begin{equation}
        G[\rho] = \int g(\rho(\mathbf{r}),\nabla\rho(\mathbf{r}),\nabla^{2}\rho(\mathbf{r}),\ldots) \mathrm{d}\mathbf{r}
\end{equation}  
where $g$ 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 gasses. 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})$. Lastly, we define an embedding energy for placing an atom 
$i$ in an constant density electron gas to be
\begin{equation}
 G_i(\bar\rho_i) = G[\rho_i^a + \bar\rho_i] - G[\rho_i^a] - G[\bar\rho_i].
\end{equation}
Using the two assumptions and this definition for the embedding energy, Daw and Baskes derived an expression for the cohesive energy in metals as

\begin{equation}
        V  =  \sum_{i} G_{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 then 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 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}\right)\left(\frac{Z_i^a(R)Z_j^a(R)}{R}
        \right)
        \label{eq:eamphi}
\end{equation}
where the $Z^a(R)$ are effective screened charges of the nuclei 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$,$\nu$ are determined empirically to give the desired screening to $Z(R)$. Likewise, $F[\rho]$ is defined by a universal function that relates the sublimation energy of a metal to the lattice constant.\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.aAs an example do to Voter, assume 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^2|\chi_j\right> \left<\chi_j|H^2|\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 functions set. One can show that for any basic shape for 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 $\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}, to obtain an interatomic-potential 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$ do 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 constant with the trend in metallic bond energies. This particular embedding function also has a negative slope and therefore is the total source of the 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} do 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 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}
$a$ 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 experimental lattice parameter, cohesive energy and bulk modulus the phonon frequencies, vacancy formation energy and surface energies.\cite{kimura-quantum} The new parameterization was named Quantum Sutton-Chen ({\sc Q-SC}) because the new parameterization takes into account zero-point energy. The functional form for {\sc SC} and {\sc Q-SC} is computationally advantages because it somewhat resembles the Lennard-Jones potential presented in Section \ref{introSec:LJPot}. Parameters for the various parameterizations of Sutton-Chen  are presented in Table \ref{introtab:qsc}.

\begin{table}
        \caption{ {\sc SC} and {\sc Q-SC } PARAMETERS\cite{kimura-quantum,Chen90} TABLE 1.}
        \centering
\begin{tabular}{llccccc}
        \toprule
        \toprule
 & & n & m & $\epsilon$(eV) & c & $\alpha$ a(\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}
\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 for a given system, 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$) for each particle in the system. 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 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}}$ 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 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}
Newtons 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{r}}^2/2$, gives
\begin{equation}
\frac{d}{dt}(m\dot{\mathbf{r}})+\frac{dV}{d\mathbf{r}}=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:verletdlm}VERLET AND DLM METHODS OF INTEGRATION}
The evolution of time in Molecular Dynamics requires an integrator that is symplectic and preserves certain properties of phase space during the transformation.\cite{McQuarrie:2000yt,Tolman:1938kl} 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{p})
\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{p} = -\frac{\partial H}{\partial q}
        \label{eq:cannon1}
\end{equation}
and
\begin{equation}
        \dot{q} = -\frac{\partial H}{\partial p}.
        \label{eq:cannon2}
\end{equation}
Liouville's theorem describes phase space density,$\rho$ at a point of interest changes as a function of time
\begin{equation}
        \left(\frac{\partial \rho}{\partial t}\right)_{p,q}=-\sum_i \left(\frac{\partial \rho}{\partial q_i}\dot{q_i}+\frac{\partial \rho}{\partial p_i}\dot{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)_{p,q}=-\sum_i \left(\frac{\partial \rho}{\partial q_i}\frac{\partial H}{\partial p_i}
        +\frac{\partial \rho}{\partial p_i}\frac{\partial H}{\partial 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 q_i}\frac{\partial N}{\partial p_i}
        -\frac{\partial N}{\partial q_i}\frac{\partial M}{\partial p_i}\right).
        \label{eq:poisson}
\end{equation}
Equation \ref{eq:liouville2} is then
\begin{equation}
        \left(\frac{\partial \rho}{\partial t}\right)_{q,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{q}},\mathbf{q},t) = e^{-iLt}\rho(\mathbf{\dot{q}},\mathbf{q},0).
\end{equation}


\subsection{LANGEVIN DYNAMICS}

\section{PARALLEL MOLECULAR DYNAMICS} 
Revision:	3484
Committed:	Tue Jan 13 14:39:51 2009 UTC (15 years, 8 months ago)
Content type:	application/x-tex
File size:	40205 byte(s)
Log Message:	This commit was manufactured by cvs2svn to create tag 'start'.