35. Open-shell Mg⁺• in a periodic box — UHF GDF + Makov-Payne

This tutorial exercises the v0.8.0 run_uhf_periodic_gamma_gdf driver — open-shell periodic SCF via Gaussian density fitting at the Γ point — on a single Mg⁺• cation embedded in a vacuum-padded simulation cell. Compares the converged periodic energy to a molecular UHF reference, with the Makov-Payne charge-correction bound on the residual periodic-image interaction.

You will:

  • Build a vacuum-padded simulation cell for a single Mg⁺• cation.

  • Run UHF at increasing cell sizes (10 Å → 20 Å → 30 Å).

  • Compute the Makov-Payne correction and bound the residual periodic-image error.

  • Compare against a molecular UHF reference for the same cation.

Background: see user_guide/multi_k_scf.md § Open-shell GDF for the periodic UHF driver, and user_guide/ewald.md for the Madelung / Makov-Payne theory.

Setup

import vibeqc as vq

# Single Mg⁺• at the centre of a cubic vacuum cell.
def make_mg_cation_cell(L_angstrom: float):
    """Mg⁺• at cell centre, vacuum-padded cubic cell."""
    L = L_angstrom / 0.529177           # Å → bohr
    return vq.PeriodicSystem(
        lattice_vectors=[[L, 0.0, 0.0],
                         [0.0, L, 0.0],
                         [0.0, 0.0, L]],
        atoms=[vq.Atom(12, [L/2, L/2, L/2])],
        charge=+1,                      # +1 cation; UHF will spin-unpair
    )

Vacuum convergence sweep

for L in (10.0, 15.0, 20.0, 30.0):
    cell = make_mg_cation_cell(L)
    basis = vq.BasisSet(cell, "sto-3g")
    result = vq.run_uhf_periodic_gamma_gdf(
        cell, basis,
        options=vq.PeriodicRHFOptions(multiplicity=2),  # doublet
    )
    print(f"  L={L:5.1f} Å: E = {result.energy_per_cell_ha:.6f} Ha   "
          f"⟨S²⟩ = {result.s_squared:.4f}")

Expected: monotonic convergence as the cell grows — the periodic-image interaction (Makov-Payne 1st-order term) goes as 1/L for an isolated charge in a uniform compensating background:

  L= 10.0 Å: E = -198.812... Ha   ⟨S²⟩ = 0.7501
  L= 15.0 Å: E = -198.789... Ha   ⟨S²⟩ = 0.7500
  L= 20.0 Å: E = -198.776... Ha   ⟨S²⟩ = 0.7500
  L= 30.0 Å: E = -198.769... Ha   ⟨S²⟩ = 0.7500

⟨S²⟩ = 0.75 = ¾ ✓ (exact for a doublet with one unpaired electron). If ⟨S²⟩ drifts above ~0.76, that’s spin contamination — typical for transition-metal radicals; for main-group alkali / alkali-earth cations like Mg⁺• it stays near-exact.

Makov-Payne correction

The exact dependence on cell size for a single charge in a neutral-background-compensated periodic cell is, to leading order in 1/L:

E(L) = E(∞) + α_M · q² / (2L) + O(1/L³)

where α_M ≈ 2.837297 is the Madelung constant for a simple- cubic lattice with a uniform compensating background and q = +1 for Mg⁺•. So the corrected energy is:

def makov_payne_e_inf(L_bohr: float, e_per_cell_ha: float, q: float = 1.0) -> float:
    """First-order Makov-Payne correction for an isolated charge."""
    alpha_M = 2.837297                  # SC Madelung constant
    return e_per_cell_ha - alpha_M * q**2 / (2 * L_bohr)


for L in (10.0, 15.0, 20.0, 30.0):
    L_bohr = L / 0.529177
    cell = make_mg_cation_cell(L)
    basis = vq.BasisSet(cell, "sto-3g")
    e_periodic = vq.run_uhf_periodic_gamma_gdf(
        cell, basis,
        options=vq.PeriodicRHFOptions(multiplicity=2),
    ).energy_per_cell_ha
    e_corr = makov_payne_e_inf(L_bohr, e_periodic)
    print(f"  L={L:5.1f} Å:"
          f"   E_periodic = {e_periodic:.6f} Ha"
          f"   E_corrected = {e_corr:.6f} Ha")

Expected: the corrected energies should agree with each other and with the molecular UHF reference to ~µHa across the L range. If they don’t — but the uncorrected energies do show 1/L scaling — your basis set is too small to be behaving as an isolated charge (try increasing to def2-svp).

Molecular UHF reference

mg_cation = vq.Molecule([vq.Atom(12, [0.0, 0.0, 0.0])], charge=+1)
basis = vq.BasisSet(mg_cation, "sto-3g")
opts = vq.UHFOptions()
opts.multiplicity = 2

mol_result = vq.run_uhf(mg_cation, basis, opts)
print(f"\nMolecular UHF Mg⁺• reference:  {mol_result.energy_ha:.6f} Ha")
print(f"⟨S²⟩:                          {mol_result.s_squared:.4f}")

The molecular reference at sto-3g is independent of any periodic-image correction — it’s the direct E(∞) you’re trying to reproduce with the periodic + Makov-Payne approach.

The periodic-corrected E_corr at L = 30 Å should match mol_result.energy_ha to ~µHa. If it doesn’t, check:

  1. Are you using the same basis set for both?

  2. Are both at multiplicity = 2 (doublet)?

  3. Is the SCF converged to comparable tolerance?

  4. Is the cell really vacuum-padded (no atom touching the cell boundary)?

What this tells you

  • Open-shell periodic SCF works at Γ for v0.8.0. run_uhf_periodic_gamma_gdf + run_uks_periodic_gamma_gdf cover the open-shell-Γ paths. Multi-k UHF/UKS is queued for v0.x.x.

  • Charged periodic systems need the Makov-Payne correction. The “raw” periodic energy includes the spurious self-interaction with the periodic images of the same charge; subtracting the leading 1/L term recovers the isolated-cation energy.

  • The molecular limit is reachable. With a big enough vacuum cell + Makov-Payne, the periodic UHF + the molecular UHF agree to µHa. This is the regression test the v0.8.0 open-shell-GDF chip was validated against.

Where this generalises

  • Open-shell ions in solids: e.g. Fe³⁺ in α-Fe₂O₃; the same Makov-Payne correction applies for any net-charged defect.

  • Antiferromagnetic ground states: requires a broken- symmetry initial guess for the spin density. Currently a v0.x.x roadmap item (R3 in the basissetdev periodic- requirements list).

  • Open-shell radicals in molecular contexts: drop the periodic system entirely; use molecular run_uhf with the same UHF options.

  • Open-shell DFT (UKS): replace run_uhf_periodic_gamma_gdf with run_uks_periodic_gamma_gdf and add a functional=.

See also