I’m working with rigid molecules for which I need to calculate the mean-squared displacement. As per the fix-rigid doc page, I have a 4 step procedure:
- From the all-atom dump file, I calculate the center of mass of each rigid molecule and generate a new dump file which contains the x,y,z positions of only the centers of mass (1 for each molecule), with ix,iy,iz all set to 0.
- From this new dump file, every timestep I analyze whether the center-of-mass has moved more than half the box length in any x,y,z direction. If so, I increase/decrease the ix,iy,iz flag accordingly. and make a new “wrapped” dump file with the proper ix,iy,iz.
- I use pizza.py to unwrap the previous dump file.
- I calculate the mean-squared displacement using the unwrapped coordinates.
My question concerns the first step. At first, I thought that the documentation meant that all atoms in a single rigid body would be forced to have the same ix,iy,iz coordinates (in the same box), even if they’re outside the box bounds. However, upon examining a sample dump file, I found this not to be the case. If it’s the case that the image flags are as they appear to be, it seems rather trivial to me to compute the center of mass of each molecule by “unwrapping” each atom in the molecule and doing the simple center of mass calculation. However, if it is truly this trivial, I can’t understand what LAMMPS is doing exactly such that it couldn’t calculate the center of mass of each molecule by itself, which the documentation explicitly states it cannot do.
In looking at the source code, I found the following comment in fix_rigid.cpp:
“remap xcm of each rigid body back into periodic simulation box done during pre_neighbor so will be after call to pbc() and after fix_deform::pre_exchange() may have flipped box
use domain-remap() in case xcm is far away from box due to 1st definition of rigid body or due to box flip
if don’t do this, then atoms of a body which drifts far away from a triclinic box will be remapped back into box with huge displacements when the box tilt changes via set_x()
adjust image flag of body and image flags of all atoms in body”
This statement explains to me the rationale for why LAMMPS does something different with the ix,iy,iz flags for rigid bodies (because it’s necessary for triclinic boxes), but it still doesn’t tell me WHAT it does differently.
I’ve read through the documentation and searched through the mailing list pretty thoroughly on this issue, and although I find lots of places in which the issue is discussed, I don’t see a very detailed explanation of what exactly is going on, and I’d like to be 100% sure that what I’m doing is correct.
Thanks in advance.
Graduate Research Assistant
University of California, Berkeley