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
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)
boundary p p p
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
variable dt_3 equal 0.000002
fix 1 all nve
compute msd all msd com yes
variable twopointmsd equal c_msd/(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
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
#fix 5 all vector 1 c_VACF
#variable diffvacf equal v_dt_3*trap(f_5)
#variable diffvacfd3 equal v_diffvacf/3
#thermo_style custom step v_diffvacfd3
Can anybody tell me, where is the problem ?