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 reduced-scaling DLPNO-CCSD solver and the local DLPNO-CCSD(T). The matching ready-to-run scripts are input-h2o-dlpno-mp2.py, input-h2o-dlpno-ccsd.py, and input-h2o-dlpno-ccsd-t.py.

One-call DLPNO-MP2

The whole method is reachable through run_job with method="dlpno-mp2". This runs DLPNO-MP2 on water in a cc-pVDZ basis and pulls the correlation energy off the result object:

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. Remaining DLPNO accuracy and reduced-scaling gates live in HANDOVER_GATED_ITEMS.md.

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: the reduced-scaling local solver

The same PNO machinery carries the coupled-cluster ansatz. method="dlpno-ccsd" runs the reduced-scaling local solver: each occupied pair’s CCSD residual is evaluated in its own PNO basis, with amplitudes projected between pair domains, no full-system tensor is formed. It is FCI-anchored, in the full-domain limit it reproduces canonical closed-shell CCSD bit-for-bit (≤ 1 µHa; on H₂ that is FCI), and default truncation recovers ≈ 99.9 % of E_corr (tests/test_dlpno_ccsd_solver.py).

result = vq.run_job(mol, basis="cc-pvdz", method="dlpno-ccsd",
                    output="h2o_cc")
cc = result.dlpno_ccsd
print(f"E_corr(CCSD) = {cc.e_corr:.8f}  total = {result.energy_total:.8f}")

Output:

E_corr(CCSD) = -0.21017221  total = -76.23448609

(15 pairs, ~14 PNOs each from the cc-pVDZ virtual space.) Pass dlpno_ccsd_options=LocalCCSDOptions(tcut_pno=1e-9) (from vibeqc.dlpno.ccsd_local_solver) to tighten the PNO truncation toward the canonical limit, or tune coupling_radius (bohr; default 12), how far each pair’s occupied coupling set reaches, the lever that bounds the occupied coupling. It is exact at a radius wider than the molecule and a controlled approximation below the PNO truncation error once fragments separate: in a chain of H₂ molecules the average coupled-occupied count saturates while the system grows, instead of tracking it. On a water this small everything is within the default radius, so it is a no-op; the radius earns its keep on extended systems (set coupling_radius=0 for exact full coupling).

DLPNO fits the correlation integrals against an RI auxiliary basis, which the runner auto-resolves from the orbital basis (the RI: line in the log). Correlation-consistent (cc-pVnZ) and def2 orbital bases have a bundled default RI aux; minimal / Pople bases (sto-3g, 6-31g*, …) do not, and a DLPNO run on one raises NotImplementedError: No default ri-aux registered …. Supply an aux explicitly through the option object to override: dlpno_ccsd_options=LocalCCSDOptions(aux_basis="cc-pvdz-ri") (closed shell) or DLPNOUCCSDPilotOptions(aux_basis="cc-pvdz-ri") (open shell). DLPNO on a minimal orbital basis is rarely physically meaningful, but the override keeps it available for smoke tests and parity checks.

CCSD(T): DLPNO-(T1)

method="dlpno-ccsd(t)" runs the same local solver and then DLPNO-(T1) on the converged amplitudes (vibeqc.dlpno.triples_local): the perturbative triples are evaluated per occupied triple in a TNO domain, with the off-diagonal localised Fock coupling restored iteratively (Guo, Riplinger et al., J. Chem. Phys. 148, 011101 (2018)) – the same exact-(T) accuracy as canonical CCSD(T) at TNO-domain scaling. At full domains it reproduces canonical CCSD(T) to machine precision (the (T1) parity ratchet). Set triples_mode="local" for the older DLPNO-(T0) (diagonal localised Fock, ~0.1 kcal/mol looser). The (T) vanishes identically for two-electron systems.

result = vq.run_job(mol, basis="cc-pvdz", method="dlpno-ccsd(t)",
                    output="h2o_cct")
cc = result.dlpno_ccsd
print(f"E_corr(CCSD) = {cc.e_corr:.8f}  (T) = {cc.e_t:.8f}  "
      f"total = {result.energy_total:.8f}")
# E_corr(CCSD) = -0.21017219  (T) = -0.00267814  total = -76.23716421

The (T) uses a spatial closed-shell kernel (no spin-orbital redundancy). The O(N⁶) correctness pilot is still available by passing a DLPNOCCSDPilotOptions, useful as an independent benchmark anchor.

Is DLPNO-CCSD(T) accurate enough?

The exactness limit above answers the question in principle, but production runs use default thresholds, not the full-domain limit. The defaults (tcut_pno=1e-7, tcut_pairs=1e-4, DLPNO-(T1)) give a mean absolute error of about 0.37 kcal/mol, within chemical accuracy. Tightening only the PNO cutoff to tcut_pno=1e-8 gives 0.17 kcal/mol, matching ORCA 6.1’s 0.16 on the same set. The benchmark benchmark-dlpno-ccsd-t.py runs vibe-qc’s local DLPNO-CCSD(T) against its own canonical DF-CCSD(T) on a seven-molecule def2-SVP set; the table below is the tightened tcut_pno=1e-8 high-accuracy recipe:

Molecule	Canonical CCSD(T) (Ha)	DLPNO error (kcal/mol)	Recovery
H₂O	-76.177021	+0.135	99.90 %
NH₃	-56.357944	+0.046	99.99 %
CH₄	-40.360454	-0.073	100.06 %
HF	-100.141044	+0.101	99.93 %
CO	-112.955086	+0.288	99.86 %
N₂	-109.180282	+0.401	99.82 %
H₂CO	-114.125245	+0.117	99.94 %

So every molecule sits comfortably inside chemical accuracy, and the tightened option matches the gold-standard DLPNO implementation’s aggregate MAE on this benchmark (ORCA 6.1, cross-validated on the identical basis / all-electron setup).

Legacy fast-scan defaults and the high-accuracy recipe

The old fast-scan pairing was tcut_pno=1e-7 + DLPNO-(T0) (MAE 0.30). The looser PNO caused CCSD over-recovery that accidentally cancelled the (T0) semicanonical error. (T1) removes the semicanonical error at the (T0) cost, and tcut_pno=1e-8 removes the CCSD over-recovery – the combined result is a reliable, cancellation-free MAE at the ORCA level.

The old defaults are still available for fast scans:

LocalCCSDOptions(tcut_pno=1e-7, triples_mode="local", compute_triples=True)

For an even tighter, O(N⁷) exact-(T) oracle:

LocalCCSDOptions(tcut_pno=1e-8, triples_mode="exact", compute_triples=True)

The numbers in this section are the measured output of benchmark-dlpno-ccsd-t.py; re-run it to reproduce them.

The point is reduced scaling

Accuracy is only half the story; DLPNO exists because its cost grows far more slowly than canonical CCSD(T)’s O(N⁷). The companion reach demonstration in the benchmark (water clusters, same basis) makes the contrast concrete: canonical CCSD(T) becomes intractable past roughly three to four waters, while DLPNO-CCSD(T) carries on. See the reach walkthrough in benchmark-dlpno-ccsd-t.py and the remaining DLPNO entries in HANDOVER_GATED_ITEMS.md.

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, the reduced-scaling local solver: coupled-cluster quality with PNO truncation control, full-domain == canonical CCSD. coupling_radius defaults to 12 bohr, a no-op on compact molecules, bounding the occupied coupling on extended systems (set 0 for the exact full-coupling reference). The per-pair residual runs in compiled C++ when available (4-8× faster, automatic numpy fallback); it carries a generous, overridable max_nbf guard.

• DLPNO-CCSD(T), the local DLPNO-(T) on the converged local amplitudes (per-triple TNO domains), exact == canonical CCSD(T) at full domains. The O(N⁶) pilot stays available as the exact-(T) oracle (DLPNOCCSDPilotOptions, 64-bf cap) for benchmark anchors.

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

See also

• DLPNO methods user guide, full threshold table, result-object reference, limitations.

• MP2 and double hybrids, the canonical methods DLPNO converges to.

• Tutorial 38, RI-MP2 auxiliary-basis verification

  • the RI fitting error underneath everything here.