Hi Axel and Steve,
Sorry to bother you again.
Following your suggestion, I update the LAMMPS to the last vesion (Stable version (9 Dec 2014) downloaded last night).
I changed the in file following the RIGID example in LAMMPS package. The temperature is zero too.
The warning information was seen when performing the Velocity Command (WARNING: Cannot count rigid body degrees-of-freedom before bodies are initialized (…/fix_rigid.cpp:1123))
I attach the in and data file I used again. Pleas help me check what’s wrong (bug or in file).
Thanks in advance.
The following is the log file.
test.data (9.28 KB)
in.test (876 Bytes)
Hi Axel and Steve
I have read the fix rigid.cpp many times, I found the error maybe due to the dipole moment in my data file.
Then, i tested my data file after setting all dipole moment zero. For this case. the temperature is, indeed, not zero, the program performs normally.
So, I think there may exist a bug in the void FixRigid::setup_bodies_static() for the dipole style.
By the way, I may misunderstand the code.
Anyway, I hope the problem will be solved
Wade
------------------ Original ------------------
Hi Axel and Steve
I have read the fix rigid.cpp many times, I found the error maybe due to the
dipole moment in my data file.
Then, i tested my data file after setting all dipole moment zero. For this
case. the temperature is, indeed, not zero, the program performs normally.
So, I think there may exist a bug in the void
FixRigid::setup_bodies_static() for the dipole style.
By the way, I may misunderstand the code.
Anyway, I hope the problem will be solved
as steve already indicated, you can monitor the temperature of the
rigid objects through printing the scalar property of the fix. this is
the temperature as computed from the individual rigid bodies as
opposed to the regular temperature compute that uses individual atoms
and requires the fixes to properly subtract out the degrees of freedom
removed by making objects rigid.
your test case is a bit extreme: two dimensions, only two rigid
objects, the objects don't interact.
please also note that time integration is not affected by this, only
operations that depend on a (atomic) temperature computation. so this
is mostly a diagnostic issue.
axel.
Hi Axel,
Thanks for your reply.
I now know if the rigid includes the extended style(dipole, sphere, ellipsoid), the temperature will be zero and ke in the thermo_style also is zero. Then, I turn to use the scalar property of the fix rigid to monitor the temperature.
the command is like:
fix 1 all rigid/nve molecule langevin 1.2 1.2 0.01 1234
thermo_style custom step f_1
However, the temperature of f_1 is about 1800, not around the 1.2 during the whole simulation. This seems ridiculous.
I do not know what’s happening. Is there any relation between scalar property and the temperature given in langevin dynamics or my misunderstanding? The system seems normal without the numerical unstability.
Thanks again.
after tracking down the other issue, i've now also seen, why the
regular temperature compute is off in your case.
it has to do with the fact that point dipole particles are treated as
particles with additional degrees of freedom.
i am wondering whether this is for historical reasons, since atom
style dipole is now a true point dipole and to model extended dipole
systems require to use it as a hybrid with atom style sphere. to get
the regular temperature output corrected you have to adjust the number
of degrees of freedom, e.g. though:
compute_modify thermo_temp extra $(-3*atoms)
which will re-add the degrees of freedom that were removed in excess.
this should provide the same temperature as the fix scalar property.
since in periodic boundaries the system is invariant to translations,
these number of extra degrees of freedom is set to 2 or 3 (for 2d or
3d systems, respectively) the fully consistent modification would be:
compute_modify thermo_temp extra $(-3*atoms+2)
or:
compute_modify thermo_temp extra $(-3*atoms+3)
axel.
after tracking down the other issue, i've now also seen, why the
regular temperature compute is off in your case.
it has to do with the fact that point dipole particles are treated as
particles with additional degrees of freedom.
i am wondering whether this is for historical reasons, since atom
style dipole is now a true point dipole and to model extended dipole
systems require to use it as a hybrid with atom style sphere. to get
when double checking the code for how extended particles are treated,
it appears that indeed (point) dipoles to not contribute a torque and
thus should not be counted as extended particles. however, in addition
there is a bug in the code that counts point particles as extended
particles, too, as soon as there is one extended particle in the
system. the following one line patch to fix_rigid.cpp corrects this
(and treats point dipoles as points, too).
[[email protected]... src]$ git diff lammps-ro/master RIGID/fix_rigid.cpp
diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp
index 75080ec..4161c30 100644
--- a/src/RIGID/fix_rigid.cpp
+++ b/src/RIGID/fix_rigid.cpp
@@ -1141,7 +1141,8 @@ int FixRigid::dof(int tgroup)
for (int i = 0; i < nlocal; i++)
if (body[i] >= 0 && mask[i] & tgroupbit) {
- if (extended && eflags[i]) mcount[body[i]]++;
+ // do not count point particles and point dipoles as extended particles
+ if (extended && (eflags[i] & ~(POINT|DIPOLE))) mcount[body[i]]++;
else ncount[body[i]]++;
}
with this and the other change in place, everything should now work as
advertised.
axel.