I understand that it is straightforward to dump the total or net force components on each atom using LAMMPS. However, I would like to know whether it is possible to separately extract the force contributions from different interaction types - namely, bond, angle, dihedral, and pairwise interactions - on a per-atom basis.
One approach I am considering is using “compute bond/local”` and “compute pair/local” to obtain forces for individual bonds and pairs, dump them, filter atom info from bond and pair info and summing them to get net bond and pairwise force contributions on each atom. However, this method does not extend to angles and dihedrals since there are no equivalent “compute angle/local” or “dihedral/local” force options available in my knowledge.
Is there a more straightforward or built-in method in LAMMPS to extract per-atom force components by interaction type?
PS: I am transitioning a polymer bed from OPLS-AA to ReaxFF, and I am observing a significant increase in the net atomic forces on each atom. I would like to identify which interaction terms are primarily responsible for this increase.
You are trying to compare apples and oranges. OPLS-AA and ReaxFF use a very different decomposition of the force components. ReaxFF has more of them and also a lot is based on the bond order, plus you have charge equilibration. In OPLS-AA charges are fixed based on an increment system and you have fewer but fixed components. Overall the balance between the components is different, so any differences are not relevant.
ReaxFF parameterizations are tricky and you can have a hard time finding one that matches your needs. You also need to expect lower accuracy and transferability compared to well balanced and tuned molecular force fields.
My apologies, I didn’t realize the force decomposition differs. That said, there are some common components such as bonded, van der Waals, and Coulombic interactions. Would it not be possible to compare these force contributions between OPLS-AA and ReaxFF?
Regardless of the differences in decomposition, I wanted to know is there a built-in method in LAMMPS to dump per-atom force components, (e.g., bond force, angle force, dihedral force for OPLS-AA)?
Technically yes, but in my opinion it has next to no value.
You could record a trajectory and then use rerun and turn off the force components you don’t want. That is straightforward with OPLS-AA but I have no idea how to do the same with ReaxFF.
Even those aren’t as straightforward as you think:
for ReaxFF to model chemical reactivity it must (and does) continuously switch any two atoms’ interactions from bonded to non-bonded. So even the difference between “bond” and “van der Waals” is not straightforward.
ReaxFF also uses charge equilibration to dynamically generate partial charges based on electronegativity differences and particle positions. So electrostatic energy can change due to charge redistribution in ReaxFF whereas it cannot in OPLS-AA.
Furthermore, most OPLS-AA simulations will use long-ranged electrostatics whereas ReaxFF uses shielded short range electrostatics.
Really they are very different models of a chemical system.
You can use compute pe/atom with various keywords to calculate various energy contributions to various groups of atoms. Any atom experiencing a large force will usually also have a larger-than-average energy contribution from the component(s) responsible. Having said that, it is limited use performing this analysis for OPLS-AA in your context since it is very difficult to think through what you’d be comparing this to in your ReaxFF system.