\documentclass[12pt]{article}
\usepackage{fullpage,doublespace,graphicx}
\setlength{\parindent}{0em}
\setlength{\parskip}{0.35em}
\setstretch{1.5}

\title{Second-Moment Potential for Metallic ${\rm ZrC}_{1-X}$}
\date{\today}
\author{Ju Li}
\begin{document}
\maketitle

\makeatletter
\def\thebibliography#1{\section{References\@mkboth
  {REFERENCES}{REFERENCES}}\list
 {[\arabic{enumi}]}{\settowidth\labelwidth{[#1]}\leftmargin\labelwidth
   \advance\leftmargin\labelsep
   \usecounter{enumi}}
    \def\newblock{\hskip .11em plus .33em minus .07em}
   \sloppy\clubpenalty4000\widowpenalty4000
     \sfcode`\.=1000\relax}
\makeatother
\tableofcontents

\section{Motivation}

Let us briefly review the main assumptions of the second-moment
potential

Under the tight-binding


${\rm ZrC}_{1-X}$ is a semi-metal.

\begin{equation}
V_{\rm tot} = \sum_i V_i,
\end{equation}

Willaime-Massobrio (WM)\cite{Willaime89,Willaime91} potential for Zr:
\begin{equation}
V_{\rm i} = \left \{ a \sum_{r_{ij} < r_c} \exp \left [ - p \left ( 
\frac{r_{ij}}{r_0} - 1 \right ) \right ] \right \} - 
\left \{ \sum_{r_{ij} < r_c} \xi^2 \exp \left [ - 2q \left (
\frac{r_{ij}}{r_0} - 1 \right ) \right ]
\right \}^{1/2}
\label{WM_form}
\end{equation}
where the parameters are listed in Table \ref{WM_parameters}.

\begin{table}[th]
\begin{center}
\caption{Parameters of WM potential (\ref{WM_form}) for Zr}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
$a[{\rm eV}]$ & $r_c$[{\rm\AA}] & $p$ & $r_0$[{\rm\AA}] & 
$\xi[{\rm eV}]$ & $q$ \\
\hline
$0.179364$ & $6.8$ &$9.3$  & $3.1744$ & $2.2014548$ & $2.1$ \\
\hline
\end{tabular}
\end{center}
\label{WM_parameters}
\end{table}

For AB binary system, one can show that
\begin{equation}
a \rightarrow a_{AB}=a_{BA},\;\; p \rightarrow p_{AB}=p_{BA}, \;\;
\xi\rightarrow \xi_{AB}, \;\; q \rightarrow q_{AB}=q_{BA},
\end{equation}
where one can show the AB commutivity except for $\xi_{AB}$.

\section{Zr/C Binary Potential}

(\ref{WM_form}) can be massaged into a slightly simpler form for
numerical calculations
\begin{equation}
V_{\rm i} = \left \{ \sum_{r_{ij} < r_c} \exp 
\left [ A ( B - r_{ij} ) \right ] \right \} - 
\left \{ \sum_{r_{ij} < r_c} \exp \left [ C ( D - r_{ij} ) \right ]
\right \}^{1/2},
\label{WM_massaged_form}
\end{equation}
where 
\begin{equation}
A \equiv \frac{p}{r_0}, \;\; 
B \equiv r_0 \left (1 + \frac{\log(a)}{p} \right ),\;\;
C \equiv \frac{2q}{r_0}, \;\; 
D \equiv r_0 \left (1 + \frac{\log(\xi)}{q} \right ).
\end{equation}

\begin{table}[th]
\begin{center}
\caption{Parameters of massaged WM potential (\ref{WM_massaged_form}) for Zr}
\begin{tabular}{|c|c|c|c|}
\hline
$A_{\rm ZrZr}$[\rm\AA$^{-1}$] & 
$B_{\rm ZrZr}$[\rm\AA] & 
$C_{\rm ZrZr}$[\rm\AA$^{-1}$] & 
$D_{\rm ZrZr}$[\rm\AA] \\
\hline
$2.9296875$ & $2.58787395638939$ & $1.32308467741935$ & $4.3672464262563$ \\
\hline
\end{tabular}
\end{center}
\label{WM_massaged_parameters}
\end{table}

And the generalization of (\ref{WM_massaged_form}) to binary Zr/C
would be
\begin{equation}
V_{\rm i} = \left \{ \sum_{r_{ij} < r_c} \exp 
\left [ A_{\alpha_i\alpha_j} ( B_{\alpha_i\alpha_j} - r_{ij} ) \right
] \right \} - 
\left \{ \sum_{r_{ij} < r_c} \exp \left [ C_{\alpha_i\alpha_j} 
( D_{\alpha_i\alpha_j} - r_{ij} ) \right ]
\right \}^{1/2},
\label{WM_massaged_binary_form}
\end{equation}
where $\alpha_i$, $\alpha_j$ denote the chemical species of atom $i$,
$j$. 

The parameters for Zr-Zr interaction is already listed in Table
\ref{WM_massaged_parameters}. By adding carbon, I have now introduced
$8$ new parameters: $A_{\rm CC}$, $B_{\rm CC}$, $C_{\rm CC}$, $D_{\rm
CC}$, $A_{\rm ZrC}$=$A_{\rm CZr}$, $B_{\rm ZrC}$=$B_{\rm CZr}$,
$D_{\rm ZrC}$, $D_{\rm CZr}$. The decay rate of the hopping integral
is taken to be the algebraic average: $C_{\rm ZrC}$=$C_{\rm
CZr}$=($C_{\rm CC}$+$C_{\rm ZrZr}$)/2. And I use a uniform
cutoff distance of $6.8$\AA.

\begin{eqnarray}
&& \forall i\in{\rm Zr}: \nonumber \\
&& \;\;\;\;\;\;
V_i = \sum_{r_{ij}<r_c}^{j\in{\rm Zr}} \exp 
\left [ A_{\rm ZrZr} ( B_{\rm ZrZr} - r_{ij} ) \right ]
+ \sum_{r_{ij}<r_c}^{j\in{\rm C}} \exp 
\left [ A_{\rm ZrC} ( B_{\rm ZrC} - r_{ij} ) \right ] 
\;\;-\;\; X_i^{1/2},
\nonumber \\
&& \;\;\;\;\;\; 
X_i \equiv
\sum_{r_{ij}<r_c}^{j\in{\rm Zr}} \exp \left [ C_{\rm ZrZr} 
( D_{\rm ZrZr} - r_{ij} ) \right ] 
+ \sum_{r_{ij}<r_c}^{j\in{\rm C}} ( c(r_{ij}) p_j S + 1)
\exp \left [ C_{\rm ZrC}( D_{\rm ZrC} - r_{ij}) \right ], 
\nonumber \\ 
\label{new_potential_Zr_part}
\\
&& \forall i\in{\rm C}: \nonumber \\
&& \;\;\;\;\;\;
V_i = \sum_{r_{ij}<r_c}^{j\in{\rm C}} A_{\rm CC}
\left \{ \exp\left [ 2B_{\rm CC} ( C_{\rm CC} - r_{ij} ) \right ] - 
2 \exp\left [ B_{\rm CC} ( C_{\rm CC} - r_{ij} ) \right ] \right \}.
\label{new_potential_C_part}
\end{eqnarray}
where $r_c$ is a universal and rather long-ranged cutoff ($6.8$\AA).
$A_{\rm ZrZr}$, $B_{\rm ZrZr}$, $C_{\rm ZrZr}$, $D_{\rm ZrZr}$,
$A_{\rm ZrC}$, $B_{\rm ZrC}$, $S$, $C_{\rm ZrC}$, $D_{\rm ZrC}$,
$A_{\rm CC}$, $B_{\rm CC}$, $C_{\rm CC}$ are all constant
parameters. $c(r)$ in $X_i$ definition is a short-ranged cutoff
function to include only nearest neighbor C-Zr bonds:
\begin{equation}
c(r) = \left \{
\begin{array}{ll}
\exp \left ( \frac{U}{r-W} \right ), \;\;\;\;& r < W \\
0, & r \ge W \end{array},
\right . 
\label{short_cutoff_function}
\end{equation}
and $S$ is a constant describing the {\em additional} strength of
nearest neighbor C-Zr bonds with strong covalent flavor ($pd\sigma$),
compared to other C-Zr interactions of metallic character but which
are longer ranged ($pd\pi$ etc.). Largeness of $S\exp[C_{\rm ZrC}(
D_{\rm ZrC} - r_{ij})]$ compared to $\exp[C_{\rm ZrZr} ( D_{\rm ZrZr}
- r_{ij})]$ from Zr-Zr $dd$ hopping is the dominant reason that
rock-salt ZrC has the exceptionally large heat of formation of $2.14$
eV/pair\footnote{To put it in perspective, heat of formation of SiC is
only $0.68$ eV/pair. The large heat of formation is the reason for
ZrC's chemical inertness.}. However whether the strong $pd\sigma$ bond
could form between $i$-$j$ largely depends on the availability of
$p$-orbital at the carbon atom $j$ pointing in $i$'s direction, which
prefers to be orthogonal in Hilbert space to other occupied local
orbitals of $j$. This can done in pure carbon without requiring the
bonds to be mutually orthogonal in {\em real space} because of $sp$
hybridization\footnote{For example in diamond, the hybridized orbitals
$|h_{[111]}\rangle=|s\rangle+|p_x\rangle+|p_y\rangle+|p_z\rangle,\;\;$
$|h_{[1\bar{1}\bar{1}]}\rangle=
|s\rangle+|p_x\rangle-|p_y\rangle-|p_z\rangle,\;\;$
$|h_{[\bar{1}1\bar{1}]}\rangle=
|s\rangle-|p_x\rangle+|p_y\rangle-|p_z\rangle,\;\;$
$|h_{[\bar{1}\bar{1}1]}\rangle=
|s\rangle-|p_x\rangle-|p_y\rangle+|p_z\rangle$ are mutually orthogonal
in Hilbert space but not in real space.}, but without $s$-orbital it
cannot be done. Therefore to receive the benefit of strong $pd\sigma$,
it is necessary for the C-Zr bonds to be mutually orthogonal around C
in real space, and I use $p_j$ to characterize that degree of mutual
orthogonality which should be $1$ for rock-salt but penalizes
$pd\sigma$ covalent strength if otherwise. Thus in the extreme case of
$p_j\rightarrow 0$ the benefit of covalent bonding should be all but
destroyed and C-Zr interactions become mostly metallic. Furthermore,
the angular dependence of $p_j$ should stiffen up $C_{44}$ greatly at
rock-salt structure, reversing the character of positive
$C_{12}-C_{44}$ of most transition metals\footnote{$(C_{12}-C_{44})/B$
is called the Cauchy discrepancy. It is more often positive in metals
and more often negative in ceramics. When modeled by pair potentials,
it is zero. Conventional EAM potentials tend to give positive Cauchy
discrepancies.}.

\begin{figure}
\centerline{\includegraphics[scale=0.5]{angular_function.eps}}
\caption{}
\label{angular_function}
\end{figure}

I propose a form for $p_j$ as follows,
\begin{equation}
p_j = \exp\left \{  - P \sum_{k<l}^{k,l\in {\rm Zr}} c(r_{jk}) c(r_{jl}) 
\left [ \cos^3(\theta_{kjl}) + \cos^2(\theta_{kjl}) \right ] \right
\}, \;\; j\in{\rm C},
\end{equation}
where $P$ is a constant, $c(r_{jk})$, $c(r_{jl})$ are cutoff functions
in (\ref{short_cutoff_function}) to select only those Zr atoms within
nearest neighbor bonding distance to $i$. The angular function is
shown in Fig. \ref{angular_function}. It is the simplest polynomial in
$\cos\theta$ that has minima $0$ at both $90^\circ$ and $180^\circ$,
the bond angles in rock-salt structure.

When $k$ (a Zr atom) moves by $dx_{k}^{\mu}$,
\begin{eqnarray}
dp_j = - p_j P \sum_{l\ne k}^{l\in {\rm Zr}} && c^\prime(r_{jk}) 
\hat{r}_{jk}^\mu dx_{k}^{\mu} c(r_{jl}) 
\left [ \cos^3(\theta_{kjl}) + \cos^2(\theta_{kjl}) \right ] 
\nonumber \\
- && c(r_{jk}) c(r_{jl}) 
\left [ 3 \cos^2(\theta_{kjl}) + 2 \cos(\theta_{kjl}) \right ]
\sin (\theta_{kjl}) d\theta_{kjl}.
\end{eqnarray}

Define vectors
\begin{equation}
{\bf m} \equiv \hat{r}_{jl}\times \hat{r}_{jk}, \;\;\;\;
{\bf w} \equiv {\bf m}\times \hat{r}_{jk}.
\end{equation}
Clearly, $|{\bf m}|$=$|{\bf w}|$=$|\sin\theta_{kjl}|$.  Furthermore,
once one picks a sense of $\theta_{kjl}$, there is
\begin{equation}
{\bf w} \cdot \frac{d{\bf x}_k}{r_{jk}} = 
w^\mu \frac{dx_{k}^{\mu}}{r_{jk}} = 
\sin (\theta_{kjl}) d\theta_{kjl}.
\end{equation}
Therefore,
\begin{eqnarray}
\frac{dp_j}{dx_{k}^{\mu}}
= P p_j \sum_{l\ne k}^{l\in {\rm Zr}} &&
- c^\prime(r_{jk}) c(r_{jl})
\left [ \cos^3(\theta_{kjl}) + \cos^2(\theta_{kjl}) \right ] 
\hat{r}_{jk}^\mu \nonumber \\
+ && c(r_{jk}) c(r_{jl})
\left [ 3 \cos^2(\theta_{kjl}) + 2 \cos(\theta_{kjl}) \right ]
\frac{w^\mu}{r_{jk}}.
\end{eqnarray}

If one plugs in (\ref{short_cutoff_function}),
\begin{equation}
c^\prime(r) = - c(r) \frac{U}{(r-W)^2}, 
\end{equation}
there is then
\begin{eqnarray}
\frac{dp_j}{dx_{k}^{\mu}}
= P p_j \sum_{l\ne k}^{l\in {\rm Zr}}
c(r_{jk}) c(r_{jl})  && \left \{
\left [ \cos^3(\theta_{kjl}) + \cos^2(\theta_{kjl}) \right ] 
\frac{U\hat{r}_{jk}^\mu}{(r_{jk}-W)^2}
\right . \nonumber \\
&& \left .  + \left [ 3 \cos^2(\theta_{kjl}) + 2 \cos(\theta_{kjl})
\right ] \frac{w^\mu}{r_{jk}} \right \}, \;\;k\in{\rm Zr}.
\label{angular_other_derivative}
\end{eqnarray}
In the above ${\bf w}$ is defined using ${\bf m}$=${\bf
l}$$\times$${\bf k}$. When we run over $k<l$ pairs computationally, we
can calculate ${\bf m}$ just once, but {\em reversing} the sign of
${\bf m}$ when calculating the ${\bf w}$ for $l$. 

Lastly, because there is overall translational invariance in $p_j$,
there is,
\begin{equation}
\frac{dp_j}{dx_{j}^{\mu}} \;\;=\;\; -\sum_{r_{jk}<W}^{k\in {\rm Zr}}
\frac{dp_j}{dx_{k}^{\mu}}.
\end{equation}
$j$ is the only carbon atom that feels bond-bending force due to
$p_j$.

Computationally one is organized as follows: one first allocates space
and computes $p_j$ for all carbon atoms by first sorting out and then
running a double loop over all nearest neighbor Zr atoms. One then
computes $X_i$ for all Zr atoms, and carries out force evaluation as
if in a common EAM procedure, treating $p_j$ as constants. During this
process one may multiply a universal smooth cutoff function
\begin{equation}
L(r) = \left \{ \begin{array}{ll}
\exp \left ( \frac{K}{r-r_c} \right ), \;\;\;\;& r < r_c \\
0, & r \ge r_c \end{array},
\right . 
\label{long_cutoff_function}
\end{equation}
to all $r < r_c$ summed terms in (\ref{new_potential_Zr_part}) and
(\ref{new_potential_C_part}), which is missing in the original WM
potential. Including (\ref{long_cutoff_function}) takes very little
effort because
\begin{equation}
 L(r) \exp \left [ A ( B - r ) \right ] \;=\;
\exp \left [ A ( B - r ) + \frac{K}{r-r_c} \right ],
\end{equation}
and if one takes the first derivative, one simply modifies the
additional coefficient
\begin{equation}
A\;\;\rightarrow\;\; A + \frac{K}{(r-r_c)^2}.
\end{equation}

Lastly, one needs to collect the angular forces. One again runs over
all carbon atoms, each having a $p_j$. Clearly, by varying a certain
$p_j$, the total energy changes by,
\begin{equation}
\frac{dV_{\rm tot}}{dp_j} \;\;=\;\;
\sum_{i}^{i \in{\rm Zr}} \frac{d(-X_i^{\frac{1}{2}})}{dp_j}
\;\;=\;\; - \sum_{r_{ij}<W}^{i \in{\rm Zr}}
\frac{ c(r_{ij}) S \exp \left [ C_{\rm ZrC}( D_{\rm ZrC} - r_{ij})
\right ]} {2X_i^{\frac{1}{2}}} \equiv -q_j.
\label{q_j_definition}
\end{equation}
One needs to first compute the above common factor for $j$, then go
over the double loop again to obtain (\ref{angular_other_derivative})
for each Zr atom $k$,
\begin{eqnarray}
&& \Delta F_k^\mu \;\;=\;\; \left ( -\frac{dV_{\rm tot}}{dp_j} \right ) \cdot 
\frac{dp_j}{dx_{k}^{\mu}} 
\;\;=\;\; P p_j q_j \sum_{l\ne k}^{l\in {\rm Zr}} c(r_{jk})
c(r_{jl}) \times
\nonumber \\
&& \left \{ \left [ \cos^3(\theta_{kjl}) +
\cos^2(\theta_{kjl}) \right ] \frac{U\hat{r}_{jk}^\mu}{(r_{jk}-W)^2}  
+ \left [ 3 \cos^2(\theta_{kjl}) + 2
\cos(\theta_{kjl}) \right ] \frac{w^\mu}{r_{jk}} \right \},
\label{angular_other_force}
\end{eqnarray}
for all $k\in{\rm Zr}$, $r_{jk}<W$. And finally, the carbon atom
itself takes angular force of
\begin{equation}
\Delta F_j^\mu \;\;=\;\; -\sum_{r_{jk}<W}^{k\in {\rm Zr}} \Delta
F_k^\mu. \label{angular_self_force}
\end{equation}

\section{Lattice Heat Current}

EAM formulation leads naturally to a partition of the system total
energy:
\begin{equation}
E_{\rm tot} = \sum_i E_i,
\end{equation}
where each particle collects its own kinetic energy and a site
potential energy,
\begin{equation}
E_i = K_i + V_i, \;\;\;\; K_i \equiv \frac{1}{2} m_i v_i^2,
\end{equation}
with $V_i$ following the definition of (\ref{new_potential_Zr_part})
and (\ref{new_potential_C_part}). There is then
\begin{equation}
{\bf J} = \sum_i E_i {\bf v}_i + \sum_i \dot{E}_i {\bf r}_i.
\end{equation}
The first summation is simple enough to evaluate, but the second
summation has to be {\em consistent} with the first. That is, how the
total energy is partitioned into $E_i$'s can be arbitrary, but once it
is partitioned, everything has to follow the book.

We notice that
\begin{equation}
\dot{E}_i = \dot{K}_i + \dot{V}_i,
\end{equation}
and because there is
\begin{equation}
\dot{K}_i = {\bf F}_i \cdot {\bf v}_i,
\label{delta_K}
\end{equation}
and
\begin{equation}
\dot{V}_i = - \sum_j \Delta {\bf F}_j(V_i) \cdot {\bf v}_j,
\label{delta_V}
\end{equation}
($\Delta {\bf F}_j(V_i)$ means force on $j$ due to $V_i$), $\dot{E}_i$
can be decomposed into a {\em linear superposition} of contributions
from subgroups of forces. That is, we should follow the subgrouping of
interatomic forces, and be aware that each subgroup will influence
{\em both} $\dot{K}_i$'s and $\dot{V}_i$'s to influence {\bf J}. Since
$G$ usually stems from a single well-defined source, it often
satisfies
\begin{equation}
\sum_i \Delta {\bf F}_i(G) = 0.
\end{equation}

Consider a repulsive site energy located at atom $i$, for example a
$V_R(r_{ij})$ contribution in
\begin{equation}
\dot{V}_i = \sum_j V_R(r_{ij}) \;+\; ...
\end{equation}
summation due to atom $j$. Let us see how a force subgroup, in this
case of just two members, $\{\Delta{\bf F}_i(G),\Delta{\bf
F}_j(G)\}$ with $\Delta{\bf F}_i(G) = -\Delta{\bf
F}_j(G)$, influences ${\bf J}$ by influencing the
$\dot{K}$'s and $\dot{V}$'s.

In this case only three energy derivatives are influenced by $G$:
$\dot{K}_i$, $\dot{V}_i$, $\dot{K}_j$. Furthermore there is {\em
microscopic} energy conservation,
\begin{equation}
\Delta \dot{K}_i(G) + \Delta \dot{V}_i(G) + \Delta \dot{K}_j(G) = 0,
\end{equation}
since
\begin{eqnarray}
\Delta \dot{K}_i(G) = && \Delta{\bf F}_i(G)\cdot {\bf v}_i,
\nonumber \\
\Delta \dot{V}_i(G) = && - \Delta{\bf F}_i(G)\cdot {\bf v}_i - 
\Delta{\bf F}_j(G)\cdot {\bf v}_j,
\nonumber \\
\Delta \dot{K}_j(G) = && \Delta{\bf F}_j(G)\cdot {\bf v}_j.
\end{eqnarray}
Therefore, the total heat current is influenced as
\begin{equation}
\left ( \Delta \dot{K}_i(G) + \Delta\dot{V}_i(G) \right ) {\bf r}_i + 
\Delta\dot{K}_j(G) {\bf r}_j = 
\Delta\dot{K}_j(G) ({\bf r}_j - {\bf r}_i) =
\left ( \Delta{\bf F}_j(G) \cdot {\bf v}_j) \right )
({\bf r}_j - {\bf r}_i).
\end{equation}
This can be easily generalized to the case of $V_i$ depending on
multiple $j$'s, and the problem of evaluating $\sum_i \dot{E}_i {\bf
r}_i$ contribution due to site potential $V_i$ is really quite simple:
{\em one just sets the coordinate origin to be at $i$} to avoid the
complex evaluation of $\left ( \Delta \dot{K}_i(G) +
\Delta\dot{V}_i(G) \right ) {\bf r}_i$, and only evaluates $\Delta
\dot{K}_j(G) {\bf r}_j$ for $j\ne i$, which is $(\Delta{\bf
F}_j(G)\cdot {\bf v}_j)({\bf r}_j - {\bf r}_i)$. One may even formally
include $i$ in the summation because ${\bf r}_i-{\bf r}_i=0$ and it
does not matter. Therefore, the heat current due to $V_i$, in addition
to the implicit contribution in $\sum_i E_i {\bf v}_i$, is,
\begin{equation}
\Delta {\bf J}(V_i) \;=\; \sum_j (\Delta{\bf F}_j(V_i)\cdot {\bf v}_j)
({\bf r}_j - {\bf r}_i),
\label{site_energy_heat_current}
\end{equation}
where $\Delta{\bf F}_j(V_i)$ is the force increment to $j$ due to
$V_i$. Note that this result is only true for {\em site} energy form
$V_i$, and only when $V_i$ is used consistently in $E_i$ also of
$\sum_i E_i {\bf v}_i$.  Nevertheless (\ref{site_energy_heat_current})
is a fundamental result because other form of potential energy can
always be decomposed into several site energies, and then recombined,
and if possible, simplified, using (\ref{site_energy_heat_current}).

For example, starting out with a pair potential $V_P(r_{ij})$,
\begin{equation}
V_{\rm tot} \;\;=\;\; \sum_{i<j} V_P(r_{ij}) \;+\; ...
\end{equation}
instead of site energies leaves the job to oneself. One must remember
to attribute $V_P(r_{ij})/2$ to $E_i$ and $V_P(r_{ij})/2$ to $E_j$,
which needs to be reflected in $\sum_i E_i {\bf v}_i$. One then gets
{\em two} contributions from the above site formulation as
\begin{equation}
\Delta {\bf J} (V_P) = 
(\Delta {\bf F}_j (G_i) \cdot {\bf v}_j)({\bf r}_j - {\bf
r}_i) + (\Delta {\bf F}_i (G_j) \cdot {\bf v}_i)({\bf r}_i - {\bf
r}_j).
\end{equation}
Because we have given exactly {\em half} of $V_P$ to $G_i$ and half to
$G_j$, there is
\begin{equation}
\Delta {\bf F}_j (G_i) \;=\; \frac{1}{2} {\bf
F}_{ij} (V_P) \;=\; -\Delta {\bf F}_i (G_j),
\end{equation}
where ${\bf F}_{ij}(V_P)$ is the {\em total} force $j$ feels due to
$V_P(r_{ij})$, we then have
\begin{equation}
\Delta {\bf J} (V_P) = \frac{1}{2}
\left [ {\bf F}_{ij}(V_P)\cdot ({\bf v}_i + 
{\bf v}_j) \right ] ({\bf r}_j - {\bf r}_i),
\label{pair_heat_current}
\end{equation}
in agreement with old results.

Because pair counting is more efficient than site counting by saving
two terms of equal value from being computed twice,
(\ref{pair_heat_current}) should be used for repulsion between like
atoms in (\ref{new_potential_Zr_part}) and
(\ref{new_potential_C_part}), just as energies and forces are.  For
Zr-C repulsion, we have a choice. As the way it is in
(\ref{new_potential_Zr_part}) and (\ref{new_potential_C_part}),
everything is attributed to Zr site energy with nothing to C: we then
should use (\ref{site_energy_heat_current}) for Zr
site. Alternatively, we may also distribute half and half to Zr and C,
respectively, and use (\ref{pair_heat_current}) instead; but that
change must {\em also} be reflected in site energies
(\ref{new_potential_Zr_part}) and (\ref{new_potential_C_part}) at
least computationally, and consequently in $\sum_i E_i {\bf v}_i$.

EAM can use pair counting saving just like pair potential does. One
considers $ij$ pair instead of sites individually. However each $ij$
pair would now concern two site energies -- two objects for
simultaneous consideration rather than a single object of
$V_P(r_{ij})$.  Specifically $r_{ij}$ now influences two embedding
energies, $-X_i^{1/2}$ and $-X_j^{1/2}$. Call the first's force
contribution to atom $i$ and $j$: $\Delta{\bf F}_i(X_i)$ =
$-\Delta{\bf F}_j(X_i)$; the second's force contribution to $i$ and
$j$: ${\bf F}_i(X_j)$ = $-\Delta{\bf F}_j(X_j)$, then the heat current
invoked by the two embedding energies is
\begin{eqnarray}
\Delta {\bf J}(r_{ij},X_i,X_j) = && 
\left (\Delta{\bf F}_j(X_i)\cdot {\bf v}_j \right ) 
({\bf r}_j - {\bf r}_i) + 
\left (\Delta{\bf F}_i(X_j)\cdot {\bf v}_i \right ) 
({\bf r}_i - {\bf r}_j) \nonumber \\
= && \left ( \Delta{\bf F}_j(X_i)\cdot {\bf v}_j
+ {\bf F}_j(X_j)\cdot {\bf v}_i \right ) 
( {\bf r}_j - {\bf r}_i)
\nonumber \\
= && \left [ {\bf r}_{ij} \cdot 
\left ( \frac{\Delta F_j(X_i)}{r_{ij}}{\bf v}_j + 
\frac{\Delta F_j(X_j)}{r_{ij}}{\bf v}_i \right ) \right ] {\bf r}_{ij}.
\end{eqnarray}
At this point there can be no further simplification because $\Delta
F_j(X_i)$ bears no relation to $\Delta F_j(X_j)$, unlike the case of a
half-half site-split pair potential.

\begin{figure}
% \centerline{\includegraphics[scale=0.5]{angular_heat_current.eps}} %
\caption{}
\label{angular_heat_current}
\end{figure}

The heat current due to the angular forces of (\ref{q_j_definition})
to (\ref{angular_self_force}) is a bit tricky to account for. See
Fig. \ref{angular_heat_current}: the angular forces that $j$ and $k$,
$k^\prime$, $k^{\prime\prime}$, ... feel in fact come from several
site energies: $i$, $i^\prime$, $i^{\prime\prime}$, ... as $p_j$
controls the covalency of several bonds that reflects in the
(\ref{q_j_definition}) summation. In order to use
(\ref{site_energy_heat_current}), one should know exactly how much is
$\Delta{\bf F}_j(X_i)$ and $\{\Delta{\bf F}_k(X_i)\}$, and not simply
the total angular forces $\Delta{\bf F}_j$ and $\{\Delta{\bf F}_k\}$
due to $p_j$. Fortunately, a simple scheme can save a lot of
computation as follows: let us define normalized weights $\{w_i\}$:
\begin{equation}
\sum_i w_i \;\;=\;\; 1, \;\;\;\; i \in{\rm Zr},\; r_{ij} < W,
\end{equation}
for each $i$ appearing in (\ref{q_j_definition}) summation:
\begin{equation}
w_i \;\equiv\; \frac{1}{q_i}\cdot \frac{ c(r_{ij}) S \exp \left [
C_{\rm ZrC}( D_{\rm ZrC} - r_{ij}) \right ]} {2X_i^{\frac{1}{2}}}.
\end{equation}
Clearly there is proportional distribution,
\begin{equation}
\Delta{\bf F}_j(X_i) = w_i \Delta{\bf F}_j, \;\; 
\Delta{\bf F}_k(X_i) = w_i \Delta{\bf F}_k, 
\;\;...
\end{equation}
Therefore the total heat current is 
\begin{eqnarray}
&&\Delta {\bf J}(p_j,\{X_i\}) \nonumber \\
= && \sum_i \left \{ \left ( \Delta{\bf
F}_j(X_i)\cdot {\bf v}_j \right ) ({\bf r}_j - {\bf r}_i) + 
\sum_k \left (\Delta{\bf
F}_k(X_i)\cdot {\bf v}_k \right ) ({\bf r}_k - {\bf r}_i) \right \} 
\nonumber \\ 
= && \sum_i \left \{ w_i \left ( \Delta{\bf F}_j \cdot {\bf v}_j
\right ) ({\bf
r}_j - {\bf r}_i) + \sum_k w_i \left ( \Delta{\bf F}_k \cdot {\bf v}_k
\right ) ({\bf
r}_k - {\bf r}_i) \right \} \nonumber \\ 
= && \sum_i \left \{ w_i \left (\Delta{\bf F}_j \cdot {\bf v}_j \right
) ({\bf r}_j - {\bf r}_i) + \sum_k w_i \left (\Delta{\bf F}_k 
\cdot {\bf v}_k \right ) ({\bf
r}_k - {\bf r}_j + {\bf r}_j- {\bf r}_i) \right \} \nonumber \\ 
= && \sum_i \left \{ w_i
\left (\Delta{\bf F}_j \cdot {\bf v}_j \right ) 
({\bf r}_j - {\bf r}_i) + w_i \sum_k
\left (\Delta{\bf F}_k \cdot {\bf v}_k \right) 
({\bf r}_j - {\bf r}_i) + w_i \sum_k
\left (\Delta{\bf F}_k \cdot {\bf v}_k \right ) ({\bf r}_k - {\bf
r}_j) \right \}
\nonumber \\ 
= && \left ( \Delta{\bf F}_j \cdot {\bf v}_j + \sum_k
\Delta{\bf F}_k \cdot {\bf v}_k \right ) \left ( \sum_i w_i ({\bf r}_j
- {\bf r}_i) \right ) + \sum_k (\Delta{\bf F}_k \cdot {\bf v}_k) ({\bf
r}_k - {\bf r}_j).
\label{total_angular_heat_current}
\end{eqnarray}
That is, we only need to evaluate the $3$-dimensional moment during
the $q_j$ summation, then use (\ref{total_angular_heat_current}) with
the {\em total} angular forces which appear in the final force
increments.

The above process are of course error-prone. One may check one's
analytical derivation and coding by running a MD system {\em without}
PBC, and compare the calculated heat current with the numerical time
differentiation of $\sum_i E_i {\bf r}_i$ using the specific $\{E_i\}$
partition, which of course is very easy to compute. If the two agree
with each other, the derived heat current can then be used for PBC
systems.

\section{Angular Summation}

\begin{eqnarray}
&& \sum_{j<k} g(r_{ij}) g(r_{ik}) \nonumber \\ = && \frac{1}{2} \left
[ \sum_{j,k} g(r_{ij}) g(r_{ik}) - \sum_{j} g(r_{ij})^2 \right ]
\nonumber \\ = && \frac{1}{2} \left [ \left ( \sum_{j} g(r_{ij})
\right )^2 - \sum_{j} g(r_{ij})^2 \right ]
\label{zero_summation}
\end{eqnarray}

\begin{eqnarray}
&& \sum_{j<k} g(r_{ij}) g(r_{ik}) \cos(\theta_{ijk}) \nonumber \\
= && \frac{1}{2} \left [ \sum_{j,k} g(r_{ij}) g(r_{ik})  
\cos(\theta_{ijk}) - \sum_{j}
g(r_{ij})^2 \right ] \nonumber \\
= && \frac{1}{2} \left [ \sum_{j,k} g(r_{ij}) g(r_{ik}) 
\sum_{\alpha=1}^3 \hat{r}_{ij}^\alpha \hat{r}_{ik}^\alpha - \sum_{j}
g(r_{ij})^2 \right ] \nonumber \\
= && \frac{1}{2} \left [ \sum_{\alpha=1}^3 
\left ( \sum_{j} g(r_{ij}) \hat{r}_{ij}^\alpha \right )^2 - 
\sum_{j} g(r_{ij})^2 \right ]
\label{single_summation}
\end{eqnarray}

\begin{eqnarray}
&& \sum_{j<k} g(r_{ij}) g(r_{ik}) \cos^2(\theta_{ijk}) \nonumber \\
= && \frac{1}{2} \left [ \sum_{j,k} g(r_{ij}) g(r_{ik})  
\cos^2(\theta_{ijk}) - \sum_{j}
g(r_{ij})^2 \right ] \nonumber \\
= && \frac{1}{2} \left [ \sum_{j,k} g(r_{ij}) g(r_{ik}) 
\sum_{\alpha=1}^3 \sum_{\beta=1}^3 
\hat{r}_{ij}^\alpha \hat{r}_{ik}^\alpha 
\hat{r}_{ij}^\beta \hat{r}_{ik}^\beta - \sum_{j}
g(r_{ij})^2 \right ] \nonumber \\
= && \frac{1}{2} \left [ \sum_{\alpha=1}^3 \sum_{\beta=1}^3 
\left ( \sum_{j} g(r_{ij}) \hat{r}_{ij}^\alpha \hat{r}_{ij}^\beta 
\right )^2 - \sum_{j} g(r_{ij})^2 \right ] 
\label{double_summation}
\end{eqnarray}

\begin{eqnarray}
&& \sum_{j<k} g(r_{ij}) g(r_{ik}) \cos^3(\theta_{ijk}) \nonumber \\
= && \frac{1}{2} \left [ \sum_{j,k} g(r_{ij}) g(r_{ik})  
\cos^3(\theta_{ijk}) - \sum_{j}
g(r_{ij})^2 \right ] \nonumber \\
= && \frac{1}{2} \left [ \sum_{j,k} g(r_{ij}) g(r_{ik}) 
\sum_{\alpha=1}^3 \sum_{\beta=1}^3 \sum_{\gamma=1}^3 
\hat{r}_{ij}^\alpha \hat{r}_{ik}^\alpha 
\hat{r}_{ij}^\beta \hat{r}_{ik}^\beta 
\hat{r}_{ij}^\gamma \hat{r}_{ik}^\gamma - \sum_{j}
g(r_{ij})^2 \right ] \nonumber \\
= && \frac{1}{2} \left [ \sum_{\alpha=1}^3 \sum_{\beta=1}^3 
\sum_{\gamma=1}^3 
\left ( \sum_{j} g(r_{ij}) \hat{r}_{ij}^\alpha \hat{r}_{ij}^\beta 
\hat{r}_{ij}^\gamma \right )^2 - \sum_{j} g(r_{ij})^2 \right ]
\label{triple_summation}
\end{eqnarray}

Therefore, the Stillinger-Weber potential\cite{SW85}

\begin{eqnarray}
&& \sum_{j<k} g(r_{ij}) g(r_{ik}) \left ( \cos(\theta_{ijk})+\frac{1}{3}
\right )^2 \nonumber \\
= && \sum_{j<k} g(r_{ij}) g(r_{ik}) \left ( \cos^2(\theta_{ijk})+
\frac{2}{3} \cos(\theta_{ijk}) + \frac{1}{9} \right ) \nonumber \\
= && \frac{1}{2} \left ( \;\;
\left [ \sum_{\alpha=1}^3 \sum_{\beta=1}^3 
\left ( 
 \sum_{j} g(r_{ij}) \hat{r}_{ij}^\alpha \hat{r}_{ij}^\beta 
\right )^2 - \sum_{j} g(r_{ij})^2 \right ]
\right . \nonumber \\
+ && \left . \frac{2}{3} \left [ \sum_{\alpha=1}^3 
\left ( \sum_{j} g(r_{ij}) \hat{r}_{ij}^\alpha \right )^2 
- \sum_{j} g(r_{ij})^2  \right ] + \frac{1}{9} \left [ 
\sum_{j} g(r_{ij}) ^2 - 
\sum_{j} g(r_{ij})^2 \right ] \;\; \right ) \nonumber \\
= && \frac{1}{2} \sum_{\alpha=1}^3 \sum_{\beta=1}^3 
\left ( \sum_{j} g(r_{ij}) \hat{r}_{ij}^\alpha \hat{r}_{ij}^\beta 
\right )^2 + \frac{1}{3} \sum_{\alpha=1}^3 
\left ( \sum_{j} g(r_{ij}) \hat{r}_{ij}^\alpha \right )^2 
\nonumber \\ + && \frac{1}{18} \sum_{j} g(r_{ij}) ^2
- \frac{8}{9} \sum_{j} g(r_{ij})^2.
\label{SW_simplified_1}
\end{eqnarray}

One can further simplify (\ref{SW_simplified_1}) by extending
$\hat{r}_{ij}$ to four dimensions. Define vector
\begin{equation}
\tilde{r}_{ij} \;\;\equiv\;\; \left 
(\hat{r}_{ij}, \;\;\frac{1}{\sqrt{3}} \right ), 
\end{equation}
where the fourth component is constant $1/\sqrt{3}$. It can be shown
that (\ref{SW_simplified_1}) is equal to
\begin{eqnarray}
\sum_{j<k} g(r_{ij}) g(r_{ik}) \left ( \cos(\theta_{ijk})+\frac{1}{3}
\right )^2 = \frac{1}{2} \sum_{\alpha=1}^4 \sum_{\beta=1}^4
\left ( \sum_{j} g(r_{ij}) \tilde{r}_{ij}^\alpha \tilde{r}_{ij}^\beta 
\right )^2 - \frac{8}{9} \sum_{j} g(r_{ij})^2.
\label{SW_simplified_2}
\end{eqnarray} 
Both (\ref{SW_simplified_1}) and (\ref{SW_simplified_2}) indicates
that the Stillinger-Weber three-body term is a tensorial extension of
the embedded atom formalism with embedding exponent $2$.

To get the force on atom $j$ due to site energy of atom $i$, let us
define ``charge density matrix''
\begin{equation}
{\bf M}:\;\;\;\; M^{\alpha\beta}_i \;\;=\;\; \sum_{j} g(r_{ij})
\tilde{r}_{ij}^\alpha \tilde{r}_{ij}^\beta, \;\;\;\;
\alpha,\beta=1,2,3,4.
\end{equation}
I then have
\begin{equation}
v_i = \frac{1}{2} {\rm Tr}{\bf M}^2 - \frac{8}{9} \sum_{j}
g(r_{ij})^2.
\end{equation}
The second part is just an ordinary pair potential with force
contributions
\begin{equation}
\Delta F^{\mu}_j = -\Delta F^{\mu}_i = \frac{16}{9}
g(r_{ij}) g^\prime(r_{ij}) \hat{r}_{ij}^\mu.
\end{equation}
The first contribution is more interesting because it must give
bond-bending forces that are similar to vector cross products. If atom
$j$ is displaced by $dx_{j}^{\mu}$, since
\begin{equation}
dr_{ij}^{\mu} = dx_{j}^{\mu} \;\rightarrow\; 
d\tilde{r}_{ij}^{\nu} = d\hat{r}_{ij}^{\nu} = d\left (
\frac{ r_{ij}^{\nu} }{ r_{ij} } \right ) = 
\sum_{\mu=1}^3 \frac{\delta^{\nu\mu} -
\hat{r}_{ij}^{\nu} \hat{r}_{ij}^{\mu}}{ r_{ij} } dx_{j}^{\mu},
\;\;\;\; \nu=1,2,3,
\end{equation}
there is,
\begin{eqnarray}
dM^{\alpha\beta}_i = && \sum_{\mu=1}^3 \left [
g^\prime(r_{ij}) \tilde{r}_{ij}^\alpha \tilde{r}_{ij}^\beta
\hat{r}_{ij}^{\mu}
+ g(r_{ij})  \frac{
\tilde{r}_{ij}^\alpha (\delta^{\beta\mu} -
\hat{r}_{ij}^{\beta} \hat{r}_{ij}^{\mu}) + 
\tilde{r}_{ij}^\beta (\delta^{\alpha\mu} -
\hat{r}_{ij}^{\alpha} \hat{r}_{ij}^{\mu})}{ r_{ij} } \right ]
dx_{j}^{\mu} \nonumber \\
= && \sum_{\mu=1}^3 \left [
\left ( g^\prime(r_{ij}) \tilde{r}_{ij}^\alpha \tilde{r}_{ij}^\beta
- \frac{g(r_{ij})}{r_{ij}} \tilde{r}_{ij}^\alpha \hat{r}_{ij}^{\beta}
- \frac{g(r_{ij})}{r_{ij}} \hat{r}_{ij}^\alpha \tilde{r}_{ij}^{\beta}
\right ) \hat{r}_{ij}^{\mu}
+ g(r_{ij})  \frac{\tilde{r}_{ij}^\alpha \delta^{\beta\mu} + 
\tilde{r}_{ij}^\beta \delta^{\alpha\mu}}{r_{ij}} \right ]
dx_{j}^{\mu}.
\label{SW_force_intermediate}
\end{eqnarray}

The first term of (\ref{SW_force_intermediate}) still gives a pair
potential like contribution, with
\begin{equation}
\Delta F^{\mu}_j = -\Delta F^{\mu}_i = \left ( \sum_{\alpha=1}^4
\sum_{\beta=1}^4 M^{\alpha\beta}_i \left ( 
\frac{g(r_{ij})}{r_{ij}} \tilde{r}_{ij}^\alpha \hat{r}_{ij}^{\beta}
\frac{g(r_{ij})}{r_{ij}} \hat{r}_{ij}^\alpha \tilde{r}_{ij}^{\beta}
- g^\prime(r_{ij}) \tilde{r}_{ij}^\alpha \tilde{r}_{ij}^\beta 
\right ) \right ) 
\hat{r}_{ij}^{\mu},
\end{equation}
whereas the second term gives
\begin{equation}
\Delta F^{\mu}_j = -\Delta F^{\mu}_i = - 2 \sum_{\alpha=1}^4
M^{\alpha\mu}_i \frac{g(r_{ij})}{r_{ij}} \tilde{r}_{ij}^\alpha.
\end{equation}


\begin{thebibliography}{00}
\begin{singlespace}

\bibitem{Willaime89} F. Willaime, C. Massobrio, Phys. Rev. Lett. {\bf
63}, 2244 (1991).

\bibitem{Willaime91} F. Willaime, C. Massobrio, Phys. Rev. B {\bf 43},
11653 (1991).

\bibitem{SW85} F. H. Stillinger, T. A. Weber, Phys. Rev. B {\bf 31},
5262 (1985).

\end{singlespace}
\end{thebibliography}

\end{document}
