Reliably obtaining thermal properties of complex systems, which often involves computing heat flux to obtain thermal conductivity via either Fourier’s law or the Green–Kubo relation, is an important task in modern molecular dynamics simulations. In our previous work [Surblys et al., Phys. Rev. E 99, 051301(R) (2019)], we have demonstrated that atomic stress could be used to efficiently compute heat flux for molecules with angle, dihedral, or improper many-body interactions, provided a newly derived “centroid” form was used. This was later successfully implemented in the LAMMPS simulation package. On the other hand, small rigid molecules, like water and partial constraints in semi-flexible molecules, are often implemented via constraint force algorithms. There has been a lack of clarification if the constraint forces that maintain geometric constraints and can also be considered as many-body forces contribute to the overall heat flux and how to compute them correctly and efficiently. To address this, we investigate how to apply the centroid atomic stress form to reliably compute heat flux for systems with constraint or rigid body dynamics. We successfully apply the centroid atomic stress form to flexible, semi-flexible, and rigid water models; decompose the computed thermal conductivity into separate components; and demonstrate that the contribution from constraint forces to the overall heat flux and thermal conductivity is small but non-negligible. We also show that while the centroid formulation produces correct heat flux values, the original “group” formulation produces incorrect and sometimes unphysical results. Finally, we provide insight into the meaning of constraint force contribution.

Thermal conductivity of a material can be comprehended as being composed of microscopic building blocks relevant to the energy transfer due to a specific microscopic process or structure. The building block is called the partial thermal conductivity (PTC). The concept of PTC is essential to evaluate the contributions of various molecular mechanisms to heat conduction and has been providing detailed knowledge of the contribution. The PTC can be evaluated by equilibrium molecular dynamics (EMD) and non-equilibrium molecular dynamics (NEMD) in different manners: the EMD evaluation utilizes the autocorrelation of spontaneous heat fluxes in an equilibrium state whereas the NEMD one is based on stationary heat fluxes in a non-equilibrium state. However, it has not been fully discussed whether the two methods give the same PTC or not. In the present study, we formulate a Green-Kubo relation, which is necessary for EMD to calculate the PTCs equivalent to those by NEMD. Unlike the existing theories, our formulation is based on the local equilibrium hypothesis to describe a clear connection between EMD and NEMD simulations. The equivalence of the two derivations of PTCs is confirmed by the numerical results for liquid methane and butane. The present establishment of the EMD–NEMD correspondence makes the MD analysis of PTCs a robust way to clarify the microscopic origins of thermal conductivity.