Dear all,

I want to verify my MD program output using lammps.

Using my program I calculated equilibrium state of N2 (as a point)

The experimental self-diffusion coefficient is 14.2 (in our units, nm^2 / ps), so in the simulation it can be slightly bigger or less.

My MD program does such calculation, and gets 13.8 by MSD and 17.3 by Green-Kubo (VACF).

The self-diffusion coefficient, calculated using MSD and VACF in LAMMPS differs a lot from the the experimental one.

Using LAMMPS with the script below I get next values for MSD:

Step twopoint fitslope

...

6999000 38250.766 36939.342

7000000 38250.995 36939.716

Step diffvacf

...

6999000 39858.095

7000000 39858.254

They are not correct, (in lammps nano units the self-diffusion coefficient should be close to 14200, and we have 36939 by MSD and 39858 by VACF)

units nano

dimension 3

boundary p p p

atom_style atomic

read_data test.data # equilibrium state of N2 (T = 273.14K)

pair_style mie/cut 0.91225 #such potential was used in my program

# eps sig ga gr

pair_coeff 1 1 0.0135 0.3649 6 11.5

mass 1 0.000046517 # mass of N2 molecula

timestep 0.000002 # the same timestep I used in my program

run_style verlet

variable dt_3 equal 0.000002

fix 1 all nve

compute msd all msd com yes

variable twopointmsd equal c_msd[4]/(3*2)/(step*v_dt_3+2.0e-6)

# the divison is from the einshtein formula for msd method, i looked through the LAMMPS src and find out that this operation do no performed by lammps.

fix 9 all vector 10 c_msd[4]

variable fitslopemsd equal slope(f_9)/6/(10*v_dt_3)

thermo_style custom step v_twopointmsd v_fitslopemsd

#for VACF it is

#compute VACF all vacf

#variable VACF4 equal c_VACF[4]

#fix 5 all vector 1 c_VACF[4]

#variable diffvacf equal v_dt_3*trap(f_5)

#variable diffvacfd3 equal v_diffvacf/3

#thermo_style custom step v_diffvacfd3

thermo $d

run 7000000

Can anybody tell me, where is the problem ?