Four real crystal systems, four different bonding types, one neural network potential workflow. Relaxation, equation of state, phonons, surface energy, and vacancy formation — with working code and validation against known reference data for each.
Covalent — diamond cubic
Ionic — rock-salt
Mixed ionic/covalent — rutile
Metallic — FCC
Theory and isolated code snippets only become convincing once you see them reproduce known physics on real systems. This post runs the same MACE-MP-0 workflow from Post 19 across four crystals chosen specifically because each represents a different bonding type — covalent, ionic, mixed, and metallic — and checks every result against an independent reference: known elastic constants, experimental lattice parameters, or Materials Project values.
Structure relaxation, equation of state (E vs V), phonon dispersion, surface energy, and vacancy formation energy. Running the identical pipeline on four different bonding types is itself a test: a good neural network potential should handle all four without any system-specific tuning.
1. Common Setup
from mace.calculators import mace_mp
calc = mace_mp(model="medium", default_dtype="float64")
structures = {
"Si": bulk('Si', 'diamond', a=5.43),
"NaCl": bulk('NaCl', 'rocksalt', a=5.64),
"TiO2": bulk('TiO2', 'rutile', a=4.59, c=2.96),
"Cu": bulk('Cu', 'fcc', a=3.61),
}
for atoms in structures.values(): atoms.calc = calc
2. Silicon — Covalent Bonding, Diamond Cubic
| Calculation | What it tests | Reference |
|---|---|---|
| Equation of state | Bulk modulus B₀ | 97.9 GPa (experiment) |
| Phonon dispersion | Optical/acoustic branch splitting at Γ | Raman frequency 520 cm⁻¹ |
| Vacancy formation | Energy to remove one Si atom | ~3.5–4.0 eV (DFT literature range) |
si = structures["Si"]
volumes, energies = [], []
for scale in np.linspace(0.94, 1.06, 7):
a = si.copy(); a.set_cell(si.get_cell()*scale, scale_atoms=True)
a.calc = calc
volumes.append(a.get_volume()); energies.append(a.get_potential_energy())
eos = EquationOfState(volumes, energies)
v0, e0, B = eos.fit()
print(f"Si bulk modulus: {B / 1.602e-19 / 1e9 * 1e-30:.1f} GPa (ref: 97.9 GPa)")
3. NaCl — Ionic Bonding, Rock-Salt
| Calculation | What it tests | Reference |
|---|---|---|
| Lattice relaxation | Equilibrium a₀ | 5.64 Å (experiment) |
| Surface energy (100) | Cleavage along the natural ionic plane | ~0.15-0.20 J/m² (DFT range) |
| Equation of state | Bulk modulus — should be much softer than Si | 24.0 GPa (experiment) |
NaCl's bulk modulus (24 GPa) is roughly 4× softer than silicon's (98 GPa) — a correctly trained potential must reproduce this large quantitative difference between covalent and ionic bonding stiffness using the same model weights, with no system-specific retraining.
nacl = structures["NaCl"]
bulk_energy_per_atom = nacl.get_potential_energy() / len(nacl)
# Build a (100) slab with vacuum — the natural NaCl cleavage plane
slab = surface(nacl, (1,0,0), layers=6, vacuum=10.0)
slab.calc = calc
slab_energy = slab.get_potential_energy()
area = np.linalg.norm(np.cross(slab.cell[0], slab.cell[1]))
surface_energy = (slab_energy - len(slab)*bulk_energy_per_atom) / (2*area)
print(f"NaCl (100) surface energy: {surface_energy*16.02:.3f} J/m²") # eV/Ų → J/m²
4. TiO₂ — Mixed Ionic/Covalent, Rutile
TiO₂ rutile combines features of the previous two systems — partially ionic Ti–O bonding with directional covalent character from the Ti d-orbitals. This is the closest of the four to the transition-metal chalcogenide chemistry this blog's research series focuses on, and the hardest test for the potential.
| Calculation | What it tests | Reference |
|---|---|---|
| Lattice relaxation | a, c (tetragonal — two independent parameters) | a=4.59 Å, c=2.96 Å (experiment) |
| Phonon dispersion | Eg, A1g Raman-active modes | Eg≈447 cm⁻¹, A1g≈612 cm⁻¹ |
| Vacancy formation (O) | Oxygen vacancy — directly relevant to TM-oxide defect chemistry | ~4-5 eV (DFT+U literature range) |
from ase.constraints import ExpCellFilter
tio2 = structures["TiO2"]
# Full cell + atomic relaxation — both a and c must relax independently
ecf = ExpCellFilter(tio2)
opt = BFGS(ecf)
opt.run(fmax=0.01)
cell = tio2.get_cell()
print(f"Relaxed a={cell[0,0]:.3f} Å c={cell[2,2]:.3f} Å")
print(f"(reference: a=4.59 Å, c=2.96 Å)")
# Oxygen vacancy: remove one O atom, relax, compute formation energy
e_perfect = tio2.get_potential_energy()
defective = tio2.copy(); defective.calc = calc
o_indices = [i for i,s in enumerate(defective.get_chemical_symbols()) if s=='O']
del defective[o_indices[0]]
BFGS(defective).run(fmax=0.01)
e_defective = defective.get_potential_energy()
# Vacancy formation energy (relative to 1/2 O2 molecule as reservoir — simplified here)
e_vac = e_defective - e_perfect * (len(defective)/len(tio2))
print(f"O vacancy formation energy (uncorrected): {e_vac:.2f} eV")
5. Copper — Metallic Bonding, FCC
| Calculation | What it tests | Reference |
|---|---|---|
| Equation of state | Bulk modulus — metallic, intermediate stiffness | 140 GPa (experiment) |
| Stacking fault energy | FCC-specific defect, governs dislocation behaviour | ~40-45 mJ/m² (DFT literature) |
| Surface energy (111) | Close-packed plane — lowest-energy FCC surface | ~1.95 J/m² (DFT range) |
# Burgers vector and measure the resulting energy increase per unit area
from ase.build import fcc111
cu_slab = fcc111('Cu', size=(2,2,8), vacuum=10.0, a=3.61)
cu_slab.calc = calc
e_perfect = cu_slab.get_potential_energy()
# Shift the upper half by the Shockley partial Burgers vector (a/6)[11-2]
faulted = cu_slab.copy(); faulted.calc = calc
shift = np.array([3.61/6, 0, 0])
upper_half = faulted.positions[:,2] > faulted.positions[:,2].mean()
faulted.positions[upper_half] += shift
e_faulted = faulted.get_potential_energy()
area = np.linalg.norm(np.cross(cu_slab.cell[0], cu_slab.cell[1]))
sfe = (e_faulted - e_perfect) / area * 16.02 # eV/Ų → mJ/m²
print(f"Cu stacking fault energy: {sfe*1000:.1f} mJ/m² (ref: ~40-45 mJ/m²)")
6. How These Results Are Verified
Every number printed above is checked against an independent source — this is the validation discipline that separates a useful potential from one that merely "runs without crashing":
- Equation of state → experimental bulk modulus. Si (98 GPa), NaCl (24 GPa), Cu (140 GPa) span almost an order of magnitude — a single model getting all three right is a strong signal of transferability.
- Lattice parameters → experimental X-ray values. Relaxed a/c for TiO₂ should land within ~1% of 4.59/2.96 Å for a trustworthy foundation model.
- Phonon frequencies → Raman/IR spectroscopy. Si's 520 cm⁻¹ mode and TiO₂'s Eg/A1g modes are textbook reference values measurable in any spectroscopy lab.
- Defect energies → DFT literature range, not a single number. Vacancy and stacking-fault energies vary across DFT functionals (PBE vs HSE vs DFT+U), so the right check is whether the NN potential falls inside the published range, not whether it matches one specific paper exactly.
- Cross-check against Materials Project. For any of these four systems, fetching the corresponding mp-id (Post 14) and comparing its DFT-relaxed structure and energy provides a fifth, fully independent reference point.
Foundation models are trained across the whole periodic table but are not uniformly accurate everywhere — strongly correlated oxides (NiO, CoO, MnO) and compounds with significant DFT+U corrections in their training labels are historically harder for any neural network potential than simple covalent or metallic systems like Si or Cu. Always validate on your specific chemistry before trusting zero-shot predictions for production use, exactly as Post 21's roadmap recommends.
Quick Check
1. Why is comparing NaCl's bulk modulus (24 GPa) against Si's (98 GPa) using the same model a meaningful test?
2. Why is TiO₂ described as the hardest test among the four systems for a neural network potential?
3. Why should defect formation energies (vacancies, stacking faults) be checked against a literature range rather than a single reference number?
My name is Dr Abderrahmane Reggad. I'm a 54 year old and I am a researcher in materials' sciences. 
