Work conjugate pair of stress and strain in molecular dynamics

ABSTRACT Certain stress and strain form a thermodynamic conjugate pair such that their strain energy equals to a scalar-valued potential energy. Different atomistic stresses and strains are analytically derived based on the work conjugate relation. It is numerically verified with both two-body and three-body potentials that the atomistic Kirchhoff stress, first-order Piola–Kirchhoff stress and second-order Piola–Kirchhoff stress are conjugates to atomistic logarithmic strain, deformation gradient and Lagrangian strain, respectively. Virial stress at 0 K based on original volume is the special form of atomistic Kirchhoff stress for pair potential. It is numerically verified that Hencky strain is not conjugate to any stress.


Introduction
Stress tensor and strain tensor are fundamental in continuum mechanics (CM). In finite strain theory, stress and strain can take different forms that can often be converted to each other. Certain stress tensor and strain tensor form a conjugate pair such that the stress integrates over the strain yields a scalar-valued potential energy [1]. Only the conjugate pair of stress and strain should be chosen as dependent and independent constitutive variables, respectively, otherwise the existence of a scalar-valued potential energy is not guaranteed. The following conjugate pairs, (1) second-order Piola-Kirchhoff stress (PK2) and Lagrangian strain, (2) Hencky strain and Kirchhoff stress [2][3][4], (3) firstorder Piola-Kirchhoff stress (PK1) and deformation gradient, have been mentioned. It is worthwhile to note that Cauchy stress tensor is not conjugate to any strain tensor, therefore it may not be directly used in constitutive equation.
As an effort to bridge the gap between molecular dynamics (MD) and CM, it is vitally important to define stress and strain at atomistic level. Large-strain tensorial definitions of stress and strain at atomistic level allow the results of MD be expressed in the form of stress-strain relation, which can be directly fed into the material model of CM [5]. This sequential multi-scale approach may potentially replace material testing in the development of material model for finite element analysis [6,7].
CONTACT James D. Lee jdlee@gwu.edu The concepts of stress and strain are well understood in CM. The atomistic definitions of stress and strain, on the contrary, are still in the developmental stage [8][9][10][11]. Notice that the definitions of the stress and strain in CM require the material points to be infinitesimally small, which is certainly not true in MD due to its discrete description and the finiteness of the atomistic mass. To formulate stress tensor in CM, the concept of the mass of a material point needs to be replaced by that of mass density. The force equilibrium leads to the conclusion that stress follows tensorial transformation rule. Now let us formulate the stress tensor not based on the concept of infinitesimal mass point. The balance of moment is applied (cf. Equation (1)) over the material of a finite size tetrahedron (cf. Figure 1) ρaΔv ¼ t ðnÞ Δa À t k Δa k þ ρfΔv (1) where ρ is the mass density at some interior point in volume v, a is the acceleration vector at some interior point of v, f is the body force, t ðnÞ and t k are the surface traction per unit area at some point on surface area Δa and Δa k , respectively. Divide both side by Δa where n k is unit normal.

Assumption 1
In the limiting case, the ratio of volume over surface area is vanishingly small, i.e. Δv=Δa ! 0. Therefore, the inertia effect and body force effect (last two terms of Equation (2)) can be neglected. This generally means the size of the system has to be relatively small.

Assumption 2
In the static equilibrium case with no body force and acceleration, the last two terms of Equation (2) vanishes. Accepting either Assumption 1 or Assumption 2 leads to The stress tensor t kl is the l-th component of the stress vector t k acting on the positive side of the k-th coordinate surface. In other words, where i l is the unit vector along l-axis. Because both t k and i l are tensorial quantities, so is stress tensor t kl [1]. Without such assumption, there is no rigorous definition for atomistic stress.
In the derivation of strain in CM, the displacement field has to be continuous and smooth to enable the calculation of deformation gradient. In MD simulation, the displacement field is discrete; without further assumption, the rigorous definition of strain does not exist. We thus assume that all atoms follow a single uniform deformation gradient [12]. Under such assumption, all forms of strains, such as Lagrangian strain and Eulerian strain, can be rigorously defined and calculated.
Many researchers have defined various atomistic stresses [10,[13][14][15][16][17]. The most commonly used atomistic stress is virial stress of statistical mechanics. The derivation process offers no clue to tell which is the equivalence of virial stress in CM. Without knowing its equivalence, it is impossible to choose its conjugate strain in the constitutive relation.
Remark: We recall that in thermoelasticity, a branch of CM, the stress-strain-temperature relation can be expressed as where σ; e and T are Cauchy stress tensor, strain tensor and temperature variation, respectively. On the other hand, the atomistic stress in MD has two parts: the potential part and the kinetic part, which are corresponding to the two parts in the abovementioned Cauchy stress. In this work, only the potential part is considered, i.e. it is a static case of absolute zero temperature. The detailed expression of atomistic stress will be presented and discussed later.
In this article, the atomistic stress is derived based on the work conjugate relation and atomistic strain. We assume there is a single uniform deformation gradient for all atoms, therefore the strain can be defined unambiguously. The work conjugate relation requires that such stress and strain, upon integration, give total potential energy of the system, which can be readily obtained through interatomic potential. The continuum counterpart of atomistic stress and strain is identified. The conjugacy of atomistic stress and strain is numerically verified with both two-body potential and three-body potential.

Atomistic strain
The calculation of atomistic strain is based on the deformation gradient followed by a group of atoms. Practically, a single deformation gradient F may be calculated for a group of atoms by minimizing the mapping error between the deformed and the undeformed configurations [12] F ¼ min where X n and x n are the original and deformed positions of atom n; ðn ¼ 1; 2; 3:::NÞ.
Notice that a single F does not always satisfy x n ¼ FX n for a group of atoms [18]. Therefore, finding the single F becomes a problem of optimization. Remark: It is noticed that the Cauchy-Born rule is commonly used as a hypothesis to relate the motion of atoms in a crystal to the overall deformation of the continuum. However, for complex crystals, such as diamond or mono-layered crystalline (graphene, MoS 2 ), the rule needs to be modified to allow for internal degrees of freedom between the sublattices [19][20][21][22]. In this work, we follow the original idea of Cauchy-Born rule, i.e. the motion of all atoms (same kind of atoms) follow a single and uniform deformation gradient.
The Lagrangian strain can be obtained as where I is the identity matrix. The atomistic logarithmic strain L can be calculated as follows [18]. Decompose the deformation gradient F as where log λ i is the natural logarithm of λ i . Define asymmetric logarithmic strainL aŝ For abbreviation, the operation described by Equations (7)-(10) are named as 'logm'.
The logarithmic strain L is defined as the symmetric part ofL Therefore Equation (12) can also be expressed as [18] L; 1 2 logmðFÞ þ logmðFÞ T All eigenvalues λ 1 ; λ 2 ; λ 3 of deformation gradient F need to be positive in the calculation of logarithmic strain. Noticed that the natural logarithm of a negative real number is not defined (cf. Equation (9)). The atomistic logarithmic strain is objective and symmetric and is conjugate to virial stress in the case of two body potential [18] (cf. Appendix 1). Notice that the atomistic logarithmic strain is not equal to the Hencky strain. In Hencky strain calculation, one may perform polar decomposition for the deformation gradient where R is the rotational matrix. Then the Hencky strain is defined as The logarithm is performed on the deformation gradient without the rigid body rotation [2,4]. Hencky strain is not conjugate to virial stress [18].

Atomistic stress
In CM, certain stress tensor and strain tensor form conjugate pair such that the strain energy, calculated by the integration of the stress over the strain, equals to the potential energy density [1]. For instance, first-order Piola-Kirchhoff stress P and deformation gradient F form a conjugate pair Second-order Piola-Kirchhoff stress K and Lagrangian strain E form a conjugate pair; Kirchhoff stress S and logarithmic strain L forms a conjugate pair [18] where V is the volume of the system in its undeformed state, U is the strain energy of the system. In MD, system potential energy is explicitly given as the interatomic potential. Since a single deformation gradient can be calculated for a group of atoms, the stress can be determined by its corresponding strain through the work conjugate relation. In MD, where U is the interatomic potential energy, f m is the interatomic force acting on atom m due to the interaction with all other atoms in the system. Consider a system of Natoms, which has total potential energy where f m is the total force acting on atom m; x m is the current position of atom m.
Assume the motion of all atoms follows a single uniform deformation gradient F ij from its original position X m to their deformed position x m , ðm ¼ 1; 2; 3:::NÞ, it leads to Recall Equation (16) as Equate Equations (23) and (24) The atomistic first-order Piola-Kirchhoff stress is thus defined as In this work, we adopt the first-order Piola-Kirchhoff stress tensor defined by Truesdell [23][24][25], which is different from the definition given by Erigen [1,26,27] T ;J F À1 Á σ where J is the Jacobian, σ is the Cauchy stress tensor. Notice that Recall the relation between first-order Piola-Kirchhoff stress P and second-order Piola-Kirchhoff stress K Thus the atomistic second-order Piola-Kirchhoff stress is obtained as The second-order Piola-Kirchhoff stress is symmetric, therefore Equation (32) can also be written as In CM, the relation between Kirchhoff stress S and the first-order Piola-Kirchhoff stress P is Recall Equation (26) Therefore the atomistic Kirchhoff stress is defined as Due to the symmetric nature of Kirchhoff stress, Equation (37) can be also written as The Cauchy Stress and Kirchhoff stress are related as Therefore the atomistic Cauchy stress becomes where v is the volume of the system in its deformed state. Because Cauchy stress is symmetric, one may write Notice that there is no assumption that the interatomic potential is a pair potential in the derivation of Equations (27), (32) and (36). In fact, it is numerically verified that, with two body potential (Coulomb-Buckingham and Lennard-Jones) and three body potential (Tersoff), the atomistic first-order Piola-Kirchhoff stress (cf. Equation (27)) is conjugate to deformation gradient; atomistic second-order Piola-Kirchhoff stress (cf. Equation (32)) is conjugate to Lagrangian strain (cf. Equation (6)); atomistic Kirchhoff stress (cf. Equation (36)) is conjugate to atomistic logarithmic strain (cf. Equation (12)). The calculated strain energy is identically equal to the interatomic potential. It is also numerically tested that the Hencky strain as defined in Equation (15) is not conjugate to any atomistic stress.

Special case: atomistic stress for pair potential
Recall the system potential energy (cf. Equation (20)) Notice that f m is the summation of all forces acting on atomm considering all influence within the system. It is not necessary that f m can be written as the summation of all contributing pair-wise forces f mn as shown in Equation (43). For instance, the forces in Tersoff potential contain the influence of the third atom [28][29][30].
In the special case of two-body potential, the interatomic force f m on atomm can be written as where f mn is the interatomic force acting on atomm due to the interaction with atomn. In case of pair potential, the atomistic first-order Piola-Kirchhoff stress becomes (cf. Equation (27)) Because m and n can be any atom within the group, Equation (44) remains valid if m is interchanged with n In three-body potential It is numerically verified that Equation (47) is not valid in three-body Tersoff potential (cf. Appendix 2) .
In two-body potential, however Add Equations (44) and (48) which givesP, the second form of atomistic first-order Piola-Kirchhoff stress in two-body potential. It is numerically verified thatP ¼ P only in the case of two-body potential. It can be further derived from Equation (49) that In two-body potential, the atomistic force vector f mn is always parallel to the relative position vector Therefore, Thus, it arrives at the third form of first-order Piola-Kirchhoff stress for two-body potential. Similarly, the atomistic second-order Piola-Kirchhoff stress also has its twobody form.
Interchange m and n Add Equations (57) and (59) which leads to the second form of atomistic second-order Piola-Kirchhoff stress for pair potentialK Because force and relative position vector are always parallel to each other in pair potential, which leads to the third form of atomistic second-order Piola-Kirchhoff stress for pair potential.K The atomistic Kirchhoff stress also has its two-body form Interchange m and n Add Equations (68) and (70) which leads to the second form of atomistic Kirchhoff stress for pair potential Because force and relative position vector are always parallel to each other in pair potential, which leads to the third form of atomistic Kirchhoff stress for pair potential Notice that the third form of atomistic Kirchhoff stress is identical to virial stress at absolute zero temperature. Thus, it is concluded that virial stress at 0 K based on original volume is the special form of atomistic Kirchhoff stress for pair potential.

Numerical verification
A numerical process is proposed to verify the conjugacy between different pairs of atomistic stress and strain [18]. Assume N atoms move from their original positions X m i to their current position x m i following a uniform deformation gradient F, which specifies where m ¼ 1; 2; 3 Á Á Á N. In numerical process, the deformation can be divided into k incremental steps þ δ ij ; k 2 0; 1; 2; 3; :::: Interatomic potential can be calculated for given positions of all atoms. Define system potential energy k P of step k as where U is the interatomic potential. The strain energy of the system can be obtained as where Vis the volume of the system, t ij and ε ij represents any conjugate pair of atomistic stress and strain, respectively. Knowing the history of all previous steps, the strain energy k E of step k can be further derived as with the understanding that 0 t ij ¼ 0, 0 ε ij ¼ 0. The conjugacy of the chosen t ij and ε ij is verified if the strain energy calculated by Equation (82) is always identical to the interatomic potential calculated by Equation (80), i.e.
It is worthwhile to emphasise that the k P refers to the summation of interatomic potential energy of all atoms and k E refers to the strain energy of the entire specimen.
The MATLAB® code used for verification can be found in Appendix 4. In Table 1, we list general atomistic stress-strain pairs whose work conjugacy have been numerically verified. The interatomic potential used for verification is three-body Tersoff potential.
In Figures 2-4, it is seen that the interatomic potential are identical to the strain energy obtained through numerical integration of atomistic stress-strain pair. In each of the three cases, the initial atomic configuration and deformation gradient are randomly created. It is numerically verified again in Figure 5 that Hencky strain is not conjugate to atomic Kirchhoff stress. The number of atoms, N, used in Figures 2-4 is 10. Actually it was Table 1. Conjugacy verification for general atomistic stress-strain pair.

General case
General atomistic stress Atomistic strain Conjugacy verification PK1/deformation gradient Kirchhoff stress/ Hencky strain Figure 2. Interatomic potential energy compared with strain energy calculated by the atomistic firstorder Piola-Kirchhoff stress P and deformation gradient F, verified by using Tersoff potential.
numerically verified that the conclusion about conjugacy is independent of N. Also, as an observation, the conjugacy is valid regardless whether the atoms are in their equilibrium position, i.e. work conjugacy is independent of initial stress or initial potential energy.
In Table 2, we summarize the conjugacy verification of atomistic stress and strain in case of pair potential. The interatomic potential adopted here are Coulomb-Buckingham, Lennard-Jones potential and an artificial created interatomic potential.
Also, in case of pair potential, it is numerically verified that P ¼P ¼P (84) Thus, in summary, it is concluded that the strain energy calculated from the conjugate pair of atomistic stress and strain (cf. Tables 1 and 2) equals to the interatomic potential energy for any given arbitrary deformation gradient and initial atomic configuration, i.e.

Conclusion and discussion
The atomistic first-order Piola-Kirchhoff stress is analytically obtained through work conjugate relation with deformation gradient. It is numerically verified that, for both pair potential and three-body potential, the atomistic first-order Piola-Kirchhoff stress P, atomistic second-order Piola-Kirchhoff stress K and atomistic Kirchhoff stress S are conjugates to deformation gradient F, Lagrangian strain E and logarithmic strain L respectively. The conjugacy for atomistic Kirchhoff stress and logarithmic strain is valid when deformation gradient has all positive eigenvalues, otherwise the logarithmic strain L is not defined. For all other pairs of atomistic stress and strain, the only requirement for conjugacy is that the deformation gradient F is real-valued. Hencky strain H is not conjugate to any atomistic

Verified
Kirchhoff stress 2nd form /logarithmic strainS Kirchhoff stress 3rd form(Virial)/logarithmic strainŜ þ logmðFÞ T Verified stress [18]. The conjugate atomistic stress and strain pair yields interatomic potential energy upon integration. It is noticed that the virial stress at 0 K based on the original volume is a special form of atomistic Kirchhoff stress for pair potential; virial stress based on the current volume is a special form of atomistic Cauchy stress. In this work, the atomistic system used to calculate stress and strain involves all existing atoms. In the real-world application, a region has to be defined in order to calculate the stress and strain of the atoms within. A system may have several regions and may be much larger than the region specified. How to treat the atoms outside of the region and how to calculate the atomistic stress rigorously remain as an open question.

Disclosure statement
No potential conflict of interest was reported by the authors.

Notes on contributors
Leyu Wang is a research assistant professor at the Center for Collision Safety and Analysis, George Mason University. He received his PhD degree in mechanical engineering from the George Washington University in 2015. His research interests include fracture mechanics, solid mechanics, molecular dynamics, multi-scale modeling, dynamic finite element analysis, material constitutive modeling, and automobile safety.
James D. Lee is a professor in the Department of Mechanical and Aerospace Engineering, George Washington University since 1990. He is a fellow of ASME. He received his PhD from Princeton University in 1967. Since then, he has been teaching and/or doing research in universities (Purdue U., West Virginia U., U. of Akron, U. of Minnesota and GWU), industry (The General Tire and Rubber Company), and US government laboratories (NIST and NASA). His research fields include liquid crystals, fracture mechanics, composite materials, optimal control theory, robotics, finite element and meshless methods, molecular dynamics, biomechanics, and microcontinuum physics. Most recently, he has been doing research on multiple time/length scale modeling from atoms to continua. He has been a principal investigator of research grants awarded by NASA, NSF, and US DOT. He has published two books and more than 120 journal articles. While at NCAC he led efforts to develop vehicle structural and occupant models, conduct evaluations of vehicle crashworthiness, assess bumper-height compatibility issues in collisions between large and small vehicles, analyze roadside hardware safety, and conduct component and full-scale crash tests to gather data for formulating models and/or validating simulation results. Dr Kan was a pioneer in research on the implementation of the nonlinear explicit FE code on highperformance parallel computer platforms. His research and development efforts in benchmarking of high-performance computing platforms for crash simulations have become standards used today in the automotive industry. Under his leadership the NCAC landed multimillion dollar contracts that involved more than 60 task orders for testing, modeling building and validations, analysis of highway barrier crashworthiness, vehicle safety, occupant risk analysis, outreach programs, and crash information management.

Appendix 1. Remark on the objectivity of logarithmic strain and Hencky strain
In this appendix, a MATLAB® Code is presented to show that Hencky strain, logarithmic strain and asymmetric logarithmic strain are objective tensor quantities with respect to Galilean transformation [31] and they are not equal to each other in general.
Because the reference Lagrangian configuration is affected by the change of the frame of the Eulerian coordinate, the Galilean objectivity of deformation gradient F as defined [31,32].
where Q is a time independent transformation matrix. Notice that the deformation gradient F should reduce to identity matrix I at the initial undeformed stage, therefore the Lagrangian and Eulerian coordinates are essentially related. The same formula is also derived by Murdoch with different methodology [33]. First, we polar decompose the deformation gradient F as [34] where V is the left stretch tensor; U is the right stretch tensor; R is the rotational matrix with Then the left stretch tensor V and the Hencky strain are calculated as where the operation 'logm' is defined in Equations (7)- (10). If the Hencky strain H is Galilean objective, the following equation should be valid.
where H Ã , V Ã and F Ã stand for the Hencky strain, left stretch tensor and deformation gradient, respectively, in the transformed Eulerian coordinate. Similarly, the asymmetric Logarithm L is Galilean objective if the following equation is valid whereL Ã is the asymmetric logarithmic strain in the transformed Eulerian coordinate and the operation 'logm' is defined in Equation (7)- (11). The logarithmic strain L as defined in Equation (13) is Galilean objective if the following equation is valid.
It is noticed that, in general LÞH;LÞH; LÞL (A4) The following is a MATLAB® code, which enables one to prove numerically that Hencky strain, logarithmic strain, asymmetric logarithmic strain are objective. The Tersoff potential may be expressed as [28] V The atomistic forces can be obtained as It is seen that ðcos θ ijk Þ ;i þ ðcos θ ijk Þ ;j þ ðcos θ ijk Þ ;k ¼ 0, which implies, among any three atoms, i; j; and k, Tersoff potential yields Thus it is proved that in general