Exercise 3: Density functional theory

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:

  • Density functional theory calculation,

  • Calculate the equilibrium lattice constant

  • Compare the results to interatomic potentials.

Reminder

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 import ____
# Create a Project object instance for a project named dft
pr = ____("dft")
# Create a cubic aluminum fcc structure
al_fcc = pr.create_ase_bulk(____, _____=True)
# Confirm the final structure has 4 atoms by calculating 
# the length of the structure object
____(al_fcc) == 4

Density Functional Theory

To provide a brief introduction to density functional theory (DFT) we calculate the equilibrium lattice constant for a cubic fcc aluminium structure. Besides the pseudo potential which defines the electron electron interaction, the DFT precision is dominated by the convergence parameters, namely the plane wave energy cutoff and the kpoint mesh. Both can be set in pyiron using the corresponding properties encut and kpoint_mesh.

# Create a DFT job with the S/PHI/nX quantum engine named spx
job_dft_template = pr.create_job(
    job_type=pr.job_type.Sphinx,
    job_name=________
)
# Print the default DFT convergence parameters energy cutoff
print(
    job_dft_template.encut, 
    job_dft_template.kpoint_mesh
)
# Increase the energy cutoff to 400eV and increase the kpoint mesh to 5x5x5
job_dft_template.set_encut(_______)
job_dft_template.set_kpoints([__, __, __])

Energy Volume Curve

The discretisation of the plane waves on the fourier mesh and the kpoint mesh in the brillouin zone cause fluctuations in the energy surface. Therefore a minimization like we used it for the interatomic potentials in the previous section is insufficient. Instead we calculate the energy for various volumes around the equlibrium volume and then determine the equilibrium volume by interpolating the minimum between these volumes.

# Import the numerical library numpy
import numpy as np
# Create 7 strains ranging from -5% (-0.05) to +5% (0.05)
strain_lst = np.linspace(_____, ______, _____)
strain_lst
# Copy the cubic fcc aluminium supercell
al_fcc_copy = _________.copy()
# Apply of -5% to the copy of the aluminium
# supercell and compare the volume 
al_fcc_copy.apply_strain(_____)
al_fcc.get_volume(), al_fcc_copy.get_volume()
# Iterate over the list of strains 
for _______ in strain_lst: 
    job_strain = job_dft_template.copy_to(
        # Define the job name based on the current strain 
        new_job_name="spx_" + str(1 + _____).replace(".", "_"), 
        new_database_entry=False
    )
    
    # Copy the cubic fcc aluminium supercell
    al_fcc_copy = ________.copy()
    
    # Apply the strain of using the for loop
    al_fcc_copy.apply_strain(_____)
    
    # Set the strained structure to the job
    job_strain.structure = al_fcc_copy 
    
    # Execute the job
    job_strain.run()
# Check the status of the calculation in the job_table
pr.job_table()

Reminder

We again use the pyiron table object to collect the simulation results and aggregate them in a pandas DataFrame.

# Create a pyiron table object
table = pr.create_table()
# Implement a filter function, which returns true 
# for finished jobs and jobs with "spx_" in the job_name
def filter_jobs(job):
    return job.status == _______ and _____ in job.job_name
# Many commonly used functions are already available using tab completion
# We select the get_volume function and the get_energy_tot function
table.add.________
table.add._________
# Assign the filter function defined above
table.filter_function = _____
# Execute the pyiron table just like a pyiron job object
table._____
# Return a pandas DataFrame with the collected results 
df_res = table.get_dataframe()
df_res

Visualise the energy volume curve

We again use the matplotlib library to visualise the calculated energy volume curve and calculate the equilibrium volume by fitting a second order polynomial and calculate the roots of the derivative.

# Fit a second order polynomial to the energy volume curve 
fit = np.polyfit(_____.volume, _____.energy_tot, ___)
# Get the polynomial of the fit 
fit_poly = np.poly1d(fit)

# Calculate the roots of the derivate of the polynomial
vol_roots = np.polyder(fit_poly).roots

# Select the root within the volume range which is
# smaller (<) than maximum volume and larger (>) than 
# the minimum volume. 
vol_eq = vol_roots[
    (vol_roots ____ df_res["volume"].max()) & 
    (vol_roots ____ df_res["volume"].min())
][0]
# Import the matplotlib library for plotting. 
import matplotlib.pyplot as plt

# Plot the volume and the total energy from the DataFrame 
# with the collected results
plt.plot(_____.volume, _____.energy_tot)

# Plot the fitted equilibrium volume
plt.axvline(______, linestyle="--", color="red")
plt.xlabel("Volume")
plt.ylabel("total Energy");

pyiron Master Jobs

While managing calculations with with for loops and aggregating calculation results in pandas DataFrame is already a very scaleable concept, we have implemented master jobs which can execute multiple calculations in series or in parallel to automate common simulation tasks. The calculation of the energy volume curve is one of those examples.

# Create a DFT job with the S/PHI/nX quantum engine named "spxjob"
job_master_template = pr.create_job(
    job_type=pr.job_type.Sphinx,
    job_name=______
)
# Assign the cubic fcc aluminium structure to the template job 
job_master_template.structure = _____
# Use the job object to create the Murnaghan object, named "murn" 
murn = job_master_template.create_job(
    job_type=pr.job_type.Murnaghan, 
    job_name=_______
)
# Execute the Murnaghan object 
murn.run()
# Check the status of the calculation in the job_table
pr.job_table()
# Finally we can use the internal functionality of the master job 
# to visualise the energy volume curve
murn.plot()

Summary

In this section you learned:

  • how to calculate equilibrium bulk material properties using density functional theory.

  • how the finite plane wave energy cutoff and the finite kpoint mesh limit the precision of DFT calculation.

  • and how to use master jobs in pyiron to automate common tasks like calculating energy volume curves.

Suggestions:

  • To learn more about convergence you can plot the convergence of an individual energy over the change in plane wave energy cutoff or kpoint mesh.

  • To visualise the discretisation of the energy volume curve you can calculate the energy volume curve at a small volume range of +/- 1% for a low kpoint mesh of 3x3x3 and an energy cutoff of 300eV with 21 points.

  • To validate the interatomic potentials from the previous section we can calculate energy volume curves for those and compare the energy differences to the DFT results.