\documentclass[17pt,landscape,dvips]{foils}

\usepackage{amsmath,amsthm,amssymb,amscd}
\usepackage{array}
\usepackage{alltt}
\usepackage{epsf}
\usepackage{latexsym}
\usepackage{times}
\usepackage{bookman}
\usepackage{graphicx}
\usepackage{color}
\usepackage{epsfig}
\usepackage{psfrag}
\usepackage[ps2pdf]{hyperref}
\usepackage{url}
\usepackage{fancyvrb}
\usepackage{multirow}
\usepackage{epic,eepic}
\usepackage{pst-all}

%\def\basswidth{4.0in}
%\def\gutter{0.2in}

%\newenvironment{comment}{\setbox0=\vbox\bgroup}{\egroup}
%\newenvironment{singlespace}{\def\baselinestretch{1.0}\large\normalsize}{\par}

\setlength{\textheight}{8in}
\setlength{\textwidth}{9in}
\voffset=-0.5in


\newcommand{\myblue}[1]   {{\color{blue} #1}}
\newcommand{\mygreen}[1]  {{\color{green} #1}}
\newcommand{\myred}[1]    {{\color{red} #1}}
\newcommand{\mybluebf}[1] {{\blue{\bf{#1}}}}
\newcommand{\mygreenbf}[1]{{\green{\bf{#1}}}}
\newcommand{\myredbf}[1]  {{\red{\bf{#1}}}}
\newcommand{\stitle}[1]{\foilhead{\myred{#1}}\vspace{-0.5in}}
\newcommand{\bigitem}[2]
           {\item{#1\begin{small}\begin{itemize}#2\end{itemize}\end{small}}}
\newcommand{\mysum} {\displaystyle \sum}

\title{Algorithmic Techniques for Atomistic Modeling}
 
\author{
  Hasan Metin Aktulga\\
  Dept of Computer Science\\
  Purdue University\\
  \texttt{haktulga@cs.purdue.edu}
}
\date{May 14, 2009}

\begin{document}
\maketitle

{\bf ayg: chage green to a darker color}

\stitle{Outline}
\begin{enumerate} \itemsep 17pt
  \item \mybluebf{Introduction: MD Methods \& Reactive Force Fields}
  
  \bigitem{\mybluebf{Our Implementation of Reactive Force Fields (ReaxFF)}}{ \itemsep 7pt
    \item General structure of a ReaxFF implementation
    \item Implementation details of major components
    \item Additional features
  }

  \bigitem{\mybluebf{Applications using the Serial Version}}{ \itemsep 7pt
    \item Hexane simulations: Validation and performance analysis
    \item Water-Silica Surface Interaction
    \item Si/Ge Nanobar Strain Rates
  }

  \item \mybluebf{Ongoing \& Future Work}
  
  \bigitem{\mybluebf{Parallelization of ReaxFF}}{ \itemsep 7pt
    \item Challenges
    \item Our solution approaches 
  }
  
  \item \mybluebf{Current Userbase}
\end{enumerate}


\stitle{Introduction: Atomistic Modeling Methods}
\begin{itemize} \itemsep 17pt
  \bigitem{\mybluebf{\emph{Ab-initio} methods:} {\small Atomistics with electronic degrees of freedom}}{ \itemsep 7pt
    \item{Hartree-Fock(HF) methods $\rightarrow$ misses $E_{corr}$}
      
      \marginpar{$\uparrow$\\
        {\tiny accuracy\\ complexity\\ detail\\ comp. demand\\}
	\vspace{1in}
               {\tiny approx's\\ sys size\\sim. frame\\}
	       $\downarrow$}
      
    \bigitem{Post-HF methods}{ \itemsep 5pt
      \item $E_{corr}$ incorporated at great computational expense
      \item semi-empirical methods to the rescue
    }
     
    \bigitem{DFT-based methods}{ \itemsep 5pt
      \item completely different approach but similar derivations to HF theory
      \item bigger systems, longer simulation times made possible
      \item CPMD is the most popular example
    }
  }

  \bigitem{\mybluebf{Classical MD methods:}}{ \itemsep 7pt
      \item many approximations: no electronic d.o.f.
      \item electronic effects are mimiced through parametrizations
      \item static bonds, no reactions!
      \item systems at nanoscale, simulation times upto hundreds of ns
   }

  \bigitem{\mybluebf{Continuum Mechanics:} Macroscale systems}{ \itemsep 7pt
      \item additional approximations, no atomistic detail!
      \item modeled using \myblue{PDE}'s.
  }
\end{itemize}


\stitle{Reactive Force Fields (ReaxFF): Bridging the Gap}
\noindent
\begin{tabular}{|m{2in}|m{3in}|m{3in}|}
\hline
      & \textbf{ReaxFF} & \textbf{Classical MD} \\ 
\hline \hline
\multirow{3}{*}{\textbf{advantages}} & \small{reactive} & \small{non-reactive} \\ \cline{2-3}
      &\small{dynamic bonds} & \small{static bonds} \\ \cline{2-3}
      &\small{polarization with QEq} & \small{fixed charges in general (except for polarizable FF)} \\
\hline \hline
\multirow{3}{*}{\textbf{challenges}}& \small{dynamic interaction lists}& \small{static lists}\\\cline{2-3}
      &\small{complex \& costly enegies,forces}&\small{much simpler energy \& force formulas}\\\cline{2-3}
      &\small{frequent update of charges (expensive!)} & \small{static charges}\\ \cline{2-3}
      &\small{shorter timesteps ($\approx0.25$~fs)} & \small{longer timesteps ($1-10$~fs)} \\
\hline
\end{tabular}
\vspace{1in}
\noindent \myblue{Interesting applications}: large system with reactions
and charge-transfer\\ 
{\small Simulation of fuel cells, silica crack propogation, corrosion of silica in water, etc.} 

\stitle{General Flowchart of Conventional MD Programs}
\noindent \small
\begin{figure}
  \centering
    \input{md_flowchart.tex}
  \label{fig:md_flowchart}
\end{figure}
\normalsize


\stitle{ReaxFF Flowchart}
\noindent \small
\begin{figure}
  \centering
    \input{reax_flowchart.tex}
  \label{fig:reax_flowchart}
\end{figure}

{\bf ayg: What about valence corrections?}

\normalsize


\stitle{Implementation: Neighbor Generation}
\begin{itemize} \itemsep 17pt
  \bigitem{3 different neighbor lists:}{ \itemsep 7pt
    \item \myblue{near\_nbrs} for bonded forces $\rightarrow$ 
      \mygreen{bond\_cut} $\approx$ 4-5~\AA, full matrix stored
    \item \myblue{hbond\_list} for hydrogen bonds $\rightarrow$ 
      \mygreen{hbond\_cut} $\approx$ 6-7.5~\AA, only for $H$ 
    \item \myblue{far\_nbrs} for non-bonded forces $\rightarrow$ 
      \mygreen{nonb\_cut} $\approx$ 10~\AA, upper-half only
  }

  \bigitem{Bin atoms into 3D grid cells}{
    \item grid cell dims $\approx$ $\frac{1}{2}$\mygreen{nonb\_cut}
  }

  \item Verlet lists with delayed re-neighboring \myred{not} implemented $\rightarrow$ 
    little benefit

  \item Compressed adjacency list representation
\end{itemize}

\noindent
\begin{figure}
  \centering
    \input{nbr_list.tex}
    \input{3D_grid.tex}
%    \includegraphics[height=1in]{3D_grid.eps}
  \label{fig:nbrs}
\end{figure}


\stitle{Implementation: Computing Forces and Potentials}
\begin{itemize} \itemsep 17pt
  \bigitem{\myblue{Bonded Interactions:}Similar to classical MD but accounts for dynamic bonds}{ \itemsep 7pt
    \item \myred{precursor:} bond orders
    \item lone-pair energy, over/undercoordination energies
    \item bond energy
    \item valence energy (with penalty \& 3-body conjugation corrections)
    \item dihedral energy (with 4-body conjugation correction)
  }

  \bigitem{\myblue{Hydrogen Bonds}}{ \itemsep 7pt
    \item \myred{precursor:} bond orders
    \item $H$ covalently bonded to $X$ and interacting with $Z$
    \item can be considered a \myblue{bonded interaction}
  }

  \bigitem{\myblue{Non-bonded Interactions}}{ \itemsep 7pt
    \item \myred{precursor:} charge equilibration (QEq)
    \item electrostatic (Coulomb) energy, van der Waals energy
  }
\end{itemize}
\noindent \small
  \begin{eqnarray}
    E_{system} &=& E_{bond} + E_{lp} + E_{over} + E_{under} + E_{val} + E_{pen} + E_{3conj} \notag \\
    & & {} + E_{tors} + E_{4conj} + E_{H-bond} + E_{vdW} + E_{Coulomb} \notag
\end{eqnarray}
\normalsize


\stitle{Bonded Interaction: Bond Orders}
\begin{itemize} \itemsep 17pt
  \bigitem{Prior to \myblue{bonded forces}, compute \myblue{bond orders} based on the new \mygreen{near\_nbrs}}{
    \itemsep 7pt
    \item \mygreen{bond\_list:} subset of \mygreen{near\_nbrs}, stored in the same way 
    \item \myred{uncorrected} bond orders and derivatives
    \item store both bo$(i,j)$ and bo$(j,i)$ $\rightarrow$ 
      efficient construction of \myblue{angles}, \myblue{dihedrals}
    \item compute bo$(i,j)$ only if $i<j$, otherwise bo$(i,j)$=bo$(j,i)$ 
  }
\end{itemize}
\tiny
  \begin{eqnarray}
      \text{BO}'_{ij}&=&\text{BO}_{ij}^{\sigma'}+\text{BO}_{ij}^{\pi'}+\text{BO}_{ij}^{\pi\pi'} \notag\\
      \text{BO}_{ij}^{\alpha'}(r_{ij})&=& 
        exp\left[a_\alpha\left(\frac{r_{ij}}{{r_0}_\alpha}\right)^{b_\alpha}\right] \notag\\
      \text{BO}_{ij}&=&\text{BO}'_{ij} \cdot f_1(\Delta'_i,\Delta'_j) 
        \cdot f_4(\Delta'_i,\text{BO}'_{ij}) \cdot f_5(\Delta'_j,\text{BO}'_{ij}) 
      \text{ where $\Delta'_i$ is the valency of atom $i$.} \notag
  \end{eqnarray}
%\framebox[5in]{
%\begin{verbatim}
%  if i < j 
%    compute bond_order(i,j)
%  else bond_order(i,j) = bond_order(j,i)
%\end{verbatim}
%}
\begin{figure} \centering
  \includegraphics[height=2.5in, width=3in]{bond_order.eps}
  \label{fig:bond_orders}
\end{figure}
\normalsize


\stitle{Bonded Interaction: Bond Energy}
\begin{eqnarray}
\text{E}_{\text{bond}} &=& -D_e^\sigma \cdot \text{BO}_{ij}^\sigma \cdot 
             exp\left\{p_{be1}\left(1-{\left(\text{BO}_{ij}^\sigma \right)}^{p_{be2}}\right)\right\} \notag\\
        & & -D_e^\pi \cdot \text{BO}_{ij}^\pi - D_e^{\pi\pi} \cdot \text{BO}_{ij}^{\pi\pi} \notag
\end{eqnarray}
\begin{itemize} \itemsep 19pt
  \item the \myred{stronger} the bond, the \myred{lower} the associated energy
  \item sweep over the \mygreen{bond\_list}
  \item compute \myblue{bond\_energy} between $i$, $j$ only if $i<j$
\end{itemize}


\stitle{Bonded Interaction: Lone-Pair \& Over/UnderCoordination}
\begin{itemize} \itemsep 19pt
  \bigitem{\myblue{Lone-pair energy}}{ \itemsep 10pt
    \item $\Delta_i^{lp} = n_{opt}^{lp} - n_i^{lp}$
    \item energy associated with \myred{unpaired} electrons of an atom $\rightarrow$
      \myred{zero} for a fully coordinated atom
    \item single-body interaction $\rightarrow$ just sweep over \mygreen{atom\_list} 
  }
  \vspace{0.5in}
  \bigitem{\myblue{Over/undercoordination energy}}{ \itemsep 10pt
    \item \myred{ideal} \# of bonds = \# of valence electrons
    \item $\Delta_i = \mysum_{j \in nbrs(i)}{\text{bo}(i,j)- Val_i}$
    \item \mygreen{actual} \# of bonds $>$ \myred{ideal} \# of bonds ($\Delta_i>0$) $\rightarrow$ 
      over-coordination
    \item \mygreen{actual} \# of bonds $<$ \myred{ideal} \# of bonds ($\Delta_i<0$) $\rightarrow$ 
      under-coordination
    \item \mygreen{actual} \# of bonds $=$ \myred{ideal} \# of bonds ($\Delta_i=0$) $\rightarrow$ 
      \myred{no} over/undercoordination energy
    \item functionals of $\Delta_i$ and $\Delta_j$'s $\rightarrow$ just sweep over \mygreen{atom\_list} 
  }
\end{itemize}


\stitle{Bonded Interaction: Valence Angle Energy}
\tiny
\begin{eqnarray}
\text{E}_{\text{val}} &=& f_7(\text{BO}_{ij}, p_{val3}, p_{val4}) \cdot 
           f_7(\text{BO}_{jk}, p_{val3}, p_{val4}) \cdot 
           f_8(\Delta_j, p_{val5}, p_{val6}, p_{val7}) \cdot  \notag \\
       & & \left( p_{val1} - p_{val1} \cdot 
           exp \left\{ -p_{val2} \cdot {\left( \Theta_0 - \Theta_{ijk} \right)}^2 \right\} \right) \notag
\end{eqnarray}
\small
\begin{itemize} \itemsep 14pt
  \bigitem{\myblue{$\Theta_0$} is the \myred{ideal} angle, 
    \myblue{$\Theta_{ijk}$} is the \myred{actual} angle}{
    \item the \myred{closer} the $\Theta_{ijk}$ to $\Theta_0$, the \myred{lower} the energy
  }
  
  \item $f_7(\text{BO}_{ij},\ldots)$ and $f_7(\text{BO}_{jk},\ldots)$ 
    ensure $\text{E}_{\text{val}}\rightarrow0$ 
    as $\text{BO}_{ij}\rightarrow0$ or $\text{BO}_{jk}\rightarrow0$

  \item $\{\forall x,y\in\text{bo\_list}_i|x<y\}$ that meet certain criteria,
    compute the energy of $<x,i,y$

  \item $\{\forall x,y\in\text{bo\_list}_i|x>y\}$ copy $<y,i,x$ into $<x,i,y$ (\myred{for dihedrals!})

  \item $\text{E}_{\text{pen}}$ and $\text{E}_{\text{3conj}}$ account for corrections in special cases
\end{itemize}
\tiny
\begin{figure} \centering
  %\includegraphics{valence_angle.eps}
  \input{3body_list.tex}
  \label{fig:valence_angle}
\end{figure}
\normalsize


\stitle{Bonded Interaction: Dihedral Energy}
\begin{eqnarray}
E_{tors} &=& \frac{1}{2} \cdot f_{10}(BO_{ij}, BO_{jk}, BO_{kl}, p_{tor2}, 1) \cdot 
            sin\Theta_{ijk} \cdot sin\Theta_{jkl} \cdot V_{123}(\omega_{ijkl}) \notag
\end{eqnarray}
\begin{itemize} \itemsep 17pt
  \item \myblue{Dihedral angle $\omega_{ijkl}$} is the angle 
    between planes defined by positions of $i,j,k$ and $j,k,l$ 

  \item $f_{10}(\text{BO}_{ij},\text{BO}_{jk},\text{BO}_{kl},\ldots)$ 
    ensure $\text{E}_{\text{tors}}$ vanishes smoothly as any of these bonds dissociate

  \item $sin\Theta_{ijk}$ and $sin\Theta_{jkl}$ ensure that $\text{E}_{\text{tors}}\rightarrow0$ 
    as $\Theta_{ijk}\rightarrow0$ or $\Theta_{jkl}\rightarrow0$
    
  \item $\{\forall i,j,k,l \in \text{atoms} | j<k, <i,j,k \in \text{3body\_list}_{jk}, 
                                              <j,k,l \in \text{3body\_list}_{kj}\}$
     compute the energy associated with $\omega_{ijkl}$

  \item \myred{weak but very important} in determining the 3D structures 

  \item \myred{no} higher order interactions $\rightarrow$ \myred{no} storage necessary
\end{itemize}


\stitle{Hydrogen Bonds}
\scriptsize
\begin{equation}
\text{E}_{XHZ} = p_{hb1} \cdot f_7(\text{BO}_{XH}, p_{hb2}, 1) \cdot 
           sin^4 \left( \frac{\Theta_{XHZ}}{2} \right) \cdot
                  exp\left\{ -p_{hb3} \cdot \left( \frac{r_{hb}^0}{r_{HZ}} + \frac{r_{HZ}}{r_{hb}^0} - 2 \right)
                 \right\} \notag
\end{equation}
\normalsize
\begin{itemize} \itemsep 17pt
  \bigitem{Constraints of a hydrogen bond:}{ \itemsep 7pt
    \item Middle atom must be $H$
    \item $X,Z$ must be one of $N,O,P,F$
    \item $X-H$ covalently bonded, $Z\in\text{hbond\_list}_H$
  }
  
  \item $f_7(\text{BO}_{XH},\ldots)$ ensure $\text{E}_{\text{val}}\rightarrow0$ 
    as the \myred{covalent bond} breaks 
    
  \item $sin^4(\frac{\Theta_{XHZ}}{2})$  maximized when $\Theta_{XHZ}=\pi$ ensures \myred{alignment on a line}

  \item crucial for accurately describing \myred{water}, \green{DNA structure}, 
    \blue{secondary structures in proteins}, etc.
\end{itemize}


\stitle{Nonbonded Interaction: Charge Equilibriation (QEq)}

Minimizing electrostatic energy by redistributing (partial) charges.
\tiny
\begin{eqnarray}
& &\text{Minimize   } 
E(Q_1 \ldots Q_N) = \mysum_A (E_{A0} + \chi_A^0Q_A + \frac{1}{2}J_{AA}^0Q_A^2) +
                    \mysum_{A<B}(J_{AB}Q_AQ_B) \notag \\ [-0.5mm]
\notag \\ [-0.5mm]
& &\text{subject to   } 
Q_{net} = \mysum_{i=1}^{N} Q_i \notag
\end{eqnarray}

\normalsize
\begin{itemize} \itemsep 17pt
  \bigitem{Solve the optimization problem using \myblue{the method of Lagrange multipliers}}{
    \item gives a sparse linear system of equations for finding charges
  }
  
  \bigitem{\myblue{GMRES} with restarts, $GMRES(50)$}{
    \item heavy diagonal $\rightarrow$ diagonal preconditioner
    \item little configurational change between steps 
      $\rightarrow$ initial guess $q_t = \text{linear\_extrapolation}(q_{t-1},q_{t-2})$  
  }

  \bigitem{Implemented \myblue{GMRES} both with $MGS$ and $Householder$ orthogonalizations}{
    \item virtually no difference
  }

  \bigitem{Implemented \mygreen{CG} for comparison}{
    \item \myblue{GMRES} takes fewer number of \myred{matvecs} and is faster than \mygreen{CG}
  }

  \bigitem{Choose \myblue{tolerance} for \myblue{the norm of the relative residual} \myred{carefully}:}{
    \item too high $\rightarrow$ wrong results!
    \item too low $\rightarrow$ \myblue{QEq} dominates the total computation time!
    \item more on this later \ldots
  }
\end{itemize}


\stitle{Nonbonded Interactions: Coulomb \& van der Waals Energy} 
\tiny
\begin{eqnarray}
E_{Coulomb} &=& C \cdot Tap(r_{ij}) \cdot \frac{q_i \cdot q_j}{\left[r_{ij}^3+\gamma_{ij}^{-3}\right]^{\frac{1}{3}}}
\notag\\
E_{vdWaals} &=& Tap(r_{ij}) \cdot D_{ij} \cdot \notag\\
           & & \,\left[ exp\left\{\alpha_{ij}\cdot\left( 1-\frac{f_{13}(r_{ij})}{r_{vdW}} \right) \right\} - 
                        2 \cdot exp\left\{\frac{1}{2} \cdot \alpha_{ij} \cdot 
                        \left(1-\frac{f_{13}(r_{ij})}{r_{vdW}} \right) 
                        \right\} \right] \notag
\end{eqnarray}
\normalsize

\begin{itemize} \itemsep 17pt
  \item \myblue{Shielding} prevents energies from increasing drastically at close distances

  \bigitem{Long range interactions with \myred{cutoffs}}{
    \item \myblue{Taper} term ensures smooth vanishing of energies after the cutoff 
  }

  \item \myred{no $1-2,1-3 \text{ or } 1-4$ exclusions} $\rightarrow$ smooth bond forming/breaking

  \bigitem{Takes up a \myred{large portion} of the total computation time}{ \itemsep 7pt
    \item \myblue{tabulate} long range energy \& forces
    \item \mygreen{linear interpolation}
    \item large table $\rightarrow$ good approximations, reasonable memory usage
    \item more on our gains later \ldots
  }
\end{itemize}


\stitle{Summing Alltogether: Net Force}
\begin{itemize} \itemsep 17pt
  \item Let $E_i$ be the sum of energies from all interactions involving atom $i$

  \item Let $r_i$ denote the position of atom $i$

  \item $F_i = \frac{\partial E_i}{\partial r_i}$

  \bigitem{\myred{Problem:}}{ \itemsep 9pt
    \item bonded energy expressions include $\text{BO}_{ij}(\text{BO}'_{ij},\Delta'_i, \Delta'_j)$ terms 
    \item $\frac{\partial \text{BO}_{ij}}{\partial r_k}$ arise in every bonded interaction 
    \item $\frac{\partial \text{BO}_{ij}}{\partial r_k} =c_1\cdot \frac{\partial \text{BO}'_{ij}}{\partial r_k}+
           c_2\cdot \frac{\partial \Delta'_i}{\partial r_k}+c_3\cdot \frac{\partial \Delta'_j}{\partial r_k}$
    \item $\frac{\partial \text{BO}_{ij}}{\partial r_k}\neq 0,\forall k \in \text{bonds}_i \cup \text{bonds}_j$,
      \myred{huge memory overhead!}
      
    \item even if we choose to store them, 
      \myred{very time consuming} to compute each single bonded interaction! 
  }

\end{itemize}


\stitle{Summing Alltogether: Net Force - Solution}
\begin{itemize} \itemsep 17pt
    \item \myblue{idea:} distribution law of multiplication over summation
    \item let $C_0,\ldots,C_n$ be the coefficients of $\frac{\partial \text{BO}_{ij}}{\partial r_k}$
      arising in different interactions
    \item re-write $\mysum_t{C_t \times \frac{\partial \text{BO}_{ij}}{\partial r_k}}$ as
      $\frac{\partial \text{BO}_{ij}}{\partial r_k} \times \mysum_t{C_t}$
    \item while computing interactions, accumulate $C_t$'s in $C_{ijk}$
    \item \myred{delay} computation of $\frac{\partial \text{BO}_{ij}}{\partial r_k} \times \mysum_t{C_t}$ and 
      \myblue{forces} due to them until $C_{ijk}$'s are determined
    \item \myred{no additional storage}, \mygreen{important savings in CPU time}  
\end{itemize}


\stitle{Implementation: Additional Features}
\begin{itemize} \itemsep 19pt
  \bigitem{Modular implementation}{
    \item a different force field can be adopted by plugging-in new interaction routines 
  }

  \item NVE, NVT and NPT ensembles
{\bf ayg: explain these terms first and what it takes to implement them}
  
  \item Compressible custom trajectory format

  \bigitem{Tools for performing common \myblue{analysis}}{ \itemsep 9pt
    \item detection of \myred{reactions} (\myblue{on-the-fly})
    \item \myred{property calculations} such as drift coefficient, dipole moment (\myblue{on-the-fly})
    \item \myred{distributions} of bond lengths, strengths, valence angles, charges, etc.
      (\myblue{over the trajectory file})
  }
\end{itemize}


\stitle{Applications}
\begin{itemize} \itemsep 19pt
  \bigitem{\myblue{Hexane Simulations: Validation and Performance Analysis}}{ \itemsep 9pt
    \item Preparation of systems
    \item Effects of QEq tolerance on accuracy and performance
    \item Effects of tabulating long range interactions on accuracy and performance
    \item Hexane structure verification
    \item Scalabilty of ReaxFF compared to \emph{ab-initio} and classical MD
  }

  \item{\myblue{Corrosion of silica surface in water} 
    \tiny(in collaboration with Dr Pandit's group)\normalsize}

  \item{\myblue{Measuring the strain tensor of $Si$/$Ge$ nanobar}
    \tiny(in collaboration with Dr Strachan's group)\normalsize}
  
\end{itemize}


\stitle{Hexane Simulations: Preparation od the Systems}
\begin{itemize} \itemsep 19pt
  \item \mygreen{Hexane:} $C_6H_{14}$, hydrocarbon, constituent of gasoline
  
  \bigitem{Initial configuration setup:}{ \itemsep 7pt
    \item \myblue{Very large box} compared to the ideal volume $\rightarrow$ \myred{reduces overlaps}!
    \item \myred{Randomly spread copies} of a model hexane molecule
    \item \myred{Rotations} of the model molecule around $x$,$y$ and $z$ axis to increase randomness
    \item Various system sizes for scalability analysis ($343$,$512$,$1000$,$1728$, $3375$ molecules)}{ 
  }

\bigitem{Energy minimization and NPT simulations using \myblue{Gromacs}}{ \itemsep 7pt
  \item brings the systems to the \myred{ideal volume} quickly
  \item \mygreen{output} to be used for \myblue{ReaxFF} studies and scalability analysis
}

\bigitem{Energy minimization and NVT equilibration using \myblue{ReaxFF}}{ \itemsep 7pt
  \item added $H$ to \myblue{Gromacs} output conf using the \myblue{Avogadro} program
  \item energy minimization for $2.5$~ps
  \item NVT equilibration at $200$~K for $2.5$~ps
  \item QEq tolerance set to $1e-8$ to be safe
}
\end{itemize}


\stitle{Hexane Simulations: Effect of QEq Tolerance on Accuracy}
\begin{itemize} \itemsep 17pt
  \item \myred{Chosen:} hexane system with $343$ molecules = $6860$ atoms

  \item \myred{Restart} from the system equilibrated at $200$~K

  \bigitem{\myred{How to determine the ``right'' tolerance:}}{ \itemsep 5pt
    \item observe how the same system evolves over time at \myblue{different QEq tolerances}
    \item pick \mygreen{the highest one} with \mygreen{reasonable accuracy}
    \item $\text{tol}_1=1e-3$,$\text{ tol}_2=1e-4$,$\text{ tol}_3=1e-8 \rightarrow \text{control run}$
  }
\end{itemize}
\begin{figure} 
  \centering
  \includegraphics[height=2in,width=7in]{hex_343.ideal.qeq_energy.eps}
  \includegraphics[height=2in,width=7in]{hex_343.ideal.qeq_temperature.eps}
\end{figure}


\stitle{Hexane Simulations: Effect of Tabulation on Accuracy}
\small
\noindent 
tol$_2=1e-4$ looks good enough, now \myblue{turn on tabulation of long range interactions}, too!
\begin{figure} 
  \centering
  \includegraphics[height=1.5in,width=7in]{hex_343.ideal.tab_energy.eps}
  \includegraphics[height=1.5in,width=7in]{hex_343.ideal.tab_temperature.eps}
\end{figure}
\noindent 
\mygreen{More in depth comparison} shows they are \myred{almost identical}:
\tiny
\begin{center}
\begin{tabular}{l||c|c|c}
\hline
\textbf{property}&\textbf{tol$_3=1e-8$}&\textbf{tol$_2=1e-4$}&\textbf{tol$_2=1e-4$ with opt.} \\
\hline
{C-H bond}   & $1.09\pm0.01$  & $1.09\pm0.01$ & $1.09\pm0.01$ \\ 
{C-C bond}   & $1.57\pm0.01$  & $1.57\pm0.01$ & $1.57\pm0.01$ \\ 
{$<$C-C-C}   & $108.0\pm2.9$  & $107.9\pm2.9$ & $108.0\pm2.9$ \\
{$<$C-C-H}   & $111.0\pm0.0$  & $111.0\pm0.0$ & $111.0\pm0.0$ \\
{$<$H-C-H}   & $106.6\pm0.0$  & $106.6\pm0.0$ & $106.6\pm0.0$ \\
{$q_{\text{C-tip}}$} & $-0.171$ & $-0.171$      & $-0.171$      \\
{$q_{\text{C-mid}}$} & $-0.080$ & $-0.080$      & $-0.080$      \\
{$q_{\text{H-tip}}$} & $0.040$  & $0.040$       & $0.040$       \\
{$q_{\text{H-mid}}$} & $0.040$  & $0.040$       & $0.040$       \\
\hline
\end{tabular}
\end{center}
\normalsize


\stitle{Hexane Simulations: Profiling Analysis \& Scalability}
\begin{itemize}
 \item Different \myblue{qeq tolerances}, with/without \mygreen{tabulation}
 \item Used the $head.cs$ cluster \myred{except for the first case}
\end{itemize}
\small

{\bf ayg: What are the numbers here.. are they times in seconds?}
\begin{center}
  \begin{tabular}{ | l | c  c  c  c  c  c  c | }
    \hline
                       & total  & neighbors & bonded & nonb    & QEq    & matvecs & QEq\%\\ \hline
    tol=$1e-4$ w/opt\S & $0.93$ & $0.22$    & $0.14$ & $0.22$  & $0.34$ & $9.9$   & 37\% \\
    tol=$1e-4$ w/opt   & $3.26$ & $0.47$    & $0.32$ & $0.65$  & $1.78$ & $9.9$   & 55\% \\
    tol=$1e-4$         & $4.41$ & $0.45$    & $0.31$ & $1.92$  & $1.67$ & $9.9$   & 38\% \\
    tol=$1e-6$ w/opt   & $3.35$ & $0.45$    & $0.31$ & $0.59$  & $1.97$ & $13.6$  & 59\% \\
    tol=$1e-6$         & $4.65$ & $0.45$    & $0.31$ & $1.92$  & $1.95$ & $13.5$  & 42\% \\
    tol=$1e-8$         & $7.56$ & $0.44$    & $0.30$ & $1.91$  & $4.89$ & $46.8$  & 65\% \\
    \hline
  \end{tabular}
\end{center}
\small
\noindent
{\bf ayg: Dont say Dell Studio.. instead, say what the processor is, etc.}
\S \myred{Architecture can make a huge difference}. Same system with same parameters 
on a Dell Studio XPS with $2.67GHz$ quad-core $i7$ processor and $1066MHz$ memory.

\normalsize
\begin{itemize}
 \bigitem{\myred{Lessons learnt:}}{ \itemsep 7pt
   \item QEq tolerance is crucial for accuracy
   \item arbitrarily large QEq tolerance might cause QEq domination
{\bf ayg: What does the above bullet mean?}
   \item QEq is just a precursor to electrostatics yet takes up at least one third of total time!
   \item QEq must be improved to make \myblue{ReaxFF} scalable.
 }
\end{itemize}

\stitle{Hexane Simulations: Validation}
\noindent Compare the structure of our hexane molecules 
to those of \myblue{experimental results} in the literature and \mygreen{\emph{ab-initio} simulation}:  
\small
\begin{center}
\begin{tabular}{l||c|c|c}
\hline
\textbf{property}& \textbf{ReaxFF}&\textbf{experimental}\dag&\textbf{\emph{ab-initio}}\ddag\\
\hline
{C-H bond}   & $1.09\pm0.01$  & $1.118\pm0.006$ & $1.100$ \\ 
{C-C bond}   & $1.57\pm0.01$  & $1.533\pm0.003$ & $1.533$ \\ 
{$<$C-C-C}   & $108.0\pm2.9$  & $111.9\pm0.4$   & $114.2$ \\
{$<$C-C-H}   & $111.0\pm0.0$  & $109.5\pm0.5$   & $109.5$ \\
{$<$H-C-H}   & $106.6\pm0.0$  & $NA$            & $106.5$ \\
{$q_{\text{C-tip}}$} & $-0.171$ & $NA$            & $-0.205$\\
{$q_{\text{C-mid}}$} & $-0.080$ & $NA$            & $0.033$ \\
{$q_{\text{H-tip}}$} & $0.040$  & $NA$            & $0.047$ \\
{$q_{\text{H-mid}}$} & $0.040$  & $NA$            & $-0.10\sim0.10$\\
\hline
\end{tabular}
\end{center}

\scriptsize
\noindent
\dag R. A. Bonham, L. S. Bartell, and D. A. Kohl. 
``The Molecular Structures of n-Pentane, n-Hexane and n-Heptane''
J. Am. Chem. Soc., $1959$, $81$ ($18$), $4765-4769$

\noindent
\ddag geometry optimization of an isolated hexane 
using CPMD v$3.13.2$ with $PBE$ Troullier-Martins pseudopotentials
\normalsize


\stitle{Hexane Simulations: Scalability Analysis}
\begin{figure} 
  \centering
  \includegraphics[height=3in,width=7in]{systemsize_speed_chart.eps}
  \includegraphics[height=3in,width=7in]{systemsize_memory_chart.eps}
  \label{fig:qeq_tol}
\end{figure}


\stitle{Silica Surface Corrosion in Water} 


\stitle{$Si$/$Ge$ Nanobar} 


\stitle{Ongoing \& Future Work}
\begin{itemize} \itemsep 17pt
  \bigitem{\myblue{ParallelReax}}{ \itemsep 7pt
    \item complete \& verify the implementation
    \item a large scale application (maybe with the PRISM device)
  }

  \bigitem{\mygreen{Integration into LAMMPS} {\small(a well-known, widely used MD package from SNL)}}{ 
    \bigitem{QEq integrated independently}{ \itemsep 7pt
      \item opens the door for polarizable-ff inside LAMMPS
      \item compatibility issues to be sorted out!
    }
  }

  \bigitem{\myred{Better solvers for QEq}}{ \itemsep 7pt
    \item block Jacobi type pre-conditioner
    \item inner-outer schemes (use a lower cutoff inner solve to precondition
    outer solve)
    \item use a fast multipole-type preconditioner.
  } 

  \bigitem{\myred{Tight relations} among items on the agenda}{ \itemsep 7pt
    \item better QEq solvers \mygreen{necessary} for scalable ParallelReax
    \item better QEq solvers \mygreen{necessary} for QEq in LAMMPS
    \item completion of ParallelReax \mygreen{necessary} for Reax in LAMMPS
  } 
\end{itemize}


\stitle{Parallelization of ReaxFF}
\begin{itemize} \itemsep 17pt
  \bigitem{A \mygreen{draft version} of \myblue{ParallelReax}}{ \itemsep 5pt
    \bigitem{domain decomposition technique}{
      \item \textbf{repeat:} get my share of atoms $\rightarrow$ communicate boundaries 
        $\rightarrow$ compute forces $\rightarrow$ move my atoms
    }
    \item \myred{not fully verified}
    \item \myred{inefficient} handling of \myblue{processor boundaries}
  }
\end{itemize}

\noindent \myred{Two big challenges:}
\begin{itemize} \itemsep 17pt
  \bigitem{Parallelization of QEq}{ \itemsep 7pt
    \item even \myblue{CG} needs at least 4 communications per iteration!
    \item \myblue{parallel GMRES} (with \emph{Householder orth.}) is even worse
    \item \myblue{QEq} will \myred{dominate even more}
    \item \myred{definitely need:} better solvers for QEq
  }

  \bigitem{Processor boundaries}{ \itemsep 7pt
    \item avoid double computation at boundaries!
    \item avoid thick boundaries, retain accuracy for any system! 
  }
\end{itemize}


\stitle{Parallelization: Our Solution Approaches}
\begin{itemize} \itemsep 15pt
  \bigitem{\myred{Parallelization of QEq:} better solvers for QEq}{ \itemsep 5pt
    \item block Jacobi type pre-conditioner
    \item inner-outer schemes
  }

  \bigitem{\mygreen{Avoid double computation}}{ \itemsep 5pt
    \bigitem{coordination through \myred{the mid-point rule}}{ \itemsep 7pt
      \item \myblue{bond$(i,j)$:} owner$(\frac{1}{2}(r_i + r_j))$
      \item \myblue{$<(i,j,k)$:} owner$(r_j)$
      \item \myblue{dihedral$(i,j,k,l)$:} owner$(\frac{1}{2}(r_j + r_k))$
      \item \myblue{hbond$(X,H,Z)$:} owner$(r_H)$
      \item \myblue{nonbonded$(i,j)$:} owner$(\frac{1}{2}(r_i + r_j))$
    }
  }
\end{itemize}
\tiny
\begin{figure} 
  \centering
  \input{boundary_coordination.tex}
\end{figure}
\normalsize

\stitle{Parallelization: Our Solution Approaches}
\begin{itemize} \itemsep 17pt
  \bigitem{\myblue{Why avoid thick boundaries?}}{ \itemsep 5pt
    \item a $20$~\AA~cubic box $\approx 1000$ atoms
    \item assume $nonb\_cut$ long boundaries $\rightarrow$ a few thousand atoms/processor 
      $\rightarrow$ upto a couple of seconds per iteration!
    \item \myred{definitely need: thin boundaries}
  }

  \bigitem{\mygreen{How to avoid thick boundaries?}}{ \itemsep 5pt
    \item \myred{Mid-point rule} to the rescue!
    \item boundary thickness: $max(\frac{3}{2}bond\_cut, hbond\_cut, \frac{1}{2}nonb\_cut)$
    \item gets better if no hydrogen bonds present: $max(\frac{3}{2}bond\_cut, \frac{1}{2}nonb\_cut)$
  }
\end{itemize}
\tiny
\begin{figure} 
  \centering
  \input{thin_boundary.tex}
\end{figure}
\normalsize


\stitle{Conclusions}


\stitle{Current Userbase}
\begin{itemize} \itemsep 19pt
  \bigitem{Dr Pandit's group at USF}{
    \item silica-water systems
  }

  \bigitem{Dr Strachan's group at Purdue}{
    \item $Si$/$Ge$ nanobar
  }

  \bigitem{Dr Buehler's group at MIT}{
    \item Silica cracking with strain
  }

  \item{Dr van Duin at PennState}

  \item{Dr Goddard's group at Caltech}

  \item{Dr Aluru's group at UIUC}
\end{itemize}

\end{document}
