[lammps-users] Bug (?): LAMMPS gives 10x too much kinetic energy for "velocity set" atoms

I have discovered unexpected, unpredictable behavior with LAMMPS's
"velocity set" command. In the example below, I've generated two helium
atoms (mass 4.002 Da) moving in opposite directions parallel to the
z-axis. The set velocity, in "metal" units, is 219.57 Angstrom/ps (and
-219.57 A/ps for the other atom), which should yield a kinetic energy of
10.0 eV per atom (20 eV total).

However, LAMMPS reports throughout the simulation that each helium atom
has a speed of 695.1 Angstroms/ps, with total kinetic energy 200 eV. I
saw this behavior in another simulation, but this very simple system
reproduces the problem nicely. I have included a minimal working example
at the bottom of this message. We are using LAMMPS release 3-18-2010.

Any suggestions? Is this a bug, or am I doing something wrong?

Karl D. Hammond
[email protected]...
[email protected]...

Working example (exhibits behavior with LAMMPS version 3-18-2010 running
   on one processor):

# LAMMPS file to test velocity units
units metal
dimension 3
boundary p p p
atom_style atomic
lattice bcc 3.1656
atom_modify map array
region simbox block 0 10 0 10 0 10
create_box 1 simbox
mass 1 4.0026032497
pair_style lj/cut 2.5
pair_coeff 1 1 5.9225e-4 1.333

create_atoms 1 single 5 5 9
group helium type 1
velocity helium set 0 0 -219.57 # 10 eV for helium-4
create_atoms 1 single 1 1 9
group allhelium type 1
group newhelium subtract allhelium helium
velocity newhelium set 0 0 219.57

variable speed equal sqrt(vz[1]^2+vy[1]^2+vx[1]^2)
variable zvel equal vz[1]
variable speed2 equal sqrt(vz[2]^2+vy[2]^2+vx[2]^2)
variable zvel2 equal vz[2]

thermo_style custom step ke etotal
thermo 1

fix 1 all print 1 "Speed of helium #1 = \{speed\} \({zvel} in z-direction)"
fix 3 all print 1 "Speed of helium #2 = \{speed2\} \({zvel2} in z-direction)"
fix 2 all nve
run 50

The default for the velocity command is units lattice, which means the
vel you set is in units of lattice spacings/ps, not Angs/ps. Since your
lattice is 3.1, you are seeing a velocity 3x larger than you expect.
Use "units box" if you don't want the lattice setting.

Steve