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.

The standard equilibrium Green-Kubo and nonequilibrium molecular dynamics (MD) methods for computing thermal transport coefficients in solids typically require relatively long simulation times and large system sizes. To this end, we revisit here the homogeneous nonequilibrium MD method by Evans [Phys. Lett. A 91, 457 (1982)] and generalize it to many-body potentials that are required for more realistic materials modeling. We also propose a method for obtaining spectral conductivity and phonon mean free path from the simulation data. This spectral decomposition method does not require lattice dynamics calculations and can find important applications in spatially complex structures. We benchmark the method by calculating thermal conductivities of three-dimensional silicon, two-dimensional graphene, and a quasi-one-dimensional carbon nanotube and show that the method is about one to two orders of magnitude more efficient than the Green-Kubo method. We apply the spectral decomposition method to examine the long-standing dispute over thermal conductivity convergence vs. divergence in carbon nanotubes.

Two-dimensional materials have unusual phonon spectra due to the presence of flexural (out-of-plane) modes. Although molecular dynamics simulations have been extensively used to study heat transport in such materials, conventional formalisms treat the phonon dynamics isotropically. Here, we decompose the microscopic heat current in atomistic simulations into in-plane and out-of-plane components, corresponding to in-plane and outof-plane phonon dynamics, respectively. This decomposition allows for direct computation of the corresponding thermal conductivity components in two-dimensional materials. We apply this decomposition to study heat transport in suspended graphene, using both equilibrium and nonequilibrium molecular dynamics simulations. We show that the flexural component is responsible for about two-thirds of the total thermal conductivity in unstrained graphene, and the acoustic flexural component is responsible for the logarithmic divergence of the conductivity when a sufficiently large tensile strain is applied.