Application of atomic stress to compute heat flux via molecular dynamics for systems with many-body interactions

Although the computation of heat flux and thermal conductivity either via Fourier's law or the Green-Kubo relation has become a common task in molecular dynamics simulation, contributions of three-body and larger many-body interactions have always proved problematic to compute. In recent years, due to the success when applying to pressure tensor computation, atomic stress approximation has been widely used to calculate heat flux, where the lammps molecular dynamics package is the most prominent propagator. We demonstrated that the atomic stress approximation, while adequate for obtaining pressure, produces erroneous results in the case of heat flux when applied to systems with many-body interactions, such as angle, torsion, or improper potentials. This also produces incorrect thermal conductivity values. To remedy this deficiency, by starting from a strict formulation of heat flux with many-body interactions, we reworked the atomic stress definition which resulted in only a simple modification. We modified the lammps package accordingly to demonstrate that the new atomic stress approximation produces excellent results close to that of a rigid formulation.

In this work, we will briefly discuss computation of the local pressure tensor inside a control volume for many-body interactions and the atomic stress approximation, most prominently used by the current versions of the LAMMPS package [16]. We will then show how the atomic stress approximation is applied to heat flux computation, resulting in incorrect heat flux values for three-body and larger many-body interactions. To remedy this, a new atomic stress definition will be derived in accordance to precise formulation of heat flux for manybody interactions. Finally, it will be demonstrated via several test simulations that the new atomic stress approximation produces excellent agreement with results obtained via rigorous computation, while the conventional atomic stress approximation substantially underestimates heat flux and overestimates thermal conductivity.
In standard molecular dynamics simulation, pressure computation is based on the virial theorem. In a periodic system of V volume that consists of N number of particles, the pressure tensor is expressed as where m i and v i represent the mass and velocity vectors of the ith atom and angle brackets indicate time or ensemble mean [47]. The component W is the virial contribution to the pressure and for a system consisting of only pairwise interactions where r i j = r i − r j is the relative position of atom i in relation to atom j and F i j is the force acting upon atom i due to interaction with atom j. For the sake of simplicity, we assume that the system is large enough so that atoms do not interact with their own images. This kind of "i j" notation is essential to produce correct results with periodic systems [48]. Although initially W was only clarified for pairwise interactions, it has been extended for many-body interactions as well. Let us say that there are K groups of many-body interactions among atoms, where a single atom can participate in multiple groups. Thompson et al. provided a concise expression that they named "group" form to obtain the virial contribution as [10] where r k i and F k i represent the position of the ith atom in the kth group and the force acting on it due to group interactions, and i ∈ k indicates the set of atoms participating in the kth many-body interaction group. In the case of periodic boundary condition, the "k" superscript in r k i indicates that the actual coordinate values are either those of the original simulation cell or one of the image cells, selected depending on the group, i.e., r k i might have different values in different groups. Because of this and the fact that total force of a many-body interaction is i∈k F k i = 0, W group in Eq. (3) works correctly with periodic boundary condition, even though it no longer uses the "i j" notation, as detailed in the original paper [10]. It is assumed that many-body potentials only depend on the relative positions of the atoms in the group, whereas long-range electrostatic interactions would require different treatment [49,50].
Unlike the system pressure tensor, several possible expressions for local Irving-Kirkwood-type pressure have been proposed, which do not necessarily produce the same results [11]. A widely used method in molecular dynamics is the method of planes (MOP) [9], where the local pressure tensor of the control plane is obtained from the momentum and force that are transferred through it. Another popular method is the volume averaging (VA) method, where the pressure tensor in the control volume is computed by adding a weighting function to the virial [8]. For pairwise interactions, the weighting function is simply a fraction of the interparticle distance that is inside the control volume. It has been demonstrated that these two methods are equivalent, as VA expression can be obtained by spatially averaging the MOP expression [51]. Although both MOP and VA methods are straightforward for pairwise interactions, many-body interactions require more involved treatment, such as central force decomposition [5,11,12].
Another approach exists, where the virial contribution is distributed to each atom, and is referred to as atomic stress [13,15]. For pairwise interactions, we define atomic stress of the ith atom as where the sum over all atoms is equal to the total system virial contribution W pair = − N i σ pair i . Strictly speaking, Eq. (4) has the dimension of energy and should be divided by some type of atomic volume to give the dimension of stress. The form in Eq. (4) was chosen to match the definition of virial contribution component W pair in Eq. (2), which was itself chosen to match that used by Thompson et al. [10], and lets us avoid dealing with atomic volume such as taking into account the relation between the atomic and system volumes. Note that in other literature the kinetic term, i.e., the first right-hand side term of Eq. (1), is sometimes included, the 1/2 coefficient can be omitted, various definitions of atomic volume are used, and stress might be exchanged to pressure. For many-body interactions, similarly to MOP and VA methods, the more involved central force decomposition can be used [5,11,12]. Although Thompson et al. admitted in their paper that there is no strict physical basis [10], it is often simpler and more convenient to equally distribute the virial from many-body interaction to the each atom in the group in accordance to Eq. (3), where K i is the set of many-body interaction groups which contain atom i, and N k is the number of atoms in the k-th many-body interaction group. This reduces to σ pair i for pairwise interactions. A somewhat crude, but convenient method to assess the local pressure inside a control volume is to simply sum the atomic stress values of the contained atoms where denotes the control volume and either σ pair i or σ group i can be used for σ i . Strictly speaking, this method is inferior to VA, and is known to produce artifacts [14], but is sufficiently accurate when the control volume is set in a homogeneous area, where molecule composition and density are mostly uniform, as the artifacts cancel out. Satisfactory accuracy can also be reached for sufficiently large systems, where the control volume to surface ratio is high, although exact estimates are difficult and would depend on the actual pressure and conditions at the control volume boundary. The expression in Eq. (6) also produces the exact pressure tensor when used over the whole system. This is the scheme used in LAMMPS with σ group i when computing pressure or heat flux. As the pressure tensor expression in Eq. (1) can be derived from momentum conservation law, so can the heat flux vector J of the whole system be derived in a similar manner from energy conservation law as where e i is the kinetic and potential energy of the ith atom, assuming equal potential energy distribution for many-body interactions [3]. We also assume no streaming velocity for simplicity. The Q component is the virial contribution to the total heat flux, and for pairwise interactions where the definition from Eq. (4) was used. This ability to reuse atomic stress values greatly eases the implementation and computational cost of heat flux calculation. The local heat flux in a control volume can be defined in an equivalent manner to Eq. (6): and while similarly to Eq. (6) this is a somewhat crude approximation, sufficient precision can be obtained if the control volume is located in a homogeneous area and only pairwise interactions are present. It might be tempting to use the same approach for many-body interactions and substitute σ pair · v i , as is currently implemented in the present versions of LAMMPS. However, as will be demonstrated later, this approximation introduces significant artifacts and produces incorrect results. Fan et al. have demonstrated these inconsistencies using the Tersoff potential and central force decomposition [7].
Our group has previously published works detailing the heat flux equivalents of the MOP and VA methods for manybody interactions, which do not require decomposition of central force [4,52]. We will use these results to derive an appropriate expression of atomic stress suitable for use in heat flux computations. As before, if we assume equal potential distribution and limit ourselves to many-body interactions only explicitly depending on the relative atom position, the virial contribution to heat flux can be reduced to where script notations correspond to those used in Eq. (3), and r k 0 is the geometric center, i.e., centroid, of the kth group r k 0 = 1 N k i∈k r k i . As can be derived from the original paper [4], this geometric center form is the result of assuming equal potential energy distribution over each atom in a many-body interaction, and although in principle this assumption is not required, it is seemingly most natural and widely used. We will refer to Eq. (10) as the "centroid" form to distinguish it from other forms. We can now define an atomic stress tensor analogous to that in Eqs. (4) and (5) as As in previous expressions, · v i and σ cntr i reduces to σ pair i for pairwise interactions. It must be noted that unlike atomic stress tensors σ pair i and σ group i , σ cntr i is not symmetric, which is necessary to obtain correct heat flux via Eq. (10). On the other hand, as total force due to many-body interaction is zero, this does not affect the virial contribution to the system pressure, i.e., N i=1 σ group i = N i=1 σ cntr i . Also note that while the virial contribution to the total system TABLE I. Conditions for the simulation systems, where N mol is the number of molecules, and L x , L y , and L z are the simulation box dimensions in the corresponding directions, T is the average temperature, and ρ is the system density. The thermal conductivity is also indicated by λ, where parentheses indicate the standard error of mean of the last digit, which was obtained identically to Fig. 1 (2) heat flux Q cntr , which would be used in EMD systems, is a rigorous implementation of Ref. [4], it is no longer equivalent to Ref. [4] when computing the local heat flux in a control volume via Eq. (9), which would be used in NEMD systems, as Q (cntr) ≈ − i∈ σ cntr i · v i is an approximation. To demonstrate the validity of the centroid form atomic stress, we conducted several computations with different systems composed of either butane, octane, or 20 monomer polystyrene molecules with both standard LAMMPS and a modified version to adhere to the centroid form of atomic stress in Eq. (11). For butane and octane, the NERD potential model was used [53]. For polystyrene, the potential model developed by Karanikas and Economou was used [54]. Both of the models are united atom and contain van der Waals, bond, angle, and torsion potentials, with polystyrene also containing an improper potential. The choice of systems was purely pragmatic, as detailed results were available for butane and octane from our group's previous works [52,55] and polystyrene had an improper potential and was available in our in-house program. For each molecule type, a NEMD system with an induced heat flux of 300 MW/m 2 along the z direction and an EMD system were created, with system sizes, molecule numbers, and other properties shown in Table I. The NEMD systems were used to directly compute heat flux inside control volumes and investigate its components, while the Green-Kubo relation for heat flux fluctuation was used in EMD systems to obtain thermal conductivity of the whole systems, which was decomposed into contributions of each potential type via a method proposed by our group [55]. The EMD and NEMD system setup for butane is detailed in Ref. [55], except that the eHEX algorithm was used for inducing thermal flux in the NEMD system [56]. For octane and polystyrene: in the case of NEMD systems, an identical procedure to that of butane was followed, with the equilibration run being at least 5 ns, and in the case of EMD systems, NEMD systems were reused, instead being equilibrated without an induced heat flux. The original intent of not using cubic systems for EMD was to eliminate a difference in thermal transport properties due to different system aspect ratios and changes in initial structure configurations when compared to NEMD systems, but proved to be mostly unnecessary. For all NEMD systems, two control volumes with a length of 0.3L z in the z direction were placed in each half of the systems between the cold and hot temperature regions. Details about temperature settings and resulting system states can be obtained from Refs. [55] and [52] for butane and octane, respectively, while the mean temperatures for all systems are indicated in Table I. Polystyrene resulted in an amorphous state and the cold and hot regions had approximate temperatures of 278 and 322 K, respectively. In the case of butane and octane, the rRESPA multitimescale integrator was used with a 1.0 fs time step for van der Waals interactions and a 0.2 fs time step for remaining potential types [57], while the velocity Verlet algorithm with a time step of 1.0 fs was used for polystyrene systems. Production runs of 150 ns were conducted for all NEMD systems. In the case of EMD systems, production runs of 150 and 90 ns were conducted for butane and octane systems, and for the polystyrene system, respectively, producing the mean of heat flux autocorrelation function over 10 ps.
The results are shown in Fig. 1. The heat flux inside the NEMD systems in Fig. 1(a) was obtained by substituting σ group i or σ cntr i into Eq. (9), and taking the average between the values in the two control volumes. The thermal conductivity of EMD systems in Fig. 1(b) was obtained by substituting Q group or Q cntr into Eq. (7) and applying the Green-Kubo relation [55]. The actual thermal conductivity values in Fig. 1 are the mean of x, y, and z components, as for even the elongated systems of octane and polystyrene the difference between z and x or y components was comparable to that between x and y, being at most 4%, i.e., thermal transport was isotropic in all systems. The reference values for butane and octane were obtained from Refs. [55] and [52], respectively, and values for polystyrene were computed with an in-house program from a 50-ns production run. Also note that due to a technical limitation, the reference of polystyrene displays the sum of torsion and improper contributions (tor+imp). All of these reference values were obtained with rigorous implementation presented in Ref. [4]. Thermal conductivity values are also displayed in Table I, where Fourier's law was used to obtain thermal conductivity for NEMD systems via the temperature gradient inside the control volume regions and the induced heat flux of 300 MW/m 2 , and for the EMD systems, the results obtained via the centroid form of atomic stress were used, corresponding to "Modified" in Fig. 1(b). It should be noted that a non-negligible difference in the measured thermal conductivity exists between the polystyrene NEMD and EMD systems. The exact reason for this is unclear, and might be due to the amorphous nature of the system, as only the average thermal conductivity of the control volume regions was measured. Indeed, for polystyrene the difference between the two regions is about 4.5%, where for butane and octane it is approximately 2.2% and 1.5%, respectively. This is not of great concern, however, as we mainly compare among either NEMD or EMD systems.
We first discuss the decomposition of the heat flux inside the control volumes of NEMD systems in Fig. 1(a). It is important to stress that if the expression of J does not violate energy conservation, the average heat flux computed when substituting the appropriate atomic stress definition σ i into Eq. (9) must be equal to the applied heat flux. Regardless of the molecule type, the total heat flux obtained via the group  [52,55]. A zero axis and reference values are also indicated by a horizontal dashed line where applicable. Horizontally shifting black error bars indicate standard error of mean of each component that was scaled by 2 for clarity, while horizontally centered thicker red error bars at the top indicate that of the total, and were computed via block average method by dividing sampled data into 30 blocks [58]. form atomic stress σ group i with unmodified LAMMPS is below the induced heat flux of 300 MW/m 2 . The reason is obvious when compared with the reference values: the contributions from three-body and four-body interactions are much smaller than the reference. Upon closer inspection, it is apparent that torsion and improper contributions even became slightly negative. On the other hand, van der Waals and bond contribution show a good match, which is not surprising, as all atomic stress definitions in this paper reduce to the same expression at pairwise interactions. Only the polystyrene system has a non-negligible angle contribution, although it is significantly smaller than the reference. The reason for this is not immediately clear, but it is obvious that a simple presence of a non-negligible multibody contribution does not necessarily indicate correctness. In contrast to the "LAMMPS" results, heat flux values, indicated by "Modified" in Fig. 1(a), and obtained via the centroid form of atomic stress σ cntr i show an excellent match with reference values. The polystyrene system does show larger uncertainty than other systems, which is most likely because of its amorphous state that limits atom migration, but a good agreement with the reference is achieved nevertheless. It is clear that the heat flux computed via the group form atomic stress σ group i approximation used by the current versions of LAMMPS significantly underestimates the contribution of three-body and larger many-body interactions, which is critical for systems with substantial contribution from nonpairwise many-body interactions, such as octane or polystyrene in Fig. 1.
Next, we discuss the decomposition of thermal conductivity in EMD systems in Fig. 1(b). Total thermal conductivity of butane and octane obtained via the centroid form of atomic stress σ cntr i and indicated by "Modified" in the figure showed an excellent match with the reference values, which is not surprising as there are no approximations in Eq. (7) when used with Q cntr . Due to this, the scale between vertical axes of heat flux in Fig. 1(a) and thermal conductivity in Fig. 1(b) was adjusted so that the reference thermal conductivity values for butane and octane, and the "Modified" value for polystyrene would align to 300 MW/m 2 . We see an excellent match between the decomposition results of "Modified" and the reference values for both EMD and NEMD systems, which is the result of decomposition equivalence discussed by Matsubara et al. in an earlier work [55]. On the other hand, in all systems where the thermal conductivity was obtained via the group form atomic stress σ group i with unmodified LAMMPS, indicated by "LAMMPS" in Fig. 1(b), the total thermal conductivity is overestimated. This happens because the group form atomic stress σ group i approximation underestimates the total heat flux, as discussed in the previous paragraph, and thus the autocorrelation function used in the Green-Kubo relation is no longer correct. As a result, even contributions from pairwise interactions that were correctly computed in NEMD systems are not guaranteed to be correct. The difference from the reference is especially clear from the bond and angle contributions in "LAMMPS," where they are vastly overestimated for all systems. Because of this, thermal conductivity was vastly overestimated in butane, even though the underestimation in NEMD "LAMMPS" was comparatively small. On the other hand, polystyrene shows an opposite tendency: underestimation of heat flux in NEMD "LAMMPS" was comparatively large, while the overestimation of thermal conductivity was less severe than in other systems. This is not universal however, as octane has significant error in "LAMMPS" for both NEMD and EMD systems. As both butane and octane molecules had identical potential parameters and only the chain length was different, not only the type of potential, but also the structure appears to have a large effect on the heat flux and thermal conductivity obtained via the group form atomic stress. In principle, in the case of heat flux it might be possible to gauge the level of error by carefully investigating how group and centroid form atomic stress values relate, and in the case of thermal conductivity, how the thermal flux correlation changes, but as demonstrated by butane and octane systems, the relation appears to be complex and lacks apparent consistency. Therefore, as long as the group form atomic stress σ group i approximation is used, neither correct heat flux, nor thermal conductivity via the Green-Kubo relation can be obtained and gauging the level of error is also a nontrivial task.
In summary, we demonstrated that the simplistic atomic stress approximation most widely used by the popular LAMMPS package has a significant deficiency when applied to three-body or large many-body interactions, which prevents obtaining correct heat flux and thermal conductivity values that satisfy energy conservation. For simple manybody interactions such as angle, torsion, or improper potentials, we demonstrated a simple modification to the atomic stress definition, which produced an excellent match with results obtained by more rigid and computationally intensive methods. Although all of our test systems had low thermal conductivities, comparably promising results were also obtained from preliminary investigation on all-atom amorphous polymer systems with higher thermal conductivity, and our work should also be applicable to systems with high thermal conductivity. It is expected that the same framework can be used to extend the scope of application to more complex many-body potentials. Modified LAMMPS code with the new atomic stress definition will be made publicly available in a timely manner. This work was supported by JST CREST Grant No. JPMJCR17I2, Japan. Computational simulations were performed on the supercomputer system "AFI-NITY" at the Advanced Fluid Information Research Center, Institute of Fluid Science, Tohoku University.