virial in systems with rigid body

Hi all,

I’m trying to understand how the virial is calculated when you have a rigid body in your system (with force and torques being off) versus a system that you just don’t include that group in time integration. I have made a simple input (attached) to grasp this calculation but I see a couple of inconsistencies that I’d appreciate it if you can help me to understand those:

  1. The reported virial by “compute pressure” seems to SUM(F.r)/V instead of SUM(F.r)/dV, where d is the dimensionality of the problem.

  2. I can get the reported virial by “compute pressure” for z direction but I can’t verify the values that are reported for x virial neither when the forces and torques are off nor when they’re on. In fact, it’s reported to be on the order of 10^-7 or -6 but what I calculate based on the forces on the constituents of the rigid body (printed in dump file) it should be something around 0.002 (if not include the dimensionality factor d for virial).

  3. For the case when rigid body is not used and I simply don’t time integrate those particle I can verify the reported virial and it’s fine.

I have attached the input file for a simulation with 3 particles which I make two of them rigid to test the aforementioned case. The dump file is attached too.

Cheers,
Kasra.

in.virial (1.09 KB)

dump.dat (2.51 KB)

Comments below.

Steve

Thank you Steve. I finally figured it out. For the first part I was confusing the total pressure virial with pressure components where the latter doesn’t divide by ‘d’.

As I understand when I use “fix rigid” the internal forces in the rigid body don’t contribute to the virial of the system and just the constraint force on each constituent is counted in the virial. So if I don’t use ‘neigh_modify exclude’ and ‘delete_bonds’ for the rigid group, I have to get the same virial value for the system as when I use them, right? and ‘neigh_modify’ and ‘delete_bonds’ are used in rigid case just for efficiency as the document says.
For the simple case that I made I can see that ‘neigh_modify exclude’ on rigid body doesn’t affect the virial whether it’s used or not, but for the complicated system that I have water with a polymer wall and I make the centre slab of the wall rigid, the virial is not the same as when I use ‘neigh_modify exclude’ on the rigid part versus when I don’t use it. so the problem should be in my setup right? as it essentially should give the same value?!!!

Cheers,

Kasra.

For the simple case that I made I can see that 'neigh_modify exclude' on
rigid body doesn't affect the virial whether it's used or not, but for the
complicated system that I have water with a polymer wall and I make the
centre slab of the wall rigid, the virial is not the same as when I use
'neigh_modify exclude' on the rigid part versus when I don't use it.

Yes (Agree again,)
In my experience, I had the best results using "neigh_modify" exclude.
I don't remember the details regarding what happened when I left
"neigh_modify" out. (I'm pretty sure that the system exploded because
in my case I had large repulsive forces between the immobile atoms.)
Anyway, it does not shock me that the virial might be different.
"neigh_modify exclude" will reduce the number of pairwise forces
calculated, and these contribute to the virial.

To calculate the virial correctly, we need to somehow exclude the
forces exerted on the immobile atoms from the virial. (..but not the
other way around. Forces that immobile atoms exert on mobile atoms
should be included in the virial.) If this is still an interesting
topic to anybody, here is my long rant on this topic from a couple
years ago:

http://lammps.sandia.gov/threads/pdfzpIPCFTvWW.pdf
http://lammps.sandia.gov/threads/msg23341.html

Perhaps this is what fix_rigid is already doing (when used with
neigh_modify exclude).

Reading the source code is the best way to figure out what is going on
with the virial. (The code confused me two years ago, but I should
take another look at this.) If I get a chance to look at the
fix_rigid code again, I'll repost with more (hopefully more useful)
comments. Cheers!

Andrew

Thank you Steve. I finally figured it out. For the first part I was
confusing the total pressure virial with pressure components where the
latter doesn't divide by 'd'.

As I understand when I use "fix rigid" the internal forces in the rigid body
don't contribute to the virial of the system and just the constraint force

One minor comment: I'm under the impression that unless you have bonds
connecting the mobile and immobile atoms, you don't need to worry
about delete_bonds. (So, yes, I agree with you.)

Hi,

I studied the fix_rigid.cpp because I’m dealing with a single rigid body so I’m just using “fix rigid” and not the “fix rigid/npt” which I assume is for a system with multiple rigid bodies and enforces NPT ensemble on those rigid bodies, right?!

In fix_rigid.cpp, virial is calculated as the constraint force on the constituents of rigid body:
i.e. to my understanding: mass of the constituent times the change in linear velocity of that constituent in time which is m[i]*( v[i] - omega x r_cm [i] ) minus of all the forces that are exerted on that constituent from other atoms except the atoms in the rigid body, basically external forces, here is the part of the code that calculates the virial:

// virial = unwrapped coords dotted into body constraint force
// body constraint force = implied force due to v change minus f external
// assume f does not include forces internal to body
// 1/2 factor b/c final_integrate contributes other half
// assume per-atom contribution is due to constraint force on that atom

if (evflag) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];

vr[0] = 0.5x0fc0;
vr[1] = 0.5x1fc1;
vr[2] = 0.5x2fc2;
vr[3] = 0.5x0fc1;
vr[4] = 0.5x0fc2;
vr[5] = 0.5x1fc2;

v_tally(1,&i,1.0,vr);

So this is what I understand, and I think that doing “neigh_modify exclude” for rigid body shouldn’t affect the total virial of the system and it’s only for efficiency but I’m not positive. The same thing for delete_bonds. Is that right?

Cheers,
Kasra.