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.

We develop a neuroevolution-potential (NEP) framework for generating neural network-based machinelearning potentials. They are trained using an evolutionary strategy for performing large-scale molecular dynamics (MD) simulations. A descriptor of the atomic environment is constructed based on Chebyshev and Legendre polynomials. The method is implemented in graphic processing units within the open-source GPUMD package, which can attain a computational speed over 10 7 atom-step per second using one Nvidia Tesla V100. Furthermore, per-atom heat current is available in NEP, which paves the way for efficient and accurate MD simulations of heat transport in materials with strong phonon anharmonicity or spatial disorder, which usually cannot be accurately treated either with traditional empirical potentials or with perturbative methods.

We extend the phase field crystal (PFC) framework to quantitative modeling of polycrystalline graphene. PFC modeling is a powerful multiscale method for finding the ground state configurations of large realistic samples that can be further used to study their mechanical, thermal or electronic properties. By fitting to quantum-mechanical density functional theory (DFT) calculations, we show that the PFC approach is able to predict realistic formation energies and defect structures of grain boundaries. We provide an in-depth comparison of the formation energies between PFC, DFT and molecular dynamics (MD) calculations. The DFT and MD calculations are initialized using atomic configurations extracted from PFC ground states. Finally, we use the PFC approach to explicitly construct large realistic polycrystalline samples and characterize their properties using MD relaxation to demonstrate their quality.

Graphics processing units have been extensively used to accelerate classical molecular dynamics simulations. However, there is much less progress on the acceleration of force evaluations for many-body potentials compared to pairwise ones. In the conventional force evaluation algorithm for many-body potentials, the force, virial stress, and heat current for a given atom are accumulated within different loops, which could result in write conflict between different threads in a CUDA kernel. In this work, we provide a new force evaluation algorithm, which is based on an explicit pairwise force expression for many-body potentials derived recently [Phys. Rev. B 92 (2015) 094301]. In our algorithm, the force, virial stress, and heat current for a given atom can be accumulated within a single thread and is free of write conflicts. We discuss the formulations and algorithms and evaluate their performance. A new open-source code, GPUMD, is developed based on the proposed formulations. For the Tersoff many-body potential, the double precision performance of GPUMD using a Tesla K40 card is equivalent to that of the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) molecular dynamics code running with about 100 CPU cores (Intel Xeon CPU X5670 @ 2.93 GHz).

Nonequilibrium molecular dynamics (NEMD) has been extensively used to study thermal transport at various length scales in many materials. In this method, two local thermostats at different temperatures are used to generate a nonequilibrium steady state with a constant heat flux. Conventionally, the thermal conductivity of a finite system is calculated as the ratio between the heat flux and the temperature gradient extracted from the linear part of the temperature profile away from the local thermostats. Here we show that, with a proper choice of the thermostat, the nonlinear part of the temperature profile should actually not be excluded in thermal transport calculations. We compare NEMD results against those from the atomistic Green's function method in the ballistic regime, and those from the homogeneous nonequilibrium molecular dynamics method in the ballistic-to-diffusive regime. These comparisons suggest that in all the transport regimes, one should directly calculate the thermal conductance from the temperature difference between the heat source and sink and, if needed, convert it to the thermal conductivity by multiplying it with the system length. Furthermore, we find that the Langevin thermostat outperforms the Nosé-Hoover (chain) thermostat in NEMD simulations because of its stochastic and local nature. We show that this is particularly important for studying asymmetric carbon-based nanostructures, for which the Nosé-Hoover thermostat can produce artifacts leading to unphysical thermal rectification. Our findings are important to obtain correct results from molecular dynamics simulations of nanoscale heat transport as the accuracy of the interatomic potentials is rapidly improving. :1905.11024v1 [cond-mat.mes-hall] arXiv

Carbon-based nanomaterials have drawn strong interest for potential applications due to their extraordinary stability and unique mechanical, electrical and thermal properties. For the minimization of microelectronics/micromechanics circuits, bridging the low dimensional microscopic structure and mesoscopic modeling is indispensable. Graphene and carbon nanotubes are suggested as ideal 'building blocks' for the bottom-up strategy, and recently the integration of both materials has stimulated research interests. In this work we investigated the thermal and mechanical performance in the pillaredgraphene -constructed by combining graphene sheets and carbon nanotubes to create a threedimensional nano network. Reverse non-equilibrium molecular dynamics simulations were carried out to analyze the thermal transport behavior. The obtained thermal conductivities are found to be possibly isotropic in two specific directions or highly anisotropic for certain structure configurations. In the mechanical performance analysis, tensile deformations are loaded along graphene plane and along tube axis. The elongation responses and stress-strain relations are observed to be nearly linear, and the calculated strength, fracture strain and Young's moduli are lower than the pristine graphene or carbon nanotubes. The alterations in the thermal and mechanical performances are ascribed to the bond conversion on the junctions.