Slabs and adsorption: the ab-initio surface workflow¶

Surface chemistry is where vibe-qc points: a slab of solid with a molecule sitting on it, and the energy released when the two meet. This tutorial runs that workflow end to end with an honest ab-initio Hartree-Fock / DFT periodic SCF, no machine-learned shortcut. We build a small MgO(100) slab, drop a water molecule on it, and compute the adsorption energy from three single-point periodic calculations. The machine-learning potentials with MACE tutorial walks the same physical system on a learned potential; this is its first-principles counterpart, slower but a genuine total-energy method with a wavefunction.

This page is deliberately concrete about what converges and what costs. A converged all-electron periodic slab is one of the harder things a quantum-chemistry code does, and the workflow below is the one that actually settles, with the cost and accuracy caveats stated plainly.

What you will build¶

The deliverable is a single number, the adsorption energy, assembled from three periodic single points that share a basis, a functional, and a Coulomb gauge so their errors cancel:

E_ads = E[slab + adsorbate] - E[slab] - E[adsorbate, gas]

A negative E_ads means binding. Every term is a periodic SCF on a three-dimensional cell with a vacuum gap, run through vibe-qc’s BIPOLE Ewald path, the route that converges on tight ionic cells where the plane-wave-free truncated-Coulomb path does not.

Setup: a slab is a 3D cell with a vacuum gap¶

The first thing to get right is the geometry. A surface is not a two-dimensional object to the SCF; it is a three-dimensional periodic cell containing a slab of atoms plus a slab of vacuum. Stack enough vacuum (10 to 15 Angstrom) above the atoms and the periodic images stop seeing each other along the surface normal, so the cell behaves like an isolated surface even though the Coulomb sum is fully three-dimensional. This is why vibeqc.build.slab and the example below both produce PeriodicSystem(dim=3, ...), not dim=2.

vibeqc.build ships closed-form surface cells for fcc/bcc/hcp metals (see Slabs and adsorbates). MgO is ionic rocksalt, not in that table, so we lay the cell out by hand, which is the general recipe for any non-metal surface. The MgO(100) face is a square net of alternating Mg and O; one atomic layer holds one Mg and one O, and successive layers shift by half the cell diagonal:

import numpy as np
import vibeqc as vq

ANG2BOHR = 1.0 / 0.529177210903

a = 4.21          # MgO lattice constant (Angstrom)
d = a / 2.0       # interlayer spacing along [100]
n_layers = 2      # a 2-layer slab is the minimum sensible thickness
vacuum = 12.0     # Angstrom of vacuum above the slab

# Third lattice vector = slab thickness + vacuum; centre the slab in z.
slab_thick = (n_layers - 1) * d
c_len = slab_thick + vacuum
z0 = (c_len - slab_thick) / 2.0

# Square in-plane cell of side a; one Mg + one O per layer, checkerboarded,
# with successive layers swapping the Mg/O sublattices.
lattice = np.zeros((3, 3))
lattice[:, 0] = [a, 0.0, 0.0]
lattice[:, 1] = [0.0, a, 0.0]
lattice[2, 2] = c_len

def to_atoms(coords):                      # coords: (Z, x, y, z) in Angstrom
    return [vq.Atom(Z, [x * ANG2BOHR, y * ANG2BOHR, z * ANG2BOHR])
            for (Z, x, y, z) in coords]

slab_coords = []
for L in range(n_layers):
    z = z0 + L * d
    if L % 2 == 0:
        slab_coords += [(12, 0.0, 0.0, z), (8, a / 2, a / 2, z)]   # Mg, O
    else:
        slab_coords += [(8, 0.0, 0.0, z), (12, a / 2, a / 2, z)]   # O, Mg

slab = vq.PeriodicSystem(3, lattice * ANG2BOHR, to_atoms(slab_coords))
print(f"slab: dim={slab.dim}, {len(slab.unit_cell)} atoms, "
      f"{slab.n_electrons()} electrons")

The slab is dim=3 with a 12 Angstrom vacuum gap. A water molecule goes on top of a surface Mg, oxygen down, about 2.2 Angstrom above the top layer; the gas-phase reference is the same molecule in the same cell so the Coulomb gauge matches:

# top-layer Mg site (the slab's top layer swaps sublattices each layer)
z_top = max(c[3] for c in slab_coords)
mg_xy = (0.0, 0.0) if (n_layers - 1) % 2 == 0 else (a / 2, a / 2)

h = 2.2  # Angstrom, O above the surface
water_on = [(8, mg_xy[0], mg_xy[1], z_top + h),
            (1, mg_xy[0] + 0.7572, mg_xy[1], z_top + h + 0.586),
            (1, mg_xy[0] - 0.7572, mg_xy[1], z_top + h + 0.586)]

zc = c_len / 2.0  # gas-phase water, same cell, centred
water_gas = [(8, 0.0, 0.0, zc),
             (1, 0.7572, 0.0, zc + 0.586),
             (1, -0.7572, 0.0, zc + 0.586)]

slab_with_water = vq.PeriodicSystem(3, lattice * ANG2BOHR,
                                    to_atoms(slab_coords + water_on))
water_box = vq.PeriodicSystem(3, lattice * ANG2BOHR, to_atoms(water_gas))

Run the three single points¶

With the geometry in hand, the calculation is three calls to the same periodic SCF driver. We use closed-shell RKS with the PBE functional and the minimal STO-3G basis, the cheapest combination that still produces a real ionic-surface energy; jk_method="bipole" selects the CRYSTAL-gauge Ewald J/K build that converges on tight ionic cells. A thin oxide slab has surface states that crowd the gap, so the plain SCF oscillates; Fermi-Dirac smearing (smearing_temperature) lets the occupations go fractional near the Fermi level and settles it. The value matters, keep it below the slab’s gap for the physical state (see the warning at the end of this section); the 0.01 Ha used here is a robust starting point that the warning then helps you refine. The same three calls share basis, functional, cell, and gauge so their errors cancel in E_ads:

def slab_energy(system, label):
    vq.run_periodic_job(
        system,
        vq.BasisSet(system.unit_cell_molecule(), "sto-3g"),
        method="RKS",
        functional="pbe",
        jk_method="bipole",          # CRYSTAL-gauge Ewald J/K
        kpoints=(1, 1, 1),           # Gamma-only here; (4, 4, 1) for production
        smearing_temperature=0.01,   # Hartree; Fermi-Dirac, settles the slab
        conv_tol_energy=1e-7,
        max_iter=120,
        output=label,                # writes <label>.out with the full trace
    )

slab_energy(slab,            "mgo_slab")
slab_energy(water_box,       "water_gas")
slab_energy(slab_with_water, "mgo_slab_water")

Each call writes a <label>.out with the banner, the SCF iteration trace, and the energy breakdown. Read the three E_total values back and take the difference:

def total_energy(out_stem):
    for line in open(f"{out_stem}.out"):
        if "E_total" in line:
            return float(line.split("=")[-1].split()[0])

e_slab  = total_energy("mgo_slab")
e_water = total_energy("water_gas")
e_sys   = total_energy("mgo_slab_water")

e_ads = e_sys - e_slab - e_water
print(f"E_ads = {e_ads:+.6f} Ha = {e_ads * 27.2114:+.3f} eV")

Each <label>.out ends with an energy summary. All three single points converge; here is the full self-consistent set (one Fermi-Dirac smearing value shared across all three so the basin and the gauge cancel in the difference):

mgo_slab        E_total = -542.74303887 Ha   (converged, 18 iter)
water_gas       E_total =  -75.22373311 Ha   (converged, 13 iter)
mgo_slab_water  E_total = -618.12112933 Ha   (converged, 27 iter)

E_ads = -618.12112933 - (-542.74303887) - (-75.22373311)
      = -0.154357 Ha = -4.200 eV

The water binds (negative E_ads), as it should: water on MgO(100) is a real, if modest, adsorption. The magnitude is another matter: -4.2 eV is far stronger than the converged literature and machine-learned value (the MACE tutorial gives ~-0.8 eV for this system). That gap is the signature of the approximations stacked into this cheap single point, the minimal STO-3G basis (large basis-set superposition error, the dominant over-binding here), no relaxation, a thin slab, and full monolayer coverage, all unpacked below. The sign and the workflow are right; the number is a starting point, not a publication value.

The slab energy is also sensitive to the smearing, which is the whole point of the warning below. The same 2-layer slab converges to -542.755 Ha with a small smearing (0.003 Ha, the physical gapped state, what the example code requests) but to -542.743 Ha with a large one (0.01 Ha, a near-metallic basin), a 12 mHa (0.3 eV) shift. Use the same, sub-gap smearing for all three single points and the shift cancels in E_ads; mix two smearings, or read the slab from the wrong basin, and the adsorption energy is corrupted.

The cross-checks that say the machinery itself is sound sit right next to the result: the same BIPOLE path converges bulk rocksalt MgO (the FCC primitive, RHF/STO-3G) to -271.752 Ha per formula unit, and the gas-phase water term, -75.224 Ha, matches a molecular PBE/STO-3G water. So the recipe is trustworthy; it is the model (basis, slab, coverage, geometry) that is deliberately cheap here. Tighten every knob in the Next section before quoting an adsorption energy in a paper, and for publication-grade screening across many sites and surfaces in seconds, let the MACE tutorial run this very system on a learned potential.

Warning

This slab is a small-gap, near-metallic case, and vibe-qc says so. A two-layer MgO/STO-3G slab has a HOMO-LUMO gap of only ~0.5 eV at the Gamma point. If the Fermi-Dirac smearing temperature is set comparable to that gap, vibe-qc warns that the SCF has settled a near-metallic basin with several fractional occupations, which can sit hundreds of mHa from the physical gapped solution. The fix it points you to is to drop the smearing well below the gap (here ~0.003 Ha, ~0.08 eV) so the SCF lands the gapped state. Heed that warning: a slab energy from the wrong basin makes E_ads meaningless. This is the periodic-SCF discipline in action, an oscillating or wrong-basin ionic slab is a signal to diagnose, not to crank the smearing until the number stops moving.

Note

Cost. An all-electron, exact-ERI periodic slab is genuinely expensive: each of these single points sums two-electron integrals over the lattice images of a large, anisotropic vacuum cell, so a 2-layer STO-3G slab single point is minutes of wall time, not seconds, and a production slab (thicker, larger basis, a surface k-mesh) is substantially more. This is the price of a true total-energy method with a wavefunction; the MACE tutorial trades it for a learned potential when you need thousands of geometries fast.

Theory: the slab model¶

The sections below unpack the three ideas the script leans on: how a surface becomes a tractable periodic cell, what the adsorption energy actually measures, and the pitfalls that separate an illustrative number from a publishable one.

Two-dimensional periodicity from a three-dimensional cell¶

A clean crystal surface is periodic in the two directions parallel to the face and aperiodic along the normal: it simply ends, with vacuum above and bulk below. A finite calculation cannot carry semi-infinite bulk, so the standard construction is the slab supercell. Take a few atomic layers thick enough to develop a bulk-like interior, place them in a cell whose third lattice vector is much longer than the slab, and fill the remainder with vacuum. Repeated periodically, this tiles space with slabs separated by vacuum gaps.

Two convergence parameters control the fidelity of the model. The vacuum gap must be wide enough that the electron density of one slab has decayed to zero before it reaches the next image; 10 to 15 Angstrom is typical for Gaussian basis sets, whose orbitals are spatially compact. The slab thickness must be large enough that the middle layers recover the bulk electronic structure and both surfaces are well separated; properties converge from below as layers are added, and a two- or three-layer slab is a starting point, not a converged answer.

vibe-qc treats the slab as a genuine dim=3 system and sums the Coulomb interaction in all three dimensions with Ewald (see Madelung constants via Ewald summation for why the three-dimensional electrostatic sum is only conditionally convergent and needs Ewald rather than naive truncation). A symmetric, charge- neutral slab carries no net dipole across the cell, so no dipole correction is required; an asymmetric slab or a polar face would.

A note on the alternative: a PeriodicSystem built with dim=2 is periodic in only two directions, and in vibe-qc that path currently falls back to a direct, non-Ewald Coulomb treatment whose energy is the isolated-cluster (molecular) limit, not a true two-dimensional surface energy. For a real surface calculation, use the dim=3 slab-plus-vacuum construction above.

What the adsorption energy measures¶

The adsorption energy is a reaction energy for

slab + molecule(gas)  ->  slab-with-molecule

E_ads = E[slab + adsorbate] - E[slab] - E[adsorbate, gas].

A negative value means the bound state is lower in energy than the separated fragments: the molecule binds. The magnitude sorts physisorption (tenths of an eV, dispersion and weak electrostatics) from chemisorption (an eV or more, a real bond). Because the same per-atom and per-cell systematic errors appear in all three terms with the same basis, functional, and cell, they cancel to a large degree in the difference, which is why a modest basis can give a meaningful E_ads even when no single total energy is near the basis-set limit.

The cancellation is only as good as the consistency. All three single points must use the same basis set, the same functional, the same k-mesh, and the same Coulomb gauge. Mixing a molecular gas-phase code with a periodic slab code, or changing the cell between terms, breaks the cancellation and corrupts E_ads.

Pitfalls: BSSE, coverage, relaxation¶

Three systematic effects separate the number above from a converged adsorption energy.

Basis-set superposition error (BSSE). With atom-centred Gaussian basis functions, the slab-plus-adsorbate complex borrows basis functions from its partner and artificially lowers its energy, over-binding the adsorbate. The counterpoise correction recomputes each fragment in the presence of the other’s (ghost) basis functions and subtracts the difference. BSSE shrinks as the basis grows; with STO-3G it is large, so the raw E_ads here is an over-bound estimate.

Coverage. A p(1x1) surface cell places one adsorbate per surface metal site, the dense full-monolayer limit, and the periodic images of the adsorbate interact laterally at the in-plane lattice spacing. The adsorption energy you compute is therefore the high-coverage value, lateral interactions included. To approach the dilute, single-molecule limit, enlarge the lateral cell to a (2x2) or (3x3) supercell so the adsorbates are far enough apart not to interact, and watch E_ads converge with cell size.

Relaxation. The energies above are single points at a fixed, hand-placed geometry. A real adsorption energy uses relaxed structures: the clean slab relaxed, the gas molecule relaxed, and the complex relaxed with the adsorbate and the top surface layers free while the bottom layers stay fixed at bulk positions. vibe-qc drives this frozen-substrate relaxation through relax_atoms with a freeze_indices list (see Slabs and adsorbates and Periodic geometry optimization); relaxation always lowers the energy of the complex and so makes E_ads more negative. On strongly basic surfaces, also check whether the molecule stayed intact or dissociated, the two are different physical states.

References¶

Solid-state basis sets (pob-TZVP). M. F. Peintinger, D. Vilela Oliveira, and T. Bredow, “Consistent Gaussian basis sets of triple-zeta valence with polarization quality for solid-state calculations,” J. Comput. Chem. 34, 451 (2013). The production basis to step up to from STO-3G; see Why solid-state calculations use pob-TZVP.
PBE exchange-correlation functional. J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett. 77, 3865 (1996).
Ewald summation for crystalline Gaussian densities. V. R. Saunders, C. Freyria-Fava, R. Dovesi, L. Salasco, and C. Roetti, “On the electrostatic potential in crystalline systems where the charge density is expanded in Gaussian functions,” Mol. Phys. 77, 629 (1992). The multipolar 3D Ewald scheme vibe-qc’s BIPOLE path implements.
Counterpoise correction for BSSE. S. F. Boys and F. Bernardi, “The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors,” Mol. Phys. 19, 553 (1970).
CRYSTAL reference. A. Erba, J. K. Desmarais, S. Casassa, et al., “CRYSTAL23: a program for computational solid state physics and chemistry,” J. Chem. Theory Comput. 19, 6891 (2023). The Gaussian-basis periodic code whose gauge vibe-qc’s BIPOLE path follows.

Next¶

Step up the basis and k-mesh. Swap STO-3G for pob-tzvp and add a surface-plane k-mesh (kmesh=(4, 4, 1)); see LiH rocksalt with pob-TZVP and LiH at multiple k-points for the multi-k machinery.
Relax the structure. Run the frozen-substrate relaxation from Slabs and adsorbates and watch E_ads drop as the complex relaxes.
Add dispersion. For physisorption, the D3-BJ correction matters; pass dispersion="pbe" (see Dispersion corrections and the periodic D3 section of Slabs and adsorbates).
Screen fast first, confirm with ab initio. The MACE tutorial screens hundreds of adsorption sites in seconds on a learned potential; use it to find candidate geometries, then confirm the winners with the first-principles workflow here.