Dear LAMMPS Developers and Users,
I am trying to calculate self-diffusion coefficient Ds of CO2 using MSD and Einstein relation. The Ds (slope of MSD(t) divided by 6) calculated using LAMMPS is too large (35 A^2/ps, should be approx. 28 A^2/ps) so, the molecules move too fast. The reference are results from 10.1021/je5009526 (calculated in GROMACS, also compared to experiment), and calculated by me in RASPA code. For CO2 I use fully flexible model (LJ and Coulomb from 10.1021/jp810871f, bond and angle same as 10.1021/je5009526).
The cubic cell I use for simulations has 65.25^3 Angst^3 and contains 500 molecules. The temperature of simulation is set using fix nvt to 450 K (for supercritical CO2). This is related to the pressure of 100 bar (from NIST, confirmed in mentioned article and by RASPA code). After equilibration in NVT ensemble, the pressure fluctuates around 200 bar (too high). After inspecting KE and virial part of pressure I figured out that KE term is (on average) approx. 325 bar (virial term is -125 bar from which pair+kspace is -2 bar and bond is -123 bar). Using ideal gas EOS (which is used in case of KE term of pressure) for 450 K pressure should be approx. 110 bar. The difference most likely caused by the fact that LAMMPS considers 1500 atoms in unit cell (2 oxygen and 1 carbon per CO2) and in ideal gas EOS I used N = 500 (number of molecules).
If the KE term of pressure is too high, then also KE is too high. Should I define temperature in a different way in my simulations (maybe based on molecules, not on individual atoms)? Is it the reason that the Ds (i.e. MSD) is too high in my case?
In case of liquid CO2 (298 K and 10 MPa) the difference in the Ds (compared with literature and RASPA) and pressure (compared with NIST) is even larger.
Definition of molecule is same as in GCMC example where only charges are changed according to force field I use.
My input (number of production steps is decreased, I tried for 5 ns simulation, including averaging MSD over different time intervals):
Units metal
atom_style full
region simbox block 0 65.25 0 65.25 0 65.25
create_box 2 simbox &
bond/types 1 &
angle/types 1 &
extra/bond/per/atom 2 &
extra/angle/per/atom 1 &
extra/special/per/atom 2
molecule CO2 CO2_new.dat
create_atoms 0 random 500 10001 simbox mol CO2 1001
pair_style lj/cut/coul/long 12 12
pair_coeff 1 1 0.002579411 2.745
pair_coeff 2 2 0.007382511 3.017
kspace_style ewald 1e-6
pair_modify mix arithmetic tail no shift yes
bond_style harmonic
angle_style harmonic
bond_coeff 1 111.305413 1.16
angle_coeff 1 12.81 180
mass 1 12.0107
mass 2 15.9994
group co2 type 1 2
group carbon type 1
group oxygen type 2
###minimization run
thermo 1
thermo_style custom step pe fmax fnorm
minimize 0.0 0.0 400 100000
parameters for dynamics
timestep 0.001
fix f2 co2 nvt temp 450 450 0.15 tchain 5
equilibration run
thermo_style custom step time temp pe ke
thermo 100
run 5000
data gathering run
compute msd_co2 co2 msd
compute temp co2 temp
thermo_style custom step time temp c_temp press c_msd_co2[4]
thermo 10
run 1000
Thank you!
Filip Formalik,
Wroclaw University of Science and Technology, Wroclaw, Poland