MSD is not always computed correctly

Dear lammps user,

My lammps version is 29 Sep 2021.

I am trying to compute the mean square displacement “MSD” using the “compute msd” command.

When I submit the below script to run 200 times in the computer cluster at my university, most of the runs generate reasonable MSD values but a few runs generate unreasonable MSD values.

I attached a plot of the MSD versus time in the log-log scale for different runs of the below script.

In the attached Figure, for example, msd.167 means the run number 167 out of 200 runs of the below script. This run, number 167, outputs unreasonable values for MSD at short time scales while the other runs generate reasonable values as it is clear from the Figure in the attached file

I run the script using :
lmp -in in.ellipsoids -screen none -var a “some integer”

Could you please tell me what could be wrong?

Thanks in advance

Here is the script

units         lj
dimension     3
boundary      p p p
atom_style    ellipsoid
lattice        sc 0.0001

region         box block 0.0 13.0 0.0 13.0 0.0 13.0 units lattice
create_box     1 box
create_atoms   1 box

set           type 1 shape 2.5 1 1
set           type 1 quat 0 0 1 0

velocity      all create 1.0 ${seed1} mom yes rot yes dist gaussian
neighbor      0.3 bin
neigh_modify  delay 0 every 1 check yes

pair_style    gayberne_rep 1.0 1.0 1.0 5.0 # gamma upsilon mu cutoff
pair_coeff    1 1 1 1 1.0 1.0 1.0 1.0 1.0 1.0

variable       t_log equal logfreq2(1,45,10)
variable       fixed_angle equal 0

timestep        0.0001

set           type 1 quat 0 0 1 ${fixed_angle}

velocity        all create 1.0 ${seed2} mom yes rot yes dist gaussian

fix            213 all brownian/asphere 1.0 ${seed3} gamma_t_eigen 1.089 1.089 0.859 gamma_r_eigen 9.351 9.351 1.667 rng gaussian

run           50000

reset_timestep  0

timestep        0.0001

set           type 1 quat 0 0 1 ${fixed_angle}

compute         msd all msd com yes
variable        twopoint equal c_msd[4]/6/(step*dt+1.0e-6)
fix             9 all vector 10 c_msd[4]
variable        fitslope equal slope(f_9)/6/(10*dt)

variable        msd_var_x equal c_msd[1]
variable        msd_var_y equal c_msd[2]
variable        msd_var_z equal c_msd[3]
variable        msd_var equal c_msd[4]

fix      msdlog all print v_t_log "$(step) $(v_msd_var_x:%6.11f) $(v_msd_var_y:%6.11f) $(v_msd_var_z:%6.11f) $(v_msd_var:%6.11f)" file msdMe.$a

dump           222 all custom 100 coordMe.$a id x y z c_orient[*]
dump_modify    222 every v_t_log first yes pbc yes sort id

run            1000000

The msd value at short time scale is irrelevant for a freely diffusing particle. According to the Einstein relation, only its long term behavior is supposed to converge to the same value and thus allows to compute the self-diffusion from its slope.