Periodic molecular dynamics with many-body potential under external stress
Gang Liu
The dynamical equations for periodic systems under constant external stress (E-print, arXiv:cond-mat/0209372 ) have been extended to many-body potential case. Rather than introducing Lagrangian dynamics, Newtonian dynamics was employed directly. The kinetic energy term of stress and the imagined force only due to particles' passing through geometric planes or boundaries were also discussed.
cond-mat0505251.tex
%\documentclass[aps]{revtex4}
%\documentclass[aps,preprint]{revtex4}
%\usepackage{epsfig}
%\documentstyle[aps,preprint]{revtex}


\documentstyle[aps]{revtex}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\addtolength{\topmargin}{0.0in}

\begin{document}
\title{Periodic molecular dynamics with many-body potential under external stress}
\author{Gang Liu}
\address{{\small (liugang@post.queensu.ca)}\\
{\small High Performance Computing Virtual Laboratory}\\
{\small Queen's University, Kingston, ON K7L 3N6, Canada}}
\date{May 10, 2005}
\maketitle

\begin{abstract}
The dynamical equations for periodic systems under constant external stress
(E-print, arXiv:cond-mat/0209372 ) have been extended to many-body potential
case. Rather than introducing Lagrangian dynamics, Newtonian dynamics was
employed directly. The kinetic energy term of stress and the imagined force
only due to particles' passing through geometric planes or boundaries were
also discussed.
\end{abstract}

\pacs{61.50.Ah, 62.20.-x, 64.10.+h.}

The molecular dynamics (MD) simulation\cite{al,fre,haile} is an important
tool and widely used in research in Condensed Matter Physics, Biophysics,
Chemical Physics, Space Physics, and even Astrophysics. The primary idea is
to simply apply Newton's Second Law onto every particle. For simplicity or
feasibility, the periodic boundary condition is often imposed, which results
in that the system is modelized as a crystal in each simulation time step.
Then only the MD particles (in the MD cell) and the periods are the degrees
of freedom. The principle physics problem becomes dynamical equations for
the degrees of freedom with the possibility of external stress/forces acting
on surfaces of the system. Obviously Newton's Second Law is absolutely
applicable directly onto the MD particles. What about the periods?

By extending Andersen's idea\cite{a} and introducing an unproved Lagrangian,
Parrinello and Rahman proposed their theory regarding the dynamical
equations for the periods in 1980\cite{pr1,pr2}. Although their theory has
been extensively employed in simulations, we recently found that the MD\
particles violate Newton's Second Law in their theory\cite{lg.eprint2}. The
reason is that the kinetic energy for the MD particles in their Lagrangian
is incomplete. Why such a Lagrangian is adopted?

The whole story is closely associated with another key problem: whether the
internal stress contains the so called kinetic energy term. As essentially
stress is the net force per unit area, it definitely contains the regular
force term. The truth is that the kinetic energy term is widely believed and
calculated as the other part of stress. Irving and Kirkwood proved the
kinetic energy term for stress with pure statistics in 1950\cite{ik}. Haile
proved it by introducing an imaged force appearing only if a particle
passing through a geometric plane in 1992\cite{haile}. Zhou recently proved
that the kinetic energy term was not physically grounded and should be
deleted from the stress expression\cite{zhou}. We made some discussion on
physics of the imaged force in the TRANSPORT OF MOMENTUM section in our
previous work\cite{lg.eprint1} and found that for a system defined with a
fixed space there should be a momentum flux term to cancel the unneeded
momentum change per unit time only due to particles' passing through the
system boundary in order to have Newton's Second Law applied correctly\cite
{lg.eprint1,lg.eprint3}. If we put this momentum flux term onto the momentum
change side of the dynamical equation, nothing new; if we move it onto the
force side, it can be regarded as an effective force and further results in
the kinetic energy term of stress in principle. When we regard fluids as
composed of independent cells (no periodicity) and separate their whole
motion into macroscopic and internal ones, the kinetic energy term of stress
can be explained as momentum flux due to particles' internal motion across
boundaries or collisions of internal motion between particles from a
statistical point of view\cite{lg.eprint3}, resulting the stress similar to
that proposed by Irving and Kirkwood\cite{ik}. The periodic system simulated
with molecular dynamics as a crystal can be regarded as a special fluid, in
which the periods represent the macroscopic motion and the motion of the MD
particles is the internal motion. Then the kinetic energy term of stress
from the MD\ particles can be suggested. Parrinello and Rahman obtained a
stress with a similar kinetic energy term by choosing the above incomplete
Lagrangian\cite{pr1,pr2}, but as a result, the MD\ particles violate
Newton's Second Law\cite{lg.eprint2}. If the complete Lagrangian were used
there, the similar kinetic energy term of the stress would have not been
generated in physics, as in our previous works\cite{lww,lwssc}.

Instead of introducing a Lagrangian, we applied Newton's Second Law onto the
MD particles directly avoiding the above problem and applied Newton's Second
Law onto halves of the system for dynamical equations of the periods (${\bf a%
}$, ${\bf b}$, and ${\bf c}$), and got the kinetic energy term of the
internal stress by considering the collisions between particles and the
boundary recently\cite{lg.eprint1}. However this work was done for
pair-potential only. Now let us extend it to many-body potential cases. The
many-body force and potential were already discussed previously\cite
{lwssc,lchuzhan} with some misprint in Reference 13, but starting from a
Lagrangian for constant-pressure molecular dynamics. Here we will start from
the forces for the constant-stress molecular dynamics. For simplicity, we
use Fig. 1 in version 6 of Reference 10 as Fig. 1 here.

For up to $M$-body interactions considered, the kernel task for this
extension is to find the net force acting on the half-line-cell bar $B_{{\bf %
h}}$ by all particles in the $L$ part (left half) of the crystal in Fig. 1 
\begin{equation}
{\bf F}_{{\bf h}}=\text{ }\sum_{m=2}^M{\bf F}_{{\bf h}}^{(m)},  \label{et01}
\end{equation}
where ${\bf F}_{{\bf h}}^{(m)}$ is the contribution of $m$-body interaction,
and any period ${\bf h=a}$, ${\bf b}$, or ${\bf c}$. For any particle $I_k$
with $k=1,2,\cdots ,$ or $M$ in the whole crystal, its position vector is $%
{\bf r}_{I_k}={\bf T}_{I_k}+{\bf r}_k$, where ${\bf r}_k$ is the position
vector for its image particle in the MD\ cell, and 
\begin{equation}
{\bf T}_{I_k}=I_{k,{\bf a}}{\bf a}+I_{k,{\bf b}}{\bf b}+I_{k,{\bf c}}{\bf c}
\label{et02}
\end{equation}
is the lattice translation vector of the cell in which it resides. For $m$%
-body forces, the interaction between the right and left half of the crystal
means all participating particles must be distributed in both halves, or all
participating particles must not be in any one of the halves. Then ${\bf F}_{%
{\bf h}}^{(m)}$ is the summation of all these possibilities devided by the
number of equivalent bars to $B_{{\bf h}}$ in the right half crystal $N_{%
{\bf h}}=\left| {\bf S}\right| /\left| {\bf \sigma }_{{\bf h}}\right| $,
where ${\bf S}$ is the cross section area vector in Fig. 1, and ${\bf \sigma 
}_{{\bf h}}=\partial \Omega /\partial {\bf h}$ with cell volume $\Omega $.
Now let us write out the force for all the possibilities with total $t$
particles ($m>t\geq 1$ ) distributed in the right half of the crystal and
the rest particles in the left half 
\begin{equation}
{\bf F}_{t,{\bf h}}^{(m)}=\text{ }\frac 1{N_{{\bf h}}}\frac 1{t!\left(
m-t\right) !}\sum_{\left\{ I_1,I_2,\cdots ,I_t\right\} }^{\left( I_{1,{\bf h}%
},I_{2,{\bf h}},\cdots ,I_{t,{\bf h}}\geq 0\right) }\sum_{\left\{
I_{t+1},I_{t+2},\cdots ,I_m\right\} }^{\left( I_{t+1,{\bf h}},I_{t+2,{\bf h}%
},\cdots ,I_{m,{\bf h}}<0\right) }\sum_{\mu =1}^t{\bf f}_{I_\mu }^{(m)}({\bf %
r}_{I_1},{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}_{I_m}),  \label{et03}
\end{equation}
where we use $\sum_{\{\cdots \}}^{(\cdots )}$ to denote the summation over
indexes in $\{\cdots \}$ with conditions in $(\cdots )$, and ${\bf f}_{I_\mu
}^{(m)}({\bf r}_{I_1},{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}_{I_m})$ is
the $m$-body force acting on particle $I_\mu $. We also require that any two
of the indexes in $\{\cdots \}$ must not refer to the same one physical
particle, while any particle and its image particle in any other cell are
different physical particles. Remembering 
\begin{equation}
\sum_{I_1}^{I_{1,{\bf h}}\geq 0}\cdots =\sum_{I_{1,{\bf h}}}^{I_{1,{\bf h}%
}\geq 0}\sum_{I_{1,{\bf h}^{\prime }}}\sum_{I_{1,{\bf h}^{\prime \prime
}}}\sum_{i_1\in Cell_{MD}}\cdots  \label{et04}
\end{equation}
with ${\bf h}^{\prime }$, ${\bf h}^{\prime \prime }{\bf =a}$, ${\bf b}$, or $%
{\bf c}$, but ${\bf h}^{\prime }\neq {\bf h}$, ${\bf h}^{\prime \prime }\neq 
{\bf h}$, and ${\bf h}^{\prime \prime }\neq {\bf h}^{\prime }$, $N_{{\bf h}%
}=\sum_{I_{1,{\bf h}^{\prime }}}\sum_{I_{1,{\bf h}^{\prime \prime }}}1$, and
the crystal translatability, we have 
\begin{eqnarray}
{\bf F}_{t,{\bf h}}^{(m)} &=&\frac 1{t!\left( m-t\right) !}\sum_{I_{1,{\bf h}%
}=0}^{+\infty }\sum_{\left\{ i_1,I_2,I_3,\cdots ,I_t\right\} }^{\left( I_{2,%
{\bf h}},I_{3,{\bf h}},\cdots ,I_{t,{\bf h}}\geq 0\right) }\sum_{\left\{
I_{t+1},I_{t+2},\cdots ,I_m\right\} }^{\left( I_{t+1,{\bf h}},I_{t+2,{\bf h}%
},\cdots ,I_{m,{\bf h}}<0\right) }\sum_{\mu =1}^t{\bf f}_{I_\mu }^{(m)}({\bf %
r}_{I_1},{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}_{I_m})  \label{et05} \\
\ &=&\frac 1{t!\left( m-t\right) !}\sum_{I_{1,{\bf h}}=0}^{+\infty
}\sum_{\left\{ i_1,I_2,I_3,\cdots ,I_t\right\} }^{\left( I_{2,{\bf h}},I_{3,%
{\bf h}},\cdots ,I_{t,{\bf h}}\geq 0\right) }\sum_{\left\{
I_{t+1},I_{t+2},\cdots ,I_m\right\} }^{\left( I_{t+1,{\bf h}},I_{t+2,{\bf h}%
},\cdots ,I_{m,{\bf h}}<0\right) }t{\bf f}_{I_1}^{(m)}({\bf r}_{I_1},{\bf r}%
_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}_{I_m}).  \label{et06}
\end{eqnarray}
With the crystal translatability and moving the system so that the cell
containing particle $I_1$ is translated to the MD cell, it becomes 
\begin{equation}
{\bf F}_{t,{\bf h}}^{(m)}=\frac t{t!\left( m-t\right) !}\sum_{l=0}^{-\infty
}\sum_{\left\{ i,I_2,I_3,\cdots ,I_t\right\} }^{\left( I_{2,{\bf h}},I_{3,%
{\bf h}},\cdots ,I_{t,{\bf h}}\geq l\right) }\sum_{\left\{
I_{t+1},I_{t+2},\cdots ,I_m\right\} }^{\left( I_{t+1,{\bf h}},I_{t+2,{\bf h}%
},\cdots ,I_{m,{\bf h}}<l\right) }{\bf f}_i^{(m)}({\bf r}_i,{\bf r}_{I_2},%
{\bf r}_{I_3},\cdots ,{\bf r}_{I_m}),  \label{et08}
\end{equation}
and Eq. (\ref{et01}) becomes 
\begin{eqnarray}
{\bf F}_{{\bf h}} &=&\text{ }\sum_{m=2}^M\sum_{t=1}^{m-1}{\bf F}_{t,{\bf h}%
}^{(m)}  \label{et09} \\
&=&\sum_{m=2}^M\sum_{t=1}^{m-1}\frac 1{\left( t-1\right) !\left( m-t\right) !%
}\sum_{l=0}^{-\infty }\sum_{\left\{ i,I_2,I_3,\cdots ,I_t\right\} }^{\left(
I_{2,{\bf h}},I_{3,{\bf h}},\cdots ,I_{t,{\bf h}}\geq l\right)
}\sum_{\left\{ I_{t+1},I_{t+2},\cdots ,I_m\right\} }^{\left( I_{t+1,{\bf h}%
},I_{t+2,{\bf h}},\cdots ,I_{m,{\bf h}}<l\right) }{\bf f}_i^{(m)}({\bf r}_i,%
{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}_{I_m}).  \label{et10}
\end{eqnarray}

Meanwhile the total up to $M$-body cell potential energy is 
\begin{equation}
E_{cell}=\sum_{m=2}^ME_{cell}^{(m)},  \label{et21}
\end{equation}
where the $m$-body cell potential energy 
\begin{equation}
E_{cell}^{(m)}=\frac 1{m!}\sum_{\left\{ i,I_2,I_3,\cdots ,I_m\right\}
}\varphi ^{(m)}({\bf r}_i,{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}_{I_m}),
\label{et22}
\end{equation}
with $m$-body potential $\varphi ^{(m)}({\bf r}_{I_1},{\bf r}_{I_2},\cdots ,%
{\bf r}_{I_m})$. This form of $E_{cell}^{(m)}$ ensures that only a fraction $%
\frac km$ of $\varphi ^{(m)}$ belong to the cell when $k$ $\left( m\geq
k\geq 0\right) $ particles are in the cell. Now let us take a derivative 
\begin{equation}
-\frac \partial {\partial {\bf h}}E_{cell}^{(m)}=\text{ }\frac 1{m\left(
m-2\right) !}\sum_{\left\{ i,I_2,I_3,\cdots ,I_m\right\} }I_{m,{\bf h}}{\bf f%
}_{I_m}^{(m)}({\bf r}_i,{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}_{I_m}),
\label{et25}
\end{equation}
with 
\begin{equation}
{\bf f}_{I_k}^{(m)}({\bf r}_{I_1},{\bf r}_{I_2},\cdots ,{\bf r}%
_{I_m})=-\nabla _{{\bf r}_{I_k}}\varphi ^{(m)}({\bf r}_{I_1},{\bf r}%
_{I_2},\cdots ,{\bf r}_{I_m}).  \label{et26}
\end{equation}
The right side of Eq. (\ref{et25}) can be split into two terms based on the
sign of $I_{m,{\bf h}}$ 
\begin{equation}
-\frac \partial {\partial {\bf h}}E_{cell}^{(m)}=\text{ }{\bf F}_{{\bf h,}%
+}^{(m)}+{\bf F}_{{\bf h,}-}^{(m)},  \label{et30}
\end{equation}
where 
\begin{eqnarray}
{\bf F}_{{\bf h,}+}^{(m)} &=&\text{ }\frac 1{m\left( m-2\right) !}%
\sum_{\left\{ i,I_2,I_3,\cdots ,I_m\right\} }^{I_{m,{\bf h}}>0}I_{m,{\bf h}}%
{\bf f}_{I_m}^{(m)}({\bf r}_i,{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}%
_{I_m})  \label{et41} \\
\ &=&\text{ }\frac 1{m\left( m-2\right) !}\sum_{l=0}^{+\infty }\sum_{\left\{
i,I_2,I_3,\cdots ,I_m\right\} }^{I_{m,{\bf h}}>l}{\bf f}_{I_m}^{(m)}({\bf r}%
_i,{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}_{I_m}),  \label{et42}
\end{eqnarray}
\begin{eqnarray}
{\bf F}_{{\bf h,}-}^{(m)} &=&\text{ }\frac 1{m\left( m-2\right) !}%
\sum_{\left\{ i,I_2,I_3,\cdots ,I_m\right\} }^{I_{m,{\bf h}}<0}I_{m,{\bf h}}%
{\bf f}_{I_m}^{(m)}({\bf r}_i,{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}%
_{I_m})  \label{et45} \\
\ &=&\text{ }\frac{-1}{m\left( m-2\right) !}\sum_{l=0}^{-\infty
}\sum_{\left\{ i,I_2,I_3,\cdots ,I_m\right\} }^{I_{m,{\bf h}}<l}{\bf f}%
_{I_m}^{(m)}({\bf r}_i,{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}_{I_m}).
\label{et46}
\end{eqnarray}
Eq. (\ref{et42}) can be expanded regarding $t^{\prime }$ , the number of
particles distributed in the right half of the crystal, of particles indexed
from $I_2$ to $I_{m-1}$ 
\begin{equation}
{\bf F}_{{\bf h,}+}^{(m)}=\text{ }\frac 1m\sum_{t^{\prime
}=0}^{m-2}\sum_{l=0}^{+\infty }\sum_{I_m}^{I_{m,{\bf h}}>l}\sum_{\left\{
i,I_2,I_3,\cdots ,I_{t^{\prime }+1}\neq I_m\right\} }^{\left( I_{2,{\bf h}%
},I_{3,{\bf h}},\cdots ,I_{t^{\prime }+1,{\bf h}}\geq 0\right)
}\sum_{\left\{ I_{t^{\prime }+2},I_{t^{\prime }+3},\cdots ,I_{m-1}\right\}
}^{\left( I_{t^{\prime }+2,{\bf h}},I_{t^{\prime }+3,{\bf h}},\cdots ,I_{m-1,%
{\bf h}}<0\right) }\frac{{\bf f}_{I_m}^{(m)}({\bf r}_i,{\bf r}_{I_2},{\bf r}%
_{I_3},\cdots ,{\bf r}_{I_m})}{t^{\prime }!\left( m-t^{\prime }-2\right) !}.
\label{et52}
\end{equation}
Then let us use the translatability, move the system so that the cell
containing particle $I_m$ is translated to the MD cell, and rename particle $%
I_m$ as particle $i$ and particle $i$ as particle $I_m$%
\begin{eqnarray}
{\bf F}_{{\bf h,}+}^{(m)} &=&\text{ }\frac 1m\sum_{t^{\prime
}=0}^{m-2}\sum_{l=0}^{-\infty }\sum_{\left\{ i,I_2,I_3,\cdots ,I_{t^{\prime
}+1}\neq I_m\right\} }^{\left( I_{2,{\bf h}},I_{3,{\bf h}},\cdots
,I_{t^{\prime }+1,{\bf h}}\geq l\right) }\sum_{\left\{ I_{t^{\prime
}+2},I_{t^{\prime }+3},\cdots ,I_{m-1},I_m\right\} }^{\left( I_{t^{\prime
}+2,{\bf h}},I_{t^{\prime }+3,{\bf h}},\cdots ,I_{m-1,{\bf h}},I_{m,{\bf h}%
}<l\right) }\frac{{\bf f}_i^{(m)}({\bf r}_i,{\bf r}_{I_2},{\bf r}%
_{I_3},\cdots ,{\bf r}_{I_m})}{t^{\prime }!\left( m-t^{\prime }-2\right) !}
\label{et55} \\
\ &=&\sum_{t=1}^{m-1}\frac{m-t}m{\bf F}_{t,{\bf h}}^{(m)}.  \label{et56}
\end{eqnarray}
Please note the $m$-body force and potential are independent on the order of
the $m$ particles' position variables in the expression of ${\bf f}%
_{I_k}^{(m)}({\bf r}_{I_1},{\bf r}_{I_2},\cdots ,{\bf r}_{I_m})$ and $%
\varphi ^{(m)}({\bf r}_{I_1},{\bf r}_{I_2},\cdots ,{\bf r}_{I_m})$
respectively.

For each given $l$ (negative) value, Eq. (\ref{et46}) can be expanded with $%
t^{\prime }$ particles of those indexed from $I_2$ to $I_{m-1}$ distributed
in the right half of the crystal defined by $I_{k,{\bf h}}\geq l$ 
\begin{equation}
{\bf F}_{{\bf h,}-}^{(m)}=\text{ }\frac{-1}m\sum_{t^{\prime
}=0}^{m-2}\sum_{l=0}^{-\infty }\sum_{\left\{ i,I_2,I_3,\cdots ,I_{t^{\prime
}+1}\right\} }^{\left( I_{2,{\bf h}},I_{3,{\bf h}},\cdots ,I_{t^{\prime }+1,%
{\bf h}}\geq l\right) }\sum_{\left\{ I_{t^{\prime }+2},I_{t^{\prime
}+3},\cdots ,I_{m-1},I_m\right\} }^{\left( I_{t^{\prime }+2,{\bf h}%
},I_{t^{\prime }+3,{\bf h}},\cdots ,I_{m-1,{\bf h}},I_{m,{\bf h}}<l\right) }%
\frac{{\bf f}_{I_m}^{(m)}({\bf r}_i,{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf %
r}_{I_m})}{t^{\prime }!\left( m-t^{\prime }-2\right) !}.  \label{et58}
\end{equation}
Making use of crystal translatability again, and 
\begin{equation}
\sum_{\mu =1}^m{\bf f}_{I_\mu }^{(m)}({\bf r}_{I_1},{\bf r}_{I_2},{\bf r}%
_{I_3},\cdots ,{\bf r}_{I_m})=0  \label{et58.2}
\end{equation}
for conservation laws of energy and momentum, we have $(l\leq 0)$ 
\begin{eqnarray}
&&\sum_{\left\{ i,I_2,I_3,\cdots ,I_{t^{\prime }+1}\right\} }^{\left( I_{2,%
{\bf h}},I_{3,{\bf h}},\cdots ,I_{t^{\prime }+1,{\bf h}}\geq l\right)
}\sum_{\left\{ I_{t^{\prime }+2},I_{t^{\prime }+3},\cdots ,I_m\right\}
}^{\left( I_{t^{\prime }+2,{\bf h}},I_{t^{\prime }+3,{\bf h}},\cdots ,I_{m,%
{\bf h}}<l\right) }-\left( m-t^{\prime }-1\right) {\bf f}_{I_m}^{(m)}({\bf r}%
_i,{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}_{I_m})  \nonumber
\label{et58} \\
&=&\sum_{\left\{ i,I_2,I_3,\cdots ,I_{t^{\prime }+1}\right\} }^{\left( I_{2,%
{\bf h}},I_{3,{\bf h}},\cdots ,I_{t^{\prime }+1,{\bf h}}\geq l\right)
}\sum_{\left\{ I_{t^{\prime }+2},I_{t^{\prime }+3},\cdots ,I_m\right\}
}^{\left( I_{t^{\prime }+2,{\bf h}},I_{t^{\prime }+3,{\bf h}},\cdots ,I_{m,%
{\bf h}}<l\right) }\left( t^{\prime }+1\right) {\bf f}_i^{(m)}({\bf r}_i,%
{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}_{I_m}),  \label{et59}
\end{eqnarray}
then 
\begin{equation}
{\bf F}_{{\bf h,}-}^{(m)}=\text{ }\sum_{t=1}^{m-1}\frac tm{\bf F}_{t,{\bf h}%
}^{(m)}.  \label{et60}
\end{equation}
As a result, we have 
\begin{equation}
{\bf F}_{{\bf h}}=\text{ }\sum_{m=2}^M\sum_{t=1}^{m-1}{\bf F}_{t,{\bf h}%
}^{(m)}=\sum_{m=2}^M\left( {\bf F}_{{\bf h,}+}^{(m)}+{\bf F}_{{\bf h,}%
-}^{(m)}\right) =-\frac \partial {\partial {\bf h}}E_{cell}.  \label{et61}
\end{equation}

With the main interaction tensor defined for up to $M$-body cases 
\begin{equation}
\stackrel{\rightharpoonup \rightharpoonup }{{\cal \varepsilon }_m}{\bf =}%
\frac{-1}\Omega \left[ \left( \frac{\partial E_{cell}}{\partial {\bf a}}%
\right) {\bf a+}\left( \frac{\partial E_{cell}}{\partial {\bf b}}\right) 
{\bf b+}\left( \frac{\partial E_{cell}}{\partial {\bf c}}\right) {\bf c}%
\right] =\sum_{m=2}^M\stackrel{\rightharpoonup \rightharpoonup }{\cal %
\varepsilon }_m^{(m)}  \label{et62}
\end{equation}
with 
\begin{equation}
\stackrel{\rightharpoonup \rightharpoonup }{\cal \varepsilon }_m^{(m)}=\frac{%
-1}\Omega \left[ \left( \frac{\partial E_{cell}^{(m)}}{\partial {\bf a}}%
\right) {\bf a+}\left( \frac{\partial E_{cell}^{(m)}}{\partial {\bf b}}%
\right) {\bf b+}\left( \frac{\partial E_{cell}^{(m)}}{\partial {\bf c}}%
\right) {\bf c}\right] ,  \label{et62.0}
\end{equation}
we get 
\begin{equation}
{\bf F}_{{\bf h}}{\bf =}\stackrel{\rightharpoonup \rightharpoonup }{{\cal %
\varepsilon }_m}\cdot {\bf \sigma }_{{\bf h}}.  \label{et62.1}
\end{equation}
The $m$-body contribution can also be written as 
\begin{eqnarray}
\stackrel{\rightharpoonup \rightharpoonup }{\cal \varepsilon }_m^{(m)} &&%
{\bf =}\frac 1{\Omega m!}\sum_{\left\{ i,I_2,I_3,\cdots ,I_m\right\}
}\sum_{\mu =2}^m{\bf f}_{I_\mu }^{(m)}({\bf r}_i,{\bf r}_{I_2},{\bf r}%
_{I_3},\cdots ,{\bf r}_{I_m}){\bf T}_{I_\mu }  \label{et62.2} \\
\ &&{\bf =}\frac{m-1}{\Omega m!}\sum_{\left\{ i,I_2,I_3,\cdots ,I_m\right\} }%
{\bf f}_{I_m}^{(m)}({\bf r}_i,{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf r}%
_{I_m}){\bf T}_{I_m}.  \label{et62.4}
\end{eqnarray}
The general interaction tensor should be 
\begin{equation}
\stackrel{\rightharpoonup \rightharpoonup }{\cal \varepsilon }=\sum_{m=2}^M%
\stackrel{\rightharpoonup \rightharpoonup }{\cal \varepsilon }^{(m)}\text{,}
\label{i01}
\end{equation}
where the $m$-body contribution can be extended as 
\begin{eqnarray}
\stackrel{\rightharpoonup \rightharpoonup }{\cal \varepsilon }^{(m)} &=&%
\frac 1{Vm!}\sum_{\left\{ I_1,I_2,I_3,\cdots ,I_m\right\} }\sum_{\mu =1}^m%
{\bf f}_{I_\mu }^{(m)}({\bf r}_{I_1},{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,%
{\bf r}_{I_m})\left( {\bf r}_{I_\mu }-{\bf r}_0\right)  \label{i01} \\
\ &=&\frac 1{Vm!}\sum_{\left\{ I_1,I_2,I_3,\cdots ,I_m\right\} }\sum_{\mu
=1}^m{\bf f}_{I_\mu }^{(m)}({\bf r}_{I_1},{\bf r}_{I_2},{\bf r}_{I_3},\cdots
,{\bf r}_{I_m}){\bf r}_{I_\mu }
\end{eqnarray}
for an arbitrary given position vector ${\bf r}_0$, with volume $V$ of the
whole system and Eq. (\ref{et58.2}). With translatability, it can then be
written as 
\begin{eqnarray}
\stackrel{\rightharpoonup \rightharpoonup }{\cal \varepsilon }^{(m)} &=&%
\frac 1{\Omega m!}\sum_{\left\{ i,I_2,\cdots ,I_m\right\} }\sum_{\mu =2}^m%
{\bf f}_{I_\mu }^{(m)}({\bf r}_i,{\bf r}_{I_2},\cdots ,{\bf r}_{I_m})\left( 
{\bf T}_{I_\mu }+{\bf r}_\mu \right) +  \nonumber \\
&&\ \ \ \ \ +\frac 1{\Omega m!}\sum_{\left\{ i,I_2,\cdots ,I_m\right\} }{\bf %
f}_i^{(m)}({\bf r}_i,{\bf r}_{I_2},\cdots ,{\bf r}_{I_m}){\bf r}_i  \nonumber
\\
\ &=&\frac 1{\Omega m!}\sum_{\left\{ i,I_2,\cdots ,I_m\right\} }\sum_{\mu
=2}^m{\bf f}_{I_\mu }^{(m)}({\bf r}_i,{\bf r}_{I_2},\cdots ,{\bf r}_{I_m})%
{\bf T}_{I_\mu }+  \nonumber \\
&&\ \ \ \ \ +\frac 1{\Omega m!}\sum_{\left\{ i,I_2,\cdots ,I_m\right\} }{\bf %
f}_i^{(m)}({\bf r}_i,{\bf r}_{I_2},\cdots ,{\bf r}_{I_m}){\bf r}_i+ 
\nonumber \\
&&\ \ \ \ \ +\frac 1{\Omega m!}\sum_{\left\{ i,I_2,\cdots ,I_m\right\}
}\sum_{\mu =2}^m{\bf f}_{I_\mu }^{(m)}({\bf r}_i,{\bf r}_{I_2},\cdots ,{\bf r%
}_{I_m}){\bf r}_\mu  \nonumber \\
\ &=&\stackrel{\rightharpoonup \rightharpoonup }{\cal \varepsilon }_m^{(m)}+%
\frac 1{\Omega m!}\sum_{\left\{ i,I_2,\cdots ,I_m\right\} }{\bf f}_i^{(m)}(%
{\bf r}_i,{\bf r}_{I_2},\cdots ,{\bf r}_{I_m}){\bf r}_i+  \nonumber \\
&&\ \ \ \ \ +\frac 1{\Omega m!}\sum_{\left\{ i,I_2,\cdots ,I_m\right\}
}\sum_{\mu =2}^m{\bf f}_\mu ^{(m)}({\bf r}_i-{\bf T}_{I_\mu },{\bf r}_{I_2}-%
{\bf T}_{I_\mu },\cdots ,{\bf r}_\mu \cdots ,{\bf r}_{I_m}-{\bf T}_{I_\mu })%
{\bf r}_\mu  \nonumber \\
\ &=&\stackrel{\rightharpoonup \rightharpoonup }{\cal \varepsilon }_m^{(m)}+%
\stackrel{\rightharpoonup \rightharpoonup }{\cal \varepsilon }_p^{(m)}\text{,%
}  \label{i04}
\end{eqnarray}
where with $n,$ the total number of particles in the MD\ cell, 
\begin{eqnarray}
\stackrel{\rightharpoonup \rightharpoonup }{\cal \varepsilon }_p^{(m)}\ &=&%
\frac 1\Omega \sum_{i=1}^n{\bf F}_i^{\left( m\right) }{\bf r}_i,
\label{et70} \\
{\bf F}_i^{\left( m\right) } &=&\text{ }\frac 1{\left( m-1\right) !}%
\sum_{\left\{ I_2,I_3,\cdots ,I_m\right\} }^{\left( I_2,I_3,\cdots ,I_m\neq
i\right) }{\bf f}_i^{(m)}({\bf r}_i,{\bf r}_{I_2},{\bf r}_{I_3},\cdots ,{\bf %
r}_{I_m})=-\nabla _{{\bf r}_i}E_{cell}^{\left( m\right) }.  \label{et70.2}
\end{eqnarray}

By applying Newton's Second Law onto the right half of the crystal and
applying statistics over all possible collisions between particles and
boundaries and over all possible cut bulks of complete cells from the same
crystal shown in Fig. 2 of version 6 of Reference 10, the same form of
dynamical equations are generated for many-body interactions as in Section
V. of version 6 of Reference 10. For an example, Eq. (56) of version 6 of
Reference 10 is still 
\begin{equation}
W\ddot {{\bf H}}=\left( \stackrel{\rightharpoonup \rightharpoonup }{\pi }+%
\stackrel{\rightharpoonup \rightharpoonup }{\Upsilon }\right) \cdot 
\stackrel{\rightharpoonup \rightharpoonup }{\bf \sigma },  \label{et63}
\end{equation}
where $W$ is the equivalent mass for periods, $\stackrel{\rightharpoonup
\rightharpoonup }{\Upsilon }$ is the external stress, ${\bf H}$ $=\left\{ 
{\bf a,b,c}\right\} $, $\stackrel{\rightharpoonup \rightharpoonup }{\bf %
\sigma }=\left\{ {\bf \sigma }_{{\bf a}},{\bf \sigma }_{{\bf b}},{\bf \sigma 
}_{{\bf c}}\right\} $, and 
\begin{equation}
\stackrel{\rightharpoonup \rightharpoonup }{\pi }=\stackrel{\rightharpoonup
\rightharpoonup }{{\cal \varepsilon }_m}+\stackrel{\rightharpoonup
\rightharpoonup }{\cal \varepsilon }_p+\frac 1{3\Omega }\sum_im_i\left| 
\dot {{\bf r}}_i\right| ^2,  \label{et63.2}
\end{equation}
with particle mass $m_i$ and $\stackrel{\rightharpoonup \rightharpoonup }%
{\cal \varepsilon }_p=\sum_{m=2}^M\stackrel{\rightharpoonup \rightharpoonup }%
{\cal \varepsilon }_p^{(m)}$. For free systems, as $\stackrel{%
\rightharpoonup \rightharpoonup }{\Upsilon }=0$ and no collision, Eq. (\ref
{et63}) becomes 
\begin{equation}
W\ddot {{\bf H}}=\left( \stackrel{\rightharpoonup \rightharpoonup }{{\cal %
\varepsilon }_m}+\frac 1\Omega \sum_{i=1}^n{\bf F}_i{\bf r}_i\right) \cdot 
\stackrel{\rightharpoonup \rightharpoonup }{\bf \sigma },  \label{et64}
\end{equation}
where 
\begin{equation}
{\bf F}_i=\text{ }\sum_{m=2}^M{\bf F}_i^{\left( m\right) }=-\nabla _{{\bf r}%
_i}E_{cell}.  \label{et70}
\end{equation}

As stated above, the MD particles obey Newton's Second Law 
\begin{equation}
m_i{\bf \ddot r}_i={\bf F}_i\text{ }(i=1,2,\cdots ,\text{ }n)\text{ .}
\label{b03}
\end{equation}

\acknowledgments{\ The author wishes to thank Prof. Ding-Sheng Wang ,
Institute of Physics, Prof. Si-Yuan Zhang, Changchun Institute of Applied
Chemistry, Prof. S. S. Wu, Jilin University, P.R. China, and Dr. Kenneth
Edgecombe and Dr. Hartmut Schmider, Queen's University, Canada, for their
helpful discussions and earnest encouragements. }

\begin{references}
\bibitem{al}  M.P. Allen,~and D.J. Tildesley, {\it Computer Simulation of
Liquids (}Oxford University, Oxford, 1987).

\bibitem{fre}  D. Frenkel,~and B. Smit,~{\it Understanding Molecular
Simulation (}Academic, New York, 1996).

\bibitem{haile}  J. M. Haile,~{\it Molecular Dynamics Simulation Elementary
Methods (}John Wiley \& Sons, New York, 1992).

\bibitem{a}  H.C. Andersen,~J. Chem. Phys. {\bf 72,} 2384 (1980).

\bibitem{pr1}  M. Parrinello,~and A. Rahman, Phys. Rev. Lett. {\bf 45,} 1196
(1980).

\bibitem{pr2}  M. Parrinello,~and A. Rahman, J. Appl. Phys. {\bf 52,} 7182
(1981).

\bibitem{lg.eprint2}  E-print, arXiv:cond-mat/0306040
(http://xxx.lanl.gov/abs/cond-mat/0306040).

\bibitem{ik}  J.H. Irving and John G. Kirkwood, J. Chem. Phys. {\bf 18}, 817
(1950).

\bibitem{zhou}  M. Zhou, Proc. R. Soc. Lond. A {\bf 459}, 2347 (2003).

\bibitem{lg.eprint1}  E-print, arXiv:cond-mat/0209372
(http://xxx.lanl.gov/abs/cond-mat/0209372).

\bibitem{lg.eprint3}  E-print, arXiv:physics/0409135
(http://xxx.lanl.gov/abs/physics/0409135).

\bibitem{lww}  G. Liu, E.G. Wang,~and D.S. Wang,~Chin. Phys. Lett. {\bf 14}
, 764 (1997).

\bibitem{lwssc}  G. Liu,~and E.G. Wang,~Solid State Commun. {\bf 105}, 671
(1998).

\bibitem{lchuzhan}  G. Liu,~{\it A New Molecular-dynamics Theory for
Constant Pressure and its Simulation Applications} (Post-doctoral Research
Report, {Institute of Physics, Chinese Academy of Sciences, Beijing, 1998).}
\end{references}

\end{document}

cond-mat0505251.extracted_bib
\bibitem{al}  M.P. Allen,~and D.J. Tildesley, {\it Computer Simulation of
Liquids (}Oxford University, Oxford, 1987).


\bibitem{fre}  D. Frenkel,~and B. Smit,~{\it Understanding Molecular
Simulation (}Academic, New York, 1996).


\bibitem{haile}  J. M. Haile,~{\it Molecular Dynamics Simulation Elementary
Methods (}John Wiley \& Sons, New York, 1992).


\bibitem{a}  H.C. Andersen,~J. Chem. Phys. {\bf 72,} 2384 (1980).


\bibitem{pr1}  M. Parrinello,~and A. Rahman, Phys. Rev. Lett. {\bf 45,} 1196
(1980).


\bibitem{pr2}  M. Parrinello,~and A. Rahman, J. Appl. Phys. {\bf 52,} 7182
(1981).


\bibitem{lg.eprint2}  E-print, arXiv:cond-mat/0306040
(http://xxx.lanl.gov/abs/cond-mat/0306040).


\bibitem{ik}  J.H. Irving and John G. Kirkwood, J. Chem. Phys. {\bf 18}, 817
(1950).


\bibitem{zhou}  M. Zhou, Proc. R. Soc. Lond. A {\bf 459}, 2347 (2003).


\bibitem{lg.eprint1}  E-print, arXiv:cond-mat/0209372
(http://xxx.lanl.gov/abs/cond-mat/0209372).


\bibitem{lg.eprint3}  E-print, arXiv:physics/0409135
(http://xxx.lanl.gov/abs/physics/0409135).


\bibitem{lww}  G. Liu, E.G. Wang,~and D.S. Wang,~Chin. Phys. Lett. {\bf 14}
, 764 (1997).


\bibitem{lwssc}  G. Liu,~and E.G. Wang,~Solid State Commun. {\bf 105}, 671
(1998).


\bibitem{lchuzhan}  G. Liu,~{\it A New Molecular-dynamics Theory for
Constant Pressure and its Simulation Applications} (Post-doctoral Research
Report, {Institute of Physics, Chinese Academy of Sciences, Beijing, 1998).}
