How can I calculate thermal expansion coefficient from cross-correlations between volume and enthalpy?

Dear all,

I am new in LAMMPS and I have to calculate the volumetric thermal expansion coefficient (AlphaP) by performing on-the-fly Analysis.

To do this, I created input script below, where the AlphaP is calculated by using equation (10) of the following articles:
https://pubs.acs.org/doi/10.1021/ct200731v

The equation (10) is <δVδH>=kBT^2**[AlphaP],

where V, H, kB, and T are volume, enthalpy, Boltzmann constant, and temperature, respectively. Here, δV and δH mean V- and H-; the fluctuations of volume and enthalpy.

I tried this analysis towards 500 water molecules at NPT ensemble during 2ns simulations after equilibrium of initial 2.552 nm^3 box.

However, I gained negative values of AlphaP (-0.0027/K) at 300K, because δV and δH have a negative correlation coefficient (-0.30).

(The version is LAMMPS 64-bit 5Sep2018.)

When I tried towards the same system in GROMACS, AlphaP is +0.00060/K and the correlation coefficient between δV and δH is +0.24.
(the experiment data of AlphaP is +0.000256/K.)

I have attached the correlation graphs of δV-δH in LAMMPS and GROMACS. The δH in LAMMPS is larger than that in GROMACS.

This is my input file in LAMMPS.

LAMMPS.png

GROMACS.png

Your simulation might be too small to properly sample the fluctuations and averages accurately.