fix rigid/small problem on temperature calculation

Dear Lammps users

I am testing a system using “fix rigid/small molecule” command. The system contains 125 ellipsoid particles and 125 LJ atoms. Each small rigid body is comprised of one ellipsoid and one LJ atom (therefore, 125 rigid bodies). Below is my input script; however, I encountered the error “ERROR: Temperature compute degrees of freedom < 0 (…/compute_temp.cpp:101)”

I have tried to use just two rigid dimers (that is two ellipsoid+LJ atom dimers), and used “fix rigid” instead of “fix rigid/small”. In this case everything works. Can anyone give me any suggestions??

Thanks in advance!

Chun-Wei

Dear Lammps users

I am testing a system using "fix rigid/small molecule" command. The system
contains 125 ellipsoid particles and 125 LJ atoms. Each small rigid body is
comprised of one ellipsoid and one LJ atom (therefore, 125 rigid bodies).
Below is my input script; however, I encountered the error "ERROR:
Temperature compute degrees of freedom < 0 (../compute_temp.cpp:101)"

I have tried to use just two rigid dimers (that is two ellipsoid+LJ atom
dimers), and used "fix rigid" instead of "fix rigid/small". In this case
everything works. Can anyone give me any suggestions??

have you searched the mailing list archives?
this sounds a lot like something that was discussed rather recently,
but for which i don't recall the outcome.

axel.

Thank you for your prompt reply. You are right, I did find one thread of similar problem, see the link below:

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

However, after Steve asked for a simpler version of input script, there was no further followups. This is the reason why I started this thread. Any suggestions??

Chun-Wei

Thank you for your prompt reply. You are right, I did find one thread of
similar problem, see the link below:

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

However, after Steve asked for a simpler version of input script, there was
no further followups. This is the reason why I started this thread.

you are not making much sense here. you would have helped your cause
much more by acknowledging that you are come across the same issue and
can offer an alternate example input. that would have certainly raised
awareness more and demonstrated your sincerity better.

Any suggestions??

i already *did* make a suggestion ("check the archives"). if i had a
better one, don't you think i would have mentioned it??

axel.

Yes, you are right, both I and the guy who started similar thread last year are encountering similar situation here. Of course, I also pasted an alternative input script, which I hope would be helpful. It would be great if you guys could help me clarify what might happen. Our goal is making a coarse-grained model of small molecule, that is, connecting ellipsoids with bonds. Therefore, we need the “ghost atoms” such that bonds can connect.

Thank you very much!

Chun-Wei

Hi Chun-Wei,

the error message you got results from the way compute temp, which is used by thermo output by default, calculates the number of DOF, which is correct for rigid bodies of point particles (see ComputeTemp::dof_compute() in src/compute_temp.cpp).

For your case, a workaround would be to calculate the temperature from rigid body translational and rotational energies, something like:

fix 1 RIGID rigid/small molecule langevin 300.0 300.0 1.0 2211

compute rigid_ke RIGID ke/rigid 1
compute rigid_rot RIGID erotate/rigid 1
variable total_ke equal “c_rigid_ke+c_rigid_rot”
variable dof equal 6125 # # each dimer has 6 dof
variable boltz equal 1.987e-3 # Boltzmann constant in units real: 1.987e-3 kcal/mol/K
variable rigid_temp equal "2.0
v_total_ke/v_dof/v_boltz"

thermo_style custom step pe v_rigid_temp v_total_ke

run 10000

Note that if you specify ke, etotal and press in the thermo_style command, the error message will show up again, because compute temp is used by thermo output to compute ke and press by default.

Hope that helps,
-Trung

Hi Chun-Wei,

the error message you got results from the way compute temp, which is used
by thermo output by default, calculates the number of DOF, which is correct
for rigid bodies of point particles (see ComputeTemp::dof_compute() in
src/compute_temp.cpp).

For your case, a workaround would be to calculate the temperature from rigid
body translational and rotational energies, something like:

fix 1 RIGID rigid/small molecule langevin 300.0 300.0 1.0 2211

compute rigid_ke RIGID ke/rigid 1
compute rigid_rot RIGID erotate/rigid 1
variable total_ke equal "c_rigid_ke+c_rigid_rot"
variable dof equal 6*125 # # each dimer has 6 dof
variable boltz equal 1.987e-3 # Boltzmann constant in units real:
1.987e-3 kcal/mol/K
variable rigid_temp equal "2.0*v_total_ke/v_dof/v_boltz"

thermo_style custom step pe v_rigid_temp v_total_ke

run 10000

why not just add the following two lines, which replace the point
particle temp compute with one suitable for aspherical particles?

compute thermo all temp/asphere
thermo_modify temp thermo

axel.

Dear Axel and Trung

Thank you very much for your suggestions! By using the two additional lines Axel suggested, there was no problem. Hence, in similar cases, we have to overrider the default temperature calculation scheme.

Cheers

Chun-Wei