DLPNO methods (MP2, CCSD, CCSD(T))¶

DLPNO-MP2 is local second-order Møller-Plesset theory in the domain-based pair-natural-orbital framework of Pinski, Riplinger, Valeev and Neese (J. Chem. Phys. 143, 034108 (2015)). Occupied orbitals are Foster-Boys localised, every occupied pair gets a compact virtual space of pair natural orbitals built inside its own projected atomic-orbital domain, and the MP2 residual equations are solved with the localised-orbital Fock coupling. The result tracks canonical RI-MP2 at a fraction of the asymptotic cost, with controllable truncation thresholds.

At a glance¶


Entry point	`run_job(method="dlpno-mp2")`
Reference	closed-shell RHF, or UHF for open shells (auto-routed to DLPNO-UMP2)
Auxiliary basis	auto-resolved correlation (“ri”) fit of the orbital basis
Options	`dlpno_options=DLPNOMP2Options(...)`
Accuracy (defaults)	≥ 99.8 % of canonical RI-MP2 E_corr (TCutPNO = 1e-8)
Exactness limit	thresholds → 0 reproduces canonical RI-MP2 to ≤ 1 µHa

Quick start¶

from vibeqc import run_job
from vibeqc._vibeqc_core import Atom, Molecule

mol = Molecule(
    [
        Atom(8, [0.000, 0.000, 0.000]),
        Atom(1, [0.000, 1.499, -1.160]),
        Atom(1, [0.000, -1.499, -1.160]),
    ],  # bohr
    charge=0,
    multiplicity=1,
)

result = run_job(mol, basis="def2-svp", method="dlpno-mp2", output="h2o")
print(result.dlpno_mp2.e_corr)   # DLPNO-MP2 correlation energy (Ha)
print(result.energy_total)       # RHF + correlation

The .out file reports the energy decomposition, iterated pair energies, the semicanonical PNO-truncation correction, and the distant-pair dipole estimate:

  DLPNO-MP2 (Pinski 2015; RI: def2-svp-rifit)
  ------------------------------------------------------------------------------
  E(RHF reference)       =   -75.9547597194 Ha
  pairs kept / screened  = 15 / 0   (frozen core: 0)
  avg PNOs per pair      =           14.9
  E(iterated pairs)      =    -0.2062668853 Ha
  E(PNO truncation corr) =    -0.0000032714 Ha
  E(distant-pair est.)   =     0.0000000000 Ha
  E(DLPNO-MP2 corr)      =    -0.2062701567 Ha
  E(DLPNO-MP2 total)     =   -76.1610298761 Ha

Thresholds and options¶

from vibeqc.dlpno.mp2 import DLPNOMP2Options

opts = DLPNOMP2Options(
    n_frozen=1,        # freeze core orbitals (before localisation)
    tcut_pno=1e-9,     # tighter PNO truncation (default 1e-8)
    tcut_pairs=1e-6,   # screen distant pairs via the dipole estimate
)
result = run_job(mol, basis="def2-svp", method="dlpno-mp2",
                 dlpno_options=opts, output="h2o_tight")

Option	Default	Meaning
`tcut_pno`	1e-8	PNO occupation cutoff (strong pairs). Tighter = more PNOs = closer to canonical. 1e-8 recovers ≈ 99.87 % of E_corr on H₂O/def2-SVP, 1e-9 ≈ 99.93 %.
`tcut_pno_weak`	1e-7	PNO cutoff for weak pairs.
`tcut_mkn`	1e-3	Mulliken threshold for domain atoms (0 = full domains).
`tcut_pairs`	1e-6	Distant-pair screening on the semicanonical dipole-dipole estimate, keeps the pair list O(N) for extended systems; screened pairs contribute through `e_distant`. Pairs closer than `dipole_r_min` (8 bohr) are never screened, so it is a no-op on compact molecules. Set `0` for the exactness configuration.
`n_frozen`	0	Frozen-core orbitals, validated against a canonical frozen-core reference.
`localise`	`"boys"`	`"none"` keeps canonical occupieds (diagnostic).
`local_df`	`False`	Domain-restricted (local) density fitting, see below.
`fit_buffer`	4.0	Fit-domain extension (bohr) when `local_df=True`.

Setting every threshold to zero is a supported validation mode: the energy then reproduces canonical RI-MP2 with the same fitting basis to ≤ 1 µHa (this is asserted in tests/test_dlpno_mp2.py).

Local density fitting (experimental, reduced-scaling)¶

With local_df=True, each pair’s exchange integrals are refit using only the auxiliary functions within fit_buffer bohr of the pair’s domain, instead of the global RI metric, the per-pair fit cost then stays bounded by the local neighbourhood as the molecule grows (the foundation of linear scaling; result.dlpno_mp2.fit_dim_per_pair records the per-pair fit dimension). A full fit domain reproduces the global RI exactly. This is the first reduced-scaling increment; it currently complements rather than replaces the global path (default local_df=False), and the end-to-end speedup arrives with the integral-direct local build (see handovers/HANDOVER_GATED_ITEMS.md, G-CORR-005).

Result object¶

result.dlpno_mp2 carries the full decomposition: e_corr = e_corr_iterated + e_pno_correction + e_distant, per-pair energies (pair_energies, absolute occupied indices), PNO counts per pair (pno_per_pair), per-pair fit dimensions (fit_dim_per_pair), screening statistics, and the iteration trace.

DLPNO-CCSD and DLPNO-CCSD(T)¶

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

from vibeqc import run_job
from vibeqc.dlpno.ccsd_local_solver import LocalCCSDOptions

result = run_job(mol, basis="def2-svp", method="dlpno-ccsd",
                 dlpno_ccsd_options=LocalCCSDOptions(tcut_pno=1e-7),
                 output="h2o_cc")
print(result.dlpno_ccsd.e_corr, result.energy_total)

Each pair couples only to occupied orbitals within coupling_radius (bohr) of either pair index, the lever that takes the occupied coupling sums from O(N⁴) to O(N²). The 12-bohr default is a bit-identical no-op on any molecule under ~12-bohr extent (the common case); on larger systems it keeps the dropped long-range coupling well below the PNO truncation error (~7 µHa on a 25-bohr H₂ chain vs the ~0.1% PNO error, the same controlled-locality bargain as the default sparse pair list). Set coupling_radius=0 for full coupling, the exact reference the full-domain ratchet pins. result.dlpno_ccsd reports avg_coupled_occ, the mean local-set size, which saturates as the system grows at fixed radius (tests/test_dlpno_ccsd_solver.py::TestSparseCoupling).

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

Larger systems: the SCF reference and memory¶

DLPNO is a reduced-scaling correlation method, but it runs on top of a mean-field SCF reference, and that reference is what sets the memory footprint on larger systems:

The DLPNO step itself is reduced-scaling (per-pair PNO domains, no dense n_occ^2 n_virt^2 or n_mo^4 tensor). n-octane / cc-pVTZ (492 basis functions) DLPNO-MP2 peaks at a few GB, not hundreds.
Use a density-fitted SCF reference for larger systems. Set density_fit=True (with an aux_basis) on the reference options. A conventional in-core 4-index SCF would materialise the n_basis^4 ERI tensor (~436 GB at 492 functions); the integral-direct SCF (the default above 200 functions) avoids it, and density fitting is faster still:
```
from vibeqc import RHFOptions, run_job
o = RHFOptions(density_fit=True, aux_basis="cc-pvtz-jkfit")
run_job(mol, basis="cc-pvtz", method="dlpno-mp2", rhf_options=o)  # ~13 GB
```

The O(N^6) pilot (DLPNOCCSDPilotOptions, opt-in) is a different animal: it is the small-system correctness oracle the local solver is validated against, and it forms dense full-virtual-space integrals (memory ~ n_mo^4, set by the basis size and independent of tcut_pno). It is hard-capped at max_nbf=64 for that reason. Do not use it for production runs; the default run_job(method="dlpno-ccsd") already uses the reduced-scaling local solver, which has no such cap.

Accuracy¶

Against canonical CCSD(T) on a 7-molecule set (def2-SVP, all-electron) the defaults (tcut_pno=1e-7 + DLPNO-(T1) + tcut_pairs=1e-4) give a mean absolute error of 0.37 kcal/mol, within chemical accuracy:

from vibeqc.dlpno.ccsd_local_solver import LocalCCSDOptions, run_local_dlpno_ccsd
r = run_local_dlpno_ccsd(mol, basis, hf, df, LocalCCSDOptions(compute_triples=True))
# MAE 0.37 kcal/mol vs canonical CCSD(T); tcut_pno=1e-8 → 0.17 (ORCA: 0.16)

Tightening to tcut_pno=1e-8 (the compiled C++ residual keeps it affordable) brings the MAE to 0.17 kcal/mol, matching ORCA 6.1’s 0.16. An earlier tcut_pno=1e-7 + DLPNO-(T0) pairing gave MAE 0.30 only via an accidental cancellation (CCSD over-recovery against the (T0) error); (T1) removes the (T0) error honestly and 1e-8 removes the over-recovery, so they target different errors and combine to the cancellation-free 0.17. triples_mode="local" selects DLPNO-(T0) for fast scans and triples_mode="exact" the O(N⁷) (T) oracle.

Pair screening (tcut_pairs=1e-4, ORCA’s TCutPairs, default-on) treats a pair whose full-virtual MP2 energy is below the threshold at MP2 level instead of CCSD: accuracy-neutral on compact molecules (a no-op on 5-occupied cases, the weak ~10-14 % tail on larger ones) and a ~2x linear-scaling speedup on extended systems; set tcut_pairs=0 to disable. The reproducible benchmark and ORCA comparison are examples/molecular/benchmark-dlpno-ccsd-t.py.

Closed-shell `(T1)`: choosing `tcut_tno`¶

LocalCCSDOptions.tcut_tno defaults to 0.0, the full TNO union span, so a default dlpno-ccsd(t) job performs no triples-domain truncation. If you opt into a positive value, these are the measured trade-offs on H2O/def2-SVP (2026-07-29; full domain e_t = -0.003009059087680666 Ha at 19.0 TNOs per triple, which is the entire virtual space):

`tcut_tno`	error vs full domain	mean TNOs	smallest domain	collapsed triples
`1e-9`	`-0.000` µHa	18.52	7	0
`1e-7`	`+0.018` µHa	18.20	4	0
`1e-6`	`+0.086` µHa	17.80	4	0
`1e-5`	`+33.39` µHa	15.68	1	1
`1e-4`	`+582.45` µHa	10.68	1	1

The closed-shell route tolerates much coarser thresholds than the open-shell engine: 1e-6 costs under 0.1 µHa here, where the open-shell (T1) costs 30 to 64 µHa at the same setting.

The error curve breaks exactly where the first triple collapses. A spatial triple excitation promotes three electrons into three distinct virtuals, so a domain retaining fewer than three cannot host one and contributes exactly zero rather than approximately. Between 1e-6 and 1e-5 one triple falls to a single virtual and the error jumps 400-fold. LocalCCSDResult reports this directly:

result = run_local_dlpno_ccsd(mol, basis, rhf, df, options)
result.avg_tno                    # mean retained domain per distinct triple
result.n_degenerate_tno_triples   # triples truncated below three virtuals
result.tno_per_triple             # per-triple domain sizes

Treat a nonzero n_degenerate_tno_triples as a signal that tcut_tno has truncated past the point where part of the correction exists, rather than as a merely coarse setting. triples_mode="exact" expands to the full virtual space and builds no TNO domain, so it reports an empty tno_per_triple.

A run_job(method="dlpno-ccsd(t)") run reports the same thing without any Python bookkeeping. The .out block carries a domain line beside the pair line, and a collapsed domain raises a warning:

  pairs / avg PNOs         = 15 / 19.0   (frozen core: 0)
  E(CCSD correlation)      =    -0.2130338714 Ha
  triples / avg TNOs       = 25 / 15.7   (smallest domain: 1)
  E((T) correction)        =    -0.0029756723 Ha
  ...
  WARNING: tcut_tno truncated 1 of 25 triples below three virtuals; those
  contribute exactly zero to (T), so part of the correction is missing rather
  than approximated. Lower tcut_tno (0 keeps the full domain).

coupling_radius is the occupied-side counterpart: it drops whole occupied triple keys rather than shrinking each domain, and it can be far more aggressive than it looks. On an H10 chain at 1.8 bohr spacing (16.2 bohr long, strongly coupled occupieds) a coupling_radius of 5 bohr drops 16 of the 25 distinct triple keys and takes 90% of the triples energy with them. That is a property of the system rather than a defect: correlation in a delocalized chain is genuinely long-ranged, so a locality cutoff is the wrong approximation there. The count is now reported so the choice is visible:

  triples / avg TNOs       = 9 / 5.0   (smallest domain: 5)   (16 screened by coupling_radius)

LocalCCSDResult.n_triple_keys and .n_triple_keys_screened carry the same numbers. Screening behaves comparably in both triples modes: measured against its own unscreened reference it costs 0.074 µHa under (T1) and 0.071 µHa under (T0) on H2O + H2 at 8 Å, and 1459 µHa against 1142 µHa on the H10 chain.

The dlpno_ccsd_done structured-log event carries n_tno_triples, avg_tno_per_triple, min_tno_per_triple, n_degenerate_tno_triples and n_triple_keys_screened alongside the existing pair fields. The open-shell route runs the UCCSD(T) pilot, which builds no TNO domain, so the line is omitted there rather than printed as zeros.

Current limitations¶

DLPNO-MP2 runs on closed-shell RHF and open-shell UHF references: run_job(method="dlpno-mp2") auto-routes a multiplicity > 1 system to UHF + DLPNO-UMP2 (spin-channel resolved αα/ββ/αβ; res.dlpno_ump2). Pass a DLPNOUMP2Options as dlpno_options to control its thresholds (localise, tcut_pno, tcut_pairs, n_frozen). DLPNO-CCSD/(T) also run on open-shell UHF references: run_job(method="dlpno-ccsd"| "dlpno-ccsd(t)") auto-routes a multiplicity > 1 system to UHF + the spin-orbital DLPNO-UCCSD(T) pilot (res.dlpno_ccsd; pass a DLPNOUCCSDPilotOptions as dlpno_ccsd_options). The open-shell path is the O(N⁶) correctness pilot (capped at max_nbf=64), not the reduced-scaling closed-shell local solver; a near-linear-scaling open-shell engine is the roadmap item. The first opt-in engine rung is available directly as run_local_dlpno_uccsd: it evaluates native UCCSD residuals in per-pair PNO/local-occupied domains and is exact against the pilot at full domains. Its default-off tcut_pairs option screens weak spin-orbital pairs using their full-virtual MP2 energies and retains those MP2 contributions; compact OH and a separated H + H2 doublet pin the exact default and the lossless distant-pair limit. It is not wired into run_job yet. A second default-off option, tcut_mkn, selects each pair’s atoms from the two occupied-orbital Mulliken populations and builds separate alpha/beta semicanonical PAOs on their union. A separated OH + H2 doublet pins the reduced pair-virtual span and exact fragment limit. The default-off tcut_pno_singles option sums each occupied orbital’s full-virtual MP2 pair-density contributions and retains same-spin singles-specific natural orbitals above the threshold, optionally inside that occupied orbital’s PAO atom domain. The zero default keeps the full spin-virtual singles space. A separated OH + H2 doublet pins both the domain reduction and its exact fragment limit. The same direct solver has an experimental, default-off open-shell DLPNO-(T0) path: set compute_triples=True and control the triple-natural-orbital occupation cutoff with tcut_tno. It streams each distinct occupied triple through the native spin-orbital kernel, reports e_t plus triple/TNO diagnostics, and reproduces the dense UCCSD(T) pilot when every threshold is zero and localise="none". This direct API is not yet a public production route. Set triples_mode="t1" to opt into spin-block semicanonical occupied coupling: it removes the (T0) localized occupied-Fock error while retaining per-triple TNO streaming. That rotation takes the occupied indices out of the localised basis, where a weak-pair label or an occupied distance has no meaning, so (T1) cannot carry a screened triple list. It is decided on the realised screened set rather than on the raw thresholds: tcut_pairs or coupling_radius values that screen no triple leave (T1) reachable and exact on the full triple list, while any setting that actually removes one fails closed with the offending count. Use triples_mode="t0" for a genuinely screened open-shell triples run. Broader radical benchmarks and public wiring remain gated.

Positive tcut_tno under (T1) is calibrated on OH/def2-SVP (2026-07-29; full domain e_t = -1.778474224055274e-03 Ha at 29.0 spin virtuals per triple). Truncation only ever loses correlation, so every entry is a positive error:

`tcut_tno`	error vs full domain	mean TNOs/triple
`1e-12`	`+11.27` µHa	26.07
`1e-9`	`+12.39` µHa	24.73
`1e-8`	`+20.05` µHa	23.49
`1e-7`	`+34.39` µHa	21.85
`1e-6`	`+63.37` µHa	19.92
`1e-5`	`+156.16` µHa	16.77
`1e-4`	`+759.06` µHa	10.55 (collapsing; see below)

1e-9, the T_CutTNO default of Guo et al. (2018), costs 12 µHa (0.008 kcal/mol) for a 15% smaller triples domain and is the recommended starting point. There is no free threshold: the smallest value tested already costs 11 µHa, so tcut_tno=0 remains the only exact setting.

The same sweep on two more radicals, including a high-spin case, gives the same picture (error vs full domain; d = triples that collapsed below three spin virtuals):

system	multiplicity	`1e-9`	`1e-6`	`1e-4`
OH	2	`+12.39` µHa	`+63.37` µHa	`+759.06` µHa, `d=13`
NH2	2	`+7.62` µHa	`+37.02` µHa	`+1177.91` µHa, `d=9`
CH2	3	`+5.73` µHa	`+30.25` µHa	`+706.18` µHa, `d=10`

Across all three, 1e-9 stays within 12 µHa and 1e-6 within 64 µHa, and the collapse sets in between 1e-6 and 1e-4 in every case. The triplet is not an outlier, so the recommendation does not depend on spin multiplicity over this (small) sample.

1e-6 is the coarsest calibrated value on every system tested. Through 1e-6 every triple keeps at least three spin virtuals (smallest domain 8), so the error is pure truncation. By 1e-4 the smallest domain is a single virtual and 13 of the 84 triples have collapsed to identically zero, so that row’s +759 µHa mixes truncation with collapse and should not be read as a truncation error.

Warning

Do not calibrate this threshold on a minimal basis. OH/STO-3G retains only 3.0 TNOs per triple at full domain, and a triple excitation needs three distinct spin virtuals, so any positive tcut_tno drops 75 of its 84 triples below three and e_t collapses from -2.29e-07 Ha to order 1e-46: the correction vanishes identically instead of degrading. The new LocalUCCSDResult.n_degenerate_tno_triples counts triples in that state; a nonzero value means the threshold has truncated past the point where the correction exists, not that it is merely small.

Energies only, no analytic gradients.
dlpno-ccsd(t) runs the local solver + local DLPNO-(T1) on the converged amplitudes (exact == canonical CCSD(T) at full domains); the (T) uses a spatial closed-shell kernel (no spin-orbital redundancy), validated to machine precision against the spin-orbital reference. The per-pair CCSD residual runs in compiled C++ when the core provides it (vibeqc::dlpno_pair_residual, bit-for-bit identical to the numpy kernel, 4-8× faster and growing with PNO size), with an automatic numpy fallback; the (T) is numpy (already BLAS-bound). The local CCSD solver couples all pairs by default (reduced-scaling) with an opt-in coupling_radius for linear scaling on extended systems. Remaining production follow-ups are tracked in handovers/HANDOVER_GATED_ITEMS.md, G-CORR-005.

Citations¶

Jobs running dlpno-mp2 emit the method papers into the .references / .bibtex outputs automatically: Møller-Plesset 1934, Feyereisen 1993 (RI), Pinski 2015 (DLPNO-MP2), and Foster-Boys 1960 (localisation).