Open-shell Green’s-function surface embedding

Experimental method

Green’s-function surface embedding is experimental and not yet validated to thick-slab GDF parity. The numbers in this tutorial are qualitative: they exercise the machinery and the open-shell plumbing, not a converged adsorption energy. See the user-guide page for the feature-status detail.

Surface embedding treats an adsorbate plus a patch of substrate as a localized perturbation on a laterally periodic, semi-infinite clean surface, reaching the dilute-limit adsorbate with no slab-thickness artifact and no adsorbate-adsorbate image interaction. This tutorial drives the open-shell path: an odd-electron cell solved with a spin-polarised (UHF) embedded density, through the ASE calculator.

You will:

• Build an odd-electron slab cell (a 3-atom H chain, a doublet).

• Run it through EmbeddedSurfaceCalculator with multiplicity=2 and scf_method="uhf".

• Get the embedded energy and numerical forces, the way an ASE optimiser would.

• See why an open-shell cell needs both the right multiplicity and a spin-polarised density route.

Background: the surface-embedding user guide for the API and the two layers; Open-shell Mg⁺• for open-shell periodic SCF in the conventional supercell picture; ASE integration for the calculator interface.

Why open-shell needs two ingredients

A cell with an odd electron count cannot be a closed-shell singlet. Two independent things have to agree:

1. Multiplicity. 2S + 1 must match the electron-count parity. A 3-electron cell is a doublet (multiplicity=2). If you leave the default multiplicity=1, the unit-cell molecule fails its parity check (“n_electrons and multiplicity are inconsistent”) before any embedding work runs.

2. A spin-polarised density. The default density method (one_shot/hf) is closed-shell. Spin-up and spin-down must be allowed to differ, so you pass scf_method="uhf".

Set one without the other and you either trip the parity error or silently restrict the density to closed-shell.

Setup

The geometry is a 3-atom H chain stacked along the surface normal in a slab cell (a long c-axis marks it as a 2D-periodic slab). Three H atoms carry three electrons, an open-shell doublet. Centring the chain in the lateral plane forces the in-plane forces to vanish by symmetry, leaving only the surface-normal (z) force.

import numpy as np
from ase import Atoms
from vibeqc.periodic_embedding.ase_calc import EmbeddedSurfaceCalculator

# 3-atom H chain (slab cell): 3 electrons -> doublet.
a, c, sp = 4.0, 18.0, 1.5            # Angstrom
z0 = (c - 2 * sp) / 2.0
atoms = Atoms(
    "H3",
    positions=[(a / 2, a / 2, z0 + i * sp) for i in range(3)],
    cell=[a, a, c], pbc=[True, True, True],
)
atoms.set_tags([0, 0, 1])            # bottom 2 = substrate, top = region I

The EmbeddedSurfaceCalculator carries the open-shell choice: multiplicity=2 (threaded into the PeriodicSystem it builds) and scf_method="uhf" (forwarded to the runner). The contour window (e_bottom, e_fermi) is fixed explicitly so the numerical-force loop holds it constant across displacements; force_atom_indices=[2] makes only the region-I atom mobile.

atoms.calc = EmbeddedSurfaceCalculator(
    basis_name="sto-3g",
    surface_k_mesh=(1, 1),
    contour_n_nodes=24,
    e_bottom=-2.2, e_fermi=-1.0,
    charge=0,
    multiplicity=2,                  # 3 electrons -> doublet
    scf_method="uhf",                # spin-polarised density
    force_atom_indices=[2],          # only the region-I atom is mobile
)

Run it

energy = atoms.get_potential_energy()   # eV
forces = atoms.get_forces()             # eV / Angstrom

print(f"Embedded energy: {energy:.4f} eV")
for i, f in enumerate(forces):
    print(f"  atom {i}: {f[2]:+.4f} eV/A (z)")
print(f"Region-I occupation: {atoms.calc.last_result.occupation:.4f} e-")

Expected (coarse 24-node contour, STO-3G, qualitative):

Embedded energy: 128.0366 eV
  atom 0: +0.0000 eV/A (z)
  atom 1: +0.0000 eV/A (z)
  atom 2: +11.5797 eV/A (z)
Region-I occupation: 0.2420 e-

The substrate atoms (outside force_atom_indices) keep exactly zero force; the in-plane components of the mobile atom vanish by the chain’s x/y symmetry; only the surface-normal force is nonzero. The force is the exact central-difference gradient of the returned energy, so as the energy model matures the force tracks it automatically.

Note

The embedded “energy” here is the placeholder Layer-A total (band_energy + nuclear_repulsion), not a converged adsorption energy, hence the un-physical magnitude. What this tutorial validates is the open-shell plumbing: the doublet cell builds, the UHF density runs, and energy plus forces come back finite and symmetry-correct.

The complete script is examples/embedding/open_shell_adsorbate.py.

What this tells you

• Open-shell embedded surfaces run via multiplicity plus scf_method="uhf". charge and multiplicity are threaded into the PeriodicSystem the calculator builds; the UHF route gives the spin-polarised density.

• Numerical forces come for free. get_forces() central-differences the embedded energy with the contour window held fixed, so any ASE optimiser or vibrational driver works unchanged.

• The method is experimental. Absolute energies are qualitative until thick-slab GDF parity lands; use it to explore the embedding workflow, not for production adsorption energies yet.

Theory

The substrate is eliminated exactly in the Gaussian basis: region I (the tagged atoms plus the surface layers) couples to the semi-infinite substrate \(S\) only through an energy-dependent embedding self-energy, the LCAO form of Inglesfield’s embedding potential (Inglesfield, J. E. (1981). *A method of embedding.* Journal of Physics C: Solid State Physics, **14**(26), 3795-3806. doi:[10.1088/0022-3719/14/26/015](https://doi.org/10.1088/0022-3719/14/26/015)). On the UHF route each spin \(\sigma\) has its own region-I Green’s function at complex energy \(z\) and surface k-point,

\[ \mathbf G^{\sigma}_I(z) = \Big[\, z\,\mathbf S_{II} - \mathbf h_{II} - \mathbf J[\mathbf D^{\alpha} + \mathbf D^{\beta}] + \mathbf K[\mathbf D^{\sigma}] - \boldsymbol\Sigma_{\mathrm{emb}}(z) \Big]^{-1}, \qquad \boldsymbol\Sigma_{\mathrm{emb}}(z) = (\mathbf H_{IS} - z\,\mathbf S_{IS})\, \mathbf g_S(z)\,(\mathbf H_{SI} - z\,\mathbf S_{SI}), \]

where \(\mathbf h_{II}\) is the one-electron block and the same \(\boldsymbol\Sigma_{\mathrm{emb}}\) serves both spins. The substrate surface Green’s function \(\mathbf g_S\) comes from principal-layer decimation (López Sancho, M. P., López Sancho, J. M., and Rubio, J. (1985). *Highly convergent schemes for the calculation of bulk and surface Green functions.* Journal of Physics F: Metal Physics, **15**(4), 851-858. doi:[10.1088/0305-4608/15/4/009](https://doi.org/10.1088/0305-4608/15/4/009)), and \(\boldsymbol\Sigma_{\mathrm{emb}}\) is linearised about a reference energy inside the occupied window (Ishida, Hiroshi (2001). *Surface-embedded Green-function method: A formulation using a linearized-augmented-plane-wave basis set.* Physical Review B, **63**(16), 165409. doi:[10.1103/PhysRevB.63.165409](https://doi.org/10.1103/PhysRevB.63.165409)). Each spin density is the contour integral up to the fixed Fermi level,

\[ \mathbf D^{\sigma}_I = -\frac{1}{\pi}\operatorname{Im} \int_{E_{\mathrm{bottom}}}^{E_F} \mathbf G^{\sigma}_I(E + i0^{+})\,\mathrm dE , \]

evaluated on a complex contour with contour_n_nodes points, and the J and K terms are iterated to self-consistency. Holding e_bottom and e_fermi fixed is what makes the central-difference force the derivative of one energy function. The Layer B adsorption energy of the closed-shell walkthrough instead takes the change in the integrated density of states from Lloyd’s formula (Lloyd, P. (1967). *Wave propagation through an assembly of spheres: II. The density of single-particle eigenstates.* Proceedings of the Physical Society, **90**(1), 207-216. doi:[10.1088/0370-1328/90/1/323](https://doi.org/10.1088/0370-1328/90/1/323)). The full derivation is in the user-guide theory section.

Where this generalises

• Closed-shell surfaces: drop multiplicity/scf_method and use the defaults (singlet, one-shot). See examples/embedding/simple_adsorbate.py for the Layer A plus Layer B (Lloyd adsorption energy) walkthrough.

• Spin-polarised DFT (UKS): a roadmap item. The open-shell route is UHF (J + K) today; UKS with a Vxc augmentation is future work (the original investigation).

• Conventional supercell slabs: for the non-embedding picture, see Slabs and adsorbates and the slab-adsorbate DFT tutorial.

See also

• Green’s-function surface embedding: the user-guide reference for the runner and the calculator.

• Open-shell Mg⁺• in a periodic box: open-shell periodic SCF with the Makov-Payne correction.

• Open-shell systems (UHF, UKS, UMP2): molecular open-shell basics.