I am simulating homogeneous particle systems in equilibrium (NVT). There are only pair & bond interactions and time integration is performed with fix nve using a langevin thermostat. In every case, there is a small but significant offset between the normal stresses as calculated by compute pressure & compute stress/mop, with the average stress always being higher in the first case.
As a sanity check I performed a similar simulation using the “melt” case from the LAMMPS example library. This example performs just the NVE time integration without the langevin thermostat. In this case both methods gave the same average normal stress.
One explanation I came up with is that the virial term of compute pressure also accounts for the forces applied on the particle by the thermostat, even though this is not clearly stated in the documentation. Could someone more familiar with the inner workings of these computes give an opinion on that?
The commands compute pressure but also compute stress/atom have support for tallying stress contribution from fixes (if a fix modifies forces it modifies the stress) and thus it is enabled by default.
The compute stress/mop is apparently lacking that functionality (and thus the results when used with fixes that impact stress have to be questioned) and thus would be only equivalent to using compute all pressure ke pair bond angle dihedral improper.
Update:
It may be of interest to know that until rather recently, this compute only supported 2-body pair styles.
So support for bonded interactions is a rather new feature.
I was going to comment on this exact topic. I ran some more case studies and I made the following observations:
When it comes to compute pressure, there seems to be no tallying of forces from fix langevin by default. More specifically, using compute pressure without keywords gives exactly the same results (overlapping points, not just averages) as compute all pressure ke pair bond angle dihedral improper.
For systems limited just to pair interactions, compute pressure gives the same averages as compute stress/mop. The difference appears only when bond interactions are introduced in the system.
I kept track of the different contituents of stress in each case and found out that the kinetic and pair contributions give the same averages for both computes. The observed difference is indeed attributed entirely to the bond contribution.
In your opinion, does this lead to the conclusion that compute stress/mop underestimates bond contributions to stress due to an error in its implementation? I would be very interested to know since I am planning to extend my studies to cases of nonhomgeneous systems.
It will include contribution from fixes. However, whether contributions are counted depend on the settings of the individual fix. There is a fix_modify command that can add this contribution, if supported and not active by default.
I cannot comment on this because I am not familiar with the details of the theory used in compute stress/mop and am not in a position to discuss this with competence. I suggest you contact the author. There should be a mention of the name of the author for counting the contributions from bonded interactions in the source code, or you can look up the pull request on github and discuss with the author there.
compute pressure calculates the scalar global pressure and the entire pressure tensor of the whole system. compute stress/mop calculates three elements of the local stress tensor on a surface normal to x- or y- or z- cartesian axes.
The way stress/mop operates is the following: define the vector connecting particle A and particle B. if the vector crosses the plane of interest then this pair will contribute to the local stress tensor, otherwise it will be discarded. Thus, not all pair of particles will contribute to the local stress. For the case of the global pressure tensor, contributions from all pairs are considered.
In the case of the force contributed from the Langevin thermostat, I think it is not possible to define a connecting vector between a particle and the thermostat in an unambiguous way. Someone needs to define the cartesian coordinates of the thermostat.
I agree with you on the conceptual differences between compute pressure and compute stress/mop (and the observation about the force contribution from the Langevin thermostat). Nevertheless, these two different methodologies should provide the same average pressure (or normal stress on a plane) when applied on a homogeneous & equilibrated system, which is characterized by an isotropic pressure tensor.
What is really weird for me, and the reason I made this post in the first place, is the fact that that I would get statistically different pressures on such a system with each methodology. As I mentioned, all the difference could be attributed to the bond contributions.
After some more exploration I think I found the culprit: Running my simulations on a single core results in the same stresses calculated with both methods. It is only when I partition the simulation box that differences become apparent. I tried the same with your input file from the pull request. Using 2 MPI tasks actually results in no stresses being calculated at all for the case of only-bonded interactions, even though both particles are within the neighbor/communication cutoff. From what I can understand, particles located on different subdomains are not taken into account during the bond contribution calculation.
Could this be a bug or is it just happening on my machine? I am looking forward to your comments.