Skip to content

Effective core potentials (ECP): replacing the core with a pseudopotential

Heavy atoms have a lot of core electrons that barely change during chemistry but are expensive to describe and — past the third row — genuinely relativistic. An effective core potential (ECP), also called a pseudopotential, throws those inert core electrons out of the explicit calculation and replaces them with an analytic potential acting on the remaining valence electrons. You reach for one whenever a molecule contains heavier elements (roughly Br and below, and anything in the 4th row or heavier): the ECP shrinks the number of electrons and basis functions and folds scalar-relativistic effects into the potential, so a nonrelativistic code gives near-relativistic valence energetics. This tutorial computes a BHHLYP energy of hydrogen bromide (HBr) with an ECP on bromine and an ordinary all-electron basis on hydrogen.

A little theory

An ECP splits the electrons of an atom into a frozen core and an active valence. The core is not solved for; instead its effect on the valence electrons — nuclear attraction screened by the core, plus the Pauli repulsion that keeps valence orbitals orthogonal to the (absent) core — is baked into a short, angular-momentum-projected radial potential:

V_ECP(r) = V_local(r) + Σ_l V_l(r) |l><l|

Each V_l is a small sum of Gaussian terms fitted to reproduce all-electron (often relativistic) reference calculations. Because the potential carries the relativistic contraction of the core, an ECP is the standard cheap route to scalar-relativistic valence chemistry for heavy atoms. In practice you never write the potential yourself: correlation-consistent "-PP" basis sets (e.g. aug-cc-pVDZ-PP) ship a matched valence basis and the ECP that goes with it. OpenQP reads that ECP from the basis-set data and evaluates its integrals through the libecpint library; the rest of the SCF is unchanged. The one thing to get right is the per-atom basis: heavy atoms get the -PP set (with its ECP), light atoms keep an all-electron set. For depth, see the OpenQP manual.

Input-file style

The runnable deck is inputs/hbr_ecp_energy.oqp — HBr, a BHHLYP single point, with aug-cc-pVDZ-PP (ECP) on Br and the all-electron aug-cc-pVDZ on H. Annotated:

rks/bhhlyp/aug-cc-pvdz-pp;aug-cc-pvdz    # model/functional/per-atom basis list
dftgrid(rad_type=becke)                  # Becke radial grid for the XC quadrature
geom="""
Br   0.000000000   0.000000000   0.000000000
H    0.000000000   0.000000000   1.407611000
"""

Key points:

  • The basis component aug-cc-pvdz-pp;aug-cc-pvdz is the whole trick. The semicolon-separated list assigns one basis per atom, in the same order as the geometry: the first atom (Br) gets aug-cc-pVDZ-PP, the second (H) gets aug-cc-pVDZ. The -PP suffix ("pseudopotential") is what makes OpenQP load bromine's ECP and treat its core implicitly; hydrogen is light, so it stays all-electron. There is no separate "turn on ECP" keyword — choosing a -PP basis for an atom is how you request its ECP, and the potential is built automatically through libecpint.
  • rks/bhhlyp/... requests the BHHLYP hybrid, so this is a DFT single point. For a bare Hartree-Fock energy, drop the functional component: rhf/lanl2dz. The ECP mechanics are identical either way.
  • The reference is a closed-shell singlet, which rks implies — HBr has an even number of electrons. Note that the electron count is over the valence: with the Br core removed by the ECP, the SCF handles far fewer electrons than an all-electron HBr calculation would.
  • dftgrid(rad_type=becke) only affects the numerical XC quadrature for the functional; it is unrelated to the ECP.
  • No driver line means energy(). Gradients work the same way — add grad.

For a basis that is not in the built-in library, the route also accepts an explicit per-label assignment (library plus a library="..." mapping) or a JSON file (file:custom_basis-set.json); see the packaged examples/ECP/ decks.

To run the same molecule all-electron (no ECP), you would drop the -PP and give bromine an all-electron relativistic-recontracted set — but that reintroduces all 35 electrons of Br and its tight core basis functions, which is exactly the cost the ECP is there to avoid.

Python style

The equivalent calculation with the OpenQP Python API is inputs/hbr_ecp_energy.py. The ECP is requested the same way — a per-atom list passed to job.settings.basis(...), in atom order — so the -PP entry lands on Br and the all-electron entry on H.

from oqp.openqp import OpenQP

job = OpenQP("hbr_ecp_energy", silent=1)

# Geometry in Angstrom, Br then H (bond length 1.407611 Angstrom).
job.molecule(
    """
Br  0.000000000  0.000000000  0.000000000
H   0.000000000  0.000000000  1.407611
""",
    charge=0,
    multiplicity=1,
)

# Closed-shell BHHLYP (DFT uses the .dft helper; method=hf under the hood).
job.theory.dft(functional="bhhlyp")

# Per-atom basis, in atom order: ECP-carrying set on Br, all-electron on H.
job.settings.basis(["aug-cc-pVDZ-PP", "aug-cc-pVDZ"])

mol = job.run()

print("SCF energy:", mol.get_scf_energy())
print("Results:", mol.get_results())

The moving parts map one-to-one onto the input file:

  • job.theory.dft(functional="bhhlyp") is the [input] functional=bhhlyp line (and the method=hf engine underneath) — a closed-shell BHHLYP single point.
  • job.settings.basis([...]) is the basis=...;... line. Passing a Python list is the per-atom form: ["aug-cc-pVDZ-PP", "aug-cc-pVDZ"] matches the two atoms of job.molecule(...) in order. A single string (e.g. job.settings.basis("aug-cc-pVDZ-PP")) would instead apply one basis to every atom. As in the input file, the ECP comes along for free with the -PP entry — no extra call is needed.

Run it

Input-file style (from the inputs/ folder):

cd effective-core-potentials/inputs
openqp hbr_ecp_energy.oqp

Python style:

cd effective-core-potentials/inputs
python hbr_ecp_energy.py

Both need OpenQP installed (pip install openqp) and compute the identical energy. OpenQP builds the ECP integrals through libecpint, which is bundled with the standard OpenQP build — no extra install step.

Reading the output

This is an energy-only run, so the number to read is the converged SCF (BHHLYP) total energy:

  • From Python, mol.get_scf_energy() returns the total energy, and mol.get_results() gives the full results dictionary (the same fields written to <project>.json).
  • In the log file (hbr_ecp_energy.log) look for the final converged SCF energy. A quick sanity check that the ECP is active is the electron count: with bromine's core removed the SCF reports far fewer electrons than the 36 a fully all-electron HBr would have — the missing electrons are exactly Br's frozen core. If instead you see the full all-electron count, the -PP basis was not picked up (usually a per-atom ordering slip in basis=).

Because the ECP folds scalar-relativistic core effects into the valence, the HBr bond energetics from this cheap valence-only run are close to a full relativistic all-electron treatment — at a fraction of the cost.

References / manual