Local correlation with DLPNO

Canonical MP2 and coupled cluster scale steeply with system size because they correlate every occupied orbital with every other one in a delocalised basis. Domain-based local pair-natural-orbital (DLPNO) methods exploit the short-sightedness of dynamic correlation instead: occupied orbitals are localised, each orbital pair gets its own compact virtual space (the pair natural orbitals), and pairs that are far apart are estimated rather than solved. The art is that all of this is controlled by a handful of thresholds — and as the thresholds go to zero, the canonical answer comes back exactly.

This tutorial runs DLPNO-MP2 on water, shows the threshold convergence against canonical RI-MP2, and finishes with the DLPNO-CCSD(T) correctness pilot. The matching ready-to-run scripts are examples/molecular/input-h2o-dlpno-mp2.py and examples/molecular/input-h2o-dlpno-ccsd-t.py.

One-call DLPNO-MP2

import vibeqc as vq

mol = vq.Molecule([
    vq.Atom(8, [0.0,  0.00,  0.00]),
    vq.Atom(1, [0.0,  1.43, -0.98]),
    vq.Atom(1, [0.0, -1.43, -0.98]),
])

result = vq.run_job(mol, basis="cc-pvdz", method="dlpno-mp2", output="h2o_dlpno")
r = result.dlpno_mp2
print(r.e_corr, result.energy_total)

The .out file carries the decomposed solver block — iterated pair energies, the PNO-truncation correction, and the distant-pair estimate:

  DLPNO-MP2 (Pinski 2015; RI: cc-pvdz-ri)
  ------------------------------------------------------------------------------
  E(RHF reference)       =   -76.0243138804 Ha
  pairs kept / screened  = 15 / 0   (frozen core: 0)
  avg PNOs per pair      =           14.9
  E(iterated pairs)      =    -0.2004258190 Ha
  E(PNO truncation corr) =    -0.0000037342 Ha
  E(distant-pair est.)   =     0.0000000000 Ha
  E(DLPNO-MP2 corr)      =    -0.2004295532 Ha
  E(DLPNO-MP2 total)     =   -76.2247434336 Ha

Water has 15 surviving orbital pairs (5 localised valence+core occupieds → 5 diagonal + 10 off-diagonal pairs) and keeps ~15 of the 19 virtuals per pair at the default TCutPNO = 1e-8. On a molecule this small nothing is screened — the machinery earns its keep on larger systems — but the threshold behaviour is already textbook.

The experiment every DLPNO user should run once

How much correlation energy does the PNO truncation cost, and does the method really converge to the canonical answer? Sweep TCutPNO against canonical RI-MP2 with the same fitting basis:

from vibeqc.density_fitting import DensityFitting
from vibeqc.dlpno.mp2 import DLPNOMP2Options, run_dlpno_mp2

basis = vq.BasisSet(mol, "cc-pvdz")
hf = vq.run_rhf(mol, basis)

mp2_opts = vq.MP2Options()
mp2_opts.density_fit = True
mp2_opts.aux_basis = "cc-pvdz-ri"
canonical = vq.run_mp2(mol, basis, hf, mp2_opts)

df = DensityFitting(basis, vq.BasisSet(mol, "cc-pvdz-ri"),
                    aux_basis_name="cc-pvdz-ri")

for label, opts in (
    ("1e-07", DLPNOMP2Options(tcut_pno=1e-7, tcut_pno_weak=1e-6)),
    ("1e-08", DLPNOMP2Options(tcut_pno=1e-8, tcut_pno_weak=1e-7)),
    ("1e-09", DLPNOMP2Options(tcut_pno=1e-9, tcut_pno_weak=1e-8)),
    ("0*",    DLPNOMP2Options(tcut_pno=0.0, tcut_pno_weak=0.0, tcut_mkn=0.0)),
):
    r = run_dlpno_mp2(mol, basis, hf, df, opts)
    avg = sum(r.pno_per_pair.values()) / len(r.pno_per_pair)
    print(f"{label:>6s}  {r.e_corr:14.8f}  "
          f"{r.e_corr / canonical.e_correlation:9.4%}  {avg:5.1f}")

Output:

 1e-07     -0.20028030   99.7706%   14.0
 1e-08     -0.20042955   99.8450%   14.9
 1e-09     -0.20053971   99.8998%   15.7
    0*     -0.20074079  100.0000%   19.0

Three things to read off this table:

  • The error is controlled by one knob. Tightening TCutPNO by 10× buys roughly a factor-of-2 smaller truncation error at the cost of ~1 extra PNO per pair.

  • The last row is the exactness limit — no PNO truncation and full domains (tcut_mkn=0, hence the *): DLPNO-MP2 reproduces canonical RI-MP2 to better than 1 µHa. This identity is asserted permanently in tests/test_dlpno_mp2.py, on canonical and Boys-localised occupieds.

  • Recovery percentages are measured, not advertised. The same numbers appear in HANDOVER_DLPNO.md with their gates.

Frozen core and distant-pair screening

Two production options worth knowing:

opts = DLPNOMP2Options(
    n_frozen=1,       # freeze the O 1s before localisation
    tcut_pairs=1e-6,  # screen far pairs via the dipole-dipole estimate
)
result = vq.run_job(mol, basis="cc-pvdz", method="dlpno-mp2",
                    dlpno_options=opts, output="h2o_fc")

Frozen-core DLPNO-MP2 is validated against a canonical frozen-core reference to ≤ 1 µHa. Pair screening only acts on pairs whose localised-orbital centroids are ≥ 8 bohr apart (nothing in a water monomer qualifies); the screened pairs contribute through the e_distant component instead of silently vanishing.

DLPNO-CCSD(T): the correctness pilot

The same PNO machinery carries the coupled-cluster ansatz. vibe-qc’s implementation is currently a correctness pilot: the DLPNO ansatz is exact-by-construction (residuals are evaluated through an FCI-anchored spin-orbital CCSD engine and projected onto the retained PNO spaces), but the cost is O(N⁶), so jobs are capped at 64 basis functions until the reduced-scaling engine lands.

from vibeqc.dlpno.ccsd import DLPNOCCSDPilotOptions

result = vq.run_job(mol, basis="cc-pvdz", method="dlpno-ccsd(t)",
                    dlpno_ccsd_options=DLPNOCCSDPilotOptions(tcut_pno=1e-8),
                    output="h2o_cc")
cc = result.dlpno_ccsd
print(f"E_corr(CCSD) = {cc.e_corr:.8f}  (T) = {cc.e_t:.8f}  "
      f"total = {result.energy_total:.8f}")

Output (about half a minute):

E_corr(CCSD) = -0.20987650  (T) = -0.00280886  total = -76.23699924

What backs these numbers: in the untruncated limit the pilot reproduces the anchored canonical CCSD to 0.000 µHa (canonical and localised orbitals), the anchor itself sits the physically required ~0.3 % above full CI on H₂O/STO-3G, and the (T) correction vanishes identically for two-electron systems. All of it is asserted in tests/test_dlpno_ccsd.py.

When to reach for DLPNO

  • DLPNO-MP2 — whenever you would run RI-MP2 and the system is big enough that you care; defaults recover ≥ 99.8 % of E_corr, and TCutPNO=1e-9 buys 99.9 %. Energies only (no gradients yet), closed-shell references.

  • DLPNO-CCSD(T) pilot — validation-scale molecules and benchmark anchors where you want coupled-cluster quality with a verifiable error budget; not production runs (64-bf cap).

Every dlpno-* job emits its method papers (Pinski 2015, Riplinger 2013, Foster-Boys 1960, …) into the .references / .bibtex outputs automatically.

See also