When I read your codes about computing rigid pressure, I found some problem . FIrst, You compute the constrained force , then compute virial term by equation : pxx = 0.5 * f.x * r.x, f.x means the constrained force in x direction, r.x means the unwrapped coordinate. But I think r.x should be the distance from partices to masscenter of the rigid. Then the computed pressure are almost zero because postive and negative are almost equal.
When I read your codes about computing rigid pressure, I found some problem
. FIrst, You compute the constrained force , then compute virial term by
equation : pxx = 0.5 * f.x * r.x, f.x means the constrained force in x
direction, r.x means the unwrapped coordinate. But I think r.x should be the
distance from partices to masscenter of the rigid. Then the computed
pressure are almost zero because postive and negative are almost equal.
are you aware of the fact that what is computed in the fix is not the
full virial, but a correction to the unconstrained virial?
Sorry, actually I don't know what you mean . When I use my own program to
what i mean is that code in fix rigid does not compute *the* pressure
but merely an additional term.
compute pressure, I can get the same results asyours if I use the unwrap
coordinate. And I don't use the fatctor 0.5.But I get nearly zero if I use
the distance from the particles to the masscener.
so what? who says that your LAMMPS input is correct? who says that
your code is correct? no code is automatically correct solely because
it produces the number you expect to see. please note that that factor
0.5 is used for a good reason that is explained in the comments.
the fact, that you get a different number in your code means nothing.
if you believe that LAMMPS is incorrect, you have to do two things.
1) you need to spend more time reading LAMMPS' source code to really
understand how computing the pressure (or rather the virial) works
2) you need to provide convincing proof where exactly LAMMPS is
computing things wrong and why and the corresponding input decks.
but please also consider that this is not a piece of code that was
written by a grad student yesterday. even though there are
occasionally bugs (like in all things produced by humans) they are
most likely to occur in new code and in features that are infrequently
used. so even though i am not an expert in the calculation of pressure
and didn't write the piece of code in question, all odds are against
you and favor the likelihood that you made a mistake somewhere or are
misunderstanding something, and thus you have to do better than "i
think you should do XXX".
The code in fix rigid is doing something conceptually simple,
although the bookkeeping is complicated. It’s the same idea
that is used in fix shake. If you apply constraint forces
to a small cluster of atoms (rigid body, shake cluster),
to keep the body rigid (or the bond lengths fixed for SHAKE), then
the contribution to the virial is simply sum (Fi dot Ri).
Fi is the constraint force on each atom and Ri is the position
of each atom. Ri must be consistent
for all the atoms in the cluster. I.e. if the cluster straddles
a periodic boundary, you must unwrap the Ri so that they are
close together for all atoms in the cluster.
I believe sum (Fi dot Ri) is translationally invariant, so it doesn’t
matter if you use Ri as displacements from the COM of the cluster
versus Ri from the origin of the simulation box.
As the comments in the code say, the 1/2 is b/c this term is
accumulated twice during the timestep.