Exercise 4: Thermodynamic Properties

Authors: Jan Janssen, Tilmann Hickel, Jörg Neugebauer (Max-Planck-Institut für Eisenforschung)

The scope of this first exercise is to become familar with:

  • with the phonopy interface in pyiron to calculate free energies,

  • the harmonic and quasi-harmonic approximation and

  • how to combine multiple pyiron jobs in one workflow.

With this last section we have all the necessary tools to calculate phase diagrams.

Update Notice

Since the workshop in 2020, there have been some updates to pyiron commands, and here in the exercises new commands are used.

The changes include:

  • Now pyiron has various modules, like pyiron_atomistics, pyiron_continuum, pyiron_experimental. So we can import Project from pyiron_atomistics directly.

  • using pr.create to create instances of various objects, such as structures and jobs. For example:

pr.create.structure.ase.bulk() ## this replaces pr.create_ase_bulk()

, or

pr.create.job.Lammps('job_name') ## this replaces pr.pr.create_job(job_type=pr.job_type.Lammps,'job_name')

Similary, one can create pyiron tables via pr.create.table().

By this approach, the user easily gets the available options via autocompeletion.

Please note that in the video tutorials, the old commands are still presented.


In the first session we learned how to create a pyiron project object and then use this pyiron project object to create atomistic structure objects.

# Import the Project object
from pyiron_atomistics import ____
# Create a Project object instance for a project named phonons
pr = ____("phonons")
# Create a cubic aluminum fcc structure
al_fcc = pr.create.structure.ase.bulk(____, _____=True)
# Confirm the final structure has 4 atoms by calculating 
# the length of the structure object
____(al_fcc) == 4

Harmonic Approximation

Calculate the phonons at the 0K equilibrium volume and then use the model of the harmonic oscillator to calculate the free energy, heat capacity and entropy contribution. These calculation can again be either executed with density functional theory (DFT) or interatomic potentials. In this workshop we use interatomic potentials to accelerate the process, but the same simulation protocol can also be executed with DFT code as quantum engine.


# Create a LAMMPS job to relax the structure - here we the regular minimization
# for DFT calculation it is recommended to calculate the energy volume curve 
# to identify the equilibrium volume as discussed in the previous section. 
# Name the minimization job "lmp_mini".
job_mini = pr.create.job.Lammps(
potential = "2009--Mendelev-M-I--Al-Mg--LAMMPS--ipr1"
job_mini.potential = potential
# Assign the the cubic fcc aluminium structure to the LAMMPS quantum engine job
job_mini._________ = _________
# Set the pressure during the minimization to zero 
# Execute the minimization

Phonon Calculation

# Create a interatomic template job with the LAMMPS quantum engine named lmp
job_lammp_template = pr.create.job.Lammps(
job_lammp_template.potential = potential
# Assign the the final structure of the minimization calculation job as input structure
job_lammp_template._________ = _________.get_structure()
# We calculate the bulk energy using the template job in a separate calculation named "lmp_bulk"
job_bulk = job_lammp_template.copy_to(
# Create a Phonopy job from the LAMMPS quantum engine named phono
phono = pr.create.job.PhonopyJob(
phono.ref_job = job_lammp_template
# Execute the Phonopy calculation and the bulk calculation

Density of States

# Plot the density of states over energy 

Free energy calculation

# To calculate the thermodynamic properties we define a temperature range
# Starting a 0K up to the 800K using 41 steps. We import numpy and define
# a linear space starting at 0 to 800 with 41 steps. 
import numpy as np
temperatures = np.linspace(______, _______, _____)
# Calculate the thermal properties for a the defined temperature range 
bulk_thermal_properties = phono.get_thermal_properties(temperatures=_________) 
# Import the matplotlib library for plotting.  
import matplotlib.pyplot as plt

# Plot the free energy over temperature by adding the inner bulk energy and the free vibrational energy 
plt.plot(__________, job_bulk.output.energy_pot[-1] + bulk_thermal_properties.free_energies)
plt.xlabel("Temperature [K]")
plt.ylabel("Free energy  ($U+F_{vib}$)  [eV]");

Heat capacity

# Plot the heat capacity over temperature
plt.plot(___________, bulk_thermal_properties.cv)
plt.xlabel("Temperature [K]")
plt.ylabel("Heat capacity");


# Plot the entropy over temperature
plt.plot(___________, bulk_thermal_properties.entropy)
plt.xlabel("Temperature [K]")

Quasi-harmonic Approximation

The quasi-harmonic approximation includes the thermal volume expansion in contrast to the harmonic approximation which only considers the 0K equilibrium volume.


# Create a interatomic template job with the LAMMPS quantum engine named lmp_strain
job_strain_template = pr.create.job.Lammps(
job_strain_template.potential = potential
# Assign the the final structure of the minimization calculation job as input structure
job_strain_template.structure = ________.get_structure()
# Create a secondary template job to calculate the bulk energy
# The second template job we name lmp_bulk_strain 
job_strain_bulk = job_strain_template.copy_to(
# Use this second template job to create a Murnaghan job named murn_strain
murn_strain = pr.create.job.Murnaghan(
murn_strain.ref_job = job_strain_bulk
# Create a Phonopy job from the original LAMMPS quantum engine named phono_strain
phono_strain = pr.create.job.PhonopyJob(
phono_strain.ref_job = job_strain_template
# Assign this Phonopy job to a QuasiHarmonicJob to calculate the phonons at multiple volumes
# This QuasiHarmonicJob is named quasi_strain
quasi_strain = pr.create.job.QuasiHarmonicJob(
quasi_strain.ref_job = phono_strain
# Set the end temperature to be 800K and set the number of temperature steps to 41 
quasi_strain.input["temperature_end"] = ________
quasi_strain.input["temperature_steps"] = _______
# Execute the QuasiHarmonicJob calculation and the bulk calculation

Free Energy

To calculate the heat capacity and entropy at constant pressure, the free energy at constant temperature is plotted over volume. The minimum at constant temperature defines the volume of constant pressure.

# Visualise the temperature dependent free energy over volume using a matplotlib color map
import matplotlib
cmap = matplotlib.cm.get_cmap('coolwarm')
# Iterate over the the output of the QuasiHarmonicJob to plot the temperature dependent free energy over volume 
for i, [t, free_energy, v] in enumerate(
    color = cmap(i/len(quasi_strain["output/temperatures"].T))
    # Add the energy of the Murnaghan Job to the vibrational free energy
    plt.plot(v, free_energy + _________["output/energy"], color=color)

# Add labels to the plot     
plt.ylabel("Free Energy")

# Add a color bar to visualise the temperature dependence 
temperatures = quasi_strain["output/temperatures"]
normalize = matplotlib.colors.Normalize(vmin=temperatures.min(), vmax=temperatures.max())
scalarmappaple = matplotlib.cm.ScalarMappable(norm=normalize, cmap=cmap)
cbar = plt.colorbar(scalarmappaple)

Pressure optimised thermodynamic properties

After calculating the volume expansion from the free energy surface over temperature and volume, the entropy and heat capacity are calculated at these optimised volumes.

# Use the optimise_volume() function of the QuasiHarmonicJob  
v0_lst, free_eng_lst, entropy_lst, cv_lst = quasi_strain.optimise_volume(
    # It requires the output energy of the energy volume curve as additional input 
temperature_output_lst = quasi_strain["output/temperatures"][0]
# Plot the finite temperature volume over temperature 
plt.plot(temperature_output_lst, __________)
# Plot the pressure optimised free energy over temperature 
plt.plot(temperature_output_lst, ___________)
plt.ylabel("Free Energy")
# Plot the pressure optimised entropy over temperature 
plt.plot(temperature_output_lst, ____________)
# Plot the pressure optimised heat capacity over temperature 
plt.plot(temperature_output_lst, ____________)
plt.ylabel("Heat Capacity")

Concentration dependence

With the temperature dependence of a unary discussed above the next step is the concentration dependence. For the concentration dependence we use the special quasi-random structures introduced in the first section to calculate mixed structures. Here we select the Ni-Cr phase diagram, which can be calculated using the 2018--Howells-C-A--Cr-Ni--LAMMPS--ipr1 potential.

Create Endmembers

Start by creating the endmembers of both phases the nickel fcc endmember and the chromium bcc endmember.

# Create the Nickel fcc endmember 
ni_fcc = pr.create.structure.ase.bulk(____, cubic=True)

# Create the Chromium bcc endmember 
cr_bcc_small = pr.create.structure.ase.bulk(____, cubic=True)

# Compare the length of both structures
len(ni_fcc), len(cr_bcc_small)
# Repeat the bcc cell to have the same number of atoms as the fcc cell
cr_bcc = cr_bcc_small.repeat([___, ___, ____])

# Compare the length of both structure
len(cr_bcc) == len(ni_fcc)

Select concentrations

One advantage of atomistic simulation is the ability to calculate free energies of unstable phases. In this case we calculate both the fcc phase and the bcc phase for all temperature ranges. Starting at 0% Cr up to 100%.

# Create 3 mixed concentrations for a 4 atom cell within the range 0.0 to 1.0 
# for example these could be 0.25, 0.5 and 0.75
concentration_lst = [____, ____, ___]
# Create an FCC Cr lattice by replacing all elements of the FCC Ni lattice with Cr
cr_fcc = ______.copy()
cr_fcc[:] = "Cr"
# Create an BCC Ni lattice by replacing all elements of the BCC Cr lattice with Ni
ni_bcc = ________.copy()
ni_bcc[:] = "Ni"
# Create a list of all concentrations 
concentration_total_lst = [0.0] + concentration_lst + [1.0]

SQS Structures

To calculate SQS structures for multiple concentrations the SQSMaster job is used, it takes an SQS job as an input in addtion to a list of concentrations. Here it is used to calculate SQS structures for both the FCC phase and the BCC phase.


# Create an SQS job named sqs_job_ni
sqs_job_ni = pr.create.job.SQSJob(
# Assign the cubic nickel fcc structure to the SQS job
sqs_job_ni.structure = ______
# Limit the number of iterations to 1000
sqs_job_ni.input['iterations'] = _____
# Create an SQSMaster named sqs_master_ni
master_ni = pr.create.job.SQSMaster( 
master_ni.ref_job = sqs_job_ni
# Assign the mixed concentration list defined above - not the total concentration list which includes the end members
master_ni.input["fraction_lst"] = __________
master_ni.input["species_one"] = "Ni"
master_ni.input["species_two"] = "Cr"
# Execute the calculation


# Create an SQS job named sqs_job_cr
sqs_job_cr = pr.create.job.SQSJob(
# Assign the cubic chromium bcc structure to the SQS job
sqs_job_cr.structure = __________
# Limit the number of iterations to 1000
sqs_job_cr.input['iterations'] = _______
# Create an SQSMaster named sqs_master_cr
master_cr = pr.create.job.SQSMaster( 
master_cr.ref_job = sqs_job_cr
# Assign the mixed concentration list defined above - not the total concentration list which includes the end members
master_cr.input["fraction_lst"] = _________
master_cr.input["species_one"] = "Ni"
master_cr.input["species_two"] = "Cr"
# Execute the calculation

Collect structures

At this point ten structures have been generated five for each phase. Now each of these structures has to be minimized to apply the quasi-harmonic approximation for each of them. In the interest of time this step is skipped here and we only compile the list of structures.

# Combine the end members with the structures from the SQS masters 
# to have a total of five structures for each phase
structure_fcc_lst = [______] + master_ni.list_structures() + [______]
structure_bcc_lst = [______] + master_cr.list_structures() + [______]

len(structure_fcc_lst) == len(structure_bcc_lst)


In this section you learned:

  • how to calculate free energies using the harmonic and the quasi-harmonic approximation.

  • how hierachical workflows can be constructed in pyiron by combining multiple master jobs.

  • and how to generate mixed phases with special quasi random structures.


  • Calculate the concentration dependent volume expansion for both the FCC and the BCC phase.

  • Calculate the harmonic approximation for a single SQS structure. As the mixed phase breaks the symmetry additional displacments are necessary. Phonopy automatically identifies the necessary displacements, but these calculation take more time.

  • Use the option job.server.run_mode.interactive=True for the LAMMPS jobs before they are assigned to the phonopy jobs to accelerate the calculations by using LAMMPS as a C-library rather than calling the executable by writing output and input files.