Free energy calculations
In this tutorial we will demonstrate how one can compute the free energy of solids accounting for the full anharmonicity using thermodynamic integration (TI). For a comprehensive discussion of the subject see reference 1 and in particular the excellent article (reference 2) by Freitas et al.
Introduction
Thermodynamic integration
Consider a Hamiltonian \(H(\lambda)\) that is parametrically dependent on a parameter \(\lambda\). One can show that the free energy difference between two thermodynamic states corresponding to \(\lambda_i\) and \(\lambda_f\) can be computed as
The integral can be interpreted as the reversible work done along the path from \(\lambda_i\) to \(\lambda_f\),
In equilibrium TI approaches this integral is discretized and one carries out separate simulations for each value of \(\lambda\).
In practice non-equilibrium approaches, in which \(\lambda\) is updated continuously throughout a simulation, are, however, usually (much) more efficient and we also adopt this approach here. In this case, one computes the irreversible work down along a path,
where \(\lambda=\lambda(t)\) is a function of the simulation time \(t\) and \(\Gamma(t)\) is the phase-space trajectory traversed during the simulation.
This basic principle can be employed in several different ways to obtain free energy differences and, given a suitable reference state with a known free energy, ultimately also total free energies. Here, we consider * the Frenkel-Ladd approach, in which the Hamiltonian is gradually switched from one Hamiltonian (e.g., an Einstein crystal) to another Hamiltonian (e.g., the system of interest) and * the reversible scaling (RS) approach, in which one integrates the work done as thermodynamic variable (typically the temperature) changes.
Frenkel-Ladd approach
We consider two Hamiltonians (interatomic potentials), the first one \(H_A\) is the interatomic potential of interest (in this example a NEP model), and the second one \(H_B\) is a reference potential for which the free energy is analytically known. We define a mixing Hamiltonian
where \(\lambda\) is the mixing parameter, such that if \(\lambda=0\) we recover \(H_A\) and if \(\lambda=1\) we obtain \(H_B\).
The free energy differences between the system represented by \(H_A\) and \(H_B\) can be otained by integration as follows
where \(\left < H_B - H_A \right >_H\) is the ensemble average of the energy difference between potential \(B\) and \(A\).
Typically it is efficient to sample the integral in the expression above via a molecular dynamics (MD) simulation, during which \(\lambda\) is continously changed between 0 and 1. The change in \(\lambda\) can be carried out in both directions, i.e., with \(\lambda\) changing from \(0 \to 1\) or from \(1 \to 0\), which are referred to as “forward” and “backward” switching, respectively. This type of integration typically incurs some a systematic error due to the simulation not being in equilibrium due to \(\lambda\) changing continously in the simulation. However, these systematic errors are the same but with opposite sign in forward and backward integration and hence the average over these two lead to a much more accurate free energy [2],
Alternatively one can also run multiple runs with fixed values of \(\lambda\), but this is typically much less efficient.
In order for the thermodynamic integration to converge and be as accurate as possible the reference Halmitonian should be chosen as “similar” to the Hamiltonian of the system of interest as possible. This means that the phase space sampled by the reference potential at the given temperature should be as close as possible to the phase space sampled by the real potential.
A good and common choice of reference system is an Einstein crystal, where each atom is tethered to a reference position with a harmonic spring, i.e.,
where \(k\) is the spring constant (corresponding to a frequency \(\omega_\mathrm{E}\)), \(r_i\) is the position of atom \(i\) and \(r_i^0\) is the reference position. A good choice for the spring constant is such that the mean square displacements of the Einstein crystal is similar to that of the reference system. Or it can be chosen such that the frequency of the Einstein solid is similar to a typical/average phonon frequency in the system.
The free energy of the Einstein crystal is given by
Note that since we are running classical MD simulations we need to use the classical harmonic free energy and not the quantum variant.
There is a small correction due to the center of mass being fixed in MD simulations and thus effectivley, three degrees of freedoms are not sampled [2]. This is, however, a very small effect and usually (i) cancels out when considering free energy difference between two phases and (ii) vanishes in the thermodynamic limit, i.e., for large system sizes.
Note that if your system exhibits any kind of diffusion, an Einstein crystal will be a poor reference.
Reversible scaling
In the reversible scaling (RS) approach each value of \(\lambda\) corresponds to a different temperature, and one considers the scaled Hamiltonian
In practice one carries out a simulation during which \(\lambda\) varies with time. Afterwards one calculates the irreversible work
and then obtains the free energy as a (continuous) function of temperature via
Application case
In this tutorial we will compute the free energy of titanium in two phases, body-centered cubic (BCC) and hexagonal close-packed HCP using a NEP potential. Titanium undergoes a phase transition from HCP to BCC at around 1150 K experimentally. Here we will use the universal NEP from [5] provided in the file nep-unep-v1.txt.
Here we only consider a single set of parameters for system size, simulation length, and, in the case of the Frenkel-Ladd approach, the value of the spring constant. In general one should, however, carefully check the convergence of the results with respect to these parameters.
References
For more detailed information see, e.g., * [1] Understanding Molecular Simulation by D. Frenkel and B. Smit * [2] Nonequilibrium free-energy calculation of solids using LAMMPS, Freitas et al., Computational Materials Science 112, 333 (2016) * [3] https://en.wikipedia.org/wiki/Thermodynamic_integration * [4] GPUMD documentation, https://gpumd.org/gpumd/input_parameters/ensemble_ti_spring.html * [5] General-purpose machine-learned potential for 16 elemental metals and their alloys, Song et al., arxiv.2311.04732 (2024); zenodo record with model
Frenkel-Ladd approach
Thermal expansion
The switching simulations with respect to the Einstein crystal are carried out in an NVT ensemble. Therefore we need to know the temperature dependence of the lattice parameters beforehand. This information can be readily obtained through NPT simulations. Here, we use a simple polynomial to the data show below.
[1]:
def get_Ti_bcc_lattice_parameter(T):
alat = 3.2700559351528953 + 1.311174893837472e-05 * T + 6.272944257392548e-09 * T ** 2
return alat
def get_Ti_hcp_lattice_parameter(T):
alat = 2.9337472912371387 + 2.9726486008937908e-05 * T + -7.026081914635316e-09 * T ** 2
clat = 4.655766212690983 + 2.8534727442310983e-05 * T + 2.2462470604534456e-08 * T ** 2
return alat, clat
Setting up the runs
The switching is carried out by first equilibrating for n_eq number of time steps, then swithcing forward for n_sample steps, followed by equilibrating again and then switching backwards. The total simulation time is thus n_total = 2 * (n_eq + n_sample).
The simulation cell is chosen to be as close to a cubic shape as possible in order to maximize the distance between periodic images.
See the documentation of the ti_spring command for more details.
[2]:
import numpy as np
# MD parameters
temperatures = np.arange(700, 1401, 50)
dt = 2.0 # timestep in fs
T_damp = 100 # thermostat damping
thermo_interval = 100 # thermo.out write interval
spring_k = 2.0 # spring constant
potential_fname = 'nep-unep-v1.txt'
# number of timesteps to run
n_eq = 20 * 1000
n_sample = 100 * 1000
n_total = 2 * (n_eq + n_sample)
[3]:
import os
from ase.build import bulk
from calorine.gpumd.io import write_xyz
# set up MD run directories and input files
for phase in ['bcc', 'hcp']:
for T in temperatures:
# set up structure
if phase == 'bcc':
alat = get_Ti_bcc_lattice_parameter(T)
prim = bulk('Ti', crystalstructure='bcc', a=alat, cubic=True)
size = 20
elif phase == 'hcp':
alat, clat = get_Ti_hcp_lattice_parameter(T)
prim = bulk('Ti', crystalstructure='hcp', a=alat, c=clat, orthorhombic=True)
size = [20, 12, 13]
supercell = prim.repeat(size)
# set up directory
dirname = f'my_freen_runs/switch_{phase}_T{T}_nat{len(supercell)}' \
f'_nsample{n_sample}_spring{spring_k:5f}'
os.makedirs(dirname, exist_ok=True)
# write model.xyz
write_xyz(os.path.join(dirname, 'model.xyz'), supercell)
# write run.in
with open(os.path.join(dirname, 'run.in'), 'w') as f:
f.write(f'potential {os.path.abspath(potential_fname)}\n')
# write equilibration
f.write(f'velocity {1.0 * T}\n')
f.write(f'ensemble nvt_lan {T} {T} {T_damp}\n')
f.write(f'time_step {dt}\n')
f.write(f'run {n_eq}\n\n\n')
# write production
f.write(f'ensemble ti_spring temp {T} tperiod {T_damp}'
f' tequil {n_eq} tswitch {n_sample} spring Ti {spring_k}\n')
f.write(f'time_step {dt}\n')
f.write(f'dump_thermo {thermo_interval}\n')
f.write(f'run {n_total}\n')
After these simulations have been run (outside this notebook), we can collect the data and extract the free energies.
Inspecting a single run
First we carry out the integration “manually” by using the data in ti_spring.csv to make sure it works as expected for one simulation.
[4]:
import pandas as pd
from ase.io import read
dirname = 'my_freen_runs/switch_bcc_T1000' \
'_nat16000_nsample100000_spring2.000000'
# read data
df = pd.read_csv(f'{dirname}/ti_spring.csv')
atoms_ideal = read(f'{dirname}/model.xyz')
n_atoms = len(atoms_ideal)
# free energy from backward and forward switching
n2 = int(len(df) / 2)
df_fwd, df_bwd = df[:n2], df[n2:]
display(df_fwd)
| lambda | dlambda | pe | espring | |
|---|---|---|---|---|
| 0 | 0.000000e+00 | 0.000000e+00 | -7.654844 | 0.287045 |
| 1 | 1.259958e-23 | 6.299748e-23 | -7.654805 | 0.286966 |
| 2 | 4.031731e-22 | 1.007919e-21 | -7.654784 | 0.286906 |
| 3 | 3.061494e-21 | 5.102388e-21 | -7.654764 | 0.286843 |
| 4 | 1.290068e-20 | 1.612542e-20 | -7.654763 | 0.286797 |
| ... | ... | ... | ... | ... |
| 99995 | 1.000000e+00 | 0.000000e+00 | -7.339000 | 0.128034 |
| 99996 | 1.000000e+00 | 0.000000e+00 | -7.339146 | 0.128034 |
| 99997 | 1.000000e+00 | -2.273464e-18 | -7.339221 | 0.128034 |
| 99998 | 1.000000e+00 | -2.273555e-18 | -7.339266 | 0.128045 |
| 99999 | 1.000000e+00 | 0.000000e+00 | -7.339302 | 0.128053 |
100000 rows × 4 columns
[5]:
from matplotlib import pyplot as plt
fig, ax = plt.subplots(1, 1, figsize=(4.0, 2.5), dpi=200)
ax.plot(df_fwd['lambda'], (df_fwd.pe - df_fwd.espring),
'-', lw=3.0, alpha=0.5, label='forward')
ax.plot(df_bwd['lambda'], (df_bwd.pe - df_bwd.espring),
'-', lw=1.5, alpha=0.5, label='backward')
ax.set_ylabel('$U_\mathrm{NEP} - U_\mathrm{Einstein}$ (eV/atom)')
ax.set_xlabel('$\lambda$')
ax.set_xlim(0, 1)
ax.legend(frameon=False)
fig.tight_layout()
We can now obtain the free energy by integrating (or rather summing) over \(\lambda\). By averaging over the forward and backward switching runs the systematic error due to the finite switching rate cancels.
The free energy obtained in this fashion is also directly available from the yaml output file. And indeed the Free energy (average) value agrees that we obtained by “manually” integrating the data from the csv file agrees with the value for F given in the yaml output file.
[6]:
import yaml
from ase.units import _amu as amu, kB, _e as eV, _hbar as hbar
# work done by switching
work_switch_fwd = np.sum(df_fwd.dlambda * (df_fwd.pe - df_fwd.espring))
work_switch_bwd = np.sum(df_bwd.dlambda * (df_bwd.espring - df_bwd.pe))
work_switch_ave = (work_switch_fwd + work_switch_bwd) / 2
# define harmonic reference system free energy
mass = atoms_ideal.get_masses()#[0]
omegas = np.sqrt(spring_k * eV / (mass * amu)) * 1e10 # in 1/s
F_oscillators = 3 * kB * T * np.log(hbar / eV * omegas / (kB * T)) # in eV
F_einstein = np.sum(F_oscillators) / n_atoms
# free energies for NEP model
F_fwd = work_switch_fwd + F_einstein
F_bwd = work_switch_bwd + F_einstein
F_avg = work_switch_ave + F_einstein
print(f'Free energy (forward) : {F_fwd:12.7f} eV/atom')
print(f'Free energy (backward) : {F_bwd:12.7f} eV/atom')
print(f'Free energy (average) : {F_avg:12.7f} eV/atom')
print('')
# compare to values in ti_spring.yaml
with open(f'{dirname}/ti_spring.yaml') as f:
row = yaml.load(f, Loader=yaml.FullLoader)
for key, val in row.items():
print(f'{key:22} : {val:12.7f}')
Free energy (forward) : -8.5193304 eV/atom
Free energy (backward) : -8.5190133 eV/atom
Free energy (average) : -8.5191718 eV/atom
E_Einstein : -0.4847130
E_diff : -7.7187970
F : -8.2035100
T : 1000.0000000
V : 17.7965640
P : 0.0000000
G : -8.2035100
Temperature dependence
Next, we analyze the free energy as a function of temperature for both BCC and HCP phases by parsing the yaml files from a series of runs.
[7]:
from glob import glob
dfs_spring = dict()
for phase in ['bcc', 'hcp']:
records = []
for fname in glob(f'my_freen_runs/switch_{phase}*/ti_spring.yaml'):
with open(fname) as f:
row = yaml.load(f, Loader=yaml.FullLoader)
records.append(row)
df = pd.DataFrame(records).sort_values('T', ignore_index=True)
dfs_spring[phase] = df
print(phase)
display(df)
bcc
| E_Einstein | E_diff | F | T | V | P | G | |
|---|---|---|---|---|---|---|---|
| 0 | -0.274753 | -7.738216 | -8.012969 | 700.0 | 17.681046 | 0.0 | -8.012969 |
| 1 | -0.307756 | -7.735010 | -8.042766 | 750.0 | 17.698996 | 0.0 | -8.042766 |
| 2 | -0.341620 | -7.731752 | -8.073372 | 800.0 | 17.717466 | 0.0 | -8.073372 |
| 3 | -0.376293 | -7.728475 | -8.104768 | 850.0 | 17.736457 | 0.0 | -8.104768 |
| 4 | -0.411727 | -7.725250 | -8.136977 | 900.0 | 17.755969 | 0.0 | -8.136977 |
| 5 | -0.447880 | -7.722060 | -8.169939 | 950.0 | 17.776004 | 0.0 | -8.169939 |
| 6 | -0.484713 | -7.718797 | -8.203510 | 1000.0 | 17.796564 | 0.0 | -8.203510 |
| 7 | -0.522192 | -7.715521 | -8.237713 | 1050.0 | 17.817649 | 0.0 | -8.237713 |
| 8 | -0.560287 | -7.712365 | -8.272653 | 1100.0 | 17.839260 | 0.0 | -8.272653 |
| 9 | -0.598970 | -7.709060 | -8.308030 | 1150.0 | 17.861399 | 0.0 | -8.308030 |
| 10 | -0.638216 | -7.705951 | -8.344166 | 1200.0 | 17.884068 | 0.0 | -8.344166 |
| 11 | -0.678000 | -7.702832 | -8.380831 | 1250.0 | 17.907266 | 0.0 | -8.380831 |
| 12 | -0.718301 | -7.699630 | -8.417930 | 1300.0 | 17.930996 | 0.0 | -8.417930 |
| 13 | -0.759099 | -7.696553 | -8.455652 | 1350.0 | 17.955259 | 0.0 | -8.455652 |
| 14 | -0.800376 | -7.693403 | -8.493780 | 1400.0 | 17.980057 | 0.0 | -8.493780 |
hcp
| E_Einstein | E_diff | F | T | V | P | G | |
|---|---|---|---|---|---|---|---|
| 0 | -0.274753 | -7.751881 | -8.026634 | 700.0 | 17.674366 | 0.0 | -8.026634 |
| 1 | -0.307756 | -7.746786 | -8.054542 | 750.0 | 17.697600 | 0.0 | -8.054542 |
| 2 | -0.341620 | -7.741753 | -8.083373 | 800.0 | 17.720854 | 0.0 | -8.083373 |
| 3 | -0.376293 | -7.736778 | -8.113071 | 850.0 | 17.744130 | 0.0 | -8.113071 |
| 4 | -0.411727 | -7.731945 | -8.143672 | 900.0 | 17.767427 | 0.0 | -8.143672 |
| 5 | -0.447880 | -7.727120 | -8.175000 | 950.0 | 17.790744 | 0.0 | -8.175000 |
| 6 | -0.484713 | -7.722391 | -8.207104 | 1000.0 | 17.814080 | 0.0 | -8.207104 |
| 7 | -0.522192 | -7.717811 | -8.240003 | 1050.0 | 17.837436 | 0.0 | -8.240003 |
| 8 | -0.560287 | -7.713271 | -8.273558 | 1100.0 | 17.860809 | 0.0 | -8.273558 |
| 9 | -0.598970 | -7.708817 | -8.307787 | 1150.0 | 17.884201 | 0.0 | -8.307787 |
| 10 | -0.638216 | -7.704417 | -8.342633 | 1200.0 | 17.907609 | 0.0 | -8.342633 |
| 11 | -0.678000 | -7.700044 | -8.378044 | 1250.0 | 17.931033 | 0.0 | -8.378044 |
| 12 | -0.718301 | -7.695847 | -8.414148 | 1300.0 | 17.954472 | 0.0 | -8.414148 |
| 13 | -0.759099 | -7.691721 | -8.450821 | 1350.0 | 17.977925 | 0.0 | -8.450821 |
| 14 | -0.800376 | -7.687757 | -8.488134 | 1400.0 | 18.001392 | 0.0 | -8.488134 |
We can now plot the free energy as a function of temperature which shows that at some temperature between 1000 and 1300 K the BCC phase becomes more stable than the HCP phase. Due to the very similar slopes it is, however, difficult to locate the exact transition temperature. To obtain a clearer picture it is therefore preferable to consider the free energy difference \(\Delta F\) between the two phases, which is shown in the bottom panel. This yields a transition temperature of about 1140 K, in good agreement with the experimentally reported value of 1155 K.
This analysis incidentally also highlights the biggest challenge in free energy calculations, namely that one the free energy differences must be very well converged in order to control the error in the transition temperature. To obtain accuracte results it is therefore a good idea to run multiple identical independent simulatons, and average the thus obtained free energies to reduce the error. This also allows one provide a standard error estimate for the free energy by computing the standard deviation of all obtained free energies.
[8]:
fig, axes = plt.subplots(figsize=(3.4, 3.5), dpi=200, nrows=2, sharex=True)
for phase in ['bcc', 'hcp']:
axes[0].plot(dfs_spring[phase]['T'], dfs_spring[phase].F,
'-o', label=phase, alpha=0.5)
ax = axes[1]
ax.plot(dfs_spring['bcc']['T'], 1e3 * (dfs_spring['bcc'].F - dfs_spring['hcp'].F),
'-o', alpha=0.5)
ax.axvline(1136, color='0.5', alpha=0.5, ls=':')
ax.axhline(0, color='0.5', alpha=0.5, ls=':')
axes[0].set_ylabel('$F$ (eV/atom)')
axes[1].set_ylabel('$\Delta F$ (meV/atom)')
axes[-1].set_xlabel('Temperature (K)')
fig.tight_layout()
fig.align_ylabels()
fig.subplots_adjust(hspace=0)
Reversible scaling approach
We now also consider the reversible scaling (RS) approach. In the following cell we set up the input files for two runs, one for BCC and one for HCP.
[9]:
initial_temperature, final_temperature = 700, 1500
n_switch = 40000
n_equil = 1000
for phase in ['bcc', 'hcp']:
# set up structure
if phase == 'bcc':
alat = get_Ti_bcc_lattice_parameter(initial_temperature)
prim = bulk('Ti', crystalstructure='bcc', a=alat, cubic=True)
size = 20
elif phase == 'hcp':
alat, clat = get_Ti_hcp_lattice_parameter(initial_temperature)
prim = bulk('Ti', crystalstructure='hcp', a=alat, c=clat, orthorhombic=True)
size = [20, 12, 13]
supercell = prim.repeat(size)
# set up directory
dirname = f'my_freen_runs/rs_{phase}_Ti{initial_temperature}' \
f'_Ti{final_temperature}_nat{len(supercell)}_nswitch{n_switch}'
os.makedirs(dirname, exist_ok=True)
# write model.xyz
write_xyz(os.path.join(dirname, 'model.xyz'), supercell)
# write run.in
with open(os.path.join(dirname, 'run.in'), 'w') as f:
f.write(f'potential {os.path.abspath(potential_fname)}\n')
# write equilibration
f.write(f'velocity {initial_temperature}\n')
f.write(f'ensemble nvt_lan {initial_temperature} {initial_temperature} {T_damp}\n')
f.write(f'time_step {dt}\n')
f.write(f'run {n_eq}\n\n\n')
# write production
f.write(f'ensemble ti_rs temp {initial_temperature} {final_temperature}'
f' aniso 0 tswitch {n_switch} tequil {n_equil}\n')
f.write(f'time_step {dt}\n')
f.write(f'dump_thermo {thermo_interval}\n')
f.write(f'run {2 * n_switch + n_equil}\n')
After the two runs have finished we can analyze the results. As apparent from the expression for the total free energy in the RS approach (see above), we need both the variation of the energy (or enthalpy) throughout the simulation and the total free energy at the reference (initial) temperature. Here, the latter is taken from the Frenkel-Ladd approach.
For convenience we define a function that carries out all the data manipulation.
[10]:
from scipy.integrate import cumulative_trapezoid
def get_free_energy(
reference_freen_yaml: str,
reversible_switching_csv: str,
):
"""
Parameters
----------
reference_freen_yaml
Path to ti_spring.yaml file with the free energy
at the reference (initial) temperature.
reversible_switching_csv
Path to ti_rs.csv file with the data from the
reversible scaling run.
"""
# One requires the total free energy at the initial temperature.
# Here it is assumed to come from a calculation using the ti_spring command.
with open(reference_freen_yaml, 'r') as f:
y = yaml.safe_load(f)
reference_temperature = y['T']
reference_freen = y['G']
# Read data from reversible scaling (RS) run and split it into forward and backward sequences.
rs = pd.read_csv(reversible_switching_csv)
n = int(len(rs)/2)
forward = rs[:n]
backward = rs[n:][::-1].reset_index()
# Compute the irreversible work done along the lambda-path.
lmb = forward['lambda']
work_fwd = cumulative_trapezoid(forward['enthalpy'], lmb, initial=0)
work_bwd = cumulative_trapezoid(backward['enthalpy'], lmb, initial=0)
work = (work_fwd + work_bwd) / 2
# Compute the free energy.
freen = (reference_freen + 3 / 2 * kB * reference_temperature * np.log(lmb) + work) / lmb
temperature = reference_temperature / lmb
df = pd.DataFrame(np.transpose([temperature, freen]), columns=['temperature', 'free_energy'])
return df
Now we can plot the results, which shows that the results from RS approach are in very good agreement with the Frenkel-Ladd results. The RS approach requires far fewer simulations and thus is often preferable. It does, however, require a reference free energy, which needs to be very accurate in order to obtain reliable results. In practice it can therefore often be beneficial to add reference calculations (using Frenkel-Ladd) also at some higher temperature(s) in order to check the convergence of the results.
[11]:
dfs_rs = dict()
dfs_rs['bcc'] = get_free_energy(
'my_freen_runs/switch_bcc_T700_nat16000_nsample100000_spring2.000000/ti_spring.yaml',
'my_freen_runs/rs_bcc_Ti700_Ti1500_nat16000_nswitch40000/ti_rs.csv')
dfs_rs['hcp'] = get_free_energy(
'my_freen_runs/switch_hcp_T700_nat12480_nsample100000_spring2.000000/ti_spring.yaml',
'my_freen_runs/rs_hcp_Ti700_Ti1500_nat12480_nswitch40000/ti_rs.csv')
fig, axes = plt.subplots(figsize=(3.4, 3.5), dpi=200, nrows=2, sharex=True)
for phase in ['bcc', 'hcp']:
ax = axes[0]
kwargs = dict(
color = 'C0' if phase == 'hcp' else 'C1',
alpha=0.5
)
ax.plot(dfs_spring[phase]['T'], dfs_spring[phase].F, 'o', label=phase.upper(), **kwargs)
ax.plot(dfs_rs[phase].temperature, dfs_rs[phase].free_energy, **kwargs)
axes[0].legend(frameon=False, handlelength=1)
ax = axes[1]
kwargs = dict(color = '0.5', alpha=0.5)
ax.plot(dfs_spring['bcc']['T'], 1e3 * (dfs_spring['bcc'].F - dfs_spring['hcp'].F),
'o', label='Frenkel-Ladd', **kwargs)
ax.plot(dfs_rs['bcc'].temperature,
1e3 * (dfs_rs['bcc'].free_energy - dfs_rs['hcp'].free_energy),
label='reversible scaling', **kwargs)
ax.axvline(1136, ymax=0.6, color='0.5', alpha=0.5, ls=':')
ax.axhline(0, color='0.5', alpha=0.5, ls=':')
axes[0].set_ylabel('$F$ (eV/atom)')
axes[1].set_ylabel('$\Delta F$ (meV/atom)')
axes[-1].set_xlabel('Temperature (K)')
axes[-1].legend(frameon=False, handlelength=1)
fig.tight_layout()
fig.align_ylabels()
fig.subplots_adjust(hspace=0)
