[lammps-users] fix_rigid usage

LAMMPS LAMMPS Mailing List Mirror
1

Hello,
I am attempting to simulate self diffusion of various gases in a flexible solid in the NVT ensemble by measuring the MSD. Eventually, I would like to simulate nitrogen (n2) as 2 atomic sites and 3 point charges with one charge on each atomic site and the third on the center of mass (represented as a massless particle with LJ eps and sigma = 0), similar to a previous thread (http://lammps.sandia.gov/threads/msg15533.html) which requires the use of fix rigid. However, using the test case of CO2, the use of “fix rigid/nvt” gives results that are inconsistent (and look incorrect) as compared to “fix nvt” . The fix nvt results are similar to what is observed for single LJ center gases and the results make sense.

The fixes used in each simulation are as follows:
fix 1 all nvt temp 293.15 293.15 100
and:
fix 1 solid nvt temp 293.15 293.15 100
fix 2 gas rigid/nvt molecule temp 293.15 293.15 100

Am I using fix rigid/nvt incorrectly? I need to use fix rigid/nvt for N2 and would like to use it for CO2 as well, to be consistent with other simulations and to reduce degrees of freedom and increase simulation speed.

Thanks,
-Greg

2

You might send an email to the fix rigid/nvt authors and ask them
about your geometry. Are you sure you are monitoring the correct
temperature for the rigid particles, to judge whether the
themodynamic state point is correct? You may want to verify
that the degrees-of-freedom contributing to the temp are correct.
LAMMPS tries to do this right, but never hurts to double check.

Also, I doubt you will get much speed-up from rigidifying a bunch
of 2-atom bodies.

Steve

3

Dear Greg,

I will elaborate on what I was trying to do. To accomplish running N2 with 2 atomic sites and 3 point charges, I commented out the three lines in atom.cpp that check for 0 mass.

I defined 2 groups, one for the fluid (which encompassed the N’s as well as the massless site) and the other which encompassed just the N’s

So if my nitrogen was type 2 and my massless site was type 3

group fluid type 2 3
group nitrogen type 2

I then gave only the nitrogen a velocity

velocity nitrogen create 300.0 234324 dist gaussian mom yes rot yes

I then “fixed” the entire fluid (both the nitrogen and the massless site)

fix 1 fluid rigid/nvt molecule temp 300.0 300.0 0.5

It seemed that the fix/rigid was defining the temperature correctly, if I used the entire fluid (nitrogen + massless site) to define it (see http://lammps.sandia.gov/doc/fix_rigid.html “If you use a temperature compute with a group that includes particles in rigid bodies, the degrees-of-freedom removed by each rigid body are accounted for in the temperature (and pressure) computation, but only if the temperature group includes all the particles in a particular rigid body.”) . Thus,

compute new fluid temp
thermo_modify temp new

I compared with DLPOLY with rigid molecules for a box of N2, and was getting reasonable agreement I thought, but I did not compare with fix nvt for flexible N2.

I was planning to do some larger systems for comparison with DLPOLY so that my temperature fluctuations weren’t so large. I was testing with only 33 N2 molecules before.

I was planning to get back to this soon, so if you find anything helpful and can share, please let me know.

Best.

Josh

4

Steve:
By adding “compute gTemp gas temp” and adding “c_gTemp” to thermo_sytle custom output, I find that the gas temperature is equal to the total system temperature, but with slightly larger fluctuations. This suggests that fix rigid/nvt is monitoring the correct temperature. Any other ways to check that the temp is being calculated correctly?

The speed up from fix rigid on CO2 and N2 will likely be nominal, but will be more consistent with my previous simulations. While CO2 can be successfully simulated without fix rigid/nvt, it provides a test case that can be extended to N2 (which contains the slightly more complicated zero mass point charge at the N2 center of mass).

Josh:
My scripts run basically the same as you have described (both for CO2 and for the commented out atom.cpp for N2). The result, however, is that when fix rigid/nvt is used for the gas (co2 or n2) and fix nvt is used for the solid matrix, the gas MSD levels off (rather than growing proportionally to time) while the MSD of the solid grows exponentially. When I use fix nvt all (gas + matrix) the matrix simply vibrates around its equilibrium position and the gas diffuses as expected. The use of fix rigid/nvt for the gas and fix nvt for the matrix, the matrix looks to be melted and flowing (reflected in the MSD). Any other test cases I can run to figure out the root of this problem?

Is there any reason that fix rigid/nvt and fix rigid might not be working well together? Is there some better way to represent the electrostatics of N2 in LAMMPS that does not require a point charge at the COM (and will thus not require the use of fix rigid/nvt)?

Thanks,
-Greg

5

greg,

My scripts run basically the same as you have described (both for CO2 and
for the commented out atom.cpp for N2). The result, however, is that when
fix rigid/nvt is used for the gas (co2 or n2) and fix nvt is used for the
solid matrix, the gas MSD levels off (rather than growing proportionally to
time) while the MSD of the solid grows exponentially. When I use fix nvt
all (gas + matrix) the matrix simply vibrates around its equilibrium
position and the gas diffuses as expected. The use of fix rigid/nvt for the
gas and fix nvt for the matrix, the matrix looks to be melted and flowing
(reflected in the MSD). Any other test cases I can run to figure out the
root of this problem?

there is one significant difference to use fix nvt for all and
fix nvt plus fix rigid/nvt: with fix nvt for all, you only monitor
and manipulate the _total_ temperature. what you should
compare on top of your tests is to use _two_ instances
of fix nvt and compare _that_ to your fix nvt plus fix rigid/nvt
setup. fix nvt does not guarantee equipartitioning of kinetic
energy between the gas and the matrix. by using two fixes
you enforce it. so you are not really comparing apples to
apples in your simulations.

cheers,
    axel.

6

Is there any reason that fix rigid/nvt and fix rigid might not be working well together?

So long as you are using them on 2 different groups of atoms, this should be
OK to do.

Steve

7

Dear all,

For the MSD issue, as he sees a leveling off of the rigid case, couldn’t this be related to the wrapping? I was having this issue before. I could never figure it out. The MSD needs to be calculated from the unwrapped coordinates of course.

The coordinates out of LAMMPS for the rigid body are wrapped, and they didn’t seem to unwrap correctly in the output, even for the COM of the rigid body.

"For example, 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. Likewise the compute msd will not compute the expected mean-squared displacement for such atoms if the body moves across periodic boundaries. It also means that if you have bonds between a pair of rigid bodies and the bond straddles a periodic boundary, you cannot use the replicate command to increase the system size. Note that this fix does define image flags for each rigid body, which are incremented when the rigid body crosses a periodic boundary in the usual way. These image flags have the same meaning as atom images (see the “dump” command) and can be accessed and output as described below. "

Before, after some discussion with Steve, it seemed that it was decided that fix rigid did store the image flags correctly for the COM (which for N2 with the massless site is exactly the massless site), but I never was able to see the flags displayed correctly ever. If the documentation is right I would think you would be able to output the COM of the rigid body as unwrapped, and I wasn’t able to do this.

http://lammps.sandia.gov/threads/msg15664.html

My solution was to use VMD and pbctools to unwrap the coordinates, rewrite the dcd file, then calculate the MSD myself. I was comparing with DLPOLY as well (which doesn’t store unwrapped coordinates either so I had to do this for this as well so it wasn’t a big deal).

Josh

8

This is a good point. The image flags of atoms in a rigid body are
changed internally so fix rigid can work the body correctly. This
"breaks" things like compute msd and replicate. The doc page
discusses this.

However, fix rigid does expose the image flage of the COM of
the rigid body in its output. Thus if you extract these quantities
(COM and the image flags) to a file you could easily compute
your own MSD of the rigid bodies themselves, which is likely
what you want. Again, see the doc page for details.

Steve

9

I agree with Steve. The image flags of the COM of the rigid body are correct.

I think what I was trying to say is that outputting (in a dump) the unwrapped coordinates of this COM didn’t seem to work, at least for me. But as he said you have the COM and the image flags, so you can unwrap them yourself in a MSD calculation.

Josh

10

Dear all,

I think this was my point (trying to remember a few months ago).

If you are just dumping the coordinates of the massless site of the 3 point N2 molecule (which is the COM of the rigid body), then the image flags always are 0 0 0.

However, Steve added an update back in October which allowed fix ave/time to output the rigid body COM (numerically equivalent to the massless site of the N2 molecule but not the same variable) and the image flags (but not via dump I believe).

This was the only place it seemed correct to me (with fix ave/time). Dumping the image flags with “dump” didn’t seem to work. Thus maybe the restart files don’t have the image flags stored either?

I don’t think the image flags which can be accessed via dump are ever updated. However, it would seem that if you have a site which is the COM of the rigid body, then this image flag should be updated. It appears not to.

If this information (image flags to reconstruct the rigid body) is not stored in the restart file, this might cause problems I would think?

Can anyone tell me if my understanding is correct?

Josh

11

If you are just dumping the coordinates of the massless site of the 3 point N2 molecule (which is the COM of the rigid body), then the >image flags always are 0 0 0.

This is not the COM and image flags I was talking about.
In your case the massless site happens to be the COM,
but that is a special case. The fix rigid itself exposes
the COM and image flags for the entire body in its
output. You cannot dump that with an atom via
the dump command. But you can access and output
it with something like the fix ave/time command.

Then you will have to post-process that file to get
the MSD.

Steve

12

Thanks Steve. This is how I understand it too.

What about the restart file as Greg was discussing? This wouldn’t have the correct image flags in it either, or does it store (and subsequently read in during the restart) the COM of the rigid body and corresponding image flags (those accessible by fix ave/time)?

If not, do you think restarting a simulation will still be ok? fix rigid will still group the bodies correctly, and any interactions between sites within rigid bodies that are not “joined” due to wrapping will be taken care of in the minimum image convention? So restarting will be ok?

Thanks again.

Josh

13

Restarting should be OK, if you re-specify the fix rigid
correctly. However, you will lose the COM/image flag
info for the rigid bodies which you might have
been using to compute the MSD (via post-processing).

So you'd need a smarter post-processor to keep track
of COM and image flags itself to compensate.

I'll think about whether it would be possible to preserve
that extra info in the restart file.

Steve

14

Axel,
Thank you for the insight as to the kinetic energy partitioning. I ran a few different test cases varying the number of fix nvt calls and velocity assignments. Whether I use 2 fix nvt (1 for gas, 1 for polymer) calls or 1 fix nvt (for polymer) and 1 fix rigid/nvt (for gas), the results are very similar. Namely, the MSD of the polymer is large; that is, comparable to that of the gas while a single fix nvt used for the whole system results in diffusing gas and the polymer fluctuating around its initial position. For the single fix, the gas was 15 degrees cooler than the polymer (from “compute group temp”) while in the 2 fix cases have both the gas and polymer with identical temperatures.

While a single fix nvt for the whole system will work in the case of CO2, CH4, Ar, etc., the 2 LJ center / 3 point charge N2 system will require the use of fix rigid. How do I maintain the correct partitioning of the kinetic energy? Is there a way I can use a single thermostat for the whole system while using 2 fixes?

As a follow up to the MSD for fix rigid: post processing via VMD seems the easiest route - I will be looking at the RMSD trajectory tool.

thanks,
-Greg

15

hi greg,

[...]

energy? Is there a way I can use a single thermostat for the whole system
while using 2 fixes?

yes. you should be able to use fix langevin together with
fix nve and fix rigid/nve, similarly fix berendsen.

since langevin is a dissipative termostat, you'll have to watch
how the kinetic energy is partitioned between the two subsystems
depending on the thermostat time constant. the shorter it is, the
more dissipation. berendsen has the exact opposite behavior,
it will emphasize differences between the two subsystems.
for a system in equilibrium, you would actually expect equipartitioning
between the two subsystems, so my gut feeling is that fix langevin
would be the preferred choice. the fact that you do see a difference
in temperature may be a finite size effect or due to lack of equilibration,
but this is speculation from my side.

As a follow up to the MSD for fix rigid: post processing via VMD seems the
easiest route - I will be looking at the RMSD trajectory tool.

RMSDtt is not what you need. RMSD monitors the root mean squared
deviation of a group of atoms from a reference configuration. that is
used to, e.g, monitor how much a protein structure fluctuates over time.
it has nothing to do with computing self-diffusion. many people confuse
MSD and RMSD. i may still have a (minimalistic) script to compute the
msd for a small group of atoms in VMD through tcl scripting somewhere,
but the clean way to do it would be to write a small compiled plugin,
since the kind of math required can be somewhat demanding and
a compiled plugin can also be parallelized (via OpenMP).

axel.

16

My suggestion for using VMD was just to rewrite the trajectory of the COM of the rigid body unwrapped using pbctools:pbcunwrap. I was just taking the wrapped coordinates of the COM that was dumped into a dcd and then unwrapping them in VMD and rewriting the dcd.

Then the MSD calculation can be done by hand on the unwrapped dcd. The MSD depends only on the displacement from the initial coordinates, it doesn’t really matter whether these initial (time =0 ) values are wrapped or not, only that the subsequent configurations are unwrapped with respect to the time=0 point.

However, now that you can access the image flags with fix ave/time and the COM coordinates you don’t really need to do this. You can just calculate the MSD from the trajectory reading in this file output from fix ave/time and moving the coordinates based on the image flags to the correct box in the MSD calculation itself.

Writing a MSD code in fortran is fairly easy, and you can find these in various places (Frenkel and Smit for example). It doesn’t really need to be parallelized either in my opinion, even if you are using multiple time origins, etc. Unless you are calculating MSDs out to microseconds with millions or billions of particles, the calculation just takes a few seconds or possibly minutes, even with multiple time origins on a single core.

Josh

17

joshua, greg,

Writing a MSD code in fortran is fairly easy, and you can find these in
various places (Frenkel and Smit for example). It doesn't really need to be
parallelized either in my opinion, even if you are using multiple time
origins, etc. Unless you are calculating MSDs out to microseconds with
millions or billions of particles, the calculation just takes a few seconds
or possibly minutes, even with multiple time origins on a single core.

there are some subtleties involved in computing a good MSD
and converging it well. if you want to extract good statistics from
your trajectory data, you need more than a few seconds worth
of compute power and you _definitely_ have to use multiple
time origins.

i have been giving a presentation on some issues with computing
self-diffusion (and other properties) from finite size systems a long
time ago. the MSD part is on pages 10 to 13 of this pdf:

http://klein-group.icms.temple.edu/akohlmey/files/talk-trieste2004-water.pdf

cheers,
   axel

18

Ok, maybe I’m underestimating time with a “few seconds” but my point was this can be done serially fine, and I still believe so. My point was that you don’t need to write a tcl script or parallel plugin in VMD to calculate a MSD.

I’ve calculated ten thousand + configurations in a MSD with a ~hundred thousand particles with multiple time steps every 50 configurations or so (calculating MSDs out to 10 ns), and I didn’t need a parallel MSD code to do this. In this case, it did take hours (or maybe longer, I don’t recall exactly) :slight_smile:

So yes, a parallel MSD implementation would be useful.

Josh

19

Thanks very much for sending along those nice slides, Axel. Even though I don’t use water, it’s instructive to see all this data side by side as a function of system size and simulation length. This motivated me to calculate the dielectric constant for my own system.

Cheers,
Lisa