Per-atom virial formulation for many-body potentials


Has anyone stumbled on this recent paper yet: ? The authors derive expressions of interatomic force and heat current for many-body potentials (Tersoff mainly, but also Stillinger-Weber and Brenner), and define an "unambiguous" stress applicable for periodic systems. Lattice thermal conductivity calculations using the Tersoff potential are then performed using: a) the newly derived formalism b) the per-atom virial stress as it is currently implemented in LAMMPS (see Thompson et al, J Chem Phys, 131, 154107 (2009)). Strong discrepancies are found, mainly for low-dimensional systems (graphene and CNTs). See section E, and more precisely figure 6.

tl;dr: EMD calculations with many-body potentials: be careful.


Maybe a few LAMMPS folks will have comments (CCd) …


Hi Arthur,

Thanks for bringing this to the attention of the LAMMPS community.
This is a very interesting paper containing several important results:

1. A completely general and very elegant definition of effective
pairwise forces for manybody potentials. It is so simple, I am
surprised that it has not been proposed before, but it is new to me:
Eq. (57) F_ij = -F_ji = dU_i/dr_ij - dU_j/dr_ji. This effective
pairwise force can be used to define an analogous per-atom virial and
corresponding per-atom heat flux.

2. A slightly different definition of the heat flux that includes only
the second term above (NOTE: for true two-body potentials, both terms
are identical). This version is shown to be equivalent to many
previously proposed heat flux expressions which were hitherto thought
to be distinct.

3. All three heat flux definitions (from item 1, item 2, and LAMMPS)
give the same results for Tersoff silicon and Tersoff carbon, but for
graphene and CNT, all three are different.

In my opinion, these results do not imply that any of these heat flux
definitions are somehow wrong. Rather it suggests that the Green-Kubo
method for thermal conductivity in graphene and CNTs is not reliable.
This idea could be tested further by doing careful NEMD simulations
using the same models.


Hi Aidan,

Thanks for your comments. I agree with you concerning the fact that it's certainly the EMD method itself which is not adapted to this kind of geometries. When I find some time, I'll try to write a version of compute_stress_atom based on this new effective pairwise force, and use it with compute heat/flux to compare the outcome with NEMD calculations. If I find anything interesting I'll post it here.


If you look at the code in compute stress/atom and compute heat/flux,

they are simple. They just use the per-atom tallying

of the virial that each pair style performed, in its inner

loops as it calculated pairwise or many-body interactions.

They know nothing about the pair style itself nor any

forces it computes.

I think to do something like is described in the paper,

you would have to modify how a particular pair style (e.g.

Tersoff) tallies it’s per-atom virial.