Dear all,
I wanted to determine the self diffusion coefficient of CO2 in supercritical phase. THe slope of msd vs time gives me the value 1.6x10^-8 (m^2/s). According to the literature it should come out of 10^8 order i.e. between 2x10^8 to 3x10^8 (m^2/s). I did minimization then used NVE and then NVE. I achieved the equilibration in 1ns and then did production run for 1 ns. I used real unit with the following lj parameters shown in the input file. Please suggest me what changes should I make to my input script to get the accurate result.
units real
atom_style full
variable T equal 323.2
variable V equal vol
variable dt equal 1
setup problem
dimension 3
boundary p p p
read_data data.1200sco2
pair_style lj/cut/coul/long 12 12
pair_coeff 1 1 0.0573209 2.7918
pair_coeff 2 2 0.164254 3.0
pair_coeff 1 2 0.097 2.9
bond_style harmonic
angle_style harmonic
kspace_style pppm 1.0e-4
bond_coeff 1 2566.68 1.16
angle_coeff 1 295.41 180.0
neighbor 2.0 bin
neigh_modify every 10 delay 20 check yes
group co2 type 1 2
min_style cg
dump 1 all custom 1000 water2_k_MSD_RDF2.lammpstrj id mol type x y z ix iy iz
thermo 1000
thermo_style custom step etotal temp vol press density
minimize 1e-5 1e-7 5000000 10000000
equilibration and thermalization
neighbor 3.0 bin
neigh_modify delay 0 every 1 check yes
neigh_modify exclude molecule co2
velocity all create $T 102486 mom yes rot yes dist gaussian
fix NVE all nve
thermo_style custom step etotal temp vol press density
thermo 1000
timestep ${dt}
run 1000000
unfix NVE
fix NVT all nvt temp $T T 1
timestep {dt}
thermo 1000
dump 2 all custom 1000 cd2000_conductivity.lammpstrj id mol type x y z vx vy vz ix iy iz
run 1000000
compute msd co2 msd com yes
fix msd co2 ave/time 10 100 1000 c_msd[4] file zz_wat.msd
compute rdf co2 rdf 1000 1 1 1 2 2 2
fix rdf co2 ave/time 10 100 1000 c_rdf file test.rdf mode vector
thermo 1000
dump 3 all xyz 100 dump.xyz
timestep ${dt}
thermo_style custom step etotal temp vol press density
run 1000000