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. The same numbers appear in handovers/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: 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).

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 + DLPNO-(T1) + tcut_pairs=1e-4) give a mean absolute error of 0.37 kcal/mol, comfortably inside chemical accuracy. 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:

Molecule	Canonical CCSD(T) (Ha)	DLPNO error (kcal/mol)	Recovery
H₂O	-76.177021	+0.083	99.94 %
NH₃	-56.357944	-0.123	100.09 %
CH₄	-40.360454	-0.670	100.56 %
HF	-100.141044	+0.128	99.90 %
CO	-112.955086	+0.014	99.99 %
N₂	-109.180282	-0.726	100.35 %
H₂CO	-114.125245	-0.841	100.39 %

Every molecule sits inside chemical accuracy at the defaults (MAE 0.369, max 0.841 kcal/mol). For gold-standard agreement, tightening to tcut_pno=1e-8 brings the MAE to 0.17 kcal/mol, matching ORCA 6.1’s 0.16 on the identical basis / all-electron setup (cross-validated out-of-process); the C++ per-pair kernel keeps that tightening affordable. The looser tcut_pno=1e-7 default trades that for speed while staying well within chemical accuracy.

Tuning the thresholds¶

tcut_pno is the main accuracy/cost knob and DLPNO-(T1) is the default triples mode (exact (T) at full domains). The history is worth knowing: an earlier tcut_pno=1e-7 + DLPNO-(T0) pairing gave MAE 0.30 only because the looser PNO caused a CCSD over-recovery that accidentally cancelled the (T0) semicanonical error. (T1) removes the semicanonical error honestly, and tcut_pno=1e-8 removes the CCSD over-recovery, so the tight + (T1) recipe is a reliable, cancellation-free MAE 0.17 at the ORCA level. The shipped default keeps tcut_pno=1e-7 for speed (MAE 0.37, still chemical accuracy).

# fast scan: looser PNO + DLPNO-(T0)
LocalCCSDOptions(tcut_pno=1e-7, triples_mode="local", compute_triples=True)
# ORCA-matching accuracy: tight PNO + DLPNO-(T1)
LocalCCSDOptions(tcut_pno=1e-8, triples_mode="t1", compute_triples=True)
# O(N⁷) exact-(T) oracle
LocalCCSDOptions(tcut_pno=1e-8, triples_mode="exact", compute_triples=True)

Pair screening (tcut_pairs, default 1e-4). A pair whose full-virtual MP2 energy is below tcut_pairs is treated at MP2 level (its exact pair energy is added directly) instead of being iterated at the CCSD level, mirroring ORCA’s TCutPairs. For a genuinely weak pair the CCSD correction is negligible and full-virtual MP2 carries no PNO-truncation error, so screening is accuracy-neutral (it is a no-op on compact 5-occupied molecules and screens only the weak ~10-14 % tail on larger ones) while removing that pair’s CCSD cost. On extended systems it is the linear-scaling lever: two waters 8 Å apart screen 45 % of pairs (all inter-fragment) for a 0.0006 kcal/mol change and a ~2x speedup. Set tcut_pairs=0 for the no-screening reference.

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 handovers/HANDOVER_DLPNO.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 RHF and open-shell UHF references (a multiplicity > 1 system auto-routes to DLPNO-UMP2, res.dlpno_ump2).
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.
Open-shell DLPNO-CCSD/(T) (multiplicity > 1) auto-routes to UHF + the spin-orbital DLPNO-UCCSD(T) pilot (res.dlpno_ccsd; pass a DLPNOUCCSDPilotOptions as dlpno_ccsd_options). This open-shell path is the O(N⁶) correctness pilot (max_nbf=64 capped), exact vs the spin-orbital run_ref_uccsd anchor at full domains; a reduced-scaling open-shell engine is a roadmap item.

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