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.

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.

  • 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).

See also