6. Local Structure Optimisation with DFT#
6.1. Overview#
Questions
How do I find the minimum energy structure?
How do I find the minimum energy lattice parameters?
How do I optimise a defect geometry?
Objectives
Use the BFGS optimiser to relax a structure
Relax a structure and the unit cell
Understand the importance of geometry relaxation for defect calculations
Key Points
Geometry optimisation finds the minimum energy atomic configuration
BFGS is the standard choice
Convergence criterion (
fmaxin ASE) is the maximum force component (typically 0.01–0.05 eV/Å for DFT)Defect supercell relaxation is essential before computing electronic or optical properties
6.2. Lecture Slides#
The slides for this tutorial are embedded below. 📥 Download slides (.pptx) | Open in full screen
6.3. Why Geometry Optimisation?#
Before computing any electronic property, we must find the minimum energy geometry. In the context of quantum optics, the emission energy of a defect depends sensitively on the relaxed geometry of both ground and excited defect states.
6.3.1. Example: optimising FCC Al#
To begin with, let’s import the libraries we need (including GPAW and ASE), and create a copper FCC structure. We rattle the atoms to take them a little out of equilibrium so that we can demonstrate geometric relaxation.
import numpy as np
import ase, gpaw
from ase.build import bulk
from gpaw import GPAW, PW
atoms = ase.build.bulk('Al','fcc',a=4.7)
atoms.rattle(0.01)
Second, we attach GPAW as a calculator.
We specify a number of parameters:
xcsets the exchange-correlation functionalkptssets the Brillouin-zone samplingmodesets the basis set; in this case a 400 eV cutoff plane-wave basis is used.output_logset where the output will be printed to. If this is removed it will be printed in the notebook output.
For more information about valid parameters, see the GPAW docs.
atoms.calc = GPAW(xc='PBE', kpts=(3, 3, 3), mode="lcao", basis='dzp',txt='output_log.txt')
Timing: incl. excl.
-----------------------------------------------------------
Forces: 0.070 0.070 0.0% |
Hamiltonian: 0.352 0.004 0.0% |
Atomic: 0.282 0.003 0.0% |
XC Correction: 0.279 0.279 0.0% |
Calculate atomic Hamiltonians: 0.003 0.003 0.0% |
Communicate: 0.000 0.000 0.0% |
Initialize Hamiltonian: 0.000 0.000 0.0% |
Poisson: 0.001 0.001 0.0% |
XC 3D grid: 0.062 0.062 0.0% |
LCAO initialization: 4.447 0.465 0.0% |
LCAO eigensolver: 1.804 0.003 0.0% |
Calculate projections: 0.003 0.003 0.0% |
DenseAtomicCorrection: 0.009 0.009 0.0% |
Distribute overlap matrix: 0.001 0.001 0.0% |
Orbital Layouts: 0.018 0.018 0.0% |
Potential matrix: 1.747 1.747 0.1% |
Sum over cells: 0.023 0.023 0.0% |
LCAO to grid: 1.576 1.576 0.1% |
Set positions (LCAO WFS): 0.602 0.085 0.0% |
Basic WFS set positions: 0.001 0.001 0.0% |
Basis functions set positions: 0.000 0.000 0.0% |
P tci: 0.117 0.117 0.0% |
ST tci: 0.311 0.311 0.0% |
mktci: 0.087 0.087 0.0% |
PWDescriptor: 0.014 0.014 0.0% |
SCF-cycle: 33.333 0.520 0.0% |
Davidson: 26.538 2.276 0.1% |
Apply H: 5.315 5.197 0.2% |
HMM T: 0.118 0.118 0.0% |
Subspace diag: 6.046 0.071 0.0% |
calc_h_matrix: 5.570 0.237 0.0% |
Apply H: 5.332 5.212 0.2% |
HMM T: 0.121 0.121 0.0% |
diagonalize: 0.261 0.261 0.0% |
rotate_psi: 0.144 0.144 0.0% |
calc. matrices: 11.917 1.270 0.1% |
Apply H: 10.647 10.406 0.4% |
HMM T: 0.240 0.240 0.0% |
diagonalize: 0.721 0.721 0.0% |
rotate_psi: 0.264 0.264 0.0% |
Density: 2.878 0.001 0.0% |
Atomic density matrices: 0.043 0.043 0.0% |
Mix: 0.138 0.138 0.0% |
Multipole moments: 0.006 0.006 0.0% |
Pseudo density: 2.690 2.690 0.1% |
Symmetrize density: 0.001 0.001 0.0% |
Hamiltonian: 3.349 0.031 0.0% |
Atomic: 2.723 0.033 0.0% |
XC Correction: 2.691 2.691 0.1% |
Calculate atomic Hamiltonians: 0.031 0.031 0.0% |
Communicate: 0.000 0.000 0.0% |
Poisson: 0.006 0.006 0.0% |
XC 3D grid: 0.558 0.558 0.0% |
Orthonormalize: 0.047 0.006 0.0% |
calc_s_matrix: 0.012 0.012 0.0% |
inverse-cholesky: 0.010 0.010 0.0% |
projections: 0.014 0.014 0.0% |
rotate_psi_s: 0.005 0.005 0.0% |
Set symmetry: 0.011 0.011 0.0% |
Stress: 0.836 0.001 0.0% |
Stress tensor: 0.835 0.835 0.0% |
Other: 2482.533 2482.533 98.5% |--------------------------------------|
-----------------------------------------------------------
Total: 2521.596 100.0%
Memory usage: 202.39 MiB
Date: Thu Jul 23 13:27:04 2026
Finally, we run the calculations: these tell us what the energy and forces are acting on the materials before they are relaxed.
print(f"Initial energy: {atoms.get_potential_energy():.6f} eV")
print(f"Initial max force: {np.max(np.abs(atoms.get_forces())):.6f} eV/Å")
Initial energy: -6.407110 eV
Initial max force: 0.005702 eV/Å
Now we relax the atoms using the BFGS optimiser in ASE. The atoms will relax until the force acting on them is less than 0.001 eV/Å.
from ase.optimize import BFGS
opt = BFGS(atoms)
opt.run(fmax=0.01)
print(f"Final energy: {atoms.get_potential_energy():.6f} eV")
print(f"Final max force: {np.max(np.abs(atoms.get_forces())):.6f} eV/Å")
print(f"Final lattice constant: {atoms.cell.cellpar()[0]:.4f} Å")
Step Time Energy fmax
BFGS: 0 13:29:34 -6.453718 0.006634
Final energy: -6.407110 eV
Final max force: 0.005702 eV/Å
Final lattice constant: 2.8707 Å
The final lattice parameter has not changed because we have relaxed the atoms only. To relax the cell we need to attach a unitcellfilter to the Atoms object. hydrostatic_strain ensures that the cell shape is kept constant. Cell relaxations require a calculation of the stress tensor, and this is only implemented for plane-wave in GPAW. Therefore, we will attach a new plane-wave calculator.
This is a more intensive calculation and may take minutes to complete.
from ase.filters import UnitCellFilter
from ase.optimize import BFGS
from ase.io.trajectory import Trajectory
atoms.calc = GPAW(xc='PBE', kpts=(3, 3, 3), mode="pw",txt="pw_output_log.txt")
ucfilter = UnitCellFilter(atoms,hydrostatic_strain=True)
opt = BFGS(ucfilter)
opt.run(fmax=0.01)
print(f"Final energy: {atoms.get_potential_energy():.6f} eV")
print(f"Final max force: {np.max(np.abs(atoms.get_forces())):.6f} eV/Å")
print(f"Final lattice constant: {atoms.cell.cellpar()[0]:.4f} Å")
Timing: incl. excl.
-----------------------------------------------------------
Forces: 0.008 0.008 0.0% |
Hamiltonian: 0.044 0.001 0.0% |
Atomic: 0.035 0.000 0.0% |
XC Correction: 0.034 0.034 0.0% |
Calculate atomic Hamiltonians: 0.000 0.000 0.0% |
Communicate: 0.000 0.000 0.0% |
Initialize Hamiltonian: 0.000 0.000 0.0% |
Poisson: 0.000 0.000 0.0% |
XC 3D grid: 0.008 0.008 0.0% |
LCAO initialization: 0.604 0.060 0.0% |
LCAO eigensolver: 0.260 0.000 0.0% |
Calculate projections: 0.000 0.000 0.0% |
DenseAtomicCorrection: 0.001 0.001 0.0% |
Distribute overlap matrix: 0.000 0.000 0.0% |
Orbital Layouts: 0.002 0.002 0.0% |
Potential matrix: 0.254 0.254 0.2% |
Sum over cells: 0.002 0.002 0.0% |
LCAO to grid: 0.221 0.221 0.2% |
Set positions (LCAO WFS): 0.062 0.014 0.0% |
Basic WFS set positions: 0.000 0.000 0.0% |
Basis functions set positions: 0.000 0.000 0.0% |
P tci: 0.014 0.014 0.0% |
ST tci: 0.027 0.027 0.0% |
mktci: 0.008 0.008 0.0% |
PWDescriptor: 0.003 0.003 0.0% |
SCF-cycle: 3.676 0.059 0.0% |
Davidson: 2.877 0.244 0.2% |
Apply H: 0.571 0.557 0.5% |
HMM T: 0.013 0.013 0.0% |
Subspace diag: 0.661 0.008 0.0% |
calc_h_matrix: 0.609 0.027 0.0% |
Apply H: 0.582 0.568 0.5% |
HMM T: 0.014 0.014 0.0% |
diagonalize: 0.027 0.027 0.0% |
rotate_psi: 0.016 0.016 0.0% |
calc. matrices: 1.290 0.141 0.1% |
Apply H: 1.148 1.122 0.9% |
HMM T: 0.026 0.026 0.0% |
diagonalize: 0.083 0.083 0.1% |
rotate_psi: 0.029 0.029 0.0% |
Density: 0.363 0.000 0.0% |
Atomic density matrices: 0.005 0.005 0.0% |
Mix: 0.016 0.016 0.0% |
Multipole moments: 0.001 0.001 0.0% |
Pseudo density: 0.342 0.342 0.3% |
Symmetrize density: 0.000 0.000 0.0% |
Hamiltonian: 0.371 0.004 0.0% |
Atomic: 0.300 0.004 0.0% |
XC Correction: 0.296 0.296 0.2% |
Calculate atomic Hamiltonians: 0.003 0.003 0.0% |
Communicate: 0.000 0.000 0.0% |
Poisson: 0.001 0.001 0.0% |
XC 3D grid: 0.063 0.063 0.1% |
Orthonormalize: 0.006 0.001 0.0% |
calc_s_matrix: 0.001 0.001 0.0% |
inverse-cholesky: 0.002 0.002 0.0% |
projections: 0.002 0.002 0.0% |
rotate_psi_s: 0.001 0.001 0.0% |
Set symmetry: 0.001 0.001 0.0% |
Stress: 0.104 0.000 0.0% |
Stress tensor: 0.104 0.104 0.1% |
Other: 116.026 116.026 96.3% |--------------------------------------|
-----------------------------------------------------------
Total: 120.466 100.0%
Memory usage: 202.39 MiB
Date: Thu Jul 23 13:31:36 2026
Step Time Energy fmax
BFGS: 0 13:31:40 -6.457164 0.000132
Final energy: -6.410955 eV
Final max force: 0.000111 eV/Å
Final lattice constant: 2.8707 Å
6.3.2. Discussion#
How does this compare to the Al lattice parameter calculated in last lab? Why might there be differences?
6.4. Trajectories#
A trajectory is a sequence of Atoms objects, typically snapshots from a molecular dynamics run or a geometry optimisation. ASE has a native .traj format for storing these efficiently. Let’s run through the relaxation routine again, this time for Ti.
from ase.io.trajectory import Trajectory
atoms = ase.build.bulk('Ti','fcc',a=4.7)
atoms.rattle(0.1)
atoms.calc = GPAW(xc='PBE', kpts=(3, 3, 3), mode="pw",txt="Ti_output_log.txt")
print(f"Initial energy: {atoms.get_potential_energy():.6f} eV")
print(f"Initial max force: {np.max(np.abs(atoms.get_forces())):.6f} eV/Å")
print(f"Initial lattice constant: {atoms.cell.cellpar()[0]:.4f} Å")
ucfilter = UnitCellFilter(atoms,hydrostatic_strain=True)
opt = BFGS(ucfilter, trajectory='Ti_opt.traj', logfile='Ti_opt.log')
opt.run(fmax=0.001)
print(f"Final energy: {atoms.get_potential_energy():.6f} eV")
print(f"Final max force: {np.max(np.abs(atoms.get_forces())):.6f} eV/Å")
print(f"Final lattice constant: {atoms.cell.cellpar()[0]:.4f} Å")
Initial energy: -5.677267 eV
Initial max force: 0.000433 eV/Å
Initial lattice constant: 3.3234 Å
Timing: incl. excl.
-----------------------------------------------------------
Forces: 0.007 0.007 0.0% |
Hamiltonian: 0.039 0.001 0.0% |
Atomic: 0.029 0.000 0.0% |
XC Correction: 0.029 0.029 0.0% |
Calculate atomic Hamiltonians: 0.000 0.000 0.0% |
Communicate: 0.000 0.000 0.0% |
Initialize Hamiltonian: 0.000 0.000 0.0% |
Poisson: 0.000 0.000 0.0% |
XC 3D grid: 0.009 0.009 0.0% |
LCAO initialization: 0.179 0.012 0.0% |
LCAO eigensolver: 0.053 0.001 0.0% |
Calculate projections: 0.001 0.001 0.0% |
DenseAtomicCorrection: 0.001 0.001 0.0% |
Distribute overlap matrix: 0.000 0.000 0.0% |
Orbital Layouts: 0.002 0.002 0.0% |
Potential matrix: 0.045 0.045 0.0% |
Sum over cells: 0.002 0.002 0.0% |
LCAO to grid: 0.035 0.035 0.0% |
Set positions (LCAO WFS): 0.079 0.011 0.0% |
Basic WFS set positions: 0.000 0.000 0.0% |
Basis functions set positions: 0.000 0.000 0.0% |
P tci: 0.014 0.014 0.0% |
ST tci: 0.037 0.037 0.0% |
mktci: 0.017 0.017 0.0% |
PWDescriptor: 0.002 0.002 0.0% |
SCF-cycle: 1.692 0.028 0.0% |
Davidson: 1.223 0.160 0.0% |
Apply H: 0.199 0.193 0.0% |
HMM T: 0.006 0.006 0.0% |
Subspace diag: 0.265 0.008 0.0% |
calc_h_matrix: 0.223 0.023 0.0% |
Apply H: 0.200 0.195 0.0% |
HMM T: 0.005 0.005 0.0% |
diagonalize: 0.022 0.022 0.0% |
rotate_psi: 0.011 0.011 0.0% |
calc. matrices: 0.535 0.125 0.0% |
Apply H: 0.410 0.399 0.0% |
HMM T: 0.011 0.011 0.0% |
diagonalize: 0.046 0.046 0.0% |
rotate_psi: 0.018 0.018 0.0% |
Density: 0.133 0.000 0.0% |
Atomic density matrices: 0.005 0.005 0.0% |
Mix: 0.015 0.015 0.0% |
Multipole moments: 0.001 0.001 0.0% |
Pseudo density: 0.112 0.112 0.0% |
Symmetrize density: 0.000 0.000 0.0% |
Hamiltonian: 0.302 0.005 0.0% |
Atomic: 0.233 0.004 0.0% |
XC Correction: 0.230 0.230 0.0% |
Calculate atomic Hamiltonians: 0.004 0.004 0.0% |
Communicate: 0.000 0.000 0.0% |
Poisson: 0.001 0.001 0.0% |
XC 3D grid: 0.059 0.059 0.0% |
Orthonormalize: 0.006 0.001 0.0% |
calc_s_matrix: 0.002 0.002 0.0% |
inverse-cholesky: 0.002 0.002 0.0% |
projections: 0.002 0.002 0.0% |
rotate_psi_s: 0.001 0.001 0.0% |
Set symmetry: 0.001 0.001 0.0% |
Stress: 0.099 0.000 0.0% |
Stress tensor: 0.099 0.099 0.0% |
Other: 2958.786 2958.786 99.9% |---------------------------------------|
-----------------------------------------------------------
Total: 2960.805 100.0%
Memory usage: 220.20 MiB
Date: Thu Jul 23 13:33:38 2026
Final energy: -6.410955 eV
Final max force: 0.000111 eV/Å
Final lattice constant: 2.8707 Å
Now we can pull data from each step of the trajectory to track energy throughout the geometry optimisation.
# Plot the energy convergence during optimisation
import matplotlib.pyplot as plt
traj = Trajectory('Ti_opt.traj')
energies = [atoms.get_potential_energy() for atoms in traj]
plt.figure(figsize=(7, 4))
plt.plot(energies, 'o-', color='steelblue')
plt.xlabel('Optimisation step')
plt.ylabel('Energy (eV)')
plt.title('Ti geometry optimisation convergence')
plt.tight_layout()
plt.show()
6.4.1. Case Study: Just Relax#
Find a material on the MaterialsProject and use DFT to relax it. This will be your personal “Case Study” material also used in future labs. Remember, the computational cost of DFT scales with \(n**3\) (where \(n\) is the number of atoms); so, you may want to find relatively simple materials (with one or two atoms in a unit cell).
6.5. K-point convergence#
DFT properties (e.g. total energy and forces) require integrating over the Brillouin zone; this is done numerically as a weighted sum over a discrete k-point mesh. k-point density scales inversely with cell size: a primitive cell needs a finer grid than a supercell For a large defect supercell, a single Γ-point (k = 0,0,0) is often sufficient. More k-points are always more accurate, but also more expensive.
Note that for any calculation meshes must be converged. An example k-point convergence: we run the same calculation for increasing amounts of k-points, and look for where there is convergence below a particular energy criterion (for example, 1meV/atom).
al = bulk('Al', 'fcc', a=4.05)
k_values = [2, 4, 6, 8]
energies_k = []
for k in k_values:
al.calc = GPAW(xc='PBE', kpts=(k, k, k), mode="pw",txt="k_point-log.txt")
energies_k.append(al.get_potential_energy())
plt.figure(figsize=(7, 4))
plt.plot(k_values, energies_k, 'o-', color='crimson')
plt.xlabel('k-point mesh (N×N×N)')
plt.ylabel('Energy (eV)')
plt.title('K-point convergence test')
plt.tight_layout()
plt.show()
Timing: incl. excl.
-----------------------------------------------------------
Hamiltonian: 0.035 0.000 0.0% |
Atomic: 0.026 0.000 0.0% |
XC Correction: 0.026 0.026 0.0% |
Calculate atomic Hamiltonians: 0.000 0.000 0.0% |
Communicate: 0.000 0.000 0.0% |
Initialize Hamiltonian: 0.000 0.000 0.0% |
Poisson: 0.000 0.000 0.0% |
XC 3D grid: 0.008 0.008 0.0% |
LCAO initialization: 0.290 0.013 0.0% |
LCAO eigensolver: 0.050 0.001 0.0% |
Calculate projections: 0.001 0.001 0.0% |
DenseAtomicCorrection: 0.003 0.003 0.0% |
Distribute overlap matrix: 0.000 0.000 0.0% |
Orbital Layouts: 0.004 0.004 0.0% |
Potential matrix: 0.036 0.036 0.0% |
Sum over cells: 0.004 0.004 0.0% |
LCAO to grid: 0.112 0.112 0.0% |
Set positions (LCAO WFS): 0.115 0.015 0.0% |
Basic WFS set positions: 0.001 0.001 0.0% |
Basis functions set positions: 0.000 0.000 0.0% |
P tci: 0.024 0.024 0.0% |
ST tci: 0.063 0.063 0.0% |
mktci: 0.012 0.012 0.0% |
PWDescriptor: 0.006 0.006 0.0% |
SCF-cycle: 4.407 0.054 0.0% |
Davidson: 3.651 0.430 0.0% |
Apply H: 0.641 0.626 0.0% |
HMM T: 0.015 0.015 0.0% |
Subspace diag: 0.811 0.023 0.0% |
calc_h_matrix: 0.703 0.062 0.0% |
Apply H: 0.641 0.626 0.0% |
HMM T: 0.015 0.015 0.0% |
diagonalize: 0.055 0.055 0.0% |
rotate_psi: 0.030 0.030 0.0% |
calc. matrices: 1.612 0.336 0.0% |
Apply H: 1.276 1.248 0.0% |
HMM T: 0.028 0.028 0.0% |
diagonalize: 0.105 0.105 0.0% |
rotate_psi: 0.051 0.051 0.0% |
Density: 0.448 0.000 0.0% |
Atomic density matrices: 0.010 0.010 0.0% |
Mix: 0.014 0.014 0.0% |
Multipole moments: 0.001 0.001 0.0% |
Pseudo density: 0.424 0.406 0.0% |
Symmetrize density: 0.017 0.017 0.0% |
Hamiltonian: 0.238 0.003 0.0% |
Atomic: 0.185 0.003 0.0% |
XC Correction: 0.182 0.182 0.0% |
Calculate atomic Hamiltonians: 0.002 0.002 0.0% |
Communicate: 0.000 0.000 0.0% |
Poisson: 0.000 0.000 0.0% |
XC 3D grid: 0.048 0.048 0.0% |
Orthonormalize: 0.016 0.002 0.0% |
calc_s_matrix: 0.004 0.004 0.0% |
inverse-cholesky: 0.004 0.004 0.0% |
projections: 0.004 0.004 0.0% |
rotate_psi_s: 0.001 0.001 0.0% |
Set symmetry: 0.035 0.035 0.0% |
Other: 2897.249 2897.249 99.8% |---------------------------------------|
-----------------------------------------------------------
Total: 2902.023 100.0%
Memory usage: 228.11 MiB
Date: Thu Jul 23 13:34:13 2026
6.5.1. Exercise: k-point convergence#
i) Take your material from the Just Relax exercise and perform a k-point convergence until the energy is converged to at least 5meV/atom;
ii) Now build a 2x2x2 supercell of the material and conduct the same k-point convergence.
Do the k-points converge more quickly or less quickly for the supercell? Why is there a difference between the primitive cell and the supercell?
6.5.2. Exercise: defect relaxation#
Build a 2×2×2 supercell of wurtzite GaN and introduce a nitrogen vacancy (V_N) by deleting one N atom. Relax the structure with BFGS using DFT. By how much does the nearest Ga atom move inward toward the vacancy? (This “breathing mode” relaxation is a characteristic signature of point defects.)