**I Really thanks for your helps.According to your suggestions and provide solutions, we successfully solve the rigid problem of calculation of rigid body temperature.**

**1. add the following two lines**

**>compute thermo all temp/asphere**

**>thermo_modify temp thermo**

**2.download the updated sources for compute_erotate_rigid.cpp and**

**compute_ke_rigid.cpp from the LAMMPS-ICMS repository here:**

http://git.lammps.org/git/?p=lammps-icms.git;a=tree;f=src/RIGID;hb=HEAD

**We use fix it rigid/NPT/small and rigid/npt to calculate but there is an error message (Temperature compute degrees of freedom < 0** **(…/compute_temp.cpp:101)****).As following,**

group RIGID molecule 1:3000

compute RIGID RIGID temp

velocity all create 300.0 333

neighbor 5.0 bin

timestep 0.5

thermo 1000

compute thermo all temp/asphere

thermo_modify temp thermo

fix 1 RIGID rigid/npt/small molecule temp 300.0 300.0 10.0 iso 1 1 1 **#not ok**

**If there is a wrong rigid body in calculation when the rigid/NPT?**

**Best Regards**

**Mark Li**

Dear Axel and Trung,

I Really thanks for your helps.According to your suggestions and provide

solutions, we successfully solve the rigid problem of calculation of rigid

body temperature.

1. add the following two lines

compute thermo all temp/asphere

thermo_modify temp thermo

2.download the updated sources for compute_erotate_rigid.cpp and

compute_ke_rigid.cpp from the LAMMPS-ICMS repository here:

http://git.lammps.org/git/?p=lammps-icms.git;a=tree;f=src/RIGID;hb=HEAD

We use fix it rigid/NPT/small and rigid/npt to calculate but there is an

error message (Temperature compute degrees of freedom < 0

(../compute_temp.cpp:101)).As following,

group RIGID molecule 1:3000

compute RIGID RIGID temp

velocity all create 300.0 333

neighbor 5.0 bin

timestep 0.5

thermo 1000

compute thermo all temp/asphere

thermo_modify temp thermo

fix 1 RIGID rigid/npt/small molecule temp 300.0 300.0 10.0 iso 1

1 1 #not ok

If there is a wrong rigid body in calculation when the rigid/NPT?

no. the error message comes from a "compute temp" command.

the whole purpose of step 1. was to switch away from compute temp to

compute temp/asphere, as compute temp only considers point particles,

why you have aspherical extended particles with 3 additional degrees

of freedom per particles that compute temp does not consider, but fix

rigid will (correctly) remove, if they are taken away due to forming a

rigid body.

you have the command right there.

axel.

Dear Axel

**Sorry,I didn’t describe the problem more accurately and cause you misunderstood.**

**We use fix it rigid/NPT/small and rigid/npt to calculate but there is an error message**

**(Temperature compute degrees of freedom < 0** **(…/compute_temp.cpp:101)****).**

**If there is a wrong rigid body in calculation when the rigid/NPT?**

**As follow**

units real

boundary p p p

atom_style hybrid **ellipsoid** molecular

pair_style hybrid/overlay gayberne ***

pair_coeff ***

group RIGID molecule 1:3000

compute RIGID RIGID temp

velocity all create 300.0 333

neighbor 5.0 bin

timestep 0.5

thermo 1000

**>compute thermo all temp/asphere**

**>thermo_modify temp thermo**

**>fix 1 RIGID rigid/nve/small molecule langevin 300.0 300.0 5.0 2211 #NVE is ok**

**>fix 1 RIGID rigid/npt/small molecule temp 300.0 300.0 10.0 iso 1 1 1 #NPT not ok**

**B.S.**

**Mark**

Dear Axel

* Sorry,I didn't describe the problem more accurately and cause you

misunderstood.*

what you have below is a *different* input from what you quoted previously.

don't blame me to misunderstand, when you are giving me *wrong*

information.

*We use fix it rigid/NPT/small and rigid/npt to calculate but there is an

error message *

*(Temperature compute degrees of freedom < 0** (../compute_temp.cpp:101)*

*).*

*If there is a wrong rigid body in calculation when the rigid/NPT? *

no. but you *must not* use compute temp for extended bodies.

*As follow*

units real

boundary p p p

atom_style hybrid *ellipsoid *molecular

pair_style hybrid/overlay gayberne ***

pair_coeff ***

>group RIGID molecule 1:3000

>compute RIGID RIGID temp

>velocity all create 300.0 333

>neighbor 5.0 bin

>timestep 0.5

>thermo 1000

*>compute thermo all temp/asphere *

*>thermo_modify temp thermo*

*>fix 1 RIGID rigid/nve/small molecule langevin 300.0 300.0 5.0

2211 #NVE is ok *

*>fix 1 RIGID rigid/npt/small molecule temp 300.0 300.0 10.0 iso

1 1 1 #NPT not ok*

the langevin option in rigid fixes and the temperature management in the

rigid/npt and rigid/nvt fixes work differently.

in the first case, the temperature for the rotational and translational

degrees of freedom are computed internally from the internal information

about the rigid objects and then adjusted. in the second case, an external

temperature compute is used. by default, same as for the thermo output,

this compute assumes point particles and thus you'll get the error

messages, when you apply the rigid fix to objects with additional degrees

of freedom.

so the obviously missing command here is:

fix_modify 1 temp thermo

axel.

Dear Axel,

**Thanks,I according to your suggestions for rigid/NPT and rigid/NPT/small calculation is available.**

**As follow, but results of temp ,ke, and ke_rot are NaN Value**

neighbor 2.0 bin

timestep 0.5

thermo 1000

compute thermo all temp/asphere

thermo_modify temp thermo

**fix 1 RIGID rigid/npt/small molecule temp 300.0 300.0 10.0 iso 1 1 1 #NPT ok**

**fix_modify 1 temp thermo ------->> Add in #**

compute rigid_ke RIGID ke/rigid 1

compute rigid_rot RIGID erotate/rigid 1

variable total_ke equal “c_rigid_ke+c_rigid_rot” #translational Ke +rotational Ke

variable dof equal 6*3000 # # 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.0v_total_ke/v_dof/v_boltz"

outputs:

thermo_style custom step c_thermo dt f_1 v_rigid_temp c_rigid_ke c_rigid_rot

0 284.90993 0.01 82869.875 243.63538 2720.3939 1636.5375

1000 -nan 0.43862013 -nan nan -nan -nan

2000 -nan 0.43862013 -nan nan -nan -nan

3000 -nan 0.43862013 -nan nan -nan -nan

4000 -nan 0.43862013 -nan nan -nan -nan

5000 -nan 0.43862013 -nan nan -nan -nan

Something wrong?

Regards

Mark