# Example 12: Calculation of specific heat#

In this example, the specific hear for Cu is calculated using a combination of free energy calculation and temperature sweep.

The calculation block gives the input conditions at which the calculation is carried out. First of all, the mode is ts, a temperature sweep over the temperature range given in the temperature keyword will be done.

The calculation can be run by,

calphy -i input.yaml


calculation. In this case, its 0 as we ran only a single loop (see n_iterations).

The variation of the free energy within the temperature range is given in the temperature_sweep.dat files in the simulation folder. The file contains three columns, temperature, free energy and the error in free energy. The files are read in and plotted below.

[1]:

import numpy as np
import matplotlib.pyplot as plt
import scipy.constants as sc


And some constants which we will use later,

[2]:

ev_to_j = sc.physical_constants["electron volt-joule relationship"][0]
A_to_m = 1E-10
kB = sc.physical_constants["Boltzmann constant"][0]

[3]:

temp, fe, err = np.loadtxt("ts-fcc-solid-400-0/temperature_sweep.dat", unpack=True)


## Calculation of $$c_p$$ from free energy#

The first step is to calculate the entropy from free energy using the following equation,

$-S = \bigg ( \frac{dF}{dT} \bigg )_V$

Then differentiating once more provides the specific heat,

$c_p = T\bigg ( \frac{dS}{dT} \bigg )_p$

We start by calculating the entropy by taking the gradient, converting it to J/mol K. But before that, we will fit the free energy to a polynomial so that taking gradients is easier. In this case, we will use a fifth order polynomial.

[63]:

fe_fit_fn = np.polyfit(temp, fe, 5)


Now we take the derivative of the free energy with respect to the temperature.

[69]:

entropy_fn = -np.polyder(fe_fit_fn, m=1)*ev_to_j/(1/Av)


Now we use this equation:

$c_p = T\bigg ( \frac{dS}{dT} \bigg )_p$

to calculate the specific heat

[70]:

cp_fn = np.polyder(entropy_fn, m=1)


Evaluating this polynomial gives the specific heat

[71]:

cp_temp = np.polyval(cp_fn, temp)


Note that we calculated:

$\frac{c_p}{T} = \bigg ( \frac{dS}{dT} \bigg )_p$

Note the division by $$T$$. The reason for this is to reduce the noise in our calculation when comparing with experimental results. Note that if a longer switching time and system size is used, we can use the actual equation. For more details see Fig. 4 in Ref.

In order to compare our results, we plot it together with some experimental values from Ref.

[74]:

exp_temp = [400, 450, 500, 550, 600, 650, 700, 800,900,1000, 1100, 1200]
exp_cp = [25.25, 25.61, 25.91, 26.21, 26.48, 26.73, 26.99, 27.48, 28.04, 28.66, 29.48, 30.53]


Finally the plot

[76]:

plt.plot(temp, cp_temp)
plt.plot(exp_temp, np.array(exp_cp)/np.array(exp_temp), 'o')
plt.xlabel("Temperature (K)", fontsize=12)
plt.ylabel(r"$C_p/T$ (J/mol)", fontsize=12)


We can see that there is good agreement between the two values.