I have a question regarding the calculation of pressure for systems composed of rigid molecules in LAMMPS.
I am simulating a CO₂ system with fix rigid/nvt/small, and I would like to compute the system pressure. However, I noticed that using special_bonds lj/coul 0 0 0 and special_bonds lj/coul 1 1 1 gives nearly idential pressure values.
This is surprising to me. According to the documentation, the pressure is computed from the virial, and I would expect the two settings to produce different pressures since the latter should include intramolecular virial contributions while the former excludes them.
if you want to turn off intramolecular forces in each rigid body, you would need to specify neigh_modify exclude molecular/intra rigid_grp (see neigh_modify command — LAMMPS documentation and examples/rigid/in.rigid.tnr`).
Your statement that “…the latter should include intramolecular virial contributions while the former excludes them.” is not what the special_bond command does: it specifies the weighting coefficients for pairwise energy and force contributions between pairs of atoms that are also permanently bonded to each other.
If you don’t have any bond between the constituent particles (atoms) in each rigid body, or don’t use any bond style for these intra-body bonds, and/or use exclude molecular/intra, then special_bonds coeffs wouldn’t change the virial contribution coming from pair and bond styles.
Thank you very much for your detailed explanation.
To clarify our system: we explicitly define all intramolecular bonds within each CO₂ molecule. Since the system includes Coulombic interactions and we use PPPM for long-range electrostatics, we have not applied neigh_modify exclude molecular/intra rigid_grp; because of concerns regarding its effect on the k-space contribution. Is this interpretation correct?
In fact, we tested your suggestion, and while the instantaneous pressure values differ, the time-averaged pressures remain nearly identical.
Regarding your statement: “the special_bond command does: it specifies the weighting coefficients for pairwise energy and force contributions between pairs of atoms that are also permanently bonded to each other”; I am wondering how it works when employing PPPM method or other long rang method.
In fact, what I am interested in is how the virial contribution is computed when using fix rigid? Beyond the usual contributions (pair, bond, real-space, k-space, etc.), does LAMMPS include additional virial terms associated with the rigid-body constraints? If it is so, which algorithm is employed in Lammps?
Depending on the physical regime of your simulation, this is not a particularly surprising result. Consider compressing a sample of CO2. While it remains a gas, the increased pressure will come primarily from kinetic collisions with the container wall; should it enter a supercritical or solid state, its volume will decrease dramatically as molecules come dramatically closer.
In either case, the pressure or volume changes will not come primarily from intra-molecular changes in the C=O distance. So it makes sense that the intramolecular pressure contribution is not particularly prominent in a realistic force field.
By contrast, suppose you were compressing a diamond in a high-pressure anvil. In that case any change in volume must be reflected in the average C-C bond length, and any realistic MD simulation would contain significant bond-by-bond contributions to the pressure.
The constrained forces are dependent on how rigid bodies are modeled in your system, and may cause artifacts in certain cases, such as by overlapping LJ spheres.
Because the constrained forces are fictitious as far as rigid body dynamics is applied, you may turn off the virial contribution of fix rigid to the system virial via
fix_modify fixrigid_id virial no
Regarding the neigh_modify exclude molecule/intra warning in the presence of PPPM, it is because the real-space Coulombic interaction between charged constituent particles on the same body is excluded while the k-space energy still includes the contribution of these charges to the charge density. The calculated elong value is therefore inconsistent with the real-space term ecoul. The body force is still correct because the intra-body forces sum up to zero.
You can try plugging factor_coul either 1.0 (full pair inclusion) or 0.0 (zero pair inclusion) into these equations to see how force_coul and ecoul are adjusted.