The heat capacity is calculated from the fluctuations of the internal energy for the NVT ensemble as:

Cv = /kB T^2 = (-()^2) / kB T^2

So, what you need from the trajectory of the energy, calculate the averages of the E^2 and E, and then divide by kB T^2.

@Hamid,

You have to make sure you have simulated for long time, especially if there is a conformational transition.

Dear Hiqmet

I did same simulation with methanol (opls) and water (spc/e) . I get correct value for methanol but I don't get correct value for water. Why?

Production time is 1ns.

@Ali,

yes, that formula is for equilibrium MD simulations in NVT.

About your system, you can adjust the density by running MD in NPT ensemble first at standard pressure. Then, after system is equilibrated. You can fix your V at the average value you get from MD in NPT, and switch to NVT ensemble.

Hi Hamid,

You can see these references, they also discuss calculation of heat capacity from md simulations. It seems that classical md simulations overestimate it, but quantum corrections may be an improvement. You can also see results from different water models.

Shiga and Shinoda, jcp, 123, 134502, 2005.

Levitt, Hirshberg, Sharon, Laidig, Daggett, jpc B, 101, 5051, 1997.