Machine-Learned Interatomic Potentials

8. Machine-Learned Interatomic Potentials#

8.1. Overview#

Questions

  • What are MLIPs and why use them?

  • How do I use the MACE foundational model?

Objectives

  • Use ASE and MACE to conduct geometry relaxation on large unit cells.

Keypoints

  • MLIPs are about ~1000× faster than DFT for forces and geometry relaxation

  • Foundational models are trained on DFT data for a large number of elements

  • Use MACE-MP-0 for the geometry stage; use GPAW (DFT) for the electronic structure stage

  • MACE has no concept of charge state, spin, or electronic structure

  • The multi-fidelity principle: match method to the property being simulated and the accuracy required.

8.2. Lecture Slides#

The slides for this tutorial are embedded below. 📥 Download slides (.pptx)  |  Open in full screen

8.3. Using MACE and ASE#

MACE integrates with ASE as a standard calculator. The key function is mace_mp() which downloads and caches the pre-trained model automatically.

  • model="medium" is the recommended default — good balance of speed and accuracy

  • dispersion=False van der Waals correction (use True for layered materials where this bonding is present)

import mace
from mace.calculators import mace_mp

calc = mace_mp(model="medium", dispersion=False)

8.4. Example: Diamond single-point calculation#

Single-point calculation refers to a calculation with a single electronic minimisation; ie, no geometry relaxation. We go throught the usual process: create the structure –> attach the calculator –> do the calculation –> report the results.

from ase.build import bulk

diamond = bulk('C', 'diamond', a=3.57)
diamond.calc = calc

energy = diamond.get_potential_energy()
forces = diamond.get_forces()

print(f"Diamond total energy: {energy:.4f} eV")
print(f"Max force (should be ~0 at equilibrium): {np.max(np.abs(forces)):.6f} eV/Å")

8.4.1. Example: Multi-Fidelity Geometry Relaxation#

A key strategy in computational materials science is multi-fidelity optimisation: use a cheap, approximate method to get close to the minimum, then refine with a more expensive method only when needed.

For defect calculations:

  • MACE for the geometric relaxation.

  • GPAW for the single-point electronic structure calculation on the pre-relaxed geometry

In the example below we relax a pristine (ie, non-defective) diamond supercell with mace and time how long the relaxation takes.

from ase.optimize import BFGS
from ase.build import bulk, make_supercell
import time

prim = bulk('C', 'diamond', a=3.57)
sc = make_supercell(prim, np.diag([2, 2, 2]))
sc.rattle(0.2)

print(f"Initial max force: {np.max(np.abs(sc.get_forces())):.4f} eV/Å")

sc.calc = mace_mp(model="medium", dispersion=False)

t0 = time.time()
opt = BFGS(sc, logfile=None)
opt.run(fmax=0.01)
t1 = time.time()

print(f"Final max force:   {np.max(np.abs(sc.get_forces())):.4f} eV/Å")
print(f"Relaxation time:   {t1-t0:.1f} seconds")

8.4.2. Exercise: Time comparison#

Do the same geometry relaxation as above but now using DFT. Make sure to time how long the relaxation takes. This tells you, for an 8-atom cell, the scaling between MACE and DFT. The scaling will also depend on the numerical parameters used for the DFT. If the DFT job is running for too long, feel free to “kill it”; you won’t need these results later.

8.4.3. Example: Multi-fidelity simulation of the NV centre in diamond.#

In Lab three we constructed diamond with a NV-centre. Here we construct the same system again, do a geometry relaxation with MACE, followed by a bandstructure calculation with DFT.

It is important that after building the point defect structure we give everything a little “rattle”. This is to break symmetry, so that we don’t relax to a saddle-point; rattling should help find the local minimum.

8.4.3.1. Step 1: build structure#

This follows code introduced in Lab 3. We visualise the end result to check it is correct.

from ase.geometry import get_distances

sc_matrix = np.diag([4, 4, 4])  
nv_cell = make_supercell(diamond_primitive, sc_matrix)

_, dist_matrix = get_distances(nv_cell.positions, nv_cell.positions,
                                cell=nv_cell.cell, pbc=True)
np.fill_diagonal(dist_matrix, np.inf)   
nearest_idx = np.argmin(dist_matrix[0])  

symbols = list(nv_cell.get_chemical_symbols())
symbols[0] = 'N'
nv_cell.set_chemical_symbols(symbols)

del nv_cell[nearest_idx]

fig, ax = plt.subplots(figsize=(6, 6))
plot_atoms(nv_cell, ax, rotation=('10x,10y,0z'), radii=0.4)
plt.show()

8.4.3.2. Step 2: geometry relaxation with MACE#

Give everything a little rattle, then relax. Double check that the forces are less than the convergence criteria set (0.1 eV/Angstrom).

nv_cell.rattle(0.2)

nv_cell.calc = mace_mp(model="medium", dispersion=False)
opt = BFGS(nv_cell, logfile=None)
opt.run(fmax=0.01)

print(f"\nMax residual force: {np.max(np.abs(nv_cell.get_forces())):.4f} eV/Å")

8.4.3.3. Step 3: bandstructure with DFT#

Now that we have a relaxed structure we can calculate the electronic structure. We have to do this with DFT (as MACE is used for structural properties only).

nv_cell.calc = GPAW(mode=PW(600),xc='PBE', kpts={'size': (6, 6, 4)}, txt='nv_scf.txt')
nv_cell.get_potential_energy()
nv_cell.write('nv_diamond.gpw')

path = nv_cell.cell.bandpath('GMKGALHA', npoints=100)
nv_cell.calc = GPAW('nv_diamond.gpw').fixed_density(kpts=path,txt='nv_bs.txt',symmetry='off')
nv_cell.get_potential_energy()

bs = nv_cell.calc.band_structure()
bs.plot(emin=0, emax=12, filename='nv_cell_bandstructure.png', show=True)

8.5. Exercise: MACE equation of state#

Use MACE to compute the equation of state for FCC Al. How does it compare to the value in lab four, where the same EOS was calculated using the EMT potential?

8.5.1. Case Study: Multi-fidelity simulation#

Apply the same multi-fidelity simulation technique to predict the bandstructure of a vacancy defect in your case study material. How does the material distort after relaxation? Can this be rationalised with chemical arguments?

Contents