Dear all,
I am attempting to model rigid CO2. Using fix shake with the angle
turned on of course crashes Lammps as pointed out in a recent post
(Axel and Steve seemed to give the explanation in terms of the
impossibility of constraining 180 degree angles):
http://sourceforge.net/mailarchive/message.php?msg_id=q2o33ea1b0a1004190731j671b7265p91c5d17800cfe748%40mail.gmail.com
Using only "b" or "b" and "t" constrains the bonds pretty well and
appears to give an identical answer for all the energies, cluster
sizes, etc.. The number of clusters (# of size 3 clusters) detected
seems to be the correct number, which I assume is correct as this is
the number of angles that should be constrained as well as the number
of CO2 molecules.
However, the energy for the angle is nonzero but not very large.
Checking a few angle measurements shows some slight deviations from
180 degrees (around half a degree or less).
In the manual there is a passage that says "The b keyword lists bond
types that will be constrained. The t keyword lists atom types. All
bonds connected to an atom of the specified type will be constrained.
The m keyword lists atom masses. All bonds connected to atoms of the
specified masses will be constrained (within a fudge factor of
MASSDELTA specified in fix_shake.cpp). The a keyword lists angle
types. If both bonds in the angle are constrained then the angle will
also be constrained if its type is in the list."
Does this mean that using "b" and "t" is the same thing as using "b"
and "a"? Or does just constraining the bond length constrain the
angle for CO2?
Is this the best I can do, or is it possible to restrain a 180 degree
angle for CO2 or other linear molecules?
Forgive me if I missed another message in the archive, but I didn't
see the answer for this before.
Best.
Josh
Dear all,
dear josh,
I am attempting to model rigid CO2. Using fix shake with the angle
seems to be a popular application these days.
turned on of course crashes Lammps as pointed out in a recent post
(Axel and Steve seemed to give the explanation in terms of the
impossibility of constraining 180 degree angles):
exactly.
a second option would be using fix rigid, but
there you have the problem that one of the three
moments of inertia would become infinity for a
perfectly linear molecule.
[...]
Does this mean that using "b" and "t" is the same thing as using "b"
and "a"? Or does just constraining the bond length constrain the
angle for CO2?
using "t" does the same as using "b". the angle is not constrained
in that case.
Is this the best I can do, or is it possible to restrain a 180 degree
angle for CO2 or other linear molecules?
one would have to program a special case for co2.
either by programming a variant of fix rigid that would
ignore one of the moments of inertia, or one would
program a bond constraint between the two outer atoms
then distribute forces on the carbon across all three
atoms as a center of mass force.
cheers,
axel.
Thanks Axel.
From your comments, would it be possible just to set a further
constraint using fix shake for the O-O distance (as a fake bond) in
addition to the O=C distance?
Since there can only be one fix shake this isn't possible is it? You
can't have multiple bond definitions in a single fix shake line can
you?
Josh
Hi Josh,
i've been using fix rigid for rods, which consist of several spherical beads, and had no problem with moments of inertia. In fact, one of the three principle moments of the rod is zero, and fix rigid will set the angular velocity with respect to that principle axis to be zero.
Cheers,
-Trung
Quoting Axel Kohlmeyer <[email protected]>:
Thanks Axel.
From your comments, would it be possible just to set a further
constraint using fix shake for the O-O distance (as a fake bond) in
addition to the O=C distance?
no. it will result in the same divergence than two bonds and an angle.
Since there can only be one fix shake this isn't possible is it? You
can't have multiple bond definitions in a single fix shake line can
you?
you can have multiple fix shakes and you can specify multiple
bond types, but this won't help. either your molecules will not
be properly constrained or the iterative shake procedure will
diverge at some point. there ain't no escape from the blues.
axel.
Dear Trung,
I can successfully use "fix rigid", but "fix rigid/nve" and "fix
rigid/nvt" result in "ERROR: Out of range atoms - cannot compute
PPPM".
Is there any reason that "fix rigid" would work and not the others?
Axel's comments about one of the moments of inertia being infinity and
blowing up the simulation seems valid doesn't?
I noticed that the other two were implemented by you sometime in a May
release of LAMMPS?
Thanks.
Josh
Hi Josh,
I doubt that there might be a bug with the fix rigid/nve and rigid/nvt. Could you please post a test script and data file (as simple as possible) for me to look into?
Thanks,
-Trung
Quoting "Joshua D. Moore" <[email protected]...>:
In the manual there is a passage that says "The b keyword lists bond
types that will be constrained. The t keyword lists atom types. All
bonds connected to an atom of the specified type will be constrained.
The m keyword lists atom masses. All bonds connected to atoms of the
specified masses will be constrained (within a fudge factor of
MASSDELTA specified in fix_shake.cpp). The a keyword lists angle
types. If both bonds in the angle are constrained then the angle will
also be constrained if its type is in the list."
Does this mean that using "b" and "t" is the same thing as using "b"
and "a"? Or does just constraining the bond length constrain the
angle for CO2?
using "t" does the same as using "b". the angle is not constrained
in that case.
The last line you quote from the manual means that the type in the
"a" keyword has to be used and both bonds have to be constrained
in order for the angle to also be constrained.
Steve
You can't do this. Fix shake won't let you specify a "triangle" of
atoms with all 3 bonds constrained.
Steve
Trung is correct - you can treat a linear CO2 as a rigid
body with fix rigid. It should also work with fix rigid/nve and rigid/nvt.
If it doesn't it is likely a bug.
Steve
Dear Trung and Steve,
Thank you both for your responses.
I've prepared an example of bulk CO2 here:
http://gubbins.ncsu.edu/users/joshua/LAMMPS_CO2_Question/co2_bulk_rigid_tests.tar.gz
(~0.5 MB)
I have three sets of results here, each run for 1000 steps as a test.
The output for all three are included as well as 3 resulting
configurations in the trajectory (*dcd).
1) fix rigid
This seems to run well, the energy remains constant and the
temperature fluctuates only slightly. The number of rigid bodies
detected is the number of CO2 molecules.
2) fix rigid/nve
This runs. The number of rigid bodies detected is correct. However,
the energy drifts as well as the temperature. The simulation
completes the 1000 steps. I should note that the simulation I
referred to previously (in my previous post) was for CO2 adsorbed in a
carbon structure. There it returned the same error as the fix
rigid/nvt (see below). The error was "ERROR: Out of range atoms -
cannot compute PPPM". Here the simulation for this bulk CO2
configuration ran, but the energy drifted.
3) fix rigid/nvt
This does not run, or at least does not complete 500 steps. It
returns the error "The error was "ERROR: Out of range atoms - cannot
compute PPPM". The number of rigid bodies detected is correct.
In the meantime, I've been using a model of CO2 using fixed bond
lengths (with fix shake) and a flexible angle.
Thanks in advance. It's quite possible, I just don't know how to use
fix rigid (/nve, /nvt).
Josh
I'll let Trung look at this.
Steve
dear joshua,
i played around a little bit with your inputs (they will be very useful
for the library of regression tests that we are currently building here
at temple) and i have to add one observation.
the fix rigid and fix rigid/nve inputs give pretty much the
same output, if i add the flags: "torque * off off off" to the fix.
this would be an indication that the center of mass integration
is working correctly, but the integration of the rotation in case
of rigid/nve is not correct.
on top of that, switching to "torque * on on on" is making almost
no difference which leads me to the suspicion that in the integration
of the rotation some scaling factor is missing.
i then converted the inputs to use lj units and then fix rigid/nve
became consistent with fix rigid, which confirms my suspicion.
also using fix rigid/nvt now works.
in short, currently fix rigid/nve and fix rigid/nvt currently only work
correctly with lj units.
cheers,
axel.
Thanks Axel!
Trung pointed out (off list) that my damping coefficient for "fix
rigid/nvt" was too low (0.001). I'm usually using that just to
equilibrate so I was never too concerned with it with "fix nvt". That
was the cause of "fix rigid/nvt" crashing, but energy/temperature
drift still existed for "fix rigid/nvt and fix rigid/nve.
Your findings seem to be the answer for the drift.
Thanks again.
Josh
Axel, thanks for your comments on the scaling factor with the rotation motion. I think I found and fixed the bug with updating the conjugate momenta conjqm, i.e. dtv should be replaced by (dtf * 2.0).
Fix rigid/nve and fix rigid are now consistent under units metal, and units lj, of course.
Steve, what is the best way for me to commit these changes to rigid/nve and rigid/nvt?
Cheers,
-Trung
Quoting "Joshua D. Moore" <[email protected]...>:
missing
I just spoke with Trung, and after a little bit of discussion, it seems clear why this is the source of the error.
In fix_rigid.cpp, the variables dtv, dtf, dtq represent the timestep information.
fix_rigid.cpp: lines 456-460:
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dtq = 0.5 * update->dt;
The variable force->ftm2v for LJ systems is 1, but not in any other system.
fix_rigid_nve.cpp: lines 104-107 (incorrect code):
conjqm[ibody][0] += dtv * fquat[0];
conjqm[ibody][1] += dtv * fquat[1];
conjqm[ibody][2] += dtv * fquat[2];
conjqm[ibody][3] += dtv * fquat[3];
The foregoign lines are where the conjugate quaternion momentum is updated. In LJ systems, dtf = 2.0 * dtv, but in no other systems is this true because force->ftm2v is not 1.0. So the correction is as Trung noted.
This error did not get caught in our tests because Trung and I work(ed) in LJ systems. Thanks for bringing this to our attention.
Tony
I just spoke with Trung, and after a little bit of discussion, it seems
clear why this is the source of the error.
In fix_rigid.cpp, the variables dtv, dtf, dtq represent the timestep
information.
fix_rigid.cpp: lines 456-460:
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dtq = 0.5 * update->dt;
The variable force->ftm2v for LJ systems is 1, but not in any other system.
fix_rigid_nve.cpp: lines 104-107 (incorrect code):
conjqm[ibody][0] += dtv * fquat[0];
conjqm[ibody][1] += dtv * fquat[1];
conjqm[ibody][2] += dtv * fquat[2];
conjqm[ibody][3] += dtv * fquat[3];
The foregoign lines are where the conjugate quaternion momentum is updated.
In LJ systems, dtf = 2.0 * dtv, but in no other systems is this true because
force->ftm2v is not 1.0. So the correction is as Trung noted.
actually, this only works for rigid/nve and not for rigid/nvt.
i would suspect that there is something similar happening
for coupling to the nose-hoover thermostat.
This error did not get caught in our tests because Trung and I work(ed) in
LJ systems. Thanks for bringing this to our attention.
you can get my adapted test collection from here.
http://git.icms.temple.edu/git/?p=lammps-icms.git;a=tree;f=tests/rigid-shake;h=357cbb0e5182527e7ffa2bbde52ea228406e377c;hb=c2694d0b19ca850e0508ae95e403d4ff609d2fdd
BTW: i made some tests with MPI vs. OpenMP+MPI and even when running
locally on my desktop, i see a little speedup with running 4 MPI tasks with
2 OpenMP threads each versus running 8 MPI tasks w/o OpenMP.
cheers,
axel.
Axel, you are correct. The fix rigid/nvt needs similar conversion when updating the quaternion momenta. Additionally, the translational and rotational kinetic energies should also be converted properly with regards to units:
akin_t *= force->mvv2e;
akin_r *= force->mvv2e;
before updating the Nose-Hoover thermostat chain (in update_nhcp()).
I've tested Josh's example scripts with fix rigid/nvt on my local machine and it seems to be working now. One note is that since fix rigid/nvt involves Maclaurin's expansion when updating the Nose-Hoover chain, a small absolute value of Tdamp can make the simulation crash because of a large scaling factor applied to the body velocity. For example, with Josh's inputs, if Tdamp = 0.004, the simulation will crash after about 2200 steps. I don't know how to fix this for now.
Cheers,
-Trung
Quoting Axel Kohlmeyer <[email protected]>:
Steve, what is the best way for me to commit these changes to rigid/nve and
rigid/nvt?
Just send me the updated files and I'll post a patch.
Steve
For example, with Josh's inputs, if Tdamp = 0.004, the simulation will crash after about 2200 steps. I don't know how to fix >this for now.
This seems like an unreasonably small Tdamp value for any units (lj,
real, metal)
whether you are using fix rigid or not. So I would call that user error and not
expect it to work.
Steve