# Calculating the heat capacity at constant pressure of a solid in its liquid state using AIMD as implemented in VASP

I am trying to calculate the heat capacity at constant pressure of Silicon in its liquid state. For the same, I am using the example given in VASPwiki https://www.vasp.at/wiki/index.php/Liquid_Si_-_Standard_MD. To start with, the NVT ensemble is employed and the example is repeated for different temperatures (TEBEG tag in INCAR). The time steps are increased till the complete melting of solid Si is achieved and it can be inferred from the pair correlation function. Further, the energy is plotted against the number of time steps till the energy convergence is ensured. Now, the converged energy at a fixed average temperature is known. Similarly, with varying temperature, many energies can be obtained. As per my understanding, the energy is the internal energy and its derivative with respect to temperature can give heat capacity at constant volume. However, the heat capacity at constant pressure is required to finally get the Gibbs energy.

Can anyone please explain a detailed way of getting this heat capacity at constant pressure? Also, please let me know if I am going wrong somewhere.

I don't work on this stuff, by I checked a book I use a lot. If you can find a copy, look at $$\S$$ 2.5 in 'Computer simulations of liquids (2nd ed.)' by Allen. Someone with experience in this topic can probably be more useful than me.

I also searched "specific heat molecular dynamics" in google and found a similar equation in this paper. Look at eq. 6! There are more refs. there if you want to look elsewhere too. You can also try googling the same thing.

A more naive approach might be to extract the internal energy $$U(T,P)$$ from the VASP calculation and average over the steps in the well converged region. Calculate $$U(T+\delta T,P)$$ and $$U(T-\delta T,P)$$ with $$\delta T$$ a small change in temperature. You can approximate the heat capacity as $$\left. \frac{\partial U(T,P)}{\partial T}\right|_{T,P} \approx \left. \frac{\Delta U(T,P)}{\Delta T}\right|_{T,P} = \frac{ U(T+\delta T,P) - U(T-\delta T,P)}{2\delta T}$$ and divide by mass or volume or number of particles to get a specific heat.

Hope this helps!

Ty