stresses from angular (3-body) forces

Hi everybody. We are trying to calculate per-atom stresses in a model polymeric system with angle_style harmonic. I’ve noticed that within LAMMPS, this is done using the ev_tally function in angle.cpp. However, we are trying to do this as a postprocessing step.

My question is – where are the formulas used in this function obtained from? I assume they follow the procedure described in Steve & Aidan’s 2009 J. Chem. Phys. paper “General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions.” Eqs. 19-20 of that paper give the virial in terms of the position of and force acting on each atom. However, I’m confused as to how to get from there to the formulas in ev_tally,

Looking at the ev_tally function, it seems that some of what’s going on there may have to do with the periodic-boundary-condition issues discussed in the 2009 paper. Interactions with periodic images are not an issue for us as our systems are fairly large and our forces short-range. However, aren’t there still issues with the fact that pairwise force decompositions of 3-body forces are nonunique? Does LAMMPS employ Tadmor’s central force decomposition? I can’t tell whether ev_tally is just calculating the angular contribution to Eqs. 19-20 of the 2009 paper, and didn’t see any further details in angle.cpp or the online documentation…

So – where are the formulas used in the current ev_tally obtained from?

Thanks,
Rob

Hi everybody. We are trying to calculate per-atom stresses in a model polymeric system with angle_style harmonic. I’ve noticed that within LAMMPS, this is done using the ev_tally function in angle.cpp. However, we are trying to do this as a postprocessing step.

My question is – where are the formulas used in this function obtained from? I assume they follow the procedure described in Steve & Aidan’s 2009 J. Chem. Phys. paper “General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions.” Eqs. 19-20 of that paper give the virial in terms of the position of and force acting on each atom. However, I’m confused as to how to get from there to the formulas in ev_tally,

Looking at the ev_tally function, it seems that some of what’s going on there may have to do with the periodic-boundary-condition issues discussed in the 2009 paper. Interactions with periodic images are not an issue for us as our systems are fairly large and our forces short-range. However, aren’t there still issues with the fact that pairwise force decompositions of 3-body forces are nonunique? Does LAMMPS employ Tadmor’s central force decomposition? I can’t tell whether ev_tally is just calculating the angular contribution to Eqs. 19-20 of the 2009 paper, and didn’t see any further details in angle.cpp or the online documentation…

So – where are the formulas used in the current ev_tally obtained from?

rob,

if you are looking at angle.cpp, then you are looking at stress contributions from explicit angle interactions in angle styles not implicit angles from (manybody) pair styles. those have to use tally functions defined in pair.cpp. note specifically, the v_tally#() functions.

axel

Hi Axel. That’s correct, we’re looking at stress contributions from explicit angles. We’re just looking for the source of the formulas used in the current ev_tally (in angle.cpp) so we can accurately calculate per-atom stresses in postprocessing. As we both know, angular forces’ contribution to the virial are nontrivial, so we want to make sure we’re handling these the same way LAMMPS does :slight_smile:

Thanks,
Rob

Hi Axel. That’s correct, we’re looking at stress contributions from explicit angles. We’re just looking for the source of the formulas used in the current ev_tally (in angle.cpp) so we can accurately calculate per-atom stresses in postprocessing. As we both know, angular forces’ contribution to the virial are nontrivial, so we want to make sure we’re handling these the same way LAMMPS does :slight_smile:

well, the way i read the comment for Angle::ev_tally(), those stress contributions are derived from F dot r.

axel.

Axel is correct, it is just F dot r summed over the 3 atoms, than 1/3 of the sum assigned

to each atom for the per-atom stress (for that angle). Not a more complex decomposition. Aidan

is our local expert on these stress Qs, so I’m CCing him. I think he will say

there is no “one right way” to decompose this, so what we are doing is

adequately “good”.

Steve

The formula used by LAMMPS for angles is provided unambiguously in the compute stress/atom doc page. No explicit derivation is provided in the 2009 paper, but it should be clear from the comments in the Discussion at the end of the paper that the LAMMPS formula is based on the group form of the total virial (Eq. (22). For each angle term, a virial contribution of:

r_{1_a} F_{1_b} + r_{2_a} F_{2_b} + r_{3_a} F_{3_b}

is divided equally between the 3 atoms. This is not equivalent to the Tadmor’s result or any of the many other “physics-inspired” procedures.

Aidan