Thermodynamic-integration and temperature-integration simulations¶
This example illustrates how to carry out Thermodynamic-Integration (TI) and Temperature-Integration (TEI) simulations with icet in the canonical ensemble. The TI method allows you to calculate the free energy at a given temperature for a system, whereas the TEI approach allows you to compute the change in free energy as a function of temperature.
Background¶
Thermodynamic integration is used to estimate the free energy between two given Hamiltonians, \(H_A\) and \(H_B\). This is carried out by integrating along a thermodynamic path between the two Hamiltonians. The free energy, F, of system B at temperature \(T_0\) is given by,
\[F_B = F_A + \int_0^1 \left\langle\frac{\mathrm{d}H(\lambda)} {\mathrm{d}\lambda}\right\rangle_{H} \mathrm{d}\lambda\]
where \(H(\lambda) = (1 - \lambda) H_A + \lambda H_B\).
Hence, if we chose \(H_A\) in such a way that its free energy is known, then we can then calculate the free energy of our system of interest, \(H_B\). A good choice is \(H_A = 0\) which corresponds to a completely disordered system with free energy directly proportional to the ideal mixing entropy, which is analytically known,
\[S_{\text{ideal}} = k_B\ln{\Omega},\]
where \(\Omega\) is the multiplicity, the total number of states.
This choice of \(H_A = 0\) also allows us to extract the temperature dependence of the free energy, \(F_B(T)\). This can be done from the following relations, for more details see [FreAstKon16],
\[F_B(T) = \frac{1}{\lambda} \left ( - T_{0} S_{\text{ideal}} + \int_0^\lambda \left\langle\frac{\mathrm{d}H(\lambda)} {\mathrm{d}\lambda}\right\rangle_{H} \mathrm{d}\lambda \right )\]
where the temperature T is
\[T = \frac{T_{0}}{\lambda}\]
We will also introduce temperature integration, which is a similar formulation to calculate the free energy. This requires a reference free energy, which can be obtained from e.g., thermodynamic integration or from the ideal mixing entropy in the high temperature limit. If this reference is known we can calculate the free energy as a function of temperature using the internal energy, U(T),
\[\frac{F(T_2)}{T_2} = \frac{F(T_1)}{T_1} - \int_{T_1}^{T_2}\frac{U(T)}{T^2}\mathrm{d}T\]
2D Ising model¶
The two-dimensional Ising model is well suited for demonstrating the TI and TEI algorithm and has been studied in the previous Wang-Landau simulation tutorial as well as references therein .
We use the same setup as in the Wang-Landau example here. The following code generates a cluster expansion that represents the 2D Ising model.
prim = Atoms('Au', positions=[[0, 0, 0]], cell=[1, 1, 10], pbc=True)
cs = ClusterSpace(prim, cutoffs=[1.01], chemical_symbols=['Ag', 'Au'])
ce = ClusterExpansion(cs, [0, 0, 2])
For computational convenience, we consider a very small system of only \(4\times4 = 16\) sites and we setup the calculator together with some parameters.
supercell = prim.repeat((4, 4, 1))
calc = ClusterExpansionCalculator(supercell, ce)
n_integration_steps = 400000
n_equilibration_steps = 1000
temperature_max = 30
temperature_min = 0.2
temperature_max_plot_limit = 4
n_Ag = 8
start_configuration = supercell.copy()
for i in range(n_Ag):
start_configuration[i].symbol = 'Ag'
Running a TI simulation¶
We run an equilibration using the CanonicalEnsemble
to get a reasonable starting structure at high temperatures for the thermodynamic integration.
mc = CanonicalEnsemble(
structure=start_configuration,
calculator=calc,
temperature=temperature_max,
boltzmann_constant=1,
trajectory_write_interval=None,
ensemble_data_write_interval=200)
mc.run(n_equilibration_steps)
We now want to sample the configurational space on the path between the two Hamiltonians,
i.e., on the path \(\lambda \in [0,1]\), this sampling is carried out by
ThermodynamicIntegrationEnsemble.
Note that we set the ensemble_data_write_interval to 1, this is to get the free energy
for as many temperatures as possible, but it is also to increase the accuracy of the integration later on.
Although, keep in mind that it is likely not necessary to have it as tightly sampled in the general case.
The integral is then calculated in get_free_energy_thermodynamic_integration.
We set forward to True to indicate that we are going from \(H_A\) to \(H_B\)
mc = ThermodynamicIntegrationEnsemble(
structure=mc.structure, calculator=calc,
temperature=temperature_min,
forward=True,
ensemble_data_write_interval=1,
boltzmann_constant=1,
n_steps=n_integration_steps)
mc.run()
data_container = mc.data_container
(_, free_energy_integration_forward) = \
get_free_energy_thermodynamic_integration(data_container, cs,
forward=True,
max_temperature=temperature_max_plot_limit,
boltzmann_constant=1)
We then proceed to do the same thing but in the opposite direction. This should decrease
the statistical error. forward set to False indicates that we are going from \(H_B\)
to \(H_A\)
mc = CanonicalEnsemble(
structure=mc.structure,
calculator=calc,
temperature=temperature_min,
boltzmann_constant=1,
trajectory_write_interval=None,
ensemble_data_write_interval=200)
mc.run(n_equilibration_steps)
mc = ThermodynamicIntegrationEnsemble(
structure=mc.structure, calculator=calc,
temperature=temperature_min,
forward=False,
ensemble_data_write_interval=1,
boltzmann_constant=1,
n_steps=n_integration_steps)
mc.run()
data_container = mc.data_container
(temperatures_integration, free_energy_integration_backward) = \
get_free_energy_thermodynamic_integration(data_container, cs,
forward=False,
max_temperature=temperature_max_plot_limit,
boltzmann_constant=1)
The average between the two runs is then our result for the free energy
free_energy_integration_average = 0.5 * (free_energy_integration_forward +
free_energy_integration_backward)
Running a TEI simulation¶
We can also compare with temperature integration. The procedure is very similar to the thermodynamic integration.
As previously state we need the free energy at one of the endpoints, in this example
we will use the high temperature limit, i.e., the ideal mixing entropy. This should be a good
approximation for a high enough reference temperature. This means that the temperature_max will
will be assumed to large enough for the system to be completely disordered and \(TS_B \ggg U(T)\).
Start by running a CanonicalAnnealing
simulation from high to low temperature. We start from an equilibrated structure obtained from a
CanonicalEnsemble simulation at a high
temperature.
Note that we set the ensemble_data_write_interval to 1 again, with the same justification here as
for the TI simulation.
We can then extract the temperature dependence of the free energy from the canonical annealing as
outlined above. The integral is calculated in get_free_energy_temperature_integration. Where forward set to
True indicates that we are going from high to low temperature.
mc = CanonicalEnsemble(
structure=start_configuration,
calculator=calc,
temperature=temperature_max,
boltzmann_constant=1,
trajectory_write_interval=None,
ensemble_data_write_interval=200)
mc.run(n_equilibration_steps)
mc = CanonicalAnnealing(
structure=mc.structure,
calculator=calc,
T_start=temperature_max,
T_stop=temperature_min,
cooling_function='linear',
n_steps=n_integration_steps,
boltzmann_constant=1,
trajectory_write_interval=None,
ensemble_data_write_interval=1)
mc.run()
data_container = mc.data_container
data_container = mc.data_container
(temperatures_temperature, free_energy_temperature_forward) = \
get_free_energy_temperature_integration(data_container,
cs,
forward=True,
temperature_reference=temperature_max,
max_temperature=temperature_max_plot_limit,
boltzmann_constant=1)
We do the same thing but in the opposite direction now and calculate the free energy as an average
between the two runs from low to high and high to low temperature. Again, forward set to
False indicates that we are going from low to high temperature.
mc = CanonicalEnsemble(
structure=mc.structure,
calculator=calc,
temperature=temperature_min,
boltzmann_constant=1,
trajectory_write_interval=None,
ensemble_data_write_interval=200)
mc.run(n_equilibration_steps)
mc = CanonicalAnnealing(
structure=mc.structure,
calculator=calc,
T_start=temperature_min,
T_stop=temperature_max,
cooling_function='linear',
n_steps=n_integration_steps,
boltzmann_constant=1,
trajectory_write_interval=None,
ensemble_data_write_interval=1)
mc.run()
data_container = mc.data_container
(temperatures_temperature, free_energy_temperature_backward) = \
get_free_energy_temperature_integration(data_container,
cs,
forward=False,
temperature_reference=temperature_max,
max_temperature=temperature_max_plot_limit,
boltzmann_constant=1)
The average between the two runs is then our result for the free energy
free_energy_temperature_average = 0.5 * (free_energy_temperature_forward +
free_energy_temperature_backward)
Comparing with the analytical solution¶
We can then compare the two methods with the analytical solution. The analytical solution is obtained by directly calculating the partition function. This is possible for a small configuration space when we can enumerate all possible structures.
Source code¶
The complete source code is available in
examples/free_energy_canonical_ensemble.py
from ase import Atoms
from ase.units import kB
import numpy as np
import matplotlib.pyplot as plt
from icet import ClusterExpansion, ClusterSpace
from mchammer.calculators import ClusterExpansionCalculator
from mchammer.ensembles import ThermodynamicIntegrationEnsemble
from mchammer.ensembles import CanonicalAnnealing
from mchammer.ensembles import CanonicalEnsemble
from mchammer.free_energy_tools import (get_free_energy_thermodynamic_integration,
get_free_energy_temperature_integration)
from itertools import combinations
def get_all_structures(struct, n_Ag):
"""
Generates all possible lattice occupations
"""
permutation = list(combinations(range(len(struct)), n_Ag))
atoms = []
for perm in permutation:
atom = struct.copy()
atom.symbols[np.array(perm)] = 'Ag'
atoms.append(atom)
return atoms
def get_free_energy_enumeration(temperatures, atoms, calc, boltzmann_constant=kB):
"""
Calculate the free energy by brute force
"""
energies = []
for atom in atoms:
energies.append(calc.calculate_total(occupations=atom.numbers))
energies = np.array(energies)
free_energy = []
for T in temperatures:
beta = 1 / (boltzmann_constant * T)
Z = np.exp(-beta * energies).sum()
free_energy.append(-boltzmann_constant * T * np.log(Z))
return np.array(free_energy)
##########################################################################################
# Setup
##########################################################################################
# Prepare cluster expansion
prim = Atoms('Au', positions=[[0, 0, 0]], cell=[1, 1, 10], pbc=True)
cs = ClusterSpace(prim, cutoffs=[1.01], chemical_symbols=['Ag', 'Au'])
ce = ClusterExpansion(cs, [0, 0, 2])
# Prepare initial configuration and calculator
supercell = prim.repeat((4, 4, 1))
calc = ClusterExpansionCalculator(supercell, ce)
n_integration_steps = 400000
n_equilibration_steps = 1000
temperature_max = 30
temperature_min = 0.2
temperature_max_plot_limit = 4
n_Ag = 8
start_configuration = supercell.copy()
for i in range(n_Ag):
start_configuration[i].symbol = 'Ag'
##########################################################################################
# Thermodynamic integration
##########################################################################################
# Prepare for the TI by running the canonical ensemble at a high temperature.
mc = CanonicalEnsemble(
structure=start_configuration,
calculator=calc,
temperature=temperature_max,
boltzmann_constant=1,
trajectory_write_interval=None,
ensemble_data_write_interval=200)
mc.run(n_equilibration_steps)
# Run the TI from the disordered to the ordered system.
mc = ThermodynamicIntegrationEnsemble(
structure=mc.structure, calculator=calc,
temperature=temperature_min,
forward=True,
ensemble_data_write_interval=1,
boltzmann_constant=1,
n_steps=n_integration_steps)
mc.run()
data_container = mc.data_container
(_, free_energy_integration_forward) = \
get_free_energy_thermodynamic_integration(data_container, cs,
forward=True,
max_temperature=temperature_max_plot_limit,
boltzmann_constant=1)
# Prepare for the TI by running the canonical ensemble at a low temperature.
mc = CanonicalEnsemble(
structure=mc.structure,
calculator=calc,
temperature=temperature_min,
boltzmann_constant=1,
trajectory_write_interval=None,
ensemble_data_write_interval=200)
mc.run(n_equilibration_steps)
mc = ThermodynamicIntegrationEnsemble(
structure=mc.structure, calculator=calc,
temperature=temperature_min,
forward=False,
ensemble_data_write_interval=1,
boltzmann_constant=1,
n_steps=n_integration_steps)
mc.run()
data_container = mc.data_container
(temperatures_integration, free_energy_integration_backward) = \
get_free_energy_thermodynamic_integration(data_container, cs,
forward=False,
max_temperature=temperature_max_plot_limit,
boltzmann_constant=1)
# The best approximation for TI is now given by the average of the two runs.
free_energy_integration_average = 0.5 * (free_energy_integration_forward +
free_energy_integration_backward)
##########################################################################################
# Temperature integration
##########################################################################################
# Prepare for the temperature integration by running the canonical ensemble at a high temperature.
mc = CanonicalEnsemble(
structure=start_configuration,
calculator=calc,
temperature=temperature_max,
boltzmann_constant=1,
trajectory_write_interval=None,
ensemble_data_write_interval=200)
mc.run(n_equilibration_steps)
mc = CanonicalAnnealing(
structure=mc.structure,
calculator=calc,
T_start=temperature_max,
T_stop=temperature_min,
cooling_function='linear',
n_steps=n_integration_steps,
boltzmann_constant=1,
trajectory_write_interval=None,
ensemble_data_write_interval=1)
mc.run()
data_container = mc.data_container
data_container = mc.data_container
(temperatures_temperature, free_energy_temperature_forward) = \
get_free_energy_temperature_integration(data_container,
cs,
forward=True,
temperature_reference=temperature_max,
max_temperature=temperature_max_plot_limit,
boltzmann_constant=1)
# Prepare for the temperature integration by running the canonical ensemble at a low temperature.
mc = CanonicalEnsemble(
structure=mc.structure,
calculator=calc,
temperature=temperature_min,
boltzmann_constant=1,
trajectory_write_interval=None,
ensemble_data_write_interval=200)
mc.run(n_equilibration_steps)
mc = CanonicalAnnealing(
structure=mc.structure,
calculator=calc,
T_start=temperature_min,
T_stop=temperature_max,
cooling_function='linear',
n_steps=n_integration_steps,
boltzmann_constant=1,
trajectory_write_interval=None,
ensemble_data_write_interval=1)
mc.run()
data_container = mc.data_container
(temperatures_temperature, free_energy_temperature_backward) = \
get_free_energy_temperature_integration(data_container,
cs,
forward=False,
temperature_reference=temperature_max,
max_temperature=temperature_max_plot_limit,
boltzmann_constant=1)
# The best approximation for TEI is now given by the average of the two runs.
free_energy_temperature_average = 0.5 * (free_energy_temperature_forward +
free_energy_temperature_backward)
##########################################################################################
# Analytical solution
##########################################################################################
# Get all possible occupations
atoms = get_all_structures(supercell, n_Ag)
# Gets the free energy from the enumerated structures
temperatures = np.linspace(temperature_min, temperature_max_plot_limit, 1000)
free_energy_enumeration = get_free_energy_enumeration(temperatures,
atoms,
calc,
boltzmann_constant=1)
##########################################################################################
# Plotting
##########################################################################################
natoms = len(supercell)
fig, ax = plt.subplots()
ax.plot(temperatures_integration, free_energy_integration_forward / natoms, color='tab:red',
label='forward', ls='-', lw=2)
ax.plot(temperatures_integration, free_energy_integration_backward / natoms, color='tab:blue',
label='backward', ls='-', lw=2)
ax.plot(temperatures_integration, free_energy_integration_average / natoms, color='tab:grey',
label='average', ls='-', lw=2)
ax.set_xlim((temperature_min, temperature_max_plot_limit))
ax.set_ylim((-45 / natoms, -30 / natoms))
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('free energy / atoms (eV)')
ax.legend()
fig.tight_layout()
fig.savefig('thermodynamic_integration_free_energy_forward_backward.svg')
fig, ax = plt.subplots()
ax.plot(temperatures_temperature, free_energy_temperature_forward / natoms, color='tab:red',
label='forward', ls='-', lw=2)
ax.plot(temperatures_temperature, free_energy_temperature_backward / natoms, color='tab:blue',
label='backward', ls='-', lw=2)
ax.plot(temperatures_temperature, free_energy_temperature_average / natoms, color='tab:grey',
label='average', ls='-', lw=2)
ax.set_xlim((temperature_min, temperature_max_plot_limit))
ax.set_ylim((-45 / natoms, -30 / natoms))
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('free energy / atoms (eV)')
ax.legend()
fig.tight_layout()
fig.savefig('temperature_integration_free_energy_forward_backward.svg')
fig, ax = plt.subplots()
ax.axhline(-2, label='Ground state energy', color='black', ls='--')
ax.plot(temperatures, free_energy_enumeration / natoms, color='tab:red',
label='enumeration', lw=2)
ax.plot(temperatures_integration, free_energy_integration_average / natoms, color='tab:blue',
label='thermodynamic integration', ls='-', lw=2)
ax.plot(temperatures_temperature, free_energy_temperature_average / natoms, color='tab:grey',
label='temperature integration', ls='-', lw=2)
ax.set_xlim((temperature_min, temperature_max_plot_limit))
ax.set_ylim((-45 / natoms, -30 / natoms))
ax.set_xlabel('Temperature (K)')
ax.set_ylabel('free energy / atoms (eV)')
ax.legend()
fig.tight_layout()
fig.savefig('free_energy_canonical_ensemble.svg')