Heat capacity

Hello Dear Lammps user

I used the fluctuation of energy for calculate the heat capacity at constant volume for ionic liquid in NVT ensemble, but the data that I get is abnormal.
The command that I used is:

compute 2 all pe
variable pesq equal c_2*c_2
fix Peave all ave/time 100 10 1000 c_2 v_pesq file PEavg.dat

I calculated the heat capacity at constant volume with the potential energy(c_2) and it’s squared (v_pesq) and this formula:

Cv *(kBT^2) = <E^2> - ^2

Is it a right way to calculate the heat capacity at constant volum?


The formula looks correct, as long as by kBT^2 you really mean kB * T^2. How are the values you get abnormal?

I did the simulation at 398 K, and start with a volume 48000 A with 1000 A gap to the 148000 A in NVT ensembles. Pressure varied from 10000 MPa to -25 MPa. Heat capacity fluctuated a lot, I fitted that with a equation of state and I get R square about 0.23.

Comparing your results to: experiment, another simulation of your same system, analytical results?
You could be facing a units issue in Lammps. Which units are you using and did you make sure your KbT is properly expressed in those units. There are many other potential sources for trouble…


"Heat capacity fluctuated a lot ... R square about 0.23." I think this
means you did not run
long enough to get accurate estimate of <E^2>. Just hoping that
everything will be ok generally does not work in MD. You have to
first learn how to get good results. I suggest you start by testing
out your approach against published results for a standard potential,
such as Lennard-Jones.

I will also mention that fluctuation formulae, while based on sound
statistical mechanics principles, don't always behave well in
practice. A more robust method is to estimate <E> at multiple
temperatures, do a curve fit, and then analytically generate d<E>/dT.
If you are really paranoid, run NVE dynamics, sample T, and do a curve
fit to get dT(E)/dE