High fluctuations in the thermodynamic properties of mixtures

I am trying to obtain thermodynamic properties of liquid propellant mixtures such as N2H4 + H2O2. My codes ran without errors, but the results showed properties with high fluctuations, inconsistent with the experimental data. I am simulating in the NPT ensemble, with OPLS- All Atom, in the temperature range of 250 to 500 K and constant pressure of 10 atm. What I’ve already tried:
I doubled the number of molecules from 500 to 1000.
I increased the equilibration time to 2000000.
The fluctuations are still high, I am calculating the thermal capacity, the coefficient of thermal expansion, the isothermal compressibility, the Joule-Thomson coefficient, and the Speed of Sound.
Could anyone tell me what could be causing this problem?

My input:

#N2H4_H2O2_NPT

units real
atom_style full
variable T equal 250 # [K]
variable P equal 10 # [atm]
variable sigma equal 3.368 # [Å]
variable rc equal 2.5*${sigma} # [Å]

pair_style lj/cut/coul/long ${rc}
pair_modify mix arithmetic tail yes
bond_style harmonic
angle_style harmonic
dihedral_style harmonic

#kspace_style pppm 1.0e-4
kspace_style ewald 1.0e-4

read_data N2H4_H2O2_50.lmp_data

velocity all create $T 261086 rot yes dist gaussian

timestep 1.0

fix NVE all nve
run 1000
unfix NVE

fix NPT all npt temp $T $T 100.0 iso $P $P 1000.0
run 2000000

reset_timestep 0
variable dens equal density*1000 # [kg/m³]

thermo_style custom step temp press vol v_dens pe ke etotal etail
thermo 300

dump dump all custom 100 N2H4_H2O2_NPT.dump element id x y z
#dump_modify dump element N H H O sort id

run 1000000


The beginning of N2H4_H2O2_50.lmp_data :

LAMMPS data file

4 atom types
4 bond types
3 angle types
4 dihedral types

5000 atoms
4000 bonds
4000 angles
6500 dihedrals

-20.5239716 20.5239716 xlo xhi
-20.5239716 20.5239716 ylo yhi
-20.5239716 20.5239716 zlo zhi

Masses

1 1.00784 # HIW
2 15.999 # OW
3 14.0067 # NW
4 1.00784 # HW

Pair Coeffs

1 0.054 2.332 # HIW
2 0.012 2.100 # OW
3 0.7113 3.3368 # NW
4 0.0000 0.000 # HW

Bond Coeffs

1 535 1.475 # OW OW
2 553 0.950 # OW HIW
3 360.00 1.402 # NW NW
4 445.00 1.045 # NW HW

Angle Coeffs

1 35 94.8 # HIW OW OW
2 43.600 107.0 # HW NW HW
3 271.96 108.0 # NW NW HW

Dihedral Coeffs

1 4 1 3 # HIW OW OW HIW
2 0.3920 1 3 # HW NW NW HW
3 0.3920 1 3 # HW NW NW HW
4 0.3920 1 3 # HW NW NW HW

Atoms

1 1 2 -0.312 17.7830944 14.8480721 11.3975306 # O1
2 1 2 -0.312 18.1244774 13.6185637 11.0717010 # O2
3 1 1 0.312 18.6429005 15.3300190 11.5039358 # HI1
4 1 1 0.312 17.9972916 13.5720959 10.0895948 # HI2

Hey Juliana,

When you say that your results showed “properties with high fluctuations”, do you mean that the error bars you calculated are large? Or do you mean that you estimated these properties in several different thermodynamic states within the range you mentioned and you observe a very abrupt change in the values of the properties from one thermodynamic state to another?

PS: are you sure that this model you are using is supposed to be able to reproduce the properties you are studying?

There are a number of possible reasons. Here are the ones that come to my mind immediately:

  • too large a timestep (try 0.5 fs or 0.25fs)
  • forcefield has unsuitable bond parameters for your molecules
  • forcefield has unsuitable non-bonded parameters for your molecules
  • charge distribution is not correct for the force field and the specific compounds