Force and heat current formulas for many-body potentials in molecular dynamics simulations with applications to thermal conductivity calculations

We derive expressions of interatomic force and heat current for many-body potentials such as the Tersoff, the Brenner, and the Stillinger-Weber potential used extensively in molecular dynamics simulations of covalently bonded materials. Although these potentials have a many-body nature, a pairwise force expression that follows Newton's third law can be found without referring to any partition of the potential. Based on this force formula, a stress applicable for periodic systems can be unambiguously defined. The force formula can then be used to derive the heat current formulas using a natural potential partitioning. Our heat current formulation is found to be equivalent to most of the seemingly different heat current formulas used in the literature, but to deviate from the stress-based formula derived from two-body potential. We validate our formulation numerically on various systems described by the Tersoff potential, namely three-dimensional silicon and diamond, two-dimensional graphene, and quasi-one-dimensional carbon nanotube. The effects of cell size and production time used in the simulation are examined.

Figure 1

Per-atom virial stresses in the $x$ direction on individual carbon atoms in a configuration generated by randomly shifting the positions of all the atoms from the perfect graphene structure by a small amount. Circles and crosses represent the results obtained by the LAMMPS code and our GPU code [using the formula Eq. (39)], respectively. The solid and dashed lines represent the mean values of the circles and crosses, respectively.

Figure 2

(a)–(e) Running thermal conductivities as a function of correlation time for silicon with different simulation cell sizes at 500 K. The thinner (and lighter) and the thicker (and darker) lines represent the results of independent simulations with different initial velocities and the ensemble average over the independent simulations, respectively. (f) Thermal conductivity as a function of the simulation cell size $N$. Markers with error bars represent the average values and the corresponding standard errors for a given $N$. The solid line indicates the average (147 W/mK) over the 5 simulation cell sizes and the dashed lines indicate the corresponding standard error ($\pm 2$ W/mK).

Figure 3

Same as Fig. 2, but for diamond at 300 K.

Figure 4

Same as Fig. 2, but for graphene at 300 K.

Figure 5

Same as Fig. 2, but for (10, 0)-CNT at 300 K.

Figure 6

Running thermal conductivities $\kappa \left(t\right)$ as a function of correlation time for (a) silicon at 500 K, (b) diamond at 300 K, (c) graphene at 300 K, and (d) (10, 0)-CNT at 300 K obtained by using the Hardy formula (solid lines), the stress formula (dashed lines), and LAMMPS (dot-dashed lines). For each material, the line and the shaded area represent the averaged $\kappa \left(t\right)$ and the standard error calculated from an ensemble of 10 independent simulations.