Hello,
I'm trying to compute the MSD for my simulation of rigid molecules.
From the manual 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.
So does this mean that I should post process the trajectory (using
compute-fix-dump-rerun) to compute the center of mass and then dump
the COM coordinates to a new file with image flags determined from the
atom image flags? If that is true, wouldn't I need a data file to
rerun his new COM trajectory and compute the MSD.
I think an example of what to do here would really help me out.
Thanks,
Patrick
The doc page details are describing how you could
write your own script to compute an MSD from
the dump file data. Conceptually, it’s simple.
Imagine you have an atom near the COM of each
rigid body. If your script monitors the position of
that atom it will detect when it “jumps” from one
side of the box to the other, crossing the periodic
boundary. Thus you can compute its displacement
from its original position.
I don’t think there is a simple way to auto-do this
with LAMMPS due to the issues with the image
flags for rigid bodies that the doc page discusses.
Steve
Patrick,
The com of rigid molecules can be dumped if you create an output file with the nve or npt info from your simulation. Then you must create a simple routine for the msd calculation.
Otto
Steve and Otto, thank you for your response. This problem gave me some
trouble so I thought I would follow up and let others know what worked
for me and what didn't.
My trajectory output from lammps is an unwrapped dcd file. Basically,
I used VMD with the following script to write the unwrapped
coordinates to an xyz file. Then computed the com for each of my
molecules and msd from the com trajectory.
package require pbctools
pbc unwrap
set nf [molinfo top get numframes]
for {set i 0 } {$i < $nf} {incr i} {
[atomselect top all frame $i] writexyz frame_$i.xyz
}
I also tried to post-process my trajectory file and output unwrapped
coordinates xu, yu, zu in lammps using the following script. I'm
guessing that others will know why this didn't work for me but the
outputted coordinates looked reasonable by eye. It wasn't until I
calculated the msd that I knew I had a problem. I confirmed the
coordinates were wrong by visualizing them in a trajectory.
units real
boundary p p p
atom_style full
read_data my.data
fix myrigid all rigid molecule langevin 300 300 5000 1234
dump myDump all custom 1 tmp.dump type xu yu zu
rerun my.dcd dump x y z box no format molfile dcd
Patrick
Steve and Otto, thank you for your response. This problem gave me some
trouble so I thought I would follow up and let others know what worked
for me and what didn't.
My trajectory output from lammps is an unwrapped dcd file. Basically,
I used VMD with the following script to write the unwrapped
coordinates to an xyz file. Then computed the com for each of my
molecules and msd from the com trajectory.
package require pbctools
pbc unwrap
set nf [molinfo top get numframes]
for {set i 0 } {$i < $nf} {incr i} {
[atomselect top all frame $i] writexyz frame_$i.xyz
}
I also tried to post-process my trajectory file and output unwrapped
coordinates xu, yu, zu in lammps using the following script. I'm
guessing that others will know why this didn't work for me but the
outputted coordinates looked reasonable by eye. It wasn't until I
calculated the msd that I knew I had a problem. I confirmed the
coordinates were wrong by visualizing them in a trajectory.
what is the point of using fix rigid on a rerun??
axel.
I'm assumed I needed it. I commented out the line and reran the
script. It turns out I don't but I get the same results either way.
Patrick
Did you try to have the output from myrigid?
You can do this:
fix myrigid all rigid molecule langevin 300 300 5000 1234 dump myDump
fix rigidDump all ave/time 1 1 1000 f_myrigid file myrigid.dat mode vector
and have the trajectories of coms of rigid molecules as well as the velocities etc. I mean outbox coords.
Otto