Chapter 5: Python Practice: Electronic Structure Calculations

In Chapter 1, we understood chemical bonding through classical potential models (the Madelung constant, the Morse potential, etc.); in Chapter 2, we understood band structure through the LCAO method and the tight-binding approximation; in Chapter 3, we understood the splitting of d orbitals in transition-metal complexes through crystal field theory; and in Chapter 4, we understood phase equilibria through the Gibbs free energy and the regular solution model. Each of these theories made the behavior of electrons and atoms tractable through some kind of approximation or modeling.

However, to quantitatively predict the properties of real solids and molecules, we must somehow solve the many-body problem, in which numerous electrons move while interacting with one another. The most successful practical solution to this problem is Density Functional Theory (DFT). Instead of the complicated many-electron wavefunction, DFT uses the far more tractable electron density $n(\mathbf{r})$ as its fundamental variable, making it possible to compute a material's total energy, electronic structure, and the forces acting on its atoms as an ab initio calculation that does not rely on empirical parameters.

As the final chapter of this series, we begin by learning the basic concepts of DFT. We then use the Python library ASE (Atomic Simulation Environment) to build and manipulate crystal structures, getting hands-on experience with the general workflow of structure optimization and electronic structure calculation. We also learn the idea of data-driven materials discovery using materials databases such as the Materials Project, and finally we summarize how the knowledge accumulated across Chapters 1 through 4 connects within modern computational materials chemistry powered by machine learning.

← Chapter 4 | Series Index

5.1 Fundamental Concepts of DFT (Density Functional Theory)

To exactly determine the ground state of a system containing $N$ electrons, one must solve the Schrödinger equation for the many-body wavefunction $\Psi(\mathbf{r}_1, \ldots, \mathbf{r}_N)$ defined on a $3N$-dimensional configuration space.

The groundbreaking idea for circumventing this difficulty is the pair of theorems presented by Hohenberg and Kohn in 1964, known as the Hohenberg-Kohn Theorems.

• First theorem: The ground-state electron density $n(\mathbf{r})$ of a system uniquely determines the external potential $V_{ext}(\mathbf{r})$ (up to an additive constant). Consequently, every physical quantity of the ground state can be expressed as a functional of the density $n(\mathbf{r})$.
• Second theorem: There exists a universal functional that gives the total energy $E[n]$ as a functional of the density, and the true ground-state density minimizes this functional (the variational principle).

This allows us to work with the electron density $n(\mathbf{r})$, a mere three-dimensional function, as the fundamental variable, instead of the $3N$-dimensional wavefunction. However, the Hohenberg-Kohn theorems only guarantee the "existence" of the energy functional; they do not give its explicit form. In 1965, Kohn and Sham provided a practical computational scheme by constructing a non-interacting auxiliary electron system (the Kohn-Sham system) that reproduces the same density as the true interacting electron system. This is the Kohn-Sham Equations.

The electron density is computed from the occupied Kohn-Sham orbitals as $n(\mathbf{r}) = \sum_i^{occ}|\psi_i(\mathbf{r})|^2$, but this density is needed to determine the effective potential $V_{eff}$ (in particular $V_H$ and $V_{xc}$), while $\psi_i$ — and hence $n(\mathbf{r})$ — can only be obtained once $V_{eff}$ has been solved. To resolve this circularity, one starts from some initial guess density and repeats the following cycle: (1) compute the effective potential, (2) diagonalize the Kohn-Sham equations, (3) compute a new density, and (4) return to (1) until convergence is reached. This is called a Self-Consistent Field calculation (SCF). The true exchange-correlation functional $E_{xc}[n]$ is, strictly speaking, unknown, and practical calculations use approximate functionals such as the Local Density Approximation (LDA) or the Generalized Gradient Approximation (GGA) (the PBE functional being a representative example).

Below, instead of performing an actual DFT calculation, we implement in Python a simplified toy model (a simplified model for educational purposes) that lets us intuitively understand the mechanics of the SCF cycle itself. We consider a two-electron system bound in a one-dimensional harmonic potential, and iteratively solve the Schrödinger equation by the finite-difference method under a simplified Hartree-like mean-field potential that depends on the electron density. Note that the Hartree term here is not a genuine integral but a mock term proportional to the local density. Even so, the structure of the SCF loop itself — density → effective potential → orbitals → new density — is the same as in an actual DFT calculation.

Python Implementation: A Toy Model for Self-Consistent Field (SCF) Calculation

import numpy as np
import matplotlib.pyplot as plt

def solve_1d_schrodinger(V, x):
    """Solve the one-dimensional Schrödinger equation by the finite-difference method (atomic units, hbar=m=1)"""
    N = len(x)
    dx = x[1] - x[0]
    kinetic = np.zeros((N, N))
    for i in range(N):
        kinetic[i, i] = 2.0
        if i > 0:
            kinetic[i, i - 1] = -1.0
        if i < N - 1:
            kinetic[i, i + 1] = -1.0
    kinetic /= (2 * dx**2)

    H = kinetic + np.diag(V)
    eigenvalues, eigenvectors = np.linalg.eigh(H)
    eigenvectors = eigenvectors / np.sqrt(dx)  # normalization: integral |psi|^2 dx = 1
    return eigenvalues, eigenvectors


def self_consistent_field(x, V_ext, U_hartree=0.3, n_electrons=2, max_iter=80, tol=1e-8, mix=0.3):
    """
    A toy model of a simplified self-consistent field (SCF) calculation.

    Just as in actual Kohn-Sham DFT, this repeats the cycle of solving the
    one-electron equation under a density-dependent effective potential and
    then updating the effective potential with the resulting density, until
    convergence is reached. Note that the Hartree term here is not a genuine
    integral but a simplified mock term proportional to the local density.
    """
    N = len(x)
    dx = x[1] - x[0]
    n = np.zeros(N)
    energies_history = []

    for iteration in range(max_iter):
        V_eff = V_ext + U_hartree * n
        eigenvalues, eigenvectors = solve_1d_schrodinger(V_eff, x)

        n_new = np.zeros(N)
        remaining = n_electrons
        idx = 0
        while remaining > 0:
            occ = min(2, remaining)
            n_new += occ * eigenvectors[:, idx]**2
            remaining -= occ
            idx += 1

        n_next = mix * n_new + (1 - mix) * n  # simple mixing to suppress SCF oscillations

        E_band = _occupied_energy_sum(eigenvalues, n_electrons)
        E_hartree_dc = 0.5 * U_hartree * np.sum(n_new**2) * dx  # double-counting correction
        energies_history.append(E_band - E_hartree_dc)

        if np.max(np.abs(n_next - n)) < tol:
            n = n_next
            break
        n = n_next

    return n, eigenvalues, energies_history


def _occupied_energy_sum(eigenvalues, n_electrons):
    total, remaining, idx = 0.0, n_electrons, 0
    while remaining > 0:
        occ = min(2, remaining)
        total += occ * eigenvalues[idx]
        remaining -= occ
        idx += 1
    return total


# Computational grid: a two-electron system in a one-dimensional harmonic potential (mimicking nuclear confinement)
L, N = 6.0, 300
x = np.linspace(-L, L, N)
V_ext = 0.5 * x**2

n_final, eigenvalues, energies_history = self_consistent_field(x, V_ext, U_hartree=0.3, n_electrons=2)
n_noninteracting, ev0, _ = self_consistent_field(x, V_ext, U_hartree=0.0, n_electrons=2, max_iter=1)

print(f"Number of iterations to SCF convergence: {len(energies_history)}")
print(f"Total energy after convergence: {energies_history[-1]:.6f} (atomic units)")
print(f"Eigenvalues of the lowest three Kohn-Sham orbitals: {np.round(eigenvalues[:3], 6)}")
print(f"Normalization check of the electron density: integral n(x) dx = {np.trapezoid(n_final, x):.4f} (should be close to 2 for a two-electron system)")
print(f"\nReference: lowest-orbital energy with no electron-electron interaction (U=0) = {ev0[0]:.6f}")
print("Exact solution for the harmonic oscillator (ground state, m=omega=hbar=1): E_0 = 0.5")
print(f"-> The interaction (U=0.3) pushes up the lowest orbital energy by {eigenvalues[0]-ev0[0]:.4f}")
print("  (the term mimicking Coulomb repulsion between electrons causes the electrons to occupy a more spread-out state)")

fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))
axes[0].plot(range(1, len(energies_history) + 1), energies_history, 'o-', color='#f5576c')
axes[0].set_xlabel('SCF iteration', fontsize=12)
axes[0].set_ylabel('Total energy (atomic units)', fontsize=12)
axes[0].set_title('Convergence of Total Energy During the SCF Cycle', fontsize=13, fontweight='bold')
axes[0].grid(True, alpha=0.3)

axes[1].plot(x, n_final, color='#f5576c', linewidth=2.5, label='Self-consistent density (U=0.3)')
axes[1].plot(x, n_noninteracting, '--', color='#2c3e50', linewidth=2, label='No interaction (U=0)')
axes[1].set_xlabel('Position x (atomic units)', fontsize=12)
axes[1].set_ylabel('Electron density n(x)', fontsize=12)
axes[1].set_title('Electron Density Distribution After Convergence', fontsize=13, fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('scf_toy_model.png', dpi=300)
plt.show()

Result: The SCF cycle converges after 48 iterations (using simple mixing with a mixing parameter of 0.3), and the electron density is correctly normalized to the two electrons (integral value 2.0000). Without interaction ($U=0$), the lowest orbital energy is 0.499950, in good agreement with the exact harmonic-oscillator solution $E_0=0.5$. Introducing the term mimicking Coulomb repulsion between electrons ($U=0.3$) pushes up the lowest orbital energy by 0.2284, and the electron density distribution also spreads out further. This qualitatively corresponds to the effect that the Hartree and exchange-correlation terms have on orbital energies and density distributions in an actual DFT calculation.

5.2 Using ASE (Atomic Simulation Environment)

ASE (Atomic Simulation Environment) is an open-source library that provides a unified Python API for building and manipulating atomic, molecular, and crystal structures, and for interfacing with various calculation engines (EMT, GPAW, VASP, etc.). At the core of ASE is the Atoms object, which collectively holds structural information such as the types of atoms, their position coordinates, whether periodic boundary conditions are applied, and the unit cell (the lattice vectors).

Building a crystal structure by hand from a unit cell requires knowledge of symmetry operations, but ASE's ase.build.bulk() function lets you build representative crystal structures in a single line simply by specifying the element symbol and the name of the crystal structure ('fcc', 'bcc', 'diamond', 'rocksalt', etc.).

Python Implementation: Building Various Crystal Structures with ASE

import numpy as np
from ase import Atoms
from ase.build import bulk

# Build a bulk crystal of copper (FCC structure)
cu = bulk('Cu', 'fcc', a=3.615)
print("=== Unit cell of copper (FCC) ===")
print(f"Chemical formula: {cu.get_chemical_formula()}")
print(f"Lattice constants: {cu.cell.cellpar()}")
print(f"Number of atoms: {len(cu)}")
print(f"Volume: {cu.get_volume():.4f} A^3")

# Silicon (diamond structure)
si = bulk('Si', 'diamond', a=5.431)
print("\n=== Unit cell of silicon (diamond structure) ===")
print(f"Chemical formula: {si.get_chemical_formula()}")
print(f"Lattice constants: {si.cell.cellpar()}")
print(f"Number of atoms: {len(si)}")

# NaCl structure (rock salt, containing two types of atoms)
nacl = bulk('NaCl', 'rocksalt', a=5.64)
print("\n=== Unit cell of NaCl (rock salt structure) ===")
print(f"Chemical formula: {nacl.get_chemical_formula()}")
print(f"Atomic positions (fractional coordinates):\n{nacl.get_scaled_positions()}")

Result: The FCC unit cell of copper (the default of bulk() is the primitive cell, containing one atom) has a volume of 11.8104 ų, with the angle between lattice vectors being 60° (a rhombohedral description). The diamond structure of silicon is built as a unit cell containing two atoms (Si2), and the NaCl structure is built as a rock-salt structure with Na and Cl placed at (0, 0, 0) and (0.5, 0.5, 0.5) respectively. We can see that the ionic crystal (NaCl-type) learned in Chapter 1 and the highly covalent diamond structure learned in Chapter 2 can both be handled uniformly with the same bulk() function.

In actual calculations, one often needs to build large supercells by repeating the unit cell, or to directly manipulate atomic positions to model lattice defects or surface structures. The repeat() method and the positions attribute of the Atoms object let you carry out these operations with the same feel as working with NumPy arrays.

Python Implementation: Creating a Supercell and Manipulating Atomic Positions

import numpy as np
from ase.build import bulk

cu = bulk('Cu', 'fcc', a=3.615)

# Create a 2x2x2 supercell
cu_super = cu.repeat((2, 2, 2))
print("=== 2x2x2 supercell of copper ===")
print(f"Number of atoms: {len(cu_super)} ({len(cu_super)//len(cu)} times the unit cell)")

# Direct manipulation of atomic positions: displace one Cu atom to model a lattice defect
cu_defect = cu_super.copy()
displacement = np.array([0.1, 0.0, 0.0])
cu_defect.positions[0] += displacement
disp_check = np.linalg.norm(cu_defect.positions[0] - cu_super.positions[0])
print(f"Displacement of the moved atom: {disp_check:.4f} A")

# Filtering by atomic number/chemical symbol (example: extracting only the Na atoms from the NaCl structure)
nacl_super = bulk('NaCl', 'rocksalt', a=5.64).repeat((2, 2, 2))
na_indices = [atom.index for atom in nacl_super if atom.symbol == 'Na']
print(f"\nNumber of Na atoms in the NaCl 2x2x2 supercell: {len(na_indices)} / total atoms: {len(nacl_super)}")

Result: repeat((2, 2, 2)) gives a supercell whose number of atoms is 8 times larger (from 1 atom to 8 atoms). Because the positions attribute is itself a NumPy array, structures can easily be deformed simply by selecting specific atoms and adding to their coordinates. An Atoms object can also be iterated atom by atom with a Python for loop, and filtering by chemical symbol can be written intuitively as well. By attaching a calculation engine (a Calculator) to a structure object built in this way, as described in the next section, actual calculations such as energy, forces, and structure optimization become possible.

5.3 Crystal Structure Optimization and the Electronic Structure Calculation Workflow

The typical workflow for electronic structure calculation and structure optimization in ASE can be organized into the following three steps.

1. Build an Atoms object (as described in the previous section)
2. Attach a calculation engine (a Calculator) with atoms.calc = Calculator()
3. Relax the atomic positions and cell shape with an optimization algorithm (such as BFGS), and obtain physical quantities with get_potential_energy() and get_forces()

The key here is the choice of Calculator in step 2. ASE is designed so that you can switch between a common interface, from fast empirical potentials such as EMT (Effective Medium Theory) to full-fledged DFT codes such as GPAW, VASP, and Quantum ESPRESSO. EMT is an extremely fast empirical potential parameterized for a handful of metals such as Ni, Cu, Pd, Ag, Pt, and Al; it is useful for learning the workflow itself and for grasping rough trends, but it does not have the quantum mechanical accuracy of DFT. Below, we first use EMT to actually carry out two kinds of optimization: lattice-constant optimization and structural relaxation of atomic positions.

Python Implementation: Equation-of-State (EOS) Fitting and BFGS Structure Optimization

import numpy as np
import matplotlib.pyplot as plt
from ase.build import bulk
from ase.calculators.emt import EMT
from ase.optimize import BFGS
from ase.eos import EquationOfState

a0 = 3.615  # experimental lattice constant of copper (A)

# Sweep the lattice constant, evaluate the energy with the EMT calculator, and use equation-of-state (EOS)
# fitting to find the most stable volume, lattice constant, and bulk modulus
volumes, energies = [], []
for s in np.linspace(0.94, 1.06, 9):
    atoms = bulk('Cu', 'fcc', a=a0 * s)
    atoms.calc = EMT()
    volumes.append(atoms.get_volume())
    energies.append(atoms.get_potential_energy())

eos = EquationOfState(volumes, energies, eos='birchmurnaghan')
v0, E0, B = eos.fit()
a_opt = (4 * v0) ** (1 / 3)          # FCC (1 atom/unit cell): V = a^3/4
B_GPa = B * 160.2176634              # conversion factor from eV/A^3 to GPa

print("=== Structure optimization via equation-of-state (EOS) fitting ===")
print(f"Optimal lattice constant from EMT: {a_opt:.4f} A (experimental value: {a0} A)")
print(f"Bulk modulus from EMT: {B_GPa:.2f} GPa (experimental value: about 140 GPa)")
print(f"Total energy at the most stable volume: {E0:.6f} eV/atom")

# Direct optimization of atomic positions via the BFGS method (a quasi-Newton method)
# Artificially displace one atom of a 2x2x2 supercell and relax it until the forces vanish
cu = bulk('Cu', 'fcc', a=a_opt).repeat((2, 2, 2))
cu.calc = EMT()
cu.positions[0] += [0.15, 0.10, 0.05]  # artificial displacement mimicking a lattice defect

print(f"\n=== BFGS structure optimization ===")
print(f"Maximum residual force before optimization: {np.max(np.abs(cu.get_forces())):.4f} eV/A")

opt = BFGS(cu, logfile=None)
opt.run(fmax=0.01)

print(f"Maximum residual force after optimization: {np.max(np.abs(cu.get_forces())):.6f} eV/A")
print(f"Number of steps required for optimization: {opt.get_number_of_steps()}")

# Visualize the E-V curve
plt.figure(figsize=(9, 6))
plt.plot(volumes, energies, 'o', color='#f093fb', markersize=9, label='EMT calculated points')
plt.axvline(v0, color='#f5576c', linestyle='--', label=f'Optimal volume V0 = {v0:.3f} A^3/atom')
plt.xlabel('Volume (A^3/atom)', fontsize=12)
plt.ylabel('Total energy (eV/atom)', fontsize=12)
plt.title('Equation-of-State (EOS) Fitting for Copper (FCC)', fontsize=14, fontweight='bold')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('eos_cu_final.png', dpi=300)
plt.show()

Result: The optimal lattice constant of copper from EMT is 3.5900 Å (a 0.7% error from the experimental value of 3.615 Å), and the bulk modulus is 133.01 GPa (about a 5% error from the experimental value of roughly 140 GPa), showing good agreement for a simple empirical potential. Meanwhile, in the BFGS structure optimization, the residual force on the artificially displaced atom converges from 1.28 eV/Å to 0.0038 eV/Å in 12 steps, confirming that the BFGS method (a quasi-Newton method) can efficiently search for an energy minimum by successively approximating the Hessian matrix from force information.

To make quantitative predictions in actual materials science research, one needs a full-fledged DFT code rather than an empirical potential such as EMT. ASE can interface with numerous DFT codes, including GPAW (a DFT code based on the PAW method, the Projector Augmented-Wave method). GPAW is a flexible code that can perform calculations using plane-wave, real-space grid, or LCAO bases, and it is well suited to band structure and density-of-states calculations for crystals under periodic boundary conditions. Below is an example implementation that uses GPAW to compute the band structure and density of states of copper. Running the code in this section requires a separate installation of GPAW (after conda install -c conda-forge gpaw or pip install gpaw, you also need to set up the pseudopotential/PAW datasets). Once GPAW itself is installed, this code can be run as is.

Python Implementation (requires GPAW installation): Band Structure and Density of States of Copper

# Running this code requires a separate installation of GPAW:
#   conda install -c conda-forge gpaw   (recommended)
#   or after pip install gpaw, fetch the PAW datasets with `gpaw install-data`
# GPAW depends on external libraries such as mpi4py and libxc, so an Anaconda environment is recommended.

from ase.build import bulk
from gpaw import GPAW, PW, FermiDirac

# Self-consistent (SCF) ground-state calculation for bulk copper (FCC)
atoms = bulk('Cu', 'fcc', a=3.615)
calc = GPAW(
    mode=PW(400),                  # plane-wave basis, cutoff energy 400 eV
    kpts={'size': (8, 8, 8), 'gamma': True},  # k-point sampling (discretization of the Brillouin zone)
    xc='PBE',                      # exchange-correlation functional: GGA (PBE)
    occupations=FermiDirac(0.05),  # finite-temperature Fermi-Dirac smearing for the metal
    txt='cu_scf.txt'
)
atoms.calc = calc
E_total = atoms.get_potential_energy()  # run the SCF calculation
print(f"Total energy of copper: {E_total:.4f} eV")
calc.write('cu_gs.gpw')  # save the ground-state wavefunctions

# Band structure calculation: a non-self-consistent (non-SCF) calculation along high-symmetry points
from ase.dft.kpoints import bandpath

path = bandpath('GXWLGK', atoms.cell, npoints=100)  # high-symmetry path in the first Brillouin zone
calc_bands = GPAW('cu_gs.gpw').fixed_density(
    kpts=path, symmetry='off', txt='cu_bands.txt'
)
bs = calc_bands.band_structure()
bs.write('cu_bandstructure.json')
# bs.plot(emin=-10, emax=10, filename='cu_bands.png')  # draw the band structure diagram

# Compute the density of states (DOS)
energies, dos = calc.get_dos(spin=0, npts=500, width=0.1)
print(f"Fermi level: {calc.get_fermi_level():.4f} eV")

This code computes, from first principles, for a real three-dimensional crystal (copper) the same kind of simple one-dimensional band structure and density of states that we handled with the tight-binding model in Chapter 2. The density of the k-point sampling specified with kpts, and the plane-wave cutoff energy specified with mode=PW(400), are important convergence parameters that govern computational accuracy and cost, and in practice a convergence test — running the calculation with several different values and confirming that the results no longer change — is indispensable.

5.4 Data-Driven Materials Discovery with the Materials Project API

Materials Project is one of the world's largest materials databases, openly publishing the crystal structures, energies, band gaps, elastic constants, and other properties of over 150,000 inorganic materials obtained from DFT calculations. Using mp-api (the Python client library), you can search for materials that meet specified conditions and retrieve their DFT calculation results. This makes possible data-driven materials discovery, in which candidate materials are narrowed down from existing calculated data before performing experiments or costly DFT calculations.

Using the Materials Project API requires registering for a free API key and an internet connection, so first we experience a similar idea with a small local dataset. Based on the theory of the Shockley-Queisser Limit, the optimal band gap for a light-absorbing layer usable in solar cells is considered to be roughly 1.0-1.6 eV for a single junction. Below, we use experimentally measured band gap values for representative semiconductors and insulators to perform screening (narrowing down) under this condition.

Python Implementation: Material Screening by Band Gap (Local Dataset)

import numpy as np
import pandas as pd

# Experimental band gap values for representative semiconductors and insulators (eV)
# Source: representative values such as those in Kittel, "Introduction to Solid State Physics," 8th Ed.
materials = pd.DataFrame({
    'formula':      ['Ge', 'Si', 'InP', 'GaAs', 'CdTe', 'CdS', 'SiC', 'ZnO', 'GaN', 'Diamond'],
    'band_gap_eV':  [0.67, 1.11, 1.34, 1.42, 1.50, 2.42, 2.36, 3.37, 3.40, 5.47],
    # Pauling electronegativity difference |chi_A - chi_B| (reusing the electronegativity concept introduced in Chapter 1)
    'delta_chi':    [0.00, 0.00, 0.42, 0.10, 0.54, 0.79, 0.65, 1.79, 1.23, 0.00],
})

# Screen for a band gap suited to a single-junction solar cell (roughly 1.0-1.6 eV), based on the Shockley-Queisser limit
target_min, target_max = 1.0, 1.6
candidates = materials[(materials['band_gap_eV'] >= target_min) & (materials['band_gap_eV'] <= target_max)]

print("=== Candidate materials with a band gap (1.0-1.6 eV) suited to solar cell applications ===")
print(candidates[['formula', 'band_gap_eV']].to_string(index=False))
print(f"\n{len(candidates)} out of {len(materials)} materials were extracted as candidates")

Result: Of the 10 materials, four — Si (1.11 eV), InP (1.34 eV), GaAs (1.42 eV), and CdTe (1.50 eV) — fall within the optimal range of the Shockley-Queisser limit. Indeed, GaAs and CdTe are used in practice as high-efficiency solar cell materials, and Si is the most widely deployed solar cell material. This shows that even simple condition-based filtering can be an effective first stage of narrowing down in materials discovery.

Applying this idea to a large-scale database such as the Materials Project allows instantaneous screening from hundreds of thousands of candidates. Actual search code using mp-api looks like the following (running this code requires registering for a free API key and an internet connection; you can register at https://materialsproject.org/api).

Python Implementation (requires an API key): Searching for Materials via the Materials Project API

# Running this code requires a free API key registration:
#   1. Register an account at https://materialsproject.org/api and obtain an API key
#   2. pip install mp-api
#   3. Set the environment variable MP_API_KEY, or specify the key directly in the code

from mp_api.client import MPRester

with MPRester("YOUR_API_KEY") as mpr:
    # Example: search for materials with a band gap of 1.0-1.6 eV, a negative formation
    # energy (thermodynamically stable), excluding oxides
    docs = mpr.materials.summary.search(
        band_gap=(1.0, 1.6),
        formation_energy=(None, 0.0),
        fields=["material_id", "formula_pretty", "band_gap",
                "formation_energy_per_atom", "is_stable"]
    )

    print(f"Number of materials matching the criteria: {len(docs)}")
    for doc in docs[:10]:
        print(f"{doc.material_id}: {doc.formula_pretty}, "
              f"Eg={doc.band_gap:.3f} eV, "
              f"Ef={doc.formation_energy_per_atom:.3f} eV/atom, "
              f"stable phase={doc.is_stable}")

We see that exactly the same narrowing-down logic (restricting the range of the band gap) used with the local mini-dataset can be executed on the Materials Project in combination with a much richer set of conditions, such as formation_energy (an indicator of thermodynamic stability related to the Gibbs free energy learned in Chapter 4) and is_stable (a phase-stability determination based on convex-hull analysis).

5.5 Integrating Computational Materials Chemistry with Machine Learning

Looking back over the four preceding sections, the computational approaches in materials chemistry build up as follows. DFT (Section 5.1) provides the electronic states and total energies of individual materials from first principles; ASE (Sections 5.2-5.3) provides the foundation for structure building and workflow to carry out those calculations; and the Materials Project (Section 5.4) functions as a massive database accumulating the results of DFT calculations performed around the world. Machine Learning is the final piece: it learns statistical patterns from this accumulated data to rapidly predict material properties without running new DFT calculations.

The performance of a machine learning model is heavily influenced by the quality of the descriptors (also called features) given as input. The Pauling electronegativity difference learned in Chapter 1 can serve as a descriptor characterizing the ionicity of a bond, but it does not necessarily explain the magnitude of the band gap on its own. Below, we compare how prediction performance changes when a descriptor representing the atomic size of the constituent elements (the row, or period number, in the periodic table) is combined with the electronegativity difference.

Python Implementation: Building and Evaluating a Band Gap Prediction Model Using Descriptors

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import LeaveOneOut
from sklearn.metrics import mean_absolute_error, r2_score

# Experimental band gap values (eV) and descriptors for representative semiconductors and insulators
materials = pd.DataFrame({
    'formula':      ['Ge', 'Si', 'InP', 'GaAs', 'CdTe', 'CdS', 'SiC', 'ZnO', 'GaN', 'Diamond'],
    'band_gap_eV':  [0.67, 1.11, 1.34, 1.42, 1.50, 2.42, 2.36, 3.37, 3.40, 5.47],
    'delta_chi':    [0.00, 0.00, 0.42, 0.10, 0.54, 0.79, 0.65, 1.79, 1.23, 0.00],
    # Average period number of the constituent elements (a descriptor capturing the tendency that elements with
    # larger atomic numbers have larger atomic radii and greater orbital overlap, generally leading to a smaller band gap)
    'avg_period':   [4.0, 3.0, 4.0, 4.0, 5.0, 4.0, 2.5, 3.0, 3.0, 2.0],
})

def evaluate_model(feature_cols, label):
    X = materials[feature_cols].values
    y = materials['band_gap_eV'].values
    model = LinearRegression().fit(X, y)
    r2_train = r2_score(y, model.predict(X))

    # With so little data, use Leave-One-Out cross-validation to rigorously evaluate generalization performance
    loo = LeaveOneOut()
    y_pred_loo = np.zeros_like(y)
    for train_idx, test_idx in loo.split(X):
        m = LinearRegression().fit(X[train_idx], y[train_idx])
        y_pred_loo[test_idx] = m.predict(X[test_idx])
    mae_loo = mean_absolute_error(y, y_pred_loo)
    r2_loo = r2_score(y, y_pred_loo)
    print(f"{label:28s}: R^2(train)={r2_train:.3f}  MAE(LOO)={mae_loo:.3f} eV  R^2(LOO)={r2_loo:.3f}")
    return y, y_pred_loo, model

print("=== Comparison of prediction performance by descriptor (feature) ===")
evaluate_model(['delta_chi'], 'delta_chi only (Chapter 1 metric)')
evaluate_model(['avg_period'], 'avg_period only')
y, y_pred_loo, model = evaluate_model(['delta_chi', 'avg_period'], 'delta_chi + avg_period')

print(f"\nCoefficients of the final model: delta_chi={model.coef_[0]:.3f}, avg_period={model.coef_[1]:.3f}, intercept={model.intercept_:.3f}")
print("\nCaution: this dataset contains only 10 materials, and this model is merely an educational example")
print("illustrating the idea of feature engineering. A single electronegativity difference alone has")
print("little explanatory power, and adding the atomic size (period number) improves the fit to the")
print("training data. However, as the LOO cross-validation results show, generalization performance is")
print("limited with only 10 data points, and actual materials discovery requires a large-scale database")
print("such as the Materials Project along with a richer set of descriptors.")

plt.figure(figsize=(7, 7))
plt.scatter(y, y_pred_loo, s=90, color='#f5576c', zorder=3)
lims = [0, 6]
plt.plot(lims, lims, '--', color='gray', zorder=1)
for f, yt, yp in zip(materials['formula'], y, y_pred_loo):
    plt.annotate(f, (yt, yp), textcoords="offset points", xytext=(6, 4), fontsize=9)
plt.xlabel('Measured band gap (eV)', fontsize=12)
plt.ylabel('LOO-predicted band gap (eV)', fontsize=12)
plt.title('Band Gap Prediction via delta_chi + avg_period (LOO Cross-Validation)', fontsize=13, fontweight='bold')
plt.xlim(lims); plt.ylim(lims)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('ml_bandgap_final.png', dpi=300)
plt.show()

Result: Using only the electronegativity difference delta_chi as the explanatory variable gives a low Coefficient of Determination (R²) of just 0.099 (because there are cases like Diamond, where delta_chi=0 yet the band gap is a very large 5.47 eV). Using avg_period, representing atomic size, alone improves this to R²=0.491, and combining both raises the fit to the training data to R²=0.557. However, when generalization performance is rigorously evaluated with Leave-One-Out Cross-Validation (LOO-CV), the mean absolute error is around 1 eV, and the coefficient of determination can even become negative — meaning that, with only 10 materials, this cannot be called a statistically reliable prediction model. This result should be taken seriously at face value, and it underscores why actual materials discovery requires data on the scale of tens of thousands of entries, as in the Materials Project, together with the rich sets of descriptors provided by dedicated libraries such as matminer.

Let us now organize the connections among everything learned in this chapter, and in the series as a whole. The chemical bonding theory of Chapter 1 explained, using classical models, the "origin of the energy" behind why atoms bond together. The molecular orbital theory and band theory of Chapter 2 described, within the common framework of the LCAO method, what quantum mechanical states the electrons responsible for that bonding occupy. The crystal field theory of Chapter 3 focused on the special but important class of transition-metal compounds and showed how the details of electronic states — the splitting of d orbitals — govern magnetism, color, and reactivity. The thermodynamics of Chapter 4 dealt not with individual bonds or electronic states but with the larger-scale behavior of which phase is stable for a macroscopic system composed of many atoms. And the DFT of this chapter provides the unified foundation — computing the quantum mechanical ground state of electrons from first principles — underlying all of these phenomenological models (bonding energy, orbital energy, d-orbital splitting, formation energy). The formation energies used in building a CALPHAD thermodynamic database come not only from experimental data but also from DFT calculations, and the training data for machine learning models likewise derives from the accumulation of DFT calculations (the Materials Project). In this way, the five chapters are not an unrelated collection of separate topics but form a single connected narrative spanning "why bonds form," to "how a stable structure or phase is determined," and finally to "how to compute this from first principles and search at scale."

Exercises

import numpy as np
from ase.build import bulk

al = bulk('Al', 'fcc', a=4.05)
fe = bulk('Fe', 'bcc', a=2.87)

print("=== Al (FCC) ===")
print(f"Chemical formula: {al.get_chemical_formula()}, atoms (unit cell): {len(al)}")
print(f"Volume (unit cell): {al.get_volume():.4f} A^3")

print("\n=== Fe (BCC) ===")
print(f"Chemical formula: {fe.get_chemical_formula()}, atoms (unit cell): {len(fe)}")
print(f"Volume (unit cell): {fe.get_volume():.4f} A^3")

# Packing fraction: an index of how much of space can be filled under the hard-sphere model
# The theoretical value is 0.7405 for FCC and 0.6802 for BCC
r_al = 4.05 * np.sqrt(2) / 4  # half the nearest-neighbor distance for FCC = a*sqrt(2)/4
r_fe = 2.87 * np.sqrt(3) / 4  # half the nearest-neighbor distance for BCC = a*sqrt(3)/4

packing_fcc = (len(al) * (4/3) * np.pi * r_al**3) / al.get_volume()
packing_bcc = (len(fe) * (4/3) * np.pi * r_fe**3) / fe.get_volume()
print(f"\nPacking fraction of Al (FCC): {packing_fcc:.4f} (theoretical value 0.7405)")
print(f"Packing fraction of Fe (BCC): {packing_bcc:.4f} (theoretical value 0.6802)")
print("\nFCC has a higher packing fraction than BCC, which is one of the reasons FCC metals")
print("(Al, Cu, Ni, etc.) generally exhibit superior ductility.")

import numpy as np
from ase.build import bulk
from ase.calculators.emt import EMT
from ase.eos import EquationOfState

a0 = 4.05  # experimental lattice constant of Al
volumes, energies = [], []
for s in np.linspace(0.94, 1.06, 9):
    atoms = bulk('Al', 'fcc', a=a0 * s)
    atoms.calc = EMT()
    volumes.append(atoms.get_volume())
    energies.append(atoms.get_potential_energy())

eos = EquationOfState(volumes, energies, eos='birchmurnaghan')
v0, E0, B = eos.fit()
a_opt = (4 * v0) ** (1 / 3)
B_GPa = B * 160.2176634

print(f"Optimal lattice constant of Al from EMT: {a_opt:.4f} A (experimental value: 4.05 A)")
print(f"Bulk modulus from EMT: {B_GPa:.2f} GPa (experimental value: about 76 GPa)")
error_a = abs(a_opt - 4.05) / 4.05 * 100
error_B = abs(B_GPa - 76) / 76 * 100
print(f"Error in lattice constant: {error_a:.2f}%")
print(f"Error in bulk modulus: {error_B:.1f}%")

import numpy as np
import pandas as pd

materials = pd.DataFrame({
    'formula':      ['Ge', 'Si', 'InP', 'GaAs', 'CdTe', 'CdS', 'SiC', 'ZnO', 'GaN', 'Diamond'],
    'band_gap_eV':  [0.67, 1.11, 1.34, 1.42, 1.50, 2.42, 2.36, 3.37, 3.40, 5.47],
    'delta_chi':    [0.00, 0.00, 0.42, 0.10, 0.54, 0.79, 0.65, 1.79, 1.23, 0.00],
    'avg_period':   [4.0, 3.0, 4.0, 4.0, 5.0, 4.0, 2.5, 3.0, 3.0, 2.0],
})

# Screen for wide-band-gap materials (band_gap > 3.0 eV) applicable to UV LEDs and similar applications
wide_gap = materials[materials['band_gap_eV'] > 3.0].sort_values('band_gap_eV', ascending=False)
print("=== Candidate wide-band-gap materials (band_gap > 3.0 eV) ===")
print(wide_gap[['formula', 'band_gap_eV']].to_string(index=False))

# Check the correlation coefficient between each descriptor and the band gap
corr_chi = np.corrcoef(materials['delta_chi'], materials['band_gap_eV'])[0, 1]
corr_period = np.corrcoef(materials['avg_period'], materials['band_gap_eV'])[0, 1]
print(f"\nCorrelation coefficient of delta_chi with the band gap: {corr_chi:.3f}")
print(f"Correlation coefficient of avg_period with the band gap: {corr_period:.3f}")
print("\navg_period shows the stronger correlation, suggesting that atomic size (orbital overlap)")
print("is the descriptor that more strongly governs the magnitude of the band gap. However, with")
print("only 10 materials, no statistically definitive conclusion can be drawn, and verification with")
print("a large-scale database such as the Materials Project is desirable.")

References

Checking Your Learning Objectives

← Chapter 4 | Series Index