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 *mask = atom->mask;

int nlocal = atom->nlocal;

double t = 0.0;

if (rmass) {

for (int i = 0; i < nlocal; i++)

if (mask[i] & groupbit)

t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *

rmass[i];

} else {

for (int i = 0; i < nlocal; i++)

if (mask[i] & groupbit)

t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *

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

Thanks in advance for your reply.

Best

Gerald Rosenthal