[lammps-users] Rigid molecule and NVT thermostat

Hello all,
          I do hope you can offer some insight with regards to the problems I am facing using
LAMMPS. I am attempting to model a system of dumbbells. These have a fixed bond length and
consist of 2 atoms. I am comparing the results from LAMMPS to the M.D code that I have
written up. In the atomic case (no bonds), the MSD results are in agreement within standard
deviation (100,000 steps). For the molecular case, the MSD from the M.D code is a linear plot
with time. In the case of LAMMPS, the MSD (post processing) reaches a plateau after 2000
steps i.e. the molecules freeze. For the first 500 steps the MSD from both sources are within

The work is preformed in bulk. I ensure the positions of the dumbbell are within the periodic
boundary and have a bond length of one. Positions and velocities feed into LAMMPS are
obtained from the M.D code. Using the fix rigid command I have maintain the bond length of
the dumbbells in LAMMPS to r* = 1.0. I use the NVT (Noose-Hoover) thermostat to maintain the
temperature of the system. I am aware the positions and velocities are updated with a
constant - energy time integration. If the fix NVT command appears after the fix rigid
command the thermostat does indeed maintain the specified temperature. In this scenario the
molecules after 2000 steps freeze. If the fix rigid command appears after fix NVT, the system
operates with energy time integration, reading in positions are velocities from a well
equilibrated M.D code. The molecules no longer freeze after 2000 steps, but with the lack of
a thermostat the system temperature deviates and as a result the MSD differ. The results are
consistent across a reduced density range of 0.15 to 0.3.

I have not worked with a fix langevin which would effectively maintain the temperature. This
will require a time integrator (NVE). Is this option recommended?

Another option is to use a fix shake, and apply the NVT thermostat. This is similar to my
work, where I use a rattle routine and the NVT thermostat (Noose-Hoover). In this case I can
use a harmonic bond type with a large K value (stiff spring). Do I need to specify the any
additional information (angle style) ?

Thank you for taking the time to read this. Any insight would be of tremendous value.

units lj
atom_style molecular
dimension 3d
boundary p p p

read_data data.ljbulk

mass 1 1.0
mass 2 1.0

pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0
pair_coeff 1 2 1.0 1.0
pair_coeff 2 2 1.0 1.0

group dumbbell molecule <> 1 50

neighbor 0.3 bin
neigh_modify delay 0 every 1

fix 1 dumbbell rigid molecule
fix 2 all nvt 2.0 2.0 0.01

dump 1 dumbbell xyz 100 dump.xyz

thermo 100

#thermo_style custom step temp epair ke press

run 5000
run_style verlet
timestep 0.005

You can't use fix nvt with fix rigid. They are both doing
time integration, so you are integrating the same atoms
twice (the code should have given you a warning about this).
I would thermostat with fix langevin for rigid bodies.