A new computational method gives finite quantum-chemical clusters constant-potential boundary conditions, allowing calculations to represent electron exchange without treating an entire surrounding system. Tests covered one water molecule, idealized hydrogen rings and a 40-atom lithium cluster. The work appears as an arXiv version 1 preprint.
The central device is a self-energy, a mathematical term that represents coupling to an environment. It is diagonal and independent of energy, with positive broadening under the wide-band approximation. The density-matrix integral is evaluated analytically, and the resulting matrix is used with conventional Hartree-Fock or density-functional theory calculations. The authors describe the localized self-energy as a compact way to represent a missing environment.
The custom grand-canonical self-consistent-field implementation interfaces with PySCF 2.13.0. It supports canonical Hartree-Fock and restricted Kohn-Sham density-functional calculations, and adapts grand-potential EDIIS, Anderson and a hybrid EDIIS plus Anderson mixing scheme.
A numerical test came first
The first check used a single water molecule with MINAO and STO-3G basis sets, alongside restricted Hartree-Fock and PBE density-functional treatments. The self-energy was applied to oxygen 2s and 2p and hydrogen 1s valence orbitals. The calculations compared the listed mixing schemes from a superposition-of-atomic-densities start and from converged canonical guesses, using up to 200 iterations. Grand-canonical convergence required density-matrix changes below 10^-6 and grand-potential changes below 10^-8.
At strong broadening, eta = 1 x 10^-1 Ha, grand-potential EDIIS, Anderson and the hybrid scheme reached the same reported fixed point: a grand potential of -75.496810 Ha, an electron number of 10.1389 and a residual of 1.19 x 10^-1. The hybrid was fastest, taking 27 iterations from the superposition-of-atomic-densities start and 43 from the canonical start. CDIIS did not converge from the former start. This was a numerical result about the calculation, not an experimental measurement.
Hydrogen rings show the boundary effect
The next test used a ring of 40 hydrogen atoms with 0.74-angstrom bonds. Fragments containing six to 20 atoms, in steps of two, were treated with RKS/PBE0/MINAO. Only the 1s orbitals on the outer hydrogens were coupled to the self-energy. L-BFGS optimization adjusted the coupling parameters, while the chemical potential was adjusted to keep the fragment neutral without an external charge.
Without that coupling, the fragment response was purely repulsive. After optimization, H6's well depth differed from the full-ring reference by about 0.57 kilocalories per mole, or 2.2%, at the minimum. The H38 comparison differed by about 2.7 kilocalories per mole, or 10.3%. Residual deviations still reached about 0.9 kilocalories per mole in intermediate and tail regions. The compact fragment therefore tracked the reference potential closely near its minimum, but not perfectly across the full range.
The electronic structure changed with fragment size. Chemical-potential asymmetry and the self-energy correction diminished as more of the ring was treated explicitly. The HOMO-LUMO gap, the frontier-orbital separation used in the calculation, was about 0.38 Ha for H6, 0.04 Ha for H40 and 0.11 Ha for H38. Fractional occupations corresponded to an effective temperature of roughly 2 to 3 x 10^4 K. The authors attribute the metal-like response mainly to broadened, fractional occupations rather than a complete collapse of the gap.
Lithium reveals the trade-off
The final model was Li40, built from a body-centered-cubic (110) surface with a 3.51-angstrom lattice constant and treated with PBE0/STO-3G. The self-energy shift was set to zero, selected lithium 2s orbitals were coupled, and center, surface-center and edge topologies were tested. Broadening ranged from 0.01 to 500 mHa. For a charge-response scan, the calculation used a positive point charge of +0.4 e and eta = 2.0 mHa.
The broadening sweep exposed a sharp sensitivity in the orbital spectrum. The canonical HOMO-LUMO gap was 24.5 mHa. As eta approached that scale, the HOMO occupation fell from 1.55 at 5 mHa to 1.07 at 10 mHa and nearly vanished at 20 mHa. Many previously unoccupied orbitals acquired small occupations, and the effective gap collapsed. As eta approached zero, the canonical spectrum was recovered, although the chemical potential stayed near the HOMO.
With center coupling, the positive point charge produced a grand-canonical well of about 21.5 kilocalories per mole at roughly 1.75 angstroms, compared with 19.4 kilocalories per mole for the canonical reference. As the charge approached 1 angstrom, the electron count rose monotonically to 0.20 e above the starting value. Overall density differences were small, but they were largest for edge coupling. In this model, the calculation captured the open-boundary polarization response the method was designed to test.
A method still in development
The paper does not establish how the method will perform at a realistic electrochemical interface. The authors present it as practical grand-canonical boundary conditions and a foundation for future embedding, while explicitly saying it is not a complete embedding framework. They identify solvation, automatically transferable self-energy parameters and correlated methods such as RPA as future steps.
The code is not described as a current public release. The paper says the PySCF modifications used for the results will be made available upon acceptance in a peer-reviewed journal. For now, the evidence remains a set of model calculations, with broader transferability still to be shown.
Paper data and sources
Original title: An embedding method with constant potential boundary conditions
Authors: Lisa Hetzel, Martin Head-Gordon, Christopher J. Stein
Journal/Repository: arXiv
Status: Preprint, not yet peer-reviewed
First online: 2026-08-25
DOI: Not available
Original paper · Full text