Extended Dynamical Equations of the Period Vectors of Crystals under Constant External Stress to Many-body Interactions

Since crystals are made of periodic structures in space, predicting their three period vectors starting from any values based on the inside interactions is a basic theoretical physics problem. For the general situation where crystals are under constant external stress, we derived dynamical equations of the period vectors in the framework of Newtonian dynamics, for pair potentials recently (doi:/10.1139/cjp-2014-0518). The derived dynamical equations show that the period vectors are driven by the imbalance between the internal and external stresses. This presents a physical process where when the external stress changes, the crystal structure changes accordingly, since the original internal stress can not balance the external stress. The internal stress has both a full kinetic energy term and a full interaction term. It is extended to many-body interactions in this paper. As a result, all conclusions in the pair-potential case also apply for many-body potentials.

where the integral is over all the surface of the bulk, and therefore the bulk has no acceleration. The external stress Γ is assumed to be symmetric, i.e., for all of its components Γ i,j = Γ j,i . This assumption ensures that the net external torque on the bulk is zero. As said previously, the dynamics for the MD ions is always Newton's second law m iri = F i (i = 1, 2, · · · , n), where r i is the position vector of the ith MD ion with mass m i , F i is the net force acting on MD ion i from all other ions of any cell (but no external force on MD ions due to distance from the crystal surface), and n is the total number of MD ions. Then we will derive the dynamical equations for the period vectors in the following. For general purposes, consider 2-body, 3-body, · · ·, up to M -body interactions among any group of ions in any possible configurations. Since these many-body interactions are independent on each other, forces and potentials can be written as a summation of individual m-body contributions. For example, the net force on MD ion i can be expanded as where F (m) i is the contribution of m-body interactions. For identifying an ion in the many-body interactions across the whole crystal effectively, a simplified form of index I k was used for it, so that its position vector can be expressed as where I k,a , I k,b , and I k,c are any values of integers representing the cell in which it resides, and i k , ranging from 1 to n, refers to its corresponding image ion in the MD cell. This means that I k represents the total four independent integer variables of (I k,a , I k,b , I k,c , i k ). A summation over I k means the nested summations over the four corresponding integers. As there are m distinct ions participating in any m-body interaction, the subscript k in I k is used to index the ions from 1 to m in such an interaction. Since no pair of ions can occupy the same physical location, for any pair of indexes I k and I k , the expression (I k,a − I k ,a ) 2 + (I k,b − I k ,b ) 2 + (I k,c − I k ,c ) 2 + (i k − i k ) 2 = 0 is always assumed inside any m-body interaction throughout this article. This also means that for MD ion i k and any other ion I k , the expression (I k,a ) 2 + (I k,b ) 2 + (I k,c ) 2 + (i k − i k ) 2 = 0 is always true, and that for any two MD ions i k and i k , the mutual exclusive relationship i k = i k is always true inside any m-body interaction. Based on Newton's third law, the net force of the m-body interaction in any given m-ion configuration should be zero m k=1 f (m) I k (r I1 , r I2 , r I3 , · · · , r Im ) = 0, (6) where f (m) I k (r I1 , r I2 , r I3 , · · · , r Im ) is the force acting on ion I k by all the rest total m − 1 ions. Further considering the periodicity of the system, the net m-body force acting on all MD ions should also be zero: where i1 (r i1 , r I2 , r I3 , · · · , r Im ). (8) Equation (7) means no internal force can push the system as a whole to accelerate. With Eqs. (3) and (7) combined, it follows that the net of all forces acting on all MD ions is zero, i.e.
where the summation indexes i and i 1 are identical. Employing the centre-of-mass coordinate system of the MD cell for all the work throughout this paper, the total momentum of the MD cell is zero. As the period vectors may change with time, the volume Ω = (a × b) · c and shape of the MD cell and those of the bulk should also change accordingly.

III. INSTANTANEOUS DYNAMICS
In order to find the dynamical equations for the period vectors, imagine a plane P h P h that cuts the model bulk into a right part and a left part, with S h as the area vector of the cross section between the two parts in the direction of pointing to the right part, as shown in Fig.1. Plane P h P h is chosen such that, for a given period vector h = a, b, or c, the right (R h ) part contains all T = T a a + T b b + T c c cells with T h ≥ 0, and the left (L h ) part contains all the rest T cells with T h < 0.
Apply Newton's second law to a "snapshot" of the right (R h ) part. Then, the net external force acting on the R h part is where the integral is over the surface of the bulk in the R h part. Let F L→R be the net force acting on the R h part by the L h part. Then the dynamical equation of the R h part is where M R is the total mass of the R h part andr RC is the acceleration of the centre of mass of the R h part. Since surface effects are neglected, F L→R should be uniformly distributed cell by cell across the section S h between the two parts. Dividing Eq. (11) by where: σ h = ∂Ω/∂h is the (right) surface area vector of a cell with respect to the period h, then arXiv:cond-mat/0505251v6: Extended dynamical equations of the period vectors of crystals ..., Aug. 08,2018 where which is the net force, by the L h part, acting on the "half-line-cell" bar B h composed of the MD cell and cells h, 2h, 3h, 4h, etc., till the surface, as shown in Fig.1. Using Eq. (9), the left hand side of Eq. (13) becomes where the total cell mass is M cell = n i=1 m i and the nested summations of T∈R h n i=1 mean all ions in the R h part are counted. Noticing thatT = T aä + T bb + T cc , Eq. (15) may be written as: where In the R h part, T h is always non-negative, but for any T h =h , it is assumed there exists another −T h that cancels it in the above summation. Therefore, all non-diagonal terms α h,h =h are neglected. Then Eq. (13) becomes Considering all many-body interactions, the net force F h in Eq. (18) can be written as: where F The m-body interaction between the right and left part of the crystal means that the participating ions must be distributed in both parts, namely that not all participating ions are in the same part. Then F (m) h is the net force on ions in the right part of all such possible configurations divided by N h . For total t ions (m > t ≥ 1) in the right part (the rest of the ions are in the left part at the same time), the corresponding net force for all possibilities is where (···) {···} denotes the nested summations over indexes in {· · ·} with conditions in (· · ·), and commutability among ions in each part is considered. If there are indexes listed in the condition (· · ·) expression, the condition applies to all of them, otherwise applies to all the corresponding indexes listed in {· · ·}. For example, {i,j,k} restricts i < 0 and j < 0, while (positive) {i,j,k} requires i > 0, j > 0, and k > 0. If there is only one index in the condition expression, the brackets may be omitted. Throughout this article, all layers of nested summations should be realized into reasonable forms even for special situations. For example, for k = 1: while for k = m: Remembering that where h , h are also period vectors with possible values (h, h , h ) = (a, b, c), or (b, c, a), or (c, a, b) only, considering the crystal translatability, employing N h = I 1,h I 1,h 1, and setting I 1,h = I 1,h = 0, Eq. (20) becomes Translating the system so that the cell containing ion I 1 , which is I 1,h h = 0h, 1h, 2h, 3h, · · ·, becomes the MD cell, Eq. (24) can be further written as: and Eq. (19) becomes: i1 (r i1 , r I2 , r I3 , · · · , r Im ). (26) Considering m-body potential ϕ (m) (r I1 , r I2 , · · · , r Im ), and supposing only s (m ≥ s ≥ 1) of the m ions are in the MD cell (all other ions are outside), where only a fraction s/m of the potential belongs to the cell, the sum of all such potential belonging to the cell is: Since the set of values s = s − 1 = 0, 1, 2, · · · , m − 1 means all possible situations where all ions, except ion i 1 (kept inside the cell), are placed inside or outside of the cell, one has: {I2,I3,···,Im} Then the m-body cell potential energy becomes: As a result, the total up to M -body cell potential energy is: Making use of Eq. (4), take the derivative: with the force: The right side of Eq. (31) can be split into two terms based on the sign of I m,h , so that: where: and By making use of the translatability to move the system so that the cell Renaming ion i m as ion I 1 and ion I 1 as ion i m , then Expand it with respect to t , the number of ions distributed in the part of the crystal defined by where t is actually equal to t + 1. Now let us make use of the translatability to move the system so that cell I m,h h + I m,h h + I m,h h , in which ion I m resides, is translated to the MD cell, then Eq. (35) becomes Expand the above equation with respect to t = t − 1 , the number of ions distributed in the part of the crystal defined by I k,h > l , of total m − 2 ions indexed from I 2 to I m−1 , then it changes into Meanwhile, employing Eq. (6), the first line of Eq. (20) can also be written as where in the last line, I m,h = I m,h = 0, which means the cell containing ion I m can be and only be I m,h h = −1h, −2h, −3h, · · ·. Translating the system so that the cell containing ion I m becomes the MD cell, Eq. (41) becomes: Combining Eqs. (40) and Eq. (42), then As a result Let us define the main interaction tensor for up to M -body interactions as then where h · σ h = Ω and h · σ h = h · σ h = 0 are used. Then Eq. (18) becomes The dynamical equation Eq. (48) is essentially the same as in previous work [27], for only constant external pressure being considered.

IV. MICROSCOPIC TRANSLATED STATES
As seen in Fig. 2, the two distinct states of the system are exactly the same in all microscopic details except being translated relative to each other. In these states all right parts from plane P h P h have the same number of cells conceptually. Since they cannot be distinguished from a macroscopic point of view, an unweighted average of Eq. pass through the MD cell in Fig. 1. For clarity, P h P h was used in its original meaning (i.e., with fixed position), but Q h Q h was employed for such a plane running from left to right. For simplicity, R h from Fig. 1 was separated into two parts: T h = 0 slab made of all cell of T = T a a + T b b + T c c with T h = 0, and R h part made of the rest of the cells.
As an example, for h = a, T h = 0 slab includes all cells of T = 0a + T b b + T c c with T b and T c being any integers.
Since it makes a difference only if Q h Q h meets some ion(s) when it runs, we will only consider the situation where there are s ≥ 1 ion(s) in the T h = 0 slab participating in the interactions. We will consider the following three cases in sequence.
The first case is that there are other t ≥ 1 ion(s) participating in the m-body interaction appearing in the L h part. When the plane Q h Q h runs from left to right passing through the MD cell, the probability for MD ion i k appearing on the left side of Q h Q h is Iµ (r I1 , r I2 , r I3 , · · · , r Im ) should be excluded from F h , as it was unconditionally included in the F h previously. The above equation can be simplified by translating the cell where ion I 1 resides, in the T h = 0 slab, to the MD cell (then I 1 becomes i 1 and N h can be reduced to 1 by removing summations over I 1,h and I 1,h ) where ion index numbers were re-assigned. In the above expression, there are t ≥ 1 ion(s) in the L h part. The situation with t = 0 ion(s) in the L h part is considered in the other two cases. There are s ≥ 1 ion(s) in the T h = 0 slab. As ion i 1 is in the MD cell, there are additional s − 1 ion(s) in the T h = 0 slab. Actually the above expression is valid for all situations with s = s − 1 = 0, 1, 2, · · · , m − t − 1 ions running anywhere in the T h = 0 slab (except r i1 position), and accordingly m − t − 1 − s ions running anywhere in the R h part. From total m − t − 1 ions numbered t + 2, t + 3, . . . , m running anywhere in the R h part except r i1 position, considering all the possible situations placing s ions into the T h = 0 slab and putting all the rest ions in the R h part, the following equality emerges Combining Eq. (51) and Eq. (52), The second case is that all ions participating in the m-body interaction reside in R h part, thus there is no ion in the L h part, and there must be at least one ion in the T h = 0 slab and at least one ion in the R h part. The situation with no ion in the R h part will be considered in the last case. Consider s ≥ 1 ion(s) in the T h = 0 slab and the remaining total m − s ≥ 1 ion(s) in the R h part. For a given plane Q h Q h cutting the MD cell, the net force acting on the ions on the right side of Q h Q h by those on the left side should be added to F h . Recalling Eq. (6), equivalently the net force acting on the ions on the left side of Q h Q h by those on the right side should be subtracted from F h . With ion probabilities appearing on the left side of Q h Q h considered, the averaged net force on them can be written as Iµ (r I1 , r I2 , r I3 , · · · , r Im ) The last case is that all ions participating in the m-body interaction are in the T h = 0 slab. For a given plane Q h Q h cutting the MD cell, the net force acting on the ions on the right side of Q h Q h by those on the left side should be added to F h . Based on Eq. (6), equivalently the net force acting on the ions on the left side of Q h Q h by those on the right side should be subtracted from F h . As a matter of fact, only when ions are distributed on both sides of plane Q h Q h , such forces should be considered. For a given configuration of the m ions, assuming ion I m is the last one to be crossed by plane Q h Q h when it runs from left to right, the probability for ion I k appearing on the left side of plane Q h Q h and ion I m on the right side is (r im − r i k ) · σ h /Ω, then the averaged net force of the given configuration to be subtracted from F h is Iµ (r I1 , r I2 , r I3 , · · · , r Im ).
Similarly, Eq. (53) is valid for t = 1, 2, · · · , m − 1 in the first case, while Eq. (59) is essentially the situation with t = 0, namely For fixed r i1 , only t = 0, 1, 2, · · · , m − 1 of the remaining m − 1 ions can appear in the L h part with the remaining ion(s) in the R h part at the same time, then As a result, where as defined in Eq.
is the net m-body force acting on MD ion i 1 by all other m − 1 ions in all possible configurations.
By using Eq. (7), Eq. (62) can be reduced as Then the averaged net force acting on the half-line-cell bar B h by the L h part in Fig. 1 is where Eq. (3) is used. Now, let us introduce another tensor which is zero when an equilibrium state is reached. This provides where the full interaction term of the internal stress is It can also be written as where DOF refers to all degrees of freedom of the system including the three period vectors a, b, c, and all MD ion position vectors r 1 , r 2 , · · ·, and r n . Then the period dynamics Eq. (48) can be updated into

V. MOMENTUM TRANSPORTATION
Consider an ideal gas in a fixed and closed container in an equilibrium state from a macroscopic point of view and imagine to cut it into a left half and a right half. Gas particles carrying their momentum can freely run between the two halves. Then we have two choices to study it.
One choice is to employ a material-based system. In such a system definition, the gas particles always belong to the same half system, which they belong to at the very beginning. Then at the very beginning, we have very clear half systems of gas particles. However very soon, some gas particles in one half may move into the other half, but still belong to the original half system. Then the half systems would no longer have a clear boundary. Definitely, Newton's second law still applies to the half systems, but not easy to use.
The other choice is to employ a space-based system definition [32], in which at any time, a particle belongs to the system if it is inside the corresponding space with a fixed and close geometric boundary, otherwise it is not. Then each half system is actually defined by the corresponding fixed half space inside the container, then always has a clear boundary. When a gas particle moves from one half into the other half, it leaves from the former system and joines the later system. The later system gets its momentum and the former system gets its momentum as well but in the opposite direction. For the dynamical process of each space-based half system, the total regular force we see is the net external force acting on the gas particles by the container during collisions between them, which is not zero at non-zero absolute temperature. However, the total momentum of each half system is always zero. Then the net momentum transported into and out of the half system per unit time, due to gas particles' crossing the boundary between the two halves, should also be considered as an external force acting on the half system in order to satisfy Newton's second law. As a matter of fact, these two forces balance each other. Let us call the later as the force associated with momentum transportation. Since a specific momentum transported from space-based system S A to its neighbour space-based system S B per unit time should be regarded as an external force acting on system S B by system S A , its opposite direction momentum movement rate from S B to S A should be regarded as another external force acting on system S A by system S B . They are actually action and re-action forces satisfying Newton's third law.
Similarly, if a material-based system is used to study each half part of the crystal above, when an ion "runs from one part into the other part of the crystal", it still belongs to the original part, then the motion of every individual ion must be traced all the time. Furthermore the corresponding components of the external stress Γ, acting on the surface of one part of the crystal, should also be identified acting on the ions, which belong to the other part. If a space-based system is employed, these are not needed, but the force associated with momentum transportation should be considered as an external force on the system. As in our previous work [30], for the L h and R h parts in Fig. 1, both as space-based systems, consider the above statistics over the indistinguishable translated states with the help of plane Q h Q h again, but of the force associated with momentum transportation. If the total amount of such indistinguishable translated states is assumed as the cell volume Ω, the amount of those where MD ion i can cross plane Q h Q h during a unit time is |ṙ i · σ h |, with momentum m iṙi being carried each. Then the additional averaged force associated with momentum transportation on R h part should be added to F h . As a result, Eq. (66) is updated to where the instantaneous kinetic-energy term of the internal stress is Defining the instantaneous internal stress as the period dynamics Eq. (69) becomes The observable period vectors showing fixed values under certain external conditions (e.g. constant external pressure and temperature) should not depend on the directions of ions' motions. A further unweighted average of Eq. (74) was performed over all moving directions of the MD ions. For this, the averaged Eq. (72) becomes: where E k,M D,ion is the total kinetic-energy of the MD ions. Also considering the motion of the valence electrons the same way, the averaged kinetic-energy term of internal stress should be: where E k,M D,ve is the total kinetic-energy of the valence electrons in the MD cell. The forces corresponding to this part of the internal stress should be balanced by the part of the external forces involved in collisions between the ions in the bulk surface and the surrounding external walls, as in the above example of an ideal gas. Accordingly, the averaged internal stress from Eq. (73) is Then the period dynamics Eq. (74) changes into α h,hḧ = (Π + Γ) · σ h (h = a, b, c).

VI. SUMMARY AND DISCUSSION
Keeping Newton's second law for MD ions and applying it to macroscopic half-systems with additional statistics over indistinguishable translated states and forces associated with momentum transportation applied, we arrived at the coupled dynamical equations, Eq. (2) for MD ions and Eq. (78) for the period vectors, of crystals of many-body interactions under constant external stress. Equation (78) shows that the system period vectors are driven by the imbalance between the internal and external stresses. Then when the system reaches an equilibrium state, the internal and external stresses balance each other. The internal stress has both full kinetic-energy and full interaction terms. As a result, the dynamical equations and associated formulas in this article for many-body interactions share the same form of those in our last work [30] for pair-potential only.
The kinetic-energy term was obtained from the statistics of forces associated with momentum transportation when the two halves of the system are recognized as space-based ones. Since the full interaction term of the internal stress Eq. (68) is valid for any-body interactions, it should also be valid for forces from electrons but calculated based on quantum mechanics involved. In such a situation, the effective interactions among ions through electrons are many-body ones, as the calculated state of electrons depends on the positions of all ions.
As a matter of fact, the external stress is required as a constant only in deriving Eqs.
(1) and (10) and this requirement only means that it is constant over the surface of the crystal throughout this paper. Then for such external stress but changing with real time, one can solve the crystal with Eqs. (2) and (78) to reach an equilibrium state iteratively for the given external stress at a given real time, then the next real time, ..., till end.
In the MD world, simulations are usually classified into various ensembles, based on applicable combinations of fixed volume, constant external pressure, and constant external temperature. For ensembles of fixed volume, Eq. (78) shows that an external stress balancing the internal stress should always be supplied or assumed. For ensembles of constant external pressure/stress, Eq. (78) can be used. Additionally, the straightforward ion speed rescaling method can always be a choice, for constant external temperature simulations.