Image Flags for Rigid Molecules

Hello,

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:

  1. 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.
  2. 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.
  3. I use pizza.py to unwrap the previous dump file.
  4. 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.

Efrem Braun

Graduate Research Assistant
University of California, Berkeley

I don't know if this helps.
I will say that how the atoms images are represented internally may
not be relevant if you are just trying to calculate displacement.

There's probably a way to do what you want using pizza.py, but anyway
if you just want the 3-dimensional coordinates of all of the atoms in
your dump file (unwrapped and sorted by ID#) you can try this:

dump2data.py -atomstyle full -raw -multi < dump_file > traj_coords.raw

This will add a blank line between frames/snapshots in the trajectory,
so file which is created will contain T*(N+1) lines, assuming the
original trajectory contains N atoms, and T snapshots.

I haven't taken the time to learn that tool yet, but I should The
"dump2data.py" script is located in moltemplate's "src" directory.
(So you probably already have it installed.) It's slow and klunky but
it works. (I attached a brief documentation file.)

Again, I don't know if this helps at all.
Cheers!

Andrew

docs_dump2data.txt (5.27 KB)

Hello,

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:
1. 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.
2. 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.
3. I use pizza.py to unwrap the previous dump file.
4. I calculate the mean-squared displacement using the unwrapped
coordinates.

​why not just use the output vector of the individual rigid bodies that is
provided through the fix directly?
from the fix rigid docs:

All of the *rigid* fixes except *rigid/small* compute a global array of
values which can be accessed by various output
commands<http://lammps.sandia.gov/doc/Section_howto.html#howto_15>.
The number of rows in the array is equal to the number of rigid bodies. The
number of columns is 15. Thus for each rigid body, 15 values are stored:
the xyz coords of the center of mass (COM), the xyz components of the COM
velocity, the xyz components of the force acting on the COM, the xyz
components of the torque acting on the COM, and the xyz image flags of the
COM, which have the same meaning as image flags for atom positions (see the
"dump" command). The force and torque values in the array are not affected
by the*force* and *torque* keywords in the fix rigid command; they reflect
values before any changes are made by those keywords.
The ordering of the rigid bodies (by row in the array) is as follows. For
the *single* keyword there is just one rigid body. For the *molecule* keyword,
the bodies are ordered by ascending molecule ID. For the *group* keyword,
the list of group IDs determines the ordering of bodies​

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.

Axel,

That’s a wonderful idea. With that in mind, might I be bold enough to suggest an update to the documentation on the “fix rigid” page? Currently, it states:
Here are details on how, you can post-process a dump file to calculate a diffusion coefficient for rigid bodies, -------using the altered per-atom image flags written to a dump file. The image flags for atoms in the same rigid body can be used to unwrap the body and calculate its center-of-mass (COM)--------. As mentioned above, this COM will always be inside the simulation box. Thus it will “jump” from one side of the box to the other when the COM crosses a periodic boundary. If you keep track of the jumps, you can effectively “unwrap” the COM and use that value to track the displacement of each rigid body, and thus the mean-squared displacement (MSD) of an ensemble of bodies, and thus a diffusion coefficient.

I think that the part I’ve put between the dashes should be changed to reflect your excellent suggestion, as the recommended procedure does not appear to be necessary.

Well, another thing I don’t quite get…Axel, when I use your suggestion, not all of the centers of mass have ix,iy,iz all equal to 0. Per the documentation:

IMPORTANT NOTE: The periodic image flags of atoms in rigid bodies are altered so that the rigid body can be reconstructed correctly when it straddles periodic boundaries. The atom image flags are not incremented/decremented as they would be for non-rigid atoms as the rigid body crosses periodic boundaries. Specifically, they are set so that the center-of-mass (COM) of the rigid body always remains inside the simulation box.

This means that if you output per-atom image flags you cannot interpret them as you normally would. I.e. the image flag values written to a dump file will be different than they would be if the atoms were not in a rigid body.

Well, another thing I don't quite get...Axel, when I use your suggestion,
not all of the centers of mass have ix,iy,iz all equal to 0. Per the
documentation:

IMPORTANT NOTE: The periodic image flags of atoms in rigid bodies are
altered so that the rigid body can be reconstructed correctly when it
straddles periodic boundaries. The atom image flags are not
incremented/decremented as they would be for non-rigid atoms as the rigid
body crosses periodic boundaries. Specifically, they are set so that the
center-of-mass (COM) of the rigid body always remains inside the simulation
box.

​this is a case of apples and oranges. the image flags for the rigid body
COMs are a completely different thing than the per atom image flags. this
section refers to the per atom image flags and nothing else.​

[...]

Does this not mean that the ix,iy,iz of the centers of mass should always
be 0? Otherwise, the following comment in the documentation is not
necessary:

Here are details on how, you can post-process a dump file to calculate a
diffusion coefficient for rigid bodies, using the altered per-atom image
flags written to a dump file. The image flags for atoms in the same rigid
body can be used to unwrap the body and calculate its center-of-mass (COM).
As mentioned above, this COM will always be inside the simulation box. Thus
it will "jump" from one side of the box to the other when the COM crosses a
periodic boundary. If you keep track of the jumps, you can effectively
"unwrap" the COM and use that value to track the displacement of each rigid
body, and thus the mean-squared displacement (MSD) of an ensemble of
bodies, and thus a diffusion coefficient.

​again, this refers to the per atom image flags computing the COM after the
fact. it has nothing to do with the per body information that fix rigid
maintains internally.​

If the ix,iy,iz flags correctly follow the center of mass as it moves to
the simulation box next door, there should be no need to "keep track of the
jumps;" one need simply unwrap the coordinates of the centers of mass as
per usual.

Sorry for the excessive questions recently...I just can't wrap my head
around what's going on with the "fix rigid" command in terms of the center
of mass.

​you should step back and do something else for a bit and then look again.
it seems your thinking is stuck in a corner because you are making
incorrect assumptions and then are misinterpreting the documentation
because of that. taking a break and then looking again with a clear mind
can help getting out of this.

axel.

Keep in mind if you are using "fix rigid" with the "infile" keyword,
it looks like this will wrap you molecular coordinates, at least by
default.

"The masstotal and center-of-mass coordinates (xcm,ycm,zcm) are
self-explanatory. The center-of-mass should be consistent with what is
calculated for the position of the rigid body with all its atoms
unwrapped by their respective image flags. If this produces a
center-of-mass that is outside the simulation box, LAMMPS wraps it
back into the box."
http://lammps.sandia.gov/doc/fix_rigid.html

Presumably, if you are trying to track the molecules over long times
as they move long distances, then you don't want this to happen.

In any case, the dump file should contain all the information you need
to obtain the unwrapped atom coordinates, regardless of what fix rigid
is doing.

So instead of writing the center-of-mass coordinates using "infile",
just get the unwrapped coordinates straight from the dump file. One
other way to do this is to tell LAMMPS not to include image flags in
the dump file:

dump 1 all custom 10000 traj.lammpstrj id mol type x y z

(instead of
dump 1 all custom 10000 traj.lammpstrj id mol type x y z ix iy iz)

Either way, you can still use dump2data.py, or pizza.py to extract
unwrapped coordinates from an existing dump file.

  my earlier email said:
I haven't taken the time to learn that tool yet, but I should The

(Sorry for the confusion, I was referring to pizza.py)

I hope this helps.
Andrew

P.D.
If you already have a dump file, you can use "rerun" to read a dump
file and "dump" to write it out again in a new format without the
image flags.
http://lammps.sandia.gov/doc/rerun.html

In any case, the dump file should contain all the information you need
to obtain the unwrapped atom coordinates, regardless of what fix rigid
is doing.

Oops. Not necessarily. That depends on the dump command you used...

So instead of writing the center-of-mass coordinates using "infile",
just get the unwrapped coordinates straight from the dump file. One
other way to do this is to tell LAMMPS not to include image flags in
the dump file:

dump 1 all custom 10000 traj.lammpstrj id mol type x y z

Correction: You probably want to use "xu yu zu" instead of "x y z".
So try this instead:
dump 1 all custom 10000 traj.lammpstrj id mol type xu yu zu

My apologies. It's been a while since I looked at the docs.
http://lammps.sandia.gov/doc/dump.html

Andrew

(I normally use:
dump 1 all custom 10000 traj.lammpstrj id mol type x y z ix iy iz)