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