Phase diagrams and how to calculate them

DPG Tutorial: Creating and Running Automated Workflows for Material Science Simulations
17 March 2024

In this notebook, we will use the interatomic potentials for calculation of thermodynamic properties such as Helmholtz and Gibbs free energies, which in turn can be used for the calculation of phase diagrams. We will discuss calphy, the tool for automated calculation of free energies, and the methology involved.

A simple phase diagram

Phase diagrams provide a wealth of information such as: coexisting lines, melting temperature, phase stability, nucleation mechanism.

Calculation of phase diagrams: the essential ingredients

Phase diagrams can be evaluated from free energy diagrams. The required input are:

• \(G(P, T)\) for unary systems

• \(G(x, T)\) for binary systems

To calculate this, we employ thermodynamic integration, in conjuction with a number of methodological improvements.

Theory behind free energy calculation

Please see the end of this notebook, or this publication for a longer discussion on the theory and algorithms.

Start with imports

from pyiron_atomistics import Project
import matplotlib.pyplot as plt
import numpy as np

pr = Project('thermodynamics') 

Due to our limited time, we will show the calculation of melting temperature, the end point on the phase diagram. We leverage a software tool for calculation of free energies, calphy along with pyiron for efficient and automated free energy calculations. For a detailed tutorial, see here.

Calculation of melting temperature

Solid free energy

al_sol = pr.create.job.Calphy("Al_fcc", 
                              delete_aborted_job=True, 
                              delete_existing_job=True)

There are a number of input the job can take. We can gain some information about the job from the docstrings.

structure = pr.create.structure.ase.bulk('Al', cubic=True).repeat(4)

structure.plot3d()

al_sol.structure = structure

al_sol.list_potentials()

al_sol.potential = '2009--Mendelev-M-I--Al-Mg--LAMMPS--ipr1'

al_sol.server.cores = 4

al_sol.input.equilibration_control = 'berendsen'

al_sol.calc_free_energy(temperature=[800, 1100], 
                     pressure=0, 
                     reference_phase="solid",
                     n_equilibration_steps=5000,
                     n_switching_steps=5000)

Before we actually run the calculation, let us discuss the various parameters. temperature keyword gives the temperature range over which the free energy is to be calculated. Since we provide 500, the free energy is calculated at this temperature. pressure denotes the pressure of the calculation, we chose 0 in this case. Since we are using a solid FCC lattice, we set reference_phase to "solid". This means that the Einstein crystal will be used as the reference system. Finally, we have n_equilibration_steps and n_switching_steps. n_equilibration_steps denotes the number of MD steps over which the system is equilibrated to the required temperature and pressure. n_switching_steps are the number of MD steps over which the system is continuously transformed between the given interatomic potential, and the reference Einstein crystal.

Now we can actually run the calculation

al_sol.run()

Liquid free energy

al_lqd = pr.create.job.Calphy("Al_lqd", delete_aborted_job=True)
al_lqd.structure = structure
al_lqd.potential = '2009--Mendelev-M-I--Al-Mg--LAMMPS--ipr1'
al_lqd.server.cores = 4
al_lqd.calc_free_energy(temperature=[800, 1100], 
                     pressure=0, 
                     reference_phase="liquid",
                     n_equilibration_steps=5000,
                     n_switching_steps=5000)
al_lqd.run()

Finally, the results

plt.plot(al_sol.output.temperature, al_sol.output.energy_free,
        label="Al solid", color='#C62828')
plt.plot(al_lqd.output.temperature, al_lqd.output.energy_free,
        label="Al liquid", color='#006899')
plt.xlabel("Temperature (K)")
plt.ylabel("Free energy (eV/K)")
plt.legend();

Please check out the sections below for more detailed discussions…

Calculation of free energies: Thermodynamic integration

• free energy of reference system are known: Einstein crystal, Uhlenbeck-Ford model

• the two systems are coupled by $\( H(\lambda) = \lambda H_f + (1-\lambda)\lambda H_i \)$

• Run calculations for each \(\lambda\) and integrate $\( G_f = G_i + \int_{\lambda=0}^1 d\lambda \bigg \langle \frac{\partial H(\lambda)}{\partial \lambda } \bigg \rangle \)$

To calculate \(F\),

• for each phase

  • for each pressure

    • for each temperature

      • for each \(\lambda\)

If we choose 100 different \(\lambda\) values; 100 calculations are needed for each temperature and pressure!

Dimensionality: (phase, \(P\), \(T\), \(\lambda\))

Speeding things up: Non-equilibrium calculations

Non-Equilibrium Hamiltonian Interpolation

In this method:

• Discrete \(\lambda\) parameter is replaced by a time dependent \(\lambda(t)\)

• Instead of running calculations at each \(\lambda\), run a single, short, non-equilibrium calculation

\[ G_f = G_i + \int_{t_i}^{t_f} dt \frac{d\lambda (t)}{dt} \frac{ H(\lambda)}{\partial \lambda } \]

As discussed:

• the coupling parameter \(\lambda\) earlier is replaced by a time dependent parameter

• The equation also no longer has an ensemble average

These aspects makes it quite easy and fast to estimate this integral.

However:

• this equation holds when the switching betwen the system of interest and reference system is carried out infinitely slowly

• Practically, this is not not possible.

Therefore we can write:

\[ \Delta G = W_{rev} = W_s - E_d \]

\[ W_s = \int_{t_i}^{t_f} dt \frac{d\lambda (t)}{dt} \frac{ H(\lambda)}{\partial \lambda } \]

• \(E_d\) is the energy dissipation

• \(E_d \to 0\) when \(t_f-t_i \to \infty\)

So far, so good, but how is this useful?

• Instead of a single transformation from system of interest to reference, we switch back too

• These are called forward \((i \to f)\) and backward \((f \to i)\) switching

• \(t_f - t_i = t_{sw}\) is the switching time in each direction

• If \(t_{sw}\) is long enough, \(E_d^{i \to f} = E_d^{f \to i}\)

• and \(\Delta G = \frac{1}{2} (W_s^{i \to f} - W_s^{f \to i})\)

Now, we have all the components required for actual calculations.

We have also managed to successfully reduce the dimensionality

• for each phase

  • for each pressure

    • for each temperature

      • ~~for each \(\lambda\)~~

Dimensionality: (phase, \(P\), \(T\))

So, how do we calculate the free energy of a system modelled with a given interatomic potential?

Free energy for solids

Task: Find free energy of Al in FCC lattice at 500 K and 0 pressure

1. Create an Al FCC lattice

2. Choose an interatomic potential

3. Run MD calculations at 500 K and 0 pressure to equilibrate the system

4. Introduce the reference system

5. Switch….

6. …..

Steps 1-3 should be fairly easy, we saw this in the last days and also in the first session. But how do we introduce a reference system?

• A reference system is simply one for which the free energy is analytically known (\(G_i\))

• We calculate the free energy difference between this and the system of interest.

In case of solids, a good choice of a reference system is the Einstein crystal. An Einstein crystal is a set of independent harmonic oscillators attached to the lattice positions.

The free energy of the Einstein crystal is:

\[ F_E = 3 k_B T \sum_{i} ln \bigg ( \frac{\hbar \omega_i}{k_B T} \bigg ) \]

We need to calculate:

• \(\omega\)

• A common way is $\( \frac{1}{2} k_i \langle (\Delta \pmb{r}_i)^2 \rangle = \frac{3}{2} k_\mathrm{B} T \)$

• \(\langle (\Delta \pmb{r}_i)^2 \rangle\) is the mean squared displacement.

Now that we know about the reference system, let’s update our pseudo workflow:

1. Create an Al fcc lattice

2. Choose an interatomic potential

3. Run MD calculations at 500 K and 0 pressure to equilibrate the system

4. Calculate the mean squared displacement, therefore spring constants

5. Switch system of interest to reference system

6. Equilibrate the system

7. Switch back to system of interest

8. Find the work done

9. Add to the free energy of the Einstein crystal

As you can see, there are a number of steps that need to be done. This is where calphy and pyiron come in. These tools automatise all of the above steps and makes it easy to perform these calculations.

Free energy as function of temperature

Gibb’s free energy via reversible scaling at a constant pressure is given by,

\( G(N,P,T_f) = G(N,P,T_i) + \dfrac{3}{3}Nk_BT_f\ln{\dfrac{T_f}{T_i}} + \dfrac{T_f}{T_i}\Delta G \),

Therefore, \(G(N,P,T_f)\) can be computed from \(G(N,P,T_i)\) via the free energy difference \(\Delta G\).

Here, \(\Delta G = \dfrac{1}{2}[W_{if}-W_{fi}\)] — (2)

The reversible work is related to the internal energy \(U\) by, \(W = \int_{1}^{\lambda_f}<U> \,d\lambda\) — (3)

Using MD \(W\) can be computed as:

• equilibrate for time \(t_{eq}\) in NPT ensemble

• switch \(\lambda\) : \(1->\dfrac{T_f}{T_i}\) over time \(t_{sw}\)

• calculate work \(W_{if}\) from (3)

• equilibrate for time \(t_{eq}\) in NPT ensemble

• switch \(\lambda\) : \(\dfrac{T_f}{T_i}->1\) over time \(t_{sw}\)

• calculate work \(W_{fi}\) from (3).

Free energy for liquids

How is the liquid prepared in this calculation?

• Start from the given structure

• This structure is heated until it melts.

• Melting of the structure is automatically detected by calphy

• Once melted, it is equilibrated to the required temperature and pressure.

What about the reference system for liquid?

The reference system for the Liquid structure is also different. In this case, the Uhlenbeck-Ford system is used as the reference system for liquid.

The Uhlenbeck-Ford model is described by,

\[ E = - \epsilon \log(1-\exp(-r^2/\sigma^2)) \]

where,

\[ \epsilon = p k_B T \]

\(\epsilon\) and \(\sigma\) are energy and length scales, respectively.

It is purely repulsive liquid model which does not undergo any phase transformations.

Comparison of melting temperature methods

Further reading

• Menon et al. “From electrons to phase diagrams with classical and machine learning potentials: automated workflows for materials science with pyiron”

• Menon, Sarath, Yury Lysogorskiy, Jutta Rogal, and Ralf Drautz. “Automated Free-Energy Calculation from Atomistic Simulations.” Physical Review Materials 5, no. 10 (October 11, 2021): 103801..

• Freitas, Rodrigo, Mark Asta, and Maurice de Koning. “Nonequilibrium Free-Energy Calculation of Solids Using LAMMPS.” Computational Materials Science 112 (February 2016): 333–41..

• Paula Leite, Rodolfo, and Maurice de Koning. “Nonequilibrium Free-Energy Calculations of Fluids Using LAMMPS.” Computational Materials Science 159 (March 2019): 316–26..

• Koning, Maurice de, A. Antonelli, and Sidney Yip. “Optimized Free-Energy Evaluation Using a Single Reversible-Scaling Simulation.” Physical Review Letters 83, no. 20 (November 15, 1999): 3973–77..

• Paula Leite, Rodolfo, Rodrigo Freitas, Rodolfo Azevedo, and Maurice de Koning. “The Uhlenbeck-Ford Model: Exact Virial Coefficients and Application as a Reference System in Fluid-Phase Free-Energy Calculations.” The Journal of Chemical Physics 145, no. 19 (November 21, 2016): 194101..

Software used in this notebook

• calphy

• pyiron_atomistics

• LAMMPS