# [lammps-users] Temperature and kinetic energy in lammps, inaccuracies

Hello,

I have started a simple 3d LJ-spheres simulation (lj-units) and figured
out, that for example

T= 0.74663
E_pair= -5.99469
E_tot= -4.87586
E_kin(ke)= 1.11883

have some inaccuracies compared to calculations by hand using the values
above:

T!= E_kin * 2/3=0.74588 , difference in 10^-3
E_tot!=E_kin+E_pair= -4.87586, correct
E_tot!=T * 3/2 +E_pair= -4.87474, difference in 10^-3

It occurs in every timestep, when I calculate the slope of E_kin * 2/3 as
function of T, the slope is 0.99900 (should be 1) for every timestep,
therefore s.th. is going wrong.

I do not understand is differences, because in the manual it is written,
that the temperature in thermo_style is by default calculated via
ComputeTemp and that the variable ke is calculated by default by
[(degrees of freedom) * temperature / 2], therefore I do not expect any
inaccuracy besides the rounding error of double, which is much less than
10^-3.

Below are the both methods in lammps (I have searched in the code). The
only possible explanation I would think about is, that force->boltz in
lj-units is not 1.00000, but set to s.th. like 0.998... unfortunately the
parameter boltz in force.h is not initialised in force.cpp and somehow
accessable from every (thanks to public), therefore I could not find it.

void Thermo::compute_ke()
{
dvalue = temperature->scalar;
dvalue *= 0.5 * temperature->dof * force->boltz;
if (normflag) dvalue /= natoms;
}

double ComputeTemp::compute_scalar()
{
invoked_scalar = update->ntimestep;

double **v = atom->v;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *type = atom->type;
int nlocal = atom->nlocal;

double t = 0.0;

if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
t += (v[i]*v[i] + v[i]*v[i] + v[i]*v[i]) *
rmass[i];
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
t += (v[i]*v[i] + v[i]*v[i] + v[i]*v[i]) *
mass[type[i]];
}

MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
if (dynamic) dof_compute();
scalar *= tfactor;
return scalar;
}

Here the script for the simulation:

# Point dipoles in a 3d box

units lj
atom_style atomic
dimension 3

lattice sc 0.84
region box block -5 5 -5 5 -5 5
create_box 1 box
create_atoms 1 box

mass 1 1.0

velocity all create 0.5 87287 mom yes

pair_style lj/cut 5.0
pair_coeff * * 1.0 1.0

neighbor 0.3 bin
neigh_modify delay 0

#fix 1 all nvt temp 0.75 0.75 100.0
fix 1 all nve

timestep 0.001

thermo_style custom step temp epair etotal press ke
thermo 500

dump 1 all custom 1000 dump2.lj id type x y z

run 100000