Szczepanik Research Group

Department of Theoretical Chemistry
Faculty of Chemistry, Jagiellonian University
Gronostajowa 2, 30-387 Krakow, Poland
Tel: (+48) 12 686 23 90
E-mail: dariusz.szczepanik@uj.edu.pl

◼ Manual (2026-Oct-01)

1. Introduction

runEDDB implements the Electron Density of Delocalized Bonds (EDDB) method and the corresponding electron-density partitioning schemes. The present release is a complete rewrite in Fortran 2008 of the original runEDDB R-script. Thanks to a modern linear-algebra backbone (BLAS-3 DGEMM, DSYGVD, DSYEVD, DPOSV, CGS2, etc.), the program now processes systems containing thousands of atoms on a single workstation in seconds to minutes. It is distributed as ready-to-run, pre-compiled binaries for all three major desktop platforms — Linux (x86_64), Windows (x86_64), and macOS (arm64). See the Download page for the current release and per-platform installation instructions.

The Bond-Order Projection (BOP) algorithm at the heart of EDDB requires the input wavefunction to be expressed in a basis of compact, valence-localized atomic orbitals with a hierarchical, L-resolved subshell orthogonalization — the Natural Atomic Orbital (NAO) basis. The NAO transformation maps the underlying Gaussian basis onto an atomic-orbital basis with a clean core (Cor) / valence (Val) / Rydberg (Ryd) hierarchy. Several implementations of the original definition are in routine use today (e.g., NBO7, JaNPA); runEDDB ships with its own internal implementation alongside two compact subsets — the Natural Minimal Basis (NMB = Cor ∪ Val) and the program-default Natural Valence Basis (NVB = Val).

2. Generating Input Files

The current version of runEDDB directly reads seven input formats — six wavefunction formats produced by external quantum-chemistry programs and one geometry-only XYZ format processed by the built-in GFN2-xTB engine:

• .fchk — Gaussian Formatted Checkpoint File generated by Gaussian, Q-Chem and Multiwfn.
• .json — JavaScript Object Notation generated by ORCA.
• .molden — Molden file generated by Turbomole, GAMESS, MOLCAS, PySCF, Multiwfn, and others.
• .tmol — native Turbomole bundle (concatenation of control coord basis mos from the Turbomole working directory).
• .gms — GAMESS-US output log file.
• .49 — archive file generated by Natural Bond Orbital 7.0.1 (NBO7) program.
• .xyz — plain XYZ geometry; the wavefunction is computed internally by the built-in GFN2-xTB engine (no external QC software required).

FCHK (.fchk)

The most common input format. Supported routes:

• Gaussian — add the formcheck keyword to the route section of the input; the .fchk is then written automatically alongside the binary checkpoint. (Equivalently, run formchk molecule.chk molecule.fchk after a regular job.)
• Q-Chem — add gui = 2 to the $rem block; launch Q-Chem with -save SAVENAME so the FCHK in the scratch directory is preserved (it is written as $QCSCRATCH/SAVENAME/Test.FChk).
• Multiwfn — can convert virtually any wavefunction format into .fchk as a universal bridge.

Most of our FCHK testing focused on Gaussian and Q-Chem outputs. For post-HF methods the relaxed correlated density block has to be requested explicitly — in Gaussian via density=mp2 (MP2) or density=current (CCSD, CI, CASSCF) on the route line; in Q-Chem via cc_ref_prop = true for the CC family. (Q-Chem does not currently emit the MP2 density to FCHK; for MP2 use Gaussian or ORCA.)

Molden (.molden)

Molden is one of the most popular formats for molecular density and orbital visualization. Many quantum-chemistry programs can produce it through different mechanisms — please consult the manual of your particular package for the exact keyword. Multiwfn is also a convenient universal converter to and from Molden, and runEDDB has been validated against Multiwfn-generated files.

ORCA (.json or .molden)

For ORCA we strongly recommend using .json rather than .molden. The JSON format produced by ORCA carries a more complete description of the wavefunction — including ECP definitions and other metadata — that the Molden specification cannot represent. We have validated runEDDB against both formats, generated by ORCA's standard orca_2mkl and orca_2json post-processors.

Turbomole native (.tmol)

For Turbomole calculations, runEDDB reads a single .tmol bundle that you assemble by concatenating the four key files from the Turbomole working directory:

$ cat control coord basis mos > molecule.tmol

The parser handles RHF, UHF, ROHF, and DFT wavefunctions in Turbomole's spherical basis representation, with angular momenta up to k (L = 7) — covering def2-{S,T,Q}ZVPP and cc-pVDZ through cc-pV8Z. The native Turbomole $ecp block is parsed automatically — ECP-bearing transition-metal calculations load directly without conversion. Cartesian d/f basis storage (forced by $pople CAO or $cao) is not supported — use Turbomole's default spherical mode.

C₁ symmetry is required for the .tmol lane. Under any non-trivial point group Turbomole writes the wavefunction in two ways that need irrep-aware bookkeeping not present in the .tmol parser:

• $scfmo / $uhfmo_* — each MO header carries an nsaos=<N> count of coefficients to read for that orbital. Under non-C₁ nsaos < n_ao and only the AOs of the MO's irrep are written; the parser unconditionally expects n_ao coefficients per MO and would build a scrambled MO matrix.
• $closed shells / $alpha shells / $beta shells — orbital ranges are irrep-relative (e.g. a1 1-10 and b2 1-5 reference different orbital slots); the parser collapses both onto the same global index range and corrupts occupations.

The fix is simple: re-run the SCF with $symmetry c1 in the control file (or remove the keyword). Population analysis with EDDB doesn't benefit from point-group symmetry being on during the SCF write-out, so this is the standard workflow. If a non-C₁ .tmol is passed, the parser stops with an explanatory ERROR; the run does not silently produce wrong densities. The other input formats (.fchk, .json, .molden, .49, .gms) write full-AO MOs and global orbital ordering regardless of point group, so they have no equivalent restriction.

GAMESS-US (.gms)

GAMESS-US output files (the standard text log produced by rungms) are read directly. Just rename or symlink the log to a .gms extension — PUNCH (.dat) files are not used. Supported run types: RUNTYP=ENERGY, OPTIMIZE (final equilibrium geometry is auto-located), and SADPOINT (transition state). Supported wavefunction types: RHF, UHF, ROHF, and MCSCF natural orbitals. Spherical-harmonic basis (ISPHER=1) is recommended; standard Cartesian basis is also supported. The default GAMESS print precision (F11.6) is sufficient; for high-precision basis-set studies, set DGRID=.T. in $CONTRL for F17.10 MO output.

GAMESS $ECP blocks are parsed automatically (handles both WITH ZCORE n AND LMAX l and the shared-definition shorthand ARE THE SAME AS ATOM N). Symmetry-truncated basis sections (any GAMESS run with COORD=UNIQUE under C_s, C_2v, D_2h, etc.) are reconstructed automatically — the parser detects the symmetry-implied gaps in the printed basis section and replicates each per-atom block for its unprinted symmetry-equivalent partners.

NBO archive (.49) — legacy compatibility

NBO7 .49 archives are still accepted in NAO-only mode (population analyses + BOP only; cube grids and orbital FCHK exports are disabled because NBO uses a different AO-basis convention than runEDDB's native pipeline). The .49 path is provided for backward compatibility and will be removed in a future release once the native NAO driver completes its acceptance suite. Users are strongly encouraged to switch to .fchk / .json / .molden / .tmol / .gms as soon as practical.

XYZ geometry (.xyz) — built-in GFN2-xTB engine

One of runEDDB’s most distinctive capabilities is that it can analyse chemical bonding straight from a bare set of atomic coordinates — for example a structure solved by X-ray diffraction, a molecular-dynamics snapshot, or a hand-built model — with no separate quantum-chemistry calculation. Give it a plain .xyz geometry and it builds the electronic structure itself, on the fly, using a fast built-in GFN2-xTB engine (a semi-empirical tight-binding model; Bannwarth, Ehlert & Grimme, JCTC 15 (2019) 1652), and then runs the complete EDDB analysis on top of it.

The practical payoff: in a single command, and in seconds rather than the hours a full DFT study would take, you obtain a map of where the electrons delocalize — a direct picture of the molecule’s resonance / conjugation pattern and bonding skeleton — without any external software, wavefunction file, or DFT step. This makes runEDDB a quick, stand-alone way to look at the bonding in a freshly determined or modelled structure; promising findings can then be confirmed with a higher-level wavefunction (DFT, post-HF, …) run through one of the other input lanes — again analysed by runEDDB — for a fully validated picture. The entire XYZ pipeline (GFN2-xTB SCC in its minimal valence vSTO-nG basis → AO→NAO localisation → BOP and EDDB analysis) runs end-to-end on the plain XYZ file, and the internal engine is validated against the reference tblite implementation across 81 systems covering 26 elements (H…W) at chemical accuracy.

About the GFN2-xTB basis (vSTO-nG). GFN2-xTB is a minimal valence tight-binding method: its atomic-orbital basis has just one Slater-type orbital per valence shell (e.g. one s and one p for C, N, O; one s for H; with a d shell added for heavier elements), and each Slater orbital is expanded as a short contraction of n Gaussians (an STO-nG fit). We abbreviate this vSTO-nG = valence STO-nG. It is therefore a compact, single-zeta-quality valence basis — not a double-zeta or polarized basis — which is exactly why an NAO/NBO-style decomposition of a GFN2 wavefunction looks minimal-basis-like. This keeps GFN2-xTB fast and is appropriate for the qualitative bonding screening the .xyz lane is meant for; for quantitative work, validate with a larger-basis DFT / wavefunction calculation through one of the other input lanes.

The XYZ file uses the standard XYZ format, with a runEDDB-specific optional extension on line 1 that carries the molecular charge and spin multiplicity:

N_atoms [ charge multiplicity ]
comment line
Sym1   x1   y1   z1
Sym2   x2   y2   z2
...

The optional charge and multiplicity tokens on line 1 (after the atom count) default to 0 and 1 respectively, so plain XYZ files (with only N_atoms on the first line) are processed as neutral closed-shell singlets by default. Coordinates are in Ångström. Total energies reported in the verbose output of the GFN2-xTB stage match Grimme’s reference xtb executable to within ~10⁻⁴ E_h, allowing direct sanity-checking against an external xtb run on the same geometry.

Trajectories and geometry alignment (--xyz-*)

A .xyz file may hold a trajectory — many structures concatenated one after another (e.g. an MD run or an optimisation path). By default runEDDB analyses the first structure and prints a note reporting how many it found; use --xyz-frame N to pick the N-th structure instead (1-based). Requesting a frame that does not exist is a clear error reporting the actual count.

Two further options rigid-translate the input geometry (bond lengths, angles and dihedrals are untouched, so the EDDB analysis is unchanged — only the absolute coordinates written to any exported cube / FCHK move). They apply to every input format, not just .xyz (legacy .49 excepted):

• --xyz-center [geom|mass] — shift the molecule so the origin sits at its geometric centroid (geom, the default if no mode is given) or its centre of mass (mass). Handy for tidying up a raw X-ray geometry. In verbose mode the applied translation vector is printed.
• --xyz-translate x,y,z — apply an explicit translation vector (in Ångström, comma-separated) to every atom. Mutually exclusive with --xyz-center.

Trajectory tip — comparable cubes across frames. Centering each frame independently with --xyz-center would re-centre every structure on its own (drifting) centroid, so the resulting cubes would not line up. Instead, centre one reference frame with --xyz-center geom -v, read off the printed translation vector, then run every frame with that same fixed vector via --xyz-translate (optionally combined with --cube-ref reference.cube to clone the grid) — the frames then share one coordinate frame and the EDDB cubes are directly comparable / subtractable.

Both closed- and open-shell systems are supported. Closed-shell molecules (the default, multiplicity = 1) are treated as restricted; any higher multiplicity is solved as a true spin-polarised unrestricted (UHF) calculation (spin-polarized GFN2-xTB / spGFN2-xTB; see the method references in §11) — just put the spin multiplicity on line 1 (e.g. 2 for a doublet radical, 3 for a triplet, 4 for a quartet, …).

For any open-shell run, runEDDB prints a spin-density topology table right after the GFN2-xTB step: the per-atom spin population, i.e. where the unpaired electrons actually sit. This is a quick way to read off the magnetic / radical centres of a structure — and a practical screen for diradical / broken-symmetry character: run the molecule as a triplet and, if the spin localises on two well-separated atoms or fragments, that is a strong hint of an open-shell singlet (diradical) ground state worth investigating. The NO analysis (--table no for the table, --export-fchk no for the orbitals) adds complementary diradical / multireference diagnostics — including the ⟨S²⟩ spin-contamination value, which counts roughly how many electron pairs are genuinely unpaired.

Diradicals / broken-symmetry singlets. GFN2-xTB’s simplified spin treatment cannot reliably converge a genuine broken-symmetry singlet — it tends to collapse back to the closed-shell solution or stop at partial open-shell character. Use the triplet spin-topology read-out above as a qualitative pointer only, and confirm any suspected diradical / broken-symmetry singlet with a proper broken-symmetry DFT calculation (e.g. CAM-B3LYP), which can then be analysed by runEDDB through one of the wavefunction input lanes.

Built-in model library (-i @NAME)

For quick tests, teaching, or a first look at a new option you do not need any input file at all: runEDDB carries a small library of molecular models compiled into the binary. A model is loaded exactly like an .xyz geometry — the wavefunction is computed by the built-in GFN2-xTB engine — by giving its name with a leading @ instead of a file name. The library holds 38 entries in five groups: arenes (benzene, naphthalene, anthracene, phenanthrene, pyrene, coronene, kekulene, corannulene, …), hetero (pyrrole, furan, thiophene, pyridine, borazine), rings (cyclobutadiene, cyclopentadienyl anion, tropylium, cyclooctatetraene and its dianion), macro (porphine, phthalocyanine, [18]annulene, cycloparaphenylenes) and gmtkn55 (the ten C₆₀ isomers of the GMTKN55 set). -i @list prints the full catalogue with atom counts and the geometry level of every entry:

$ runEDDB -i @benzene -q                         # bare name, resolved across the library
$ runEDDB -i @arenes/benzene --nics-rings 6      # full path (group/name)
$ runEDDB -i @pyrrole -p 1-2-3-4-5-1 -a          # N1, C2–C5: ring atoms are consecutive
$ runEDDB -i @gmtkn55/c60iso/5 -r 5,6 -q         # dataset entry, all 5- and 6-rings
$ runEDDB -i @list                               # print the catalogue

  Built-in molecular model library  (load with: runEDDB -i @NAME)

  --------------------------------------------------------------------------
  arenes
  --------------------------------------------------------------------------
    @arenes/anthracene         (24 atoms)
        MMFF94 conformer (starting-quality geometry)
    @arenes/azulene            (18 atoms)
        MMFF94 conformer (starting-quality geometry)
    @arenes/benzene            (12 atoms)
        MMFF94 conformer (starting-quality geometry)
    ...

Names are resolved in two ways: the full path @group/name is always unique, and the bare name (@benzene) works whenever it occurs only once in the library. With -v the header of the run reports the model name and its geometry level. Charge and multiplicity are part of the model (e.g. @rings/cp-anion is loaded as a singlet anion), so no --xyz-charge is needed. Windows PowerShell users: quote the name, -i '@benzene', because a bare @ is PowerShell’s splatting operator; in cmd.exe no quotes are needed.

Atom numbering. In every monocyclic and polycyclic model the ring atoms are numbered consecutively around the ring, starting at the heteroatom when there is one (pyrrole: N1, C2–C5; pyridine: N1, C2–C6), and in fused systems along the perimeter (naphthalene: 1–8 around the rim, 9 and 10 the fusion atoms); hydrogens come last, in the order of the heavy atoms they are bonded to. A pathway such as -p 1-2-3-4-5-6-1 therefore means what it says, and for fused rings -r writes the exact strings to the .rings file.

MDL Molfile (.mol) and SDF (.sdf)

The MDL Molfile (.mol) and its multi-record container SDF (.sdf) are a second geometry-only lane: like .xyz they carry just a structure (and, optionally, formal charges / radical flags), so runEDDB builds the electronic structure with the same built-in GFN2-xTB engine and runs the full EDDB analysis on top — everything said above about the GFN2-xTB lane (the vSTO-nG basis, the qualitative-screening caveat, open-shell handling, the geometry options below) applies unchanged. .mol / .sdf are ubiquitous in cheminformatics (PubChem, ChEMBL, RDKit, drawing programs), so this lets you analyse a structure straight from those sources.

runEDDB reads the standard V2000 connection table. Coordinates are in Ångström; the bond block is ignored (GFN2-xTB rebuilds connectivity itself). Molecular charge and spin are taken from the property block: total charge is summed from the M  CHG entries (falling back to the legacy per-atom charge column if absent), and multiplicity from M  RAD radical flags (doublet → 1, triplet → 2 unpaired electrons); when no radical flag is present the multiplicity follows the electron-count parity (even → singlet, odd → doublet). Two cases are refused with a clear message: the free-form V3000 format (re-export as V2000), and files explicitly flagged 2D (the dimensional code on line 2) — GFN2-xTB needs real 3D coordinates.

An .sdf file may hold many structures (conformer sets, compound libraries), each record separated by a $$$$ line. As with multi-frame .xyz (see above), runEDDB analyses the first record by default and prints a note with the count; --xyz-frame N selects the N-th record. The --xyz-center / --xyz-translate geometry options apply here too.

Tripos MOL2 (.mol2)

The Tripos .mol2 format is another geometry-only lane feeding the same GFN2-xTB engine. runEDDB reads the @<TRIPOS>MOLECULE atom count and the @<TRIPOS>ATOM records; each element is taken from the SYBYL atom type (C.3 → C, N.ar → N, or a bare symbol such as Li), and the bond block is ignored. Coordinates are in Ångström. Multiple @<TRIPOS>MOLECULE records are treated as a multi-structure file (--xyz-frame N selects one). MOL2 carries only per-atom partial charges, not a formal total charge or spin, so runEDDB defaults to a neutral closed-shell singlet; set a different state with --xyz-charge / --xyz-mult (below).

Protein Data Bank (.pdb / .ent)

The PDB lane reads ATOM / HETATM records (coordinates in Ångström). The element is taken from the dedicated element columns (77–78) when present — as written by Open Babel, PyMOL, etc. — and otherwise guessed from the atom name. Only the first alternate location (altLoc blank or A) of any disordered atom is kept. A multi-MODEL file (NMR ensembles, MD frames) is multi-structure: the first model is used by default, --xyz-frame N selects another. PDB carries no charge or spin, so the default is a neutral closed-shell singlet (override with --xyz-charge / --xyz-mult).

Hydrogens. Crystallographic PDB structures often omit hydrogens. GFN2-xTB needs a complete all-atom geometry, so a carbon-containing PDB with no hydrogens triggers a warning — add hydrogens first (e.g. obabel in.pdb -O out.pdb -h, or in PyMOL / Avogadro). Pure inorganic / metal structures that legitimately have no hydrogens are not flagged.

Overriding charge and multiplicity (--xyz-charge / --xyz-mult)

For any geometry input (.xyz, .mol, .sdf, .mol2, .pdb) the molecular charge and spin multiplicity can be set on the command line with --xyz-charge N (a signed integer) and --xyz-mult M (2S+1, a positive integer). When given they take precedence over whatever the file specifies (the .xyz header, or the M  CHG / M  RAD records of a molfile). For .mol2 and .pdb, which carry no usable charge / spin information (default: neutral singlet), these flags are the way to study an ion or an open-shell state — e.g. --xyz-charge -1 --xyz-mult 2 for a radical anion. They are rejected for the wavefunction lanes (.fchk / .json / .molden / .tmol / .gms / .49), where the charge and spin come from the wavefunction itself.

Crystallographic data (.cif) — convert to .xyz first

Crystallographic Information Files (.cif) are periodic data — a unit cell, an asymmetric unit (often only a fragment of a molecule, sitting on a symmetry element), and symmetry operators — not a ready molecular geometry. Turning that into a single, complete, chemically-sensible molecule (apply symmetry, grow the molecule across cell boundaries, resolve disorder and partial occupancies) is exactly what dedicated crystallography programs already do well. So rather than re-implement that, the recommended route is to export a complete molecule to .xyz from the tool you already use, then analyse that .xyz with runEDDB:

• Mercury (CCDC), Olex2, VESTA or Diamond: open the CIF, “grow” / complete the molecule, then File → Save/Export as .xyz (or .mol / .pdb).
• Open Babel for simple molecular cells: obabel structure.cif -O structure.xyz --fillUC strict — note this fills the unit cell rather than extracting a single molecule, so it suits discrete molecular crystals but not coordination polymers, framework solids or ionic lattices.

Once you have the molecule as .xyz, the full GFN2-xTB → EDDB pipeline runs on it directly (add hydrogens first if the structure lacks them, and set the charge / multiplicity with --xyz-charge / --xyz-mult if needed). This keeps the resonance / delocalization analysis one short step away from raw crystallographic data without duplicating the symmetry-and-disorder machinery your crystallography software already provides.

Effective Core Potentials (ECP)

Heavy elements are routinely treated with Effective Core Potentials (ECPs), which replace the chemically inert inner shells with an analytical pseudopotential. runEDDB consumes already-converged MOs and never evaluates ECP integrals itself; what each parser must do is read the per-atom effective nuclear charge Z_eff = Z_atomic − N_{core,replaced} so that atomic charges, total molecular charge, and the per-atom population tables balance. The NAO classification (Core / Val / Rydberg) also uses each atom's Z to decide how many core/valence shells to expect, so getting Z_eff right is essential — a wrong value silently misclassifies orbitals and can drop the heavy atom from the per-atom tables.

How each parser obtains Z_eff:

• Gaussian fchk: from the "Nuclear charges" record (correct since day one).
• Q-Chem fchk: Q-Chem's GUI=2 writer puts the full atomic numbers in the "Nuclear charges" record even when an ECP is in use, in violation of Gaussian's convention. runEDDB detects this inconsistency automatically (the identity sum(Nuclear charges) = N_electrons + Charge fails) and stops with an ERROR, showing the exact wrong/correct line on a worked PtF₄ example. The simplest fix is to hand-edit the "Nuclear charges" record so that each ECP'd entry holds Z_eff; alternatively rerun the calculation in Gaussian or ORCA.
• ORCA json: from the "NuclearCharge" field per atom.
• Molden: from the ORCA-extension [Pseudo] block (ORCA's molden writer emits this natively whenever an ECP is in use). When the file lacks [Pseudo] but the per-atom Mulliken charge of any atom exceeds the element's max formal oxidation state (so an ECP must be in use), runEDDB stops with an ERROR showing the offending atom(s) and the exact [Pseudo] block syntax to add by hand — insert it anywhere between [Atoms] and [GTO]. If the system really is genuinely highly-ionized and has no ECP, pass --no-ecp to skip the check. Format:
   [Pseudo] block format

  [Pseudo]
  W    1  14         ! atom 1 (W),  Z_eff = 14 (def2-ECP60)
  Pt   5  18         ! atom 5 (Pt), Z_eff = 18 (def2-ECP60)

• Turbomole tmol: from the $ecp section (per-element ncore = N lines).
• GAMESS gms: from the ECP POTENTIALS echo block (per-atom ZCORE values, with ARE THE SAME AS ATOM N shorthand resolved).
• NBO7 .49: NBO7 doesn't write ECP metadata, so the parser infers N_{core,replaced} by comparing the surviving NMB shells per (atom, L) against the aufbau prediction for the bare Z.

Once Z_eff is known, the NAO machinery automatically labels the surviving valence shells with the chemically correct principal quantum numbers. For example, tungsten under def2-ECP60 shows its valence as 6s, 5d₁..5d₅ (the actually-present W valence shells) rather than the misleading 1s, 3d₁..3d₅ labels that an ECP-blind labelling would produce.

Supported wavefunctions:

• RHF, UHF, ROHF, DFT — all standard methods are supported. The density matrix is idempotent and the standard BOP algorithm is used.
• Post-HF (MP2, CCSD, CASSCF, CI, …) — supported via the layer-by-layer BOP algorithm. The density matrix has fractional natural-orbital occupations, which are decomposed into cumulative idempotent layers. This recovers the standard HF/DFT result exactly when the input density is single-determinant. Requires the QC program to write the relaxed correlated density. Note that a relaxed (response) MP2 density is not strictly N-representable: a few natural occupations fall slightly below 0 or above 2 (benzene/cc-pVTZ: −0.00001 for the singlet, down to −0.02 per spin for the spin-contaminated UMP2 triplet). runEDDB repairs them by default (since the 2026-Oct-01 compilation): the natural-occupation spectrum is replaced by its Euclidean projection onto the admissible set {0 ≤ n_k ≤ 2 (1 per spin), Σn_k = N}, i.e. n_k = clip(n_k − μ, 0, 2) with one uniform shift μ fixed by the electron count (the nearest N-representable density in the Frobenius norm; the natural orbitals are unchanged, and N is the nearest integer to the clamped trace). The same projection is applied to HF/DFT densities that fail the check, so one rule holds throughout the program; under -v the NAO block prints the corrected population, the population and occupation range before the correction, μ and the largest change, and the BOP block prints the layer-grouping residual and the occupation of the dropped natural orbitals. To analyse the density exactly as read, add --no-nrep: the layer scheme then discards the tail of orbitals whose occupation is below 10⁻⁶ (negative ones included) and keeps the rest as they are, so the effective electron count can drift from the nominal one, and every affected run is flagged by a counted WARNING (the NBA stage prints the count of violations, the extreme occupation and the input → effective trace per spin channel); --force-nrep is kept as a legacy alias of the default. The correction is a preprocessing convention, not a recovery of the exact correlated density. State the correction in the methods when you use it. Alternatively request an unrelaxed MP2 density (ORCA) or a wavefunction-based correlated density (CASSCF, CI, CC), which do not have this problem.

3. Analysis Modes

runEDDB offers five scope modes, two NOBD-dissection options, and three analysis-basis subsets.

Scope modes (mutually exclusive):

• --global, -g — EDDB_G: entire molecular system, all atom pairs above the Wiberg bond-order threshold (default).
• --heavy-atoms, -h — EDDB_H: entire molecular system as in EDDB_G, but the connectivity / BOP mask is built and then every row and column belonging to a hydrogen is zeroed — i.e. all bonds to H are excluded. The full wavefunction is still used; the per-atom table, the summary totals (Total ED_H / EDDB_H) and the cube integral cover the heavy atoms only. Speeds up the delocalization step on typical π-conjugated systems and, by removing C–H hyperconjugation, makes π-NOBD selection cleaner (less σ/π mixing). The natural sibling of EDDB_G — and the scope recommended as default in the legacy runEDDB.R. Note: -h now selects this scope, not help — help is --help / --full-help (or a no-argument run).
• --fragment, -f LIST — EDDB_F: only bonds within the listed atoms (e.g., -f 1:6 or -f 1,2,3,4,5,6). The full molecular wavefunction is still used; the per-atom table, the summary totals (Total ED_F / EDDB_F), and the cube integral all restrict to the listed atoms.
• --pathway, -p CHAIN — EDDB_P: only the explicitly listed bonds (e.g., -p 1-2-3-4-5-6-1); repeat the first atom to close a cycle.
• --ring-detect, -r LIST — EDDB_R: automatically detect every n-membered ring of the requested size(s) in the Wiberg bond-order graph (e.g., -r 6 or -r 5,6:8) and run a Cyclic Delocalization Index (CDI) on each — an automatic -p sweep over the whole molecule. The complete ring list is written to <stem>.rings (one -p string per line, ready to feed back in); the Wiberg ring-bond cutoff is tunable with --ring-thresh (default 0.30).

Under a restricted scope (-f / -p / -r) every scope-dependent density — EDDB, EDLB and BDF — is built from density matrices with exactly the atomic support of the selected atoms, so the per-atom tables, the summary totals and the cube integrals are mutually consistent (see the note in §5).

NOBD dissection (can combine with any scope):

• --auto-pi, -a — auto-select π-type NOBDs to recover the π component of EDDB.
• --select-nobd, -s [LIST] — manual selection (e.g., -s 1:3,7). Without a LIST, runEDDB only prints the NOBD-resolution table so you can pick indices for a follow-up run.

When dissection is active, the atom table and summary show only the density of the selected NOBDs; the star marker (*) in the NOBD table flags them.

Analysis basis subset, selected by --basis VALUE (or -b VALUE). VALUE is one of:

• nvb — Natural Valence Basis (Val only). Default and fastest for native input (FCHK / JSON / Molden); matches modern aromaticity-study conventions.
• nmb — Natural Minimal Basis (Cor + Val). Legacy default of the NBO-based runEDDB; also the effective default for .49 input (NBO7 archive), since the archive carries no Cor/Val split.
• nao — Full NAO basis (Cor + Val + Ryd). Most complete and slowest; useful as a completeness check — for normal valence bases nvb / nmb agree with nao to within a few percent.

VALUE is fully case-insensitive (Nao, NMB, nVb, etc. all work).

.49 input note. The NBO7 .49 archive flags every NAO as either NMB or Rydberg, with no further Cor/Val partition. -b nvb therefore cannot be honoured on .49 input and silently falls back to -b nmb. When -b nvb is requested explicitly, a one-line clarifier — (.49 archives carry no Cor / Val split) — is printed directly under the Target NAO representation: entry in the Stage 2 detail block, aligned with the basis-label value. The clarifier is verbose-only (-v), keeping the default-mode output clean. To get a real Cor/Val split, use the matching native input (FCHK / JSON / Molden).

Basis-independent connectivity (always on). The atom-pair connectivity that BOP analyses is always derived from Wiberg bond orders evaluated in the full NAO representation (Cor + Val + Ryd), whatever -b selects; that one mask is then handed unchanged to a BOP run that uses the -b working density. Consequently a nao / nmb / nvb series differs only in the density BOP integrates, on a fixed set of atom pairs. Why it matters: the Rydberg shell carries the weak through-space and hyperconjugative couplings that sit near the 0.001 gate, so deriving the mask in a truncated basis silently prunes them — and a change in EDDB population would then mix a genuine basis effect with a change in which bonds were counted at all. The bond-order threshold, the BOP equations and the working density matrix are untouched; only the pair selection is basis-fixed. Notes: with -f the NAO connectivity is built first and the fragment restriction applied to it afterwards; -p supplies the connectivity explicitly and therefore overrides it entirely, while -r auto-detects its rings from the NAO bond-order graph. Under -v the Stage 3 detail block reports the basis actually used, as Connectivity mask basis: full NAO. Cost is negligible — one extra Wiberg pass over the NAO density; BOP itself, which dominates the runtime, still runs at nmb / nvb size. (This was briefly the opt-in --nao-mask flag in July 2026; it is now unconditional and the flag has been removed.)

Post-HF layer basis (always full NAO). For correlated (post-HF) densities the BOP analysis uses the cumulative-idempotent layer scheme: the 1-RDM is partitioned into quasi-degenerate (QD) groups of natural orbitals, and a standard BOP is run on each cumulative idempotent projector. That QD partition is always derived from the full NAO density, never from the -b sub-block. This is essential rather than cosmetic: the sub-block of an idempotent density is not idempotent (Rydberg projection removes a fraction of each MO), so grouping the sub-block reads basis truncation as correlation and shatters what should be a single layer into many, inflating EDDB. BOP, the working density and every threshold remain in the -b basis — only the layer partition is derived in full NAO. Under -v the Stage 3 block reports QD layer basis: full NAO. For an idempotent reference all occupied natural orbitals share one occupation, so there is exactly one layer and the scheme reduces to a single standard BOP call — recovering the HF/DFT answer exactly, as it must.

Bond-order thresholds. runEDDB applies two distinct Wiberg cutoffs, and only the second is user-tunable:

• The BOP pair gate — --wbo-thresh X (or -w X), default 0.001. An atom pair enters the BOP analysis when its Wiberg bond order W_AB exceeds this cutoff. The default is deliberately very loose and is not a covalent-connectivity criterion: at 0.001 the mask captures essentially every chemically meaningful interaction, including hyperconjugative donor–acceptor pairs, weak hydrogen bonds, π-stacking contacts and through-space couplings. Treat this as a study / diagnostic knob, not a tuning parameter. EDDB populations are only comparable between runs computed at the same gate, so a value other than 0.001 makes your numbers incommensurable with default-threshold runs and with every published EDDB value; under -v a non-default gate is flagged explicitly as (--wbo-thresh; NOT the default 0.00100) next to the resulting Atoms / Pairs included in BOP counts. X ≥ 0 (values above 2.0 are rejected, as no pair would survive). -w 0 is legal and keeps every pair with any non-zero coupling — i.e. the complete graph, N(N−1)/2 pairs. That is a useful completeness check on small systems, but note the cost grows as N² and BOP dominates the runtime: benzene goes from 51 pairs to 66, but a 59-atom RNA fragment goes from 260 to 1711. On anything sizeable, pair it with -f / -p / -r. The --table bo / NBA gate tracks this cutoff at ten times its value, per the method. Note that the mask is already basis-independent by construction (see Basis-independent connectivity above), so there is no need to reach for -w 0 to make a -b series comparable.
• The ring-graph cutoff — --ring-thresh X, default 0.30. Used only by --ring-detect / -r, to decide which bonds form the graph that rings are searched in. Here a covalent-grade cutoff is wanted — at the 0.001 BOP gate almost every atom pair would look “bonded” and ring perception would be meaningless. Raise it to suppress spurious rings closed through weak contacts, lower it (e.g. 0.15) to pick up rings closed by a genuinely weak bond. It requires -r and does not affect the BOP pair gate, the EDDB values, or anything else.

NAO construction algorithm, selected by --nao-algorithm ALG (native input only; .49 archives carry NBO7’s own NAOs and ignore the flag with a printed note). Background: in strongly polarised basis sets (cc-pVQZ and beyond) the textbook NAO definition suffers from High-L Redundancy Leakage (HLRL) — the molecule-wide occupancy-weighted orthogonalisation stores each NAO as a huge, near-cancelling coefficient expansion routed through the chemically insignificant high-angular-momentum shells and neighbouring centres (a benzene “carbon 2p” at cc-pV5Z carries a coefficient norm of ~12 with ~30% of it on the bonded hydrogen, cancelling down to a perfectly atom-localised orbital). The real-space orbital stays localised and integrated populations stay basis-stable, but the ill-conditioned representation has two practical consequences: most molecular viewers (Avogadro2, IQMol — verified Jul 2026; Multiwfn is correct) misrender the exported orbitals as delocalised artefacts, because their approximate high-L handling breaks the coefficient cancellation; and the NPA charges themselves lose numerical conditioning at cardinal ≥ 5. ALG is one of:

• rigorous — Standard NBO-style pipeline. Numbers are directly comparable with NBO7 and the published NPA literature. --export-fchk ao writes the coefficients raw (as for every algorithm since Jul 2026); at large basis sets the raw standard NAOs trigger the HLRL viewer artefact (see §9), so for visualisation prefer orbital cubes (--export-cube ao[LIST]) or request the legacy display projection explicitly with --export-lproj.
• frozen-nmb — Frozen-NMB two-cycle (default since Jul 2026). The Cor+Val (NMB) orbitals are built by the same rigorous machinery on a per-atom angular-momentum-truncated, Galerkin-projected density, then frozen — with exactly zero high-L coefficients by construction — and the Rydberg space is rebuilt around them. Exported orbitals are the analysis orbitals (no display projection needed). Atomic populations agree with rigorous to < 0.001 e through cc-pV5Z, and at cardinal ≥ 5 the frozen-NMB charges are actually the better-conditioned ones — at cc-pV6Z the rigorous pipeline’s near-degenerate basis breaks the symmetry-equivalence of NPA charges by up to ~0.09 e, while frozen-NMB keeps symmetry-equivalent atoms equal to ~10⁻⁴ e and converges smoothly to the basis-set limit. Since the 2026-Oct-01 compilation the Rydberg space keeps marginal near-dependencies (only genuinely singular directions are removed), so the NAO set is full rank and the NAO populations add up exactly to the electron count; NMB and NVB are unaffected, and full-NAO EDDB values may shift by up to ~0.3 % in diffuse bases.
• cc-layers (aliases cclayers / cc) — correlation-consistent layered construction (ccNAO), new in Aug 2026. Cycle 1 builds the same frozen NMB as frozen-nmb; the Rydberg space is then constructed shell-layer by shell-layer in the progression of Dunning’s correlation-consistent basis sets (carbon: NMB(1s 2s 2p) → 3s 3p 3d → 4s 4p 4d 4f → …; hydrogen: 1s → 2s 2p → 3s 3p 3d → …), each layer taking the single most-occupied ANO-like radial function per (atom, L) channel from the residual space and being frozen before the next cycle. Because each layer is selected from a subspace of the previous layer’s complement, the per-channel radial occupations are monotonically decreasing (3s ≥ 4s ≥ 5s ≥ …) for any density — including post-HF natural-orbital densities, where the single-shot Rydberg construction of frozen-nmb can reorder near-degenerate radial functions. Layers are incomplete by definition whenever the computational basis lacks the corresponding functions, and the per-channel labels ascend with the layer index, so --table subshell displays the layer hierarchy directly (and --table shell condenses each layer to a single row). Designed for shell-structured virtual-space truncation in correlation-energy calculations; full method description in docs/ccNAO_construction.txt of the source distribution.
• nbo7legacy — alias of rigorous, kept so published commands and scripts that reproduce NBO7-style raw-NAO exports keep working. (Before Jul 2026 the two differed in how --export-fchk ao coefficients were written; both now write them raw.)

The exported .NAO/.NMB/.NVB.fchk header records which variant produced it (| frozen-NMB, | cc-layers, | raw NAO for rigorous / nbo7legacy, or | L-projected when --export-lproj is given), so orbital pictures remain traceable to the algorithm that generated them.

Reproducibility across compilations. The methods are refined continuously, and a later compilation may change numbers slightly; every such change is listed in the changelog with its scale, and none so far alters the qualitative delocalization picture. State the compilation date (printed in the banner of every run) in the methods section; the compilations used in published work are listed on the download page.

4. Command-Line Reference

Basic usage:

$ runEDDB --input FILE.ext   [ Options ]

  ext = {mol | mol2 | pdb | sdf | xyz}              geometry (built-in GFN2-xTB)
        {fchk | json | molden | tmol | gms | 49}    wavefunction

The input file is the only required argument — the recognised extensions are listed right under the syntax line above (geometry formats run through the built-in GFN2-xTB engine; the rest are read as wavefunctions). Everything else has a sensible default. Output is requested with verbs — --export-cube, --export-fchk, --export-json, --export-map and --table — each taking a comma-separated list of content tokens.

Two built-in help screens. --help (or a no-argument run) prints a short, curated list — just the options most analyses reach for: the five scope modes (-g / -h / -f / -p / -r), NOBD selection (-a / -s), the three output verbs, and the common controls (-b, -c, --cube-ref, the --xtb-* engine controls, --xyz-charge / --xyz-mult, -n, -q, -v). To keep that screen uncluttered, the more specialised flags are deliberately left out of it. Run --full-help for the complete reference, which breaks the --xtb-* controls out flag-by-flag and additionally documents the less-common geometry options (--xyz-center, --xyz-translate, --xyz-frame), citation export (--cite), output naming (--out-prefix / --out-dir), the ASCII toggle (--ascii), the ring-graph cutoff (--ring-thresh), the NAO algorithm selector (--nao-algorithm) and the developer / diagnostic flags (--test-*, --no-nrep / --force-nrep, --no-ecp, --nao-debug*, --nao-atdiag-form, --density-source, --dump-element-table).

The reference table below documents every user-facing option — the curated --help set and the specialised flags that only --full-help shows on screen (the low-level diagnostic flags are listed by --full-help itself):

    --input, -i FILE          Input file (.fchk / .json / .molden / .tmol / .gms / .49 / .xyz / .mol / .sdf / .mol2 / .pdb)

    --xyz-frame N             Select structure N from a multi-frame .xyz/.sdf/.mol2/.pdb (default: 1)
    --xyz-center [geom|mass]  Recenter geometry at origin: geometric centroid (geom, default) or mass
    --xyz-translate x,y,z     Apply a fixed translation vector (Angstrom) to every atom
    --xyz-charge N            Override molecular charge      (geometry input only)
    --xyz-mult M              Override spin multiplicity 2S+1 (geometry input only)

    --global, -g              EDDB_G - entire molecular system (default)
    --heavy-atoms, -h         EDDB_H - entire system, hydrogens excluded from connectivity
    --fragment, -f ...        EDDB_F - molecular fragment, e.g., -f 1:6,10,12
    --pathway, -p ...         EDDB_P - specific bond pathway, e.g., -p 1-2-3-4-5-6-1
    --ring-detect, -r ...     EDDB_R - auto-detect n-membered rings, e.g., -r 5,6:8
    --ring-thresh X           Wiberg ring-bond cutoff for -r (default: 0.30)
    --wbo-thresh, -w X        Wiberg cutoff for the BOP pair gate (default: 0.001; 0 = keep all)

    --select-nobd, -s ...     Select NOBDs for dissection, e.g., -s 1:3,4,5
    --auto-pi, -a             Auto-select pi-type NOBDs to get pi-density component

    --basis, -b ...           NAO basis subset: nao | nmb | nvb (default: nvb; .49 input: nmb)

    --nao-algorithm ALG       NAO construction (native input only):
                                frozen-nmb : two-cycle frozen-NMB construction (DEFAULT);
                                             well-conditioned NAOs at any basis-set size
                                rigorous   : standard NBO-style pipeline (NBO7-equivalent)
                                cc-layers  : correlation-consistent layered NAOs (ccNAO);
                                             Rydberg space frozen shell-layer by shell-layer
                                             in the Dunning progression -> monotonic radial
                                             occupations (aliases: cclayers, cc)
                                nbo7legacy : alias of rigorous (kept for reproducing
                                             published raw-NAO exports)
    --export-lproj            Apply the legacy Galerkin L-projection to the coefficients
                              written by --export-fchk ao (retired former default; cubes
                              are never projected - they are exact grid evaluations)

    Output is requested with three verbs, each taking a comma-separated LIST of
    content tokens (order preserved, duplicates ignored, case-insensitive):

    --export-cube LIST        Volumetric fields to Gaussian CUBE. Tokens:
                                ed | eddb | edlb | bdf
                                ded | deddb | dedlb | dbdf   (differentials; need --cube-ref)
                                ao | no | nobd | bo | xtb    (ORBITAL cubes, one file per
                                orbital; optional 1-based [LIST], e.g. ao[2,4,11-15];
                                no list = every orbital of the lane, guarded at 64)
                                icss | icss-edlb | icss-eddb (ICSS shielding maps: sigma_iso on the
                                grid and its EDLB_G / EDDB_G shares; implies the NMR stage;
                                lattice 0.5 A step, box margin 10 A unless -c / --cube-margin
                                are given; see Section 7)
                              e.g.  --export-cube ed,eddb,deddb --cube-ref baseline.cube
                                    --export-cube ao[2,4] -b nmb     (NMB orbital cubes)
    --cube-ref FILE           Reference cube. Voxel-aligns EVERY exported cube to FILE's grid,
                              and supplies the grid + values baseline for the differential
                              tokens (required when a d-token is present; differentials: not .49)
    --cube-align [FRAGMENT]   Rigidly superpose the current geometry onto the --cube-ref frame
                              (Kabsch fit of a conserved fragment), then map the exported
                              field(s) into that frame -- differential maps AND orbital cubes
                              (ao/no/bo/nobd/xtb). Removes an arbitrary translation + rotation
                              only, reports the post-fit RMSD + max residual. No FRAGMENT =>
                              reuse --fragment. Opt-in; needs --cube-ref.
    --cube-orb2dens           Collapse each orbital lane (ao/no/bo/nobd/xtb[LIST]) into ONE
                              occupation-weighted density cube rho = sum occ_i*|psi_i|^2
                              (<stem>.<LANE>_density.cube) instead of one cube per orbital;
                              composes with --cube-ref / --cube-align. Opt-in.

    --export-fchk LIST        Orbitals / wavefunction to Gaussian FCHK. Tokens:
                                no    natural orbitals of the total density
                                nobd  natural orbitals of bond delocalization
                                ao    AO (NAO-basis) vectors, in the subset set by -b
                                bo    bond-order orbitals (one MO set; alias boo|nba)
                                xtb   GFN2-xTB wavefunction (.xyz input)
                              e.g.  --export-fchk nobd,ao   or   --export-fchk no,bo
                              (the AO subset is chosen ONLY by -b; default nvb)

    --table LIST               Analysis tables to the terminal. Tokens:
                                mo | no | bo | nobd | ao | shell | npa | nmr | nmr-nobd | nmr-mo |
                                nmr-homo[N] | nmr-somo | afm | all
                              mo  = per-MO ED/EDLB/EDDB decomposition (needs canonical MOs).
                              no  = NO occupancy table (the MR / spin-contamination
                                    indices print in the NO stage's -v detail block).
                              bo  = bond-order orbitals (Theta_cov + sigma/pi/delta/phi,
                                    in the chosen NAO|NMB|NVB basis; alias: boo|nba).
                              ao  = per-orbital (atomic-orbital) populations in the -b basis.
                              shell = subshell-resolution populations: (atom, n, L) sums of
                                    the ao rows, e.g. 2p = 2p1+2p2+2p3 (alias: subshell).
                              npa = condensed-atomic populations led by the natural charge q_N.
                              nmr = GIAO / CTOCD-DZ isotropic shielding of every nucleus split into
                                    EDLB_G (total|dia|para) + EDDB_G (total|dia|para) + other
                                    (CTOCD-DZ: total|trans|rot); with -a/-s every
                                    column is restricted to the selected NOBDs; never implied by
                                    "all" (see Section 7). Geometry input (.xyz, .mol, ...): NICS
                                    only, from the uncoupled GFN2-xTB response in the CTOCD-DZ
                                    gauge (no nuclei table, no .nmrref; a --nics* probe is required);
                                    see --xtb-coupled for the coupled response.
                              nmr-nobd = the nmr tables + .NMR_NOBD.txt: the EDDB_G column of
                                    every nucleus / NICS probe resolved by NOBD (all non-zero NOBDs,
                                    or only the -a/-s selected ones -- recommended, the file can be huge).
                              nmr-mo = the nmr tables + .NMR_MO.txt: the full shielding of every
                                    probe resolved by occupied canonical MO (alpha/beta for UHF);
                                    -a/-s do not apply (they select NOBDs, not MOs).
                              nmr-homo[N] = the nmr tables + one table with the shielding carried by
                                    HOMO-N..HOMO (N = 0..9; nmr-homo = the HOMO alone) split into
                                    EDLB_G / EDDB_G / other, plus rest = all other occupied MOs, headed
                                    by the ED/EDLB/EDDB populations of the listed MOs (as --table mo).
                              nmr-somo = the same table for all singly occupied MOs (UHF / ROHF).
                              (--nmr-mos LIST: the same table for an explicit MO set, see below.)
                              All requested tables print together in one captioned
                              "> Printing analysis details:" section before the summary,
                              in the order mo -> no -> bo -> nobd -> ao -> shell -> npa.

                              afm = nc-AFM model contrast d2V/dz2 (eV/A^2) above every bond midpoint
                                    (WBO >= 0.3): total, EDLB_G, EDDB_G with R and the Wiberg order;
                                    model, tip height and slab as --export-map afm (not in 'all')
                              stm = STM window density (10^-6 e/bohr^3, Tersoff-Hamann, --stm-bias)
                                    above every bond midpoint: total, EDLB_G, EDDB_G with R and
                                    the Wiberg order; tip height as --export-map stm (not in 'all')
    --nics-rings LIST         NICS at rings detected by size (-r syntax, e.g. 6 or 5,6:8):
                              NICS(0), NICS(1) = mean of the +1/-1 faces and their zz components
                              (implies --table nmr)
    --nics-ring CHAIN         NICS at explicit atom cycles (-p syntax, e.g. 1-2-3-4-5-6-1,7-8-...)
    --nics-scan R[,step]      NICS scan along the ring normal from -R to +R Angstrom (step 0 or
                              omitted: only -R and +R); isotropic and zz rows
    --nics x,y,z[,...]        NICS at explicit points (Angstrom, coordinate frame of the input file)
    --ring-current [LIST]     Ring current strengths I_cyc (nA/T): Ampere line integral of the
    --rcs [LIST]              induced field around a rectangle through a bond = the GIMIC plane
                              integral, split into EDLB/EDDB (and -a/-s); I_cyc = min over the
                              planes from the centre to the vacuum = the mesh current of the
                              cycle. No LIST: every NICS ring. LIST (-p syntax): closed cycles
                              i-j-k-...-i only (a bond or an open chain is an error)
    --bond-current [LIST]     Bond current strengths I_bond (nA/T), the same integral bond by
    --bcs [LIST]              bond: every bond of every NICS ring (no LIST) or of the listed
                              cycles (plus their I_cyc row), chains i-j-k-... and bonds i-j
    --bcs-box L,h[,n]         Box beyond the bond midpoint L and half-height h (Angstrom, default 5,5)
                              and Gauss points per edge n (default 12)
    --nmr-ref SPEC            Chemical shifts delta = sigma(ref) - sigma (implies --table nmr). Two
                              steps: (1) runEDDB -i tms.fchk --table nmr writes tms.nmrref (element-
                              averaged sigma AND its EDLB/EDDB partition columns); (2) runEDDB -i
                              mol.fchk --nmr-ref tms.nmrref adds a CHEMICAL SHIFTS table, delta column
                              by column (positive = deshielded vs. the reference). SPEC = FILE.nmrref,
                              or an explicit list H=31.7,C=185.4 (delta_iso only). Same method, basis,
                              gauge (GIAO/CTOCD-DZ) and -b in both runs; elements absent from the
                              reference are skipped; under -a/-s no .nmrref is written and --nmr-ref
                              is ignored (run the molecule twice for shifts AND the pi dissection).
    --nmr-kernel [SPEC]       A wavefunction file WITHOUT NMR response data still gets the NMR/NICS/
                              ring-current tables, as the APPROXIMATE CTOCD-DZ response of its MOs
                              (stage line, every caption and a WARNING say so). SPEC = the method that
                              made the MOs (b3lyp, pbe0, hf, pbe, ...; @list prints the table); without
                              it the method is read from the .fchk header, the Gaussian route, the ORCA
                              .out, the Turbomole control or the GAMESS log, else the DFT default
                              a_x = 0.20 is used and the first table says so. Ignored when the input
                              carries NMR data.
    --nmr-mos LIST            Frontier-MO shielding table (as nmr-homo[N]) for an explicit set of
                              occupied MOs, numbered as in --table mo: 5:12,20,25. UHF: a/b prefixes
                              select one spin (a5:12,b5:11); no prefix = both spins of a level.
    --nmr-xdens FILE          Turbomole: XDENS written by mpshift with $gimic in control (default:
                              XDENS in the directory of the .tmol file)
    --nmr-log FILE            Gaussian .log written with IOp(10/33=2): read the CPHF response densities
                              from its P1 blocks (required for the CTOCD-DZ lane, NMR=CSGT); default
                              probe .log / .out when the .fchk carries no P1 section

    --cite LIST               References to cite for THIS run -> file. Tokens:
                                bib | txt | both    (-> .bib / .refs.txt)
                              The reference set depends on the methods actually used;
                              it is also shown per-stage on screen under -v.

    --out-prefix NAME         Override the output filename stem for the whole batch
    --out-dir DIR             Write all generated files into DIR (created if absent)

    (The legacy --output-* flags were removed in v2.0 -- see the migration
     table in section 6 for the old -> new mapping.)

    GFN2-xTB engine (.xyz input only):
    --xtb-etemp K             Electronic (Fermi smearing) temperature in K (default ~300; 0 = off)
    --xtb-conv X              SCC convergence threshold on max|dq| (default: 1e-4)
    --xtb-maxiter N           Maximum number of SCC iterations (default: 1500)
    --xtb-noaes               Disable anisotropic electrostatics (faster; not EDDB-grade)
    --xtb-coupled [a_x]       NICS from the COUPLED tight-binding magnetic response:
                              the sum over states is iterated with an sTDA-type ZDO
                              exchange kernel (Ohno-Klopman gamma from the GFN2
                              hardnesses) scaled by ONE parameter a_x (bare flag =
                              0.20, the B3LYP exchange fraction, fitted on 16 neutral
                              six-membered rings: EDDB loop current within 7 %, NICS(1)zz
                              within 3 ppm of B3LYP/def2-TZVPP). Not for rings with a
                              pi-donor N-H or O (GFN2 underestimates them ~1.5x with or
                              without coupling) nor for 4n rings (the coupling map is
                              unstable at a_x = 0.2: WARNING + uncoupled fallback; use
                              a_x ~0.1). -v prints iterations, spectral radius and the
                              amplification. --xtb-uncoupled = the default lane.
    --xtb-frac-occ            EXPERIMENTAL: keep fractional (smeared) occupations
    --xtb-nospin              Open-shell: use ROHF instead of the default spin-polarised UHF

    --cube-size, -c N         Grid points per axis, ~|N|^3 total (default: 100, ICSS maps 50; N=0 skips)
    --cube-margin R           Box margin around the atoms of the scope, Angstrom (default 3.2;
                              ICSS maps 10); ignored with --cube-ref
    --icss-axis SPEC          Also write the nn component of the ICSS maps (NICS_zz analogue):
                              z | x | y (files ICSS_zz...), ring N (normal of NICS ring N) or
                              x,y,z (files ICSS_nn...)
    --export-map LIST         2D maps on a lattice in a plane -> .map.dat (u v x y z value per
                              line, blocks per u: gnuplot pm3d, numpy): icss | icss-edlb |
                              icss-eddb (ICSS, ppm; --icss-axis adds ICSS_zz/ICSS_nn),
                              ed | eddb | edlb | bdf (density slices, e/bohr^3, same scope
                              and dissection as the cubes), afm | afm-edlb | afm-eddb
                              (rigid-tip nc-AFM model image d2V/dz2, eV/A^2, tip 3.5 A) and
                              stm | stm-edlb | stm-eddb (constant-height STM image, Tersoff-
                              Hamann: density of the MOs in the bias window, e/bohr^3, tip 3.5 A)
    --stm-bias V              Bias window of the stm maps / table in volts (default -2.0):
                              V < 0 = occupied MOs from the HOMO down to HOMO - |V| eV, V > 0 =
                              virtual MOs from the LUMO up to LUMO + V eV (total map only)
    --afm-alpha B             Exponent of the Pauli repulsion V = A rho^B (default 0.78)
    --afm-prefactor A         Prefactor A in eV with rho in e/A^3 (default 692; Moll 2012)
    --afm-tip SIGMA[,s|pz]    Research variant: Pauli energy integrated over a Gaussian tip
                              orbital of width SIGMA (A), s or p_z along the normal (5x5x5
                              Gauss-Hermite); default 0 = point tip (Moll's fit)
    --map-plane SPEC          ring N[,h] (h A above NICS ring N) | fit[,h] (best-fit plane of
                              the heavy atoms; default, h = 1.0) | atoms i,j,k[,h] | x=V|y=V|z=V |
                              radial[,h] (--table afm only: per-bond outward normals, cages)
    --map-step H              Lattice step of the 2D maps, Angstrom (default 0.25)
    --map-slab h1,h2[,n]      Map = mean of n planes (5) between h1 and h2 above the base plane;
                              'none' = one plane. Density slices default to 0.5,1.5 unless
                              --map-plane gives a height; ICSS maps use one plane
    --icss-checkpoint [MIN]   Save the 3D ICSS pass to .ICSS*.tmp every MIN min (10) and
                              resume from it when the file exists; removed on completion
    --ncores,    -n ...       Number of CPU cores, 0 = all available (default: 1)
    --quiet,     -q           True-quiet: banner + EDDB summary block only
    --verbose,   -v           Add per-stage detail blocks (tables are opt-in via --table)
    --ascii                   Force plain 7-bit ASCII output (default: UTF-8)
    --help                    Print the short, curated end-user help (also shown with no args;
                              NOTE: -h is the EDDB_H scope, not a help alias)
    --full-help               Print the complete reference (adds --cite, --out-*, --ascii,
                              --ring-thresh, --xyz-center/-translate/-frame, per-flag --xtb-*
                              detail + diagnostics)

runEDDB has three verbosity tiers:

• Default (no flag) — banner, stage messages with per-stage timings, and the EDDB summary block. No analysis tables are printed by default — request them explicitly with --table.
• True-quiet (-q / --quiet) — banner, the EDDB summary block, and its two bookend stage lines (> Printing summary of EDDB calculations: and > Calculation completed on YYYY-MM-DD. Total time: Ns.). Per-stage timings and all detail blocks are suppressed. Useful for batch jobs where only the numerical result and a "DONE" confirmation are needed.
• Verbose (-v / --verbose) — everything the default tier prints, plus the Stage 1/2/3 detail blocks and a few one-line clarifiers (e.g. the .49-fallback note under Target NAO representation:). --verbose does not control tables — those are governed solely by --table. For the NMR family (--table nmr, --nics*, --rcs, --bcs) -v also writes .NMR_DIAG.txt: every table again with the trans/rot (dia/para) split at full precision plus the diagnostic blocks (both NICS(±1) faces, the probe-resolved sum rule, the loop-current fit, the bond-current edge terms); the terminal keeps the four-column result tables. The file is not bound by the 88-column rule.

Analysis tables are printed only when requested via --table (independent of the verbosity tier): --table mo for the per-MO ED / EDLB / EDDB decomposition (when canonical MOs are available), --table no for the NO-occupancy table (its multireference / spin-contamination indices print in the NO stage's -v detail block), --table bo for the bond-order orbitals (aliases boo / nba), --table nobd for the NOBD-resolution table, --table ao for the per-orbital (atomic-orbital) populations, --table subshell for the subshell-resolution populations (per-(atom, n, L) sums of the AO rows), --table shell for the shell-resolution populations (per-(atom, n) sums over L), --table npa for the condensed-atomic table (led by the natural charge q_N), or --table all for every table. All requested tables are emitted together in one captioned > Printing analysis details: section immediately before the EDDB summary, in the order mo → no → bo → nobd → ao → subshell → shell → npa (the orbital-based tables MO / NO / BOO come before the NOBD-driven dissection; AO / SUBSHELL / SHELL / NPA come after it). Using -s without a list still shows the NOBD-resolution table (so you can pick indices for a follow-up run). In true-quiet mode the EDDB summary still prints, and --table tables requested alongside -q are honoured.

Citing the methods you used (--cite)

runEDDB can tell you exactly which papers to cite for a given run — the set is method- and option-dependent, emitted by the stage that uses each reference:

• Always — the NAO / NPA paper (Reed et al.; the NAO basis underpins every analysis), the original BOP paper, and the core EDDB method paper.
• Only when the matching option is used — the NOBD dissection paper with -s / -a; the EDDB_P paper with -p / -r or the EDDB_F paper with -f; the AO-resolved EDDB paper with --table ao / --export-fchk ao; the Bond-Delocalization-Field paper when you export a bdf / dbdf cube.
• .xyz lane only — GFN2-xTB + DFT-D4 (GFN2-xTB swapped for spGFN2-xTB on spin-polarised runs).

The remaining algorithmic ingredients (the Wiberg index, the per-pair NBA / BOO decomposition, OWSO / Gram-Schmidt orthogonalisation, the multireference indices, ring detection, …) are not in the end-user list — they are attributed in the source-code comments instead. Two ways to get the references:

• On screen, under -v — each stage’s verbose detail block ends with a short Papers to be cited (DOI): list, formatted as <DOI> (method), so you see the relevant references right next to the step that used them. (Nothing citation-related prints in the default or quiet tiers — it is opt-in via -v.)
• To a file, with --cite LIST — --cite bib writes a BibTeX <stem>.bib, --cite txt a plain numbered <stem>.refs.txt, and --cite both writes both. The exported list is the de-duplicated set for the whole run, independent of the verbosity tier. See §11 Citations for the full reference list.

GFN2-xTB engine controls (--xtb-*)

These flags tune the built-in GFN2-xTB SCC and apply only to .xyz input (they are ignored for the wavefunction-file lanes). The defaults are chosen to be robust and EDDB-appropriate, so most runs need none of them; they exist for the occasional stiff or unusual system.

• --xtb-etemp K — electronic (Fermi-smearing) temperature in Kelvin; default ~300 K (matching tblite). A little smearing damps charge-sloshing and helps the SCC settle. If a system oscillates and will not converge, raising it (e.g. --xtb-etemp 1000) is the first thing to try. A value of 0 turns smearing off entirely (strict Aufbau occupations).
• --xtb-conv X — SCC convergence threshold on the maximum shell-charge change max|Δq|; default 1×10⁻⁴. This is deliberately loose: the bonding pattern that EDDB resolves is locked in well before max|Δq| drops below 10⁻⁴, and tighter thresholds cost extra iterations without changing the analysis. Tighten it (e.g. --xtb-conv 1e-6) only if you need converged total energies; loosen it (e.g. 1e-3) to coax a difficult case to finish.
• --xtb-maxiter N — maximum number of SCC iterations; default 1500. Raise it if a large or polarisable system reports non-convergence and you would rather let it run longer than relax the threshold.
• --xtb-noaes — disable GFN2’s anisotropic electrostatics (the atomic dipole / quadrupole, “AES”, terms). The SCC then runs faster and more easily, but the result is not EDDB-grade: AES carries the lone-pair / heteroatom anisotropy that the delocalization picture depends on. Use only as a diagnostic or last-resort convergence aid, not for production analysis.
• --xtb-frac-occ — EXPERIMENTAL / diagnostic only. By default the GFN2-xTB analysis density is rebuilt from integer orbital occupations — a true single determinant, so each spin block is idempotent — and the Fermi smearing used inside the SCC is treated purely as a convergence device that never reaches the analysis. This flag keeps the converged fractional occupations instead, building D_s = C diag(n_s) C^T from them. The resulting density is deliberately non-idempotent, is classified post-HF by the usual occupation test, and therefore takes the layer-by-layer (post-HF) BOP path. Its purpose is to manufacture cheap non-idempotent densities for exercising and validating that post-HF engine — not to produce EDDB values. Results are not comparable with any published EDDB data.
  
  Pair it with a high electronic temperature, e.g. --xtb-frac-occ --xtb-etemp 5000. The flag is a deliberate no-op when smearing is effectively off: if every converged occupation sits within 10⁻³ of an integer (the same tolerance that classifies post-HF input), runEDDB reports that it had no effect and falls back to the standard integer path rather than silently producing a density that would fail to exercise the post-HF code. When it does fire it prints a prominent warning together with the largest deviation from integer occupation — the direct measure of how non-idempotent the density is. Representative values for Li₈: 300 K → no-op; 1000 K → 0.0016; 3000 K → 0.072; 5000 K → 0.171; 10000 K → 0.321.
  
  Interpret the output with care: these fractional occupations reflect an artificial electronic temperature, not physical static correlation, and the layer decomposition was designed for CASSCF/CI natural occupations. Note also that the SCC convergence ladder may raise the temperature above your request on difficult systems; runEDDB says so explicitly when that happens, and you can pin the value with --xtb-etemp.
  
  --table mo remains available under this flag and is the most direct way to inspect what it did: the Occupancy column shows the fractional values themselves, orbital by orbital. This is a genuine exception to the usual rule that a post-HF density suppresses the MO table — there the density is not built from canonical MOs at all, whereas here the orbitals are canonical GFN2-xTB MOs and only their occupations are fractional, so the per-MO decomposition remains meaningful.
• --xtb-nospin — for open-shell systems (multiplicity > 1), solve the spin-restricted open-shell (ROHF-style) problem instead of the default spin-polarised unrestricted (UHF / spGFN2-xTB) treatment. Has no effect on closed-shell singlets.

If the SCC does not converge, runEDDB stops and prints a short reminder suggesting exactly these knobs (a higher --xtb-etemp, a larger --xtb-maxiter, or a looser --xtb-conv), so you rarely need to remember them in advance.

5. Output Files

Outputs are generated only when explicitly requested with --export-cube (cube files) or --export-fchk (FCHK files); --table prints analysis tables to the terminal without writing any file.

Gaussian cube files (.cube) — via --export-cube LIST; 3D electron density on a PCA-aligned grid; viewable in VMD, Chemcraft, GaussView, Avogadro2, IQmol, or any cube-compatible program. Tokens:

• eddb — delocalized electron density (EDDB)
• edlb — localized density (D_tot − D_del)
• ed — full 1-electron density
• bdf — Bond Delocalization Field (2·D_del − D_tot)
• ded, deddb, dedlb, dbdf — differential maps Δf = f(this) − f(REF), subtracted point-by-point for the total density (ED), the delocalized (EDDB), the localized (EDLB) or the bond-delocalization field (BDF). All differential tokens take a single shared reference cube via --cube-ref FILE. See §6 for the full workflow.
• icss, icss-edlb, icss-eddb — ICSS shielding maps: the isotropic magnetic shielding σ_iso (ppm, shielding sign, NICS = −value) at every grid point, and its EDLB_G / EDDB_G shares (under -a / -s the EDDB_G map is that of the selected NOBDs). Files <stem>.ICSS.cube, <stem>.ICSS_EDLB_G.cube, <stem>.ICSS_EDDB_G.cube; the tokens imply the NMR stage of §7 (any lane) and use a 0.5 Å lattice with a 10 Å margin unless -c / --cube-margin are given; --icss-axis adds the nn component files (<stem>.ICSS_zz.cube, …).

2D maps (.map.dat) — via --export-map LIST. The tokens ed, eddb, edlb, bdf write density slices (e/bohr³; by default the mean of five planes between 0.5 and 1.5 Å above the base plane, see --map-slab in §7; the same matrices, scope and dissection as the cubes: <stem>.EDDB_G.map.dat, <stem>.BDF_P.map.dat, …; BDF = EDDB − EDLB node by node); the tokens icss, icss-edlb, icss-eddb write the ICSS maps: the ICSS values on a lattice in a plane chosen with --map-plane (a ring at height h, the best-fit molecular plane, three atoms, or a frame plane; default the best-fit plane of the heavy atoms 1 Å above it), step --map-step (0.25 Å), extent = the projected atoms plus the box margin. Files <stem>.ICSS.map.dat, <stem>.ICSS_EDLB_G.map.dat, <stem>.ICSS_EDDB_G.map.dat (and ICSS_zz / ICSS_nn with --icss-axis): six # header lines (command, plane origin / normal / in-plane axes, lattice, columns), then one line per node with the in-plane coordinates u, v, the Cartesian x, y, z (Å) and the value (ppm), the rows of constant u separated by a blank line — ready for gnuplot (splot 'f.map.dat' using 1:2:6 with pm3d) or numpy.loadtxt. The plane origin is always a node, so the value at u = v = 0 of a ring map is the NICS(h) row of the table with the sign reversed.

nc-AFM model images (--export-map afm, afm-edlb, afm-eddb) — a simulated high-resolution AFM image with a CO tip (Gross et al. 2009) from the rigid-tip Pauli-repulsion model of Moll et al. (2012): in the short-range regime the tip–sample repulsion is a power law of the sample density alone, V = A ρ^α (their CO/PTCDA fit A = 692 eV, α = 0.78 with ρ in e/ų; --afm-prefactor, --afm-alpha), and the frequency shift of a small-amplitude sensor is Δf = (f₀/2k) ∂²V/∂z². The map holds ∂²V/∂z² in eV/Ų (central differences, 0.05 Å) on the --map-plane lattice at the tip height, by default 3.5 Å above the best-fit molecular plane (Moll’s atomic-resolution window is 3.0–4.5 Å; --map-plane fit,h or ring N,h sets it). The EDLB_G / EDDB_G shares are the proportional split V_X = V ρ_X/ρ = A ρ^α−1ρ_X, which is additive at every point for any α (it is the plain density split for α = 1); afm − afm-edlb − afm-eddb is the core/Rydberg remainder outside the analysis block. A research variant, --afm-tip SIGMA[,s|pz], replaces the point evaluation by the integral of V over a Gaussian tip orbital of width σ (s-type, or p_z along the normal; 5×5×5 Gauss–Hermite quadrature, exact for polynomials up to degree 9 per axis), E(R) = ∫|ψ_tip(r−R)|²V(r) d³r, with the same additive shares inside the integral; because the fitted A and α already absorb the CO tip, the variant is meant for exploring the tip dependence of the contrast, not as a better default (for naphthalene an s orbital of 0.5 Å doubles the values and reduces the bond discrimination, a p_z orbital triples them). The model neglects the van der Waals and electrostatic background and the CO-tip relaxation, so it reproduces the intramolecular contrast and not absolute frequency shifts. Validation (2026-09-20): on naphthalene, anthracene, phenanthrene, pyrene, triphenylene, coronene, hexa-peri-hexabenzocoronene and C₆₀ (B3LYP/def2-TZVPP) the per-bond contrast of --table afm correlates with the Wiberg bond order at r = 0.94–0.98 and with the Pauling order from the bond length at 0.74–0.96, i.e. it reproduces the bond-order discrimination of Gross et al. (Science 337, 1326, 2012): in hexabenzocoronene the five symmetry classes of C–C bonds are ordered monotonically from the peripheral benzo bonds (1.387 Å, brightest) to the inter-ring bonds (1.461 Å, ratio 1.31), and in C₆₀ the 6:6 bonds are brighter than the 6:5 bonds by a factor 1.20 (--map-plane radial gives every bond its own outward normal for cages). The EDDB_G share carries 50–98 % of the contrast, the differences between bond classes coming largely from the EDLB_G part. The apparent shortening of high-order bonds seen experimentally is a CO-tilt effect and is outside the rigid model; files <stem>.AFM.map.dat, <stem>.AFM_EDLB_G.map.dat, <stem>.AFM_EDDB_G.map.dat (scope letter as for the cubes; -a/-s give the share of the selected NOBDs). --table afm prints the same quantity above the midpoint of every bond (Wiberg order ≥ 0.3) next to the bond length and the Wiberg bond order — the per-bond form of the bond-order contrast of Gross et al. (Science 337, 1326, 2012) — with the same model, tip height and --map-slab handling. The GFN2-xTB lane carries a minimal basis without the diffuse density tail that the model probes at 3.5 Å, so its images are two orders of magnitude weaker; for a qualitative picture from a geometry lower the tip (--map-plane fit,2.5) and compare contrasts, not values.

STM images (--export-map stm, stm-edlb, stm-eddb; --stm-bias V; --table stm) — a constant-height scanning-tunnelling image in the Tersoff–Hamann approximation (Phys. Rev. B 31, 805, 1985): for an s-wave tip the tunnelling current at the tip position is proportional to the local density of the sample states between the Fermi level and the bias, which for an isolated molecule is the density of the MOs inside the bias window. --stm-bias V (default −2.0 V) selects that window with the molecular gap as the zero of the scale: V < 0 images the occupied MOs from the HOMO down to HOMO − |V| eV, V > 0 the virtual MOs from the LUMO up to LUMO + V eV. The map holds the window density in e/bohr³ on the --map-plane lattice at the tip height (default 3.5 Å above the best-fit plane, one plane unless --map-slab is given), so that it is directly comparable with the AFM maps of the same run. For an occupied window the density is partitioned as in --table mo: the EDDB_G part is the symmetrised product ½{D_EDDB, P_V} of the EDDB density matrix with the projector on the window MOs in the analysis block (EDLB_G likewise with D − D_EDDB), which is additive, integrates to the --table mo sums over the window and reduces to the density slices when the window holds all occupied MOs; stm − stm-edlb − stm-eddb is the core/Rydberg remainder outside the analysis block. A virtual window carries no bonding partition, so a positive bias writes the total map only. --table stm prints the same three quantities above every bond midpoint (Wiberg order ≥ 0.3, in 10⁻⁶ e/bohr³) next to the bond length and the Wiberg order, with the bias and the window (number of MOs, energy range) in the legend; sweeping the bias is the physical way of “reading down” the MO spectrum from the frontier: in coronene the per-bond STM density already follows the AFM contrast at −0.5 V (r = 0.98), in benzo[k]fluoranthene it does so only from about −5 V (r = 0.54 at −1 V, 0.83 at −2 V, 0.95 at −5 V), the frontier pair alone painting a different pattern. Requires an SCF (HF/DFT) wavefunction with MO coefficients and energies (FCHK, ORCA JSON, Molden, GAMESS, Turbomole; not the NBO .49 lane, not post-HF densities); the model neglects the tip structure and the substrate, and images the LDOS of the free molecule.

Scoped runs (-f / -p / -r) — exact atomic support. All three scope-dependent densities share exactly the atomic support of the selected atoms: after the congruence transformation the EDDB matrix is purified (rows and columns of out-of-scope atoms zeroed), and — as of Jul 2026 — the auxiliary π-electron density entering the EDLB difference (D_tot − D_del) and the BDF combination (2·D_del − D_tot) receives the same purification (it is reconstructed from the full molecular density via the selected π-NOBDs, so it would otherwise retain contributions on atoms outside the scope). Consequently the integrated cube density, the printed basis-trace reference and the per-atom table sums agree exactly — e.g. the EDLB_F cube integrates to ED_F − EDDB_F. The same purified matrices feed the differential tokens (dedlb / dbdf).

Grid quality is controlled by --cube-size N (default: 100; the ICSS maps use instead an exact 0.5 Å step on every axis with a node on the box centre, unless -c is given) and the box by --cube-margin R (the margin in Å added around the atoms of the active scope on every side, default 3.2 Å = 6 Bohr; the ICSS maps use 10 Å unless the option is given; a cloned grid from --cube-ref ignores it), with the grid covering the molecular bounding box (or a fragment subset, when -f / -p / -r is in effect) plus a small margin. The sign and zero of N select between three modes:

• N > 0 — ~N³ grid points covering the active scope (full molecule under -g; fragment / pathway / ring atoms under -f / -p / -r).
• N < 0 — ~|N|³ grid points but always covering the entire molecule, even when the analysis is fragment-, pathway- or ring-restricted. Useful when you want the cube rendered on the whole system for visualization while the EDDB analysis itself stays scoped.
• N = 0 — skip cube writes entirely. Every --export-cube token becomes a no-op (FCHK writes are not affected).
• --cube-ref FILE (separate flag) — clone the exact grid (origin, point counts, step vectors) from an existing .cube file. Every exported cube is then voxel-aligned with FILE so they can be subtracted, added or otherwise combined point-for-point in any cube-arithmetic tool (Multiwfn, VMD, etc.) — or, more simply, with the built-in differential tokens (§6), which use the same --cube-ref.
• --cube-orb2dens (separate flag) — for the orbital lanes (ao / no / nobd / bo / xtb[LIST]), collapse the whole selection into a single occupation-weighted density cube ρ(r) = Σ_i occ_i·|ψ_i(r)|² (written as <stem>.<LANE>_density.cube) instead of one cube per orbital — e.g. --export-cube no[449,450] --cube-orb2dens gives the density of those two (near-singly-occupied) natural orbitals in one field, integrating to their summed occupation. It reuses the same grid evaluation, so it composes directly with --cube-ref and --cube-align. As a pure modifier of the orbital lanes it is silently ignored when no orbital cube is requested (e.g. only --export-fchk, or only scalar cube tokens), so it is safe to carry in a batch script.

Gaussian FCHK files (.fchk) — via --export-fchk LIST; viewable in GaussView, Avogadro2, IQmol, or Chemcraft. Tokens:

• nobd — NOBDs as Alpha MO coefficients + EDDB density as Total SCF Density.
• no — natural orbitals (NOs) of the spinless 1-electron density (eigenvectors of the NAO-basis density matrix). Useful for inspecting occupied / partially-occupied orbitals and identifying open-shell / multireference character. Runs as a dedicated stage right after the NAO basis is built, so the export is available for any input format. (Pair with --table no to also print the NO-occupancy table and the multireference diagnostics.)
• bo — bond-order orbitals (BOOs): every selected atom pair’s Wiberg bond order resolved into orthonormal σ/π/δ/φ components, written as one MO set (the per-orbital occupation in the energies field). Pairs with --table bo; aliases boo / nba are accepted; the file is named <stem>.BO.fchk. The pair set follows the active scope (--fragment / --pathway / --ring-detect), so on large systems scoping the run keeps this file small — see §8.
• ao — AO (natural-atomic-orbital) basis vectors. The subset is chosen solely by -b (nao / nmb / nvb; default NVB) and the file is named for the resolved basis (.NAO/.NMB/.NVB.fchk). Per-orbital populations stored in the file are read from the post-dissection D_tot diagonal, so when --auto-pi or --select-nobd is in effect they exactly match the values printed by --table ao (orbitals that were not part of any selected NOBD show as 0.0). Coefficients are written raw for every algorithm by default — the former HLRL display projection is opt-in as of Jul 2026 (--export-lproj): the default frozen-NMB orbitals are intrinsically well-conditioned, and projected (polarisation-stripped) pictures invite the false reading that they are “better localised” than the honestly polarised frozen-NMB NAOs. --nao-algorithm rigorous (alias nbo7legacy) writes the raw standard NAOs. For guaranteed-faithful pictures in any viewer, prefer --export-cube ao[LIST] (see the NAO construction algorithm paragraph in Section 4).
• xtb (.xyz input only) — the converged GFN2-xTB wavefunction itself: the minimal valence STO-nG (vSTO-nG) basis, the molecular orbitals and orbital energies (α and β for open-shell), the SCC electron density (written as Total SCF Density) and, for open-shell systems, the spin density (Spin SCF Density). Written in seconds, right after the GFN2-xTB step. Two typical uses: (1) a ready-made initial guess for a follow-up DFT calculation on the same geometry (many DFT codes can read an .fchk density / orbitals as a starting point, which can shorten the DFT SCF and helps seed open-shell / broken-symmetry guesses); and (2) visualization — open it in Avogadro2, GaussView, IQmol, Chemcraft or Multiwfn to view the GFN2 density, spin density and orbitals.
  Note: this file describes the molecule in GFN2’s own minimal valence vSTO-nG basis with a parametrised effective core — it is a faithful, self-contained GFN2 wavefunction for guessing and viewing, but it is not an all-electron wavefunction and should not be treated as a production ab-initio result. (Requested for a non-.xyz input, the flag is ignored with a warning, since there is no GFN2 calculation to export.)

For all three FCHK exports the per-orbital occupation (NAO population, NOBD occupation, or NO occupation number) is stored in the Alpha Orbital Energies field, pre-scaled so that the value Avogadro2 shows in its eV-labelled energy column reads back numerically as the unscaled occupation. So a doubly-occupied orbital appears as "2.00 eV", a singly-occupied one as "1.00 eV", and so on — the unit label is technically eV but the number is the actual occupation. Other viewers (GaussView, IQmol, Chemcraft) that read FCHK energies in hartree will show the raw small values; if you need the actual occupations there, multiply by 27.2114.

Terminal tables (--table) — no file written. Every table requested with --table is collected into a single captioned > Printing analysis details: section, emitted just before the always-on EDDB summary, in the order mo → no → bo → nobd → ao → subshell → shell → npa. Tokens: mo (per-MO ED / EDLB / EDDB decomposition, when canonical MOs are present), no (NO occupancy + multireference diagnostics), bo (bond-order orbitals; σ/π/δ/φ resolution; aliases boo/nba), nobd (NOBD-resolution table), ao (per-orbital atomic-orbital populations in the -b basis), subshell (subshell-resolution populations: per-(atom, n, L) sums of the AO rows), shell (shell-resolution populations: per-(atom, n) sums over L), npa (condensed-atomic populations led by the natural charge q_N), and all (every table). The orbital-based tables (MO / NO / BOO) print before the NOBD-driven dissection; ao / subshell / shell / npa print after it. Each table is introduced by a capitalised caption; see §8 for what every column means.

NMR reference (--table nmr) — <stem>.nmrref. Written by every --table nmr run: the element-averaged isotropic shieldings and their EDLB/EDDB (dia|para or trans|rot) partition columns, plus one line per atom. Feed it back with --nmr-ref FILE.nmrref to obtain partitioned chemical shifts relative to that molecule (TMS, solvent); see §7.

Reference list (--cite) — <stem>.bib / <stem>.refs.txt. --cite bib writes a BibTeX file, --cite txt a plain numbered list, and --cite both writes both. The contents are the exact, de-duplicated set of papers to cite for the run — method- and input-dependent — and honour --out-prefix / --out-dir like every other output. The same references are shown per-stage on screen under -v (see §4); the full reference list is in §11.

Filename rules. Filenames are standardized and deterministic — each token names its own file as <stem>.<TOKEN>_<SCOPE>.<ext>, so a single verb can emit many files coherently. The <stem> is the input basename (override the whole batch with --out-prefix); <SCOPE> is the scope letter G/F/P/R (from -g/-f/-p/-r), baked into the name for scope-dependent tokens (eddb, edlb, bdf, nobd, and the differentials) and omitted for scope-independent ones (ed, no, xtb, AO). The AO export names itself by the resolved basis (NAO/NMB/NVB):

Input:   molecule.fchk          (or .json / .molden / .49 / .xyz)
ed       -> molecule.ED.cube
eddb     -> molecule.EDDB_G.cube       (_H / _F / _P / _R under -h / -f / -p / -r)
edlb     -> molecule.EDLB_G.cube
bdf      -> molecule.BDF_G.cube
ded      -> molecule.DED_G.cube        (differentials; need --cube-ref)
deddb    -> molecule.DEDDB_G.cube
dedlb    -> molecule.DEDLB_G.cube
dbdf     -> molecule.DBDF_G.cube
ao[2,4]  -> molecule.NVB_orb002.cube, molecule.NVB_orb004.cube
            (orbital cubes: one file per orbital, <stem>.<LANE>_orbNNN.cube;
             the AO lane is named for -b, the others NO / NOBD / BO / XTB)
nobd     -> molecule.NOBD_G.fchk
no       -> molecule.NO.fchk
ao       -> molecule.NVB.fchk          (named for -b: NAO / NMB / NVB; default nvb)
bo       -> molecule.BO.fchk           (bond-order orbitals)
xtb      -> molecule.XTB.fchk          (from .xyz input)

The scope (-g / -f / -p / -r) is encoded in the _G / _F / _P / _R suffix so two cubes that differ only by scope never collide. --auto-pi / --select-nobd dissection is reflected in the file contents, not the name. The export verbs deliberately drop per-file naming — use --out-prefix NAME to set the stem for the whole batch and --out-dir DIR to choose the location.

Migration from the legacy --output-* flags

The old one-flag-per-thing --output-* family (from the released runEDDB_nbo49) was removed in v2.0 — the three verbs replace it. Update old commands as follows:

--output / -o                  ->  --export-cube eddb --export-fchk nobd
--output-ed                    ->  --export-cube ed
--output-eddb                  ->  --export-cube eddb
--output-edlb                  ->  --export-cube edlb
--output-bdf                   ->  --export-cube bdf
--output-ded   REF [OUT]       ->  --export-cube ded   --cube-ref REF [--out-prefix OUT]
--output-deddb REF [OUT]       ->  --export-cube deddb --cube-ref REF [--out-prefix OUT]
--output-dedlb REF [OUT]       ->  --export-cube dedlb --cube-ref REF [--out-prefix OUT]
--output-dbdf  REF [OUT]       ->  --export-cube dbdf  --cube-ref REF [--out-prefix OUT]
--output-nobd                  ->  --export-fchk nobd
--output-no                    ->  --export-fchk no
--output-ao                    ->  --export-fchk ao           (basis via -b)
--output-xtb  [FILE]           ->  --export-fchk xtb [--out-prefix FILE]
--cube-size 0 FILE             ->  --cube-ref FILE

Lists let you combine these in one call, e.g. the old --output-eddb --output-edlb --output-bdf becomes --export-cube eddb,edlb,bdf.

6. Differential Maps (ΔED / ΔEDDB / ΔEDLB / ΔBDF)

A differential map subtracts one runEDDB density field from another, point-by-point on a shared grid, to reveal how the electron structure reorganizes between a reference and a perturbed system — a substituent change, a hydrogen / halogen bond, a π···cation or π···HF contact, solvation, or a redox step. The differential of the bond-delocalization field, ΔBDF_π, maps π-electron (de)localization changes directly in the language of resonance; the same operation applies to every density in the EDDB decomposition, so there is one differential token per density (all sharing a single --cube-ref):

• --export-cube ded — ΔED = ED − ED(REF): total electron-density difference (a classic deformation-density map).
• --export-cube deddb — ΔEDDB: change in the delocalized density.
• --export-cube dedlb — ΔEDLB: change in the localized density.
• --export-cube dbdf — ΔBDF: change in the bond-delocalization field (2·D_del − D_tot) — the sharpest π-reorganization picture.

The reference is a cube of the same density written by an earlier run, supplied via --cube-ref FILE: it provides both the grid (origin / points / steps, cloned exactly) and the reference values, so the subtraction is point-for-point and no external cube-arithmetic tool (Multiwfn, VMD, …) is needed. Outputs follow the standard naming scheme (<stem>.DED_G.cube, .DEDDB_G.cube, …; set the stem with --out-prefix). The differential maps run on the native lane (not the legacy .49).

π-resolved (ΔBDF) versus full (ΔEDDB / ΔEDLB). ΔBDF is normally a π-only quantity, so it should be built from the π-delocalization (π-NOBDs; see the workflow below). When a clean π-dissection is awkward, ΔEDDB and ΔEDLB on the full delocalized / localized densities need no π-selection at all and still show where delocalization is gained or lost — frequently the easier and more robust route.

Recommended ΔBDF_π workflow. The map is only as good as the set of NOBDs included, and the BDF cubes for both systems must be built the same way. (1) compute and save the reference π-BDF cube, e.g. runEDDB -i ref.fchk --auto-pi --export-cube bdf --out-prefix ref (writes ref.BDF_G.cube); (2) for the perturbed system run runEDDB -i pert.fchk --auto-pi --export-cube dbdf --cube-ref ref.BDF_G.cube. Capture the π-system consistently in both molecules: --auto-pi selects π-symmetry NOBDs automatically, but a perturbation can change how many NOBDs carry the π-delocalization (and the classifier may miss one in a given system). Inspect the NOBD-resolution table for each molecule (printed with --select-nobd / --auto-pi / --verbose) and, if needed, replace --auto-pi with an explicit --select-nobd LIST so that all NOBDs contributing to the π-(de)localization are included in both runs — otherwise ΔBDF shows a spurious feature wherever a real contributor was dropped from one side. Typical ΔBDF_π isosurfaces are around ±0.010.

Same coordinate frame, same scope. A point-by-point subtraction is only meaningful if the two systems share one frame. The most robust route is to align them before the wavefunction calculation — keep the common fragment’s Cartesian coordinates fixed. Alternatively, apply an after-the-fact rigid superposition at export time with --cube-align (see the box below), which absorbs an arbitrary translation and rotation. runEDDB warns if the reference and current geometries do not overlap. Both runs must also use the same scope (--global / the same --fragment / --pathway / --ring-detect) and the same density flavour, so that their grids, NOBD sets and densities are comparable.

Large, flexibly-reorganizing systems — focus on a fragment. When a perturbation reshapes a remote reactive site so the whole molecule no longer overlaps, restrict the analysis to the conserved region of interest with --fragment (or --pathway). This focuses the EDDB density, the cube grid and the geometry-overlap check on the selected atoms only, so the differential reports how that fragment’s π-system responds to a change happening elsewhere — without being penalised for motion in the part you are not studying. The overlap check is purely spatial and element-based, so the reference and perturbed structures may number their atoms differently: as long as the studied fragment is physically superimposed, the overlap is recognised. (Keep that fragment’s coordinates frozen between the two systems, so the differential reflects electronic reorganization rather than atomic displacement.)

Automatic rigid alignment — --cube-align [FRAGMENT]. When the reference and perturbed geometries are not in the same frame (e.g. two optimisations, two solvation environments, or coordinates from different sources), you do not have to pre-align them by hand. Add --cube-align to a differential or orbital-cube run: runEDDB performs a rigid-body Kabsch superposition of a conserved fragment onto its counterpart in the --cube-ref cube (finding the rotation R and translation t with x_aligned = R·x + t), then evaluates the exported field — and writes the shown geometry — in that reference frame. Because the density is a scalar field, this is done by sampling it at the back-mapped grid points; no rotation of the wavefunction or density matrix is involved, and the result is exact. The very same back-map carries the per-orbital cubes (ao / no / bo / nobd / xtb) into the reference frame — each orbital is a scalar field too — so orbital cubes emitted alongside a differential map are aligned identically; the superposition is fitted and reported once and shared by every lane. Without a FRAGMENT list it reuses the --fragment selection. The alignment atoms must appear at the same index and carry the same element in both the perturbed input and the reference cube (a mismatch is a hard error).

# ΔBDFπ of the same chromophore in two solvation environments.
# The porphyrin core (atoms 1-24) is conserved; the two geometries
# differ by an arbitrary translation + rotation (and local relaxation).

$ runEDDB -i gas.fchk    --auto-pi --export-cube bdf  --out-prefix gas
                                          # writes gas.BDF_G.cube (reference)

$ runEDDB -i water.fchk  --auto-pi --export-cube dbdf \
          --cube-ref gas.BDF_G.cube --cube-align 1:24    # aligns the core, then subtracts

# console report emitted for the aligned differential:
    Alignment atoms (--cube-align):    24
    RMSD after rigid alignment:           0.03812 Angstrom
    Max residual displacement:            0.07105 Angstrom
    (arbitrary translation + rotation removed; residual = internal deformation)

What --cube-align does and does not remove. It cancels only an arbitrary rigid translation and rotation of the perturbed system relative to the chosen fragment. It does not remove genuine internal (non-rigid) deformation: any such change of the alignment fragment survives the fit and is reported as the post-fit RMSD and the maximum residual displacement (the two lines above). When those residuals are small, the map is a clean electronic-reorganization picture; when they are large, the perturbation has deformed the fragment beyond the perturbative regime and ΔBDF_π is no longer well defined — the alignment does not, and cannot, repair that. Treat the reported RMSD as the validity gauge for the difference map.

7. NMR Shielding and NICS Partition (--table nmr)

--table nmr takes the magnetic response of a GIAO (or CSGT) NMR calculation (or, for a wavefunction file without NMR data, an approximate response built from its MOs, see below) and splits the isotropic shielding of every nucleus — and the NICS at ring centroids, along ring normals or at arbitrary points — into contributions of the localized valence density (EDLB_G) and the delocalized valence density (EDDB_G), each further separated into the part carried by the unperturbed density and the part carried by the induced response. With -a / -s every column is restricted to the selected NOBDs, so the π ring-current contribution to NICS is read directly from the π-type NOBDs instead of from canonical MOs; --table nmr-nobd and --table nmr-mo additionally write the NOBD- and MO-resolved contributions to files (see below). The one-electron shielding integrals are computed by runEDDB itself; only the CPHF response densities come from the quantum-chemistry program — or, for a geometry input, from the built-in GFN2-xTB lane, which yields approximate NICS (not nuclear shieldings) without any external program (see the end of this section).

When do the π-NOBD columns equal the π-MO contributions? The selected-NOBD value is a trace with the projector on the selected NOBDs, so it does not depend on how the orbitals inside that subspace are mixed. It therefore equals the sum of the canonical π-orbital contributions (the quantity the Steiner–Fowler rules are formulated for) only when the selected subspace is the occupied π space: a planar system in which σ and π do not mix, the full NAO representation (-b nao), and a selection — automatic (-a) or an explicit -s LIST — that captures every occupied π-NOBD and nothing else; then the π-NOBDs are just a unitary image of the π-MOs. For non-planar or σ/π-mixed systems, or when the selection is incomplete, the two numbers differ: the NOBD columns report the contribution of the selected delocalization channels, not of a symmetry block. Check the NOBD table (--table nobd, star markers) before quoting a π-restricted NICS as “the π-MO ring current”.

How the partition works. The shielding is bilinear in density-like matrices, σ = Tr[D·H⁽¹¹⁾] + Tr[D^(B)·H⁽⁰¹⁾], so any additive split of the density D and of its field derivative D^(B) splits σ exactly. runEDDB uses the EDDB fraction operator Ω, the symmetric solution of D·Ω + Ω·D = 2·D_del in the analysis block, which reproduces the EDDB density exactly (½(DΩ + ΩD) = D_del) and is diagonal in the NOBD basis (eigenvalues n_k/⟨v_k|D|v_k⟩). The response density is antisymmetric and connects occupied with virtual orbitals only, so its trace with any operator living in the occupied space (Ω, a set of NOBDs, a set of MOs) is taken with the occupied-side rule: row and column of every occupied–virtual term are assigned to the occupied side and the occupied–occupied block is counted once (G'' = 2G − G_oo). With P the projector on the occupied space this gives Tr[P·G''] = σ exactly, and the table reads EDDB_G = Tr[Ω·G''], EDLB_G = Tr[(P_V − Ω)·G''] and other = σ − Tr[P_V·G''] (the core and Rydberg share of the occupied space, chosen with -b). This is the convention of natural chemical shielding analysis and of the ipsocentric model: the whole induced current belongs to the occupied orbitals that carry it.

Input — GIAO lane (Gaussian, ORCA or Turbomole).

• Gaussian: run with the NMR keyword (GIAO is the default) and formchk the checkpoint. The .fchk then contains the section Magnetic Field P1 (GIAO) with the three field-derivative density matrices; nothing else is needed (int=nobasistransform is recommended for generally contracted bases).
• ORCA 6: request the shieldings and keep the response densities, then export the JSON with densities:

  ! B3LYP def2-TZVP NMR KeepDens
  %eprnmr NMRShielding 2  Origin GIAO end
  
  # BaseName.json.conf  (same directory and base name as the .gbw; orca_2json has no -config flag)
  { "MOCoefficients": true, "Basisset": true, "Densities": ["all"] }
  $ orca_2json BaseName.gbw

  KeepDens makes ORCA store the magnetic-field response densities (pbfield_el_0/1/2) in its .densities container; without the .json.conf file orca_2json exports no densities at all. runEDDB reads them from the .json together with the MOs and the basis set and produces exactly the same table as from a Gaussian .fchk.
• Turbomole 8: add $gimic to the control file and run mpshift after the SCF; mpshift then writes the file XDENS (the density and the three field-derivative densities in the Cartesian AO basis) next to the usual job files. Build the runEDDB input as cat control coord basis mos > mol.tmol and keep XDENS in the same directory (or pass --nmr-xdens FILE). Every functional that mpshift supports is available, including the local hybrids (e.g. $dft functional sclh22t gridsize 4; mpshift refuses the m-grids). runEDDB checks the XDENS density against the density it builds from the mos file and refuses inconsistent files. Closed-shell references only (mpshift NMR is a closed-shell module).

Restrictions: HF and DFT references (RHF; UHF with the spin-summed response from Gaussian or ORCA), Gaussian .fchk, ORCA .json or Turbomole .tmol+XDENS input, global EDDB_G only. With -b nao (full NAO basis) the total reproduces the program’s own isotropic shielding to all printed digits — a built-in sanity check.

Input — CTOCD-DZ lane (Gaussian only). Run Gaussian with NMR=CSGT IOp(10/33=2) (IGAIM or SingleOrigin work as well) and keep the .log: runEDDB reads the six CPHF response densities (three momentum and three field responses) from <stem>.log (or --nmr-log FILE) and evaluates the analytic CTOCD-DZ shielding (Lazzeretti–Zanasi; the ipsocentric model of Steiner and Fowler): the gauge origin moves with the point of observation, d(r) = r, so there is no diamagnetic term at all — the whole shielding is an induced current. The gauge in use is shown in the table title. Note that this is not the CSGT number Gaussian prints (Keith and Bader move the gauge origin to the nuclei near the atoms and integrate numerically), and that DZ shieldings at nuclei converge slowly with the basis set. The table therefore prints the fulfilment of the translational sum rule (100 % = complete basis) as a diagnostic: NICS away from nuclei are reliable from about 80 % on, nuclear shieldings need ≥ 90 %. Any large, polarised basis will do — spherical or Cartesian alike (Cartesian shells are handled exactly through their spherical + r² decomposition); the sum rule tells you whether it was enough.

Syntax.

$ runEDDB -i benzene.fchk --table nmr                        # shielding table of all nuclei
$ runEDDB -i benzene.fchk --nics-rings 6 -a                    # + NICS(0), NICS(1), NICS(1)zz of every 6-ring,
                                                              #   EDDB column = the π-type NOBDs only
$ runEDDB -i benzene.fchk --nics-ring 1-2-3-4-5-6-1 --nics-scan 2.0,0.25
$ runEDDB -i benzene.fchk --nics 0,0,0,0,0,1.5                 # explicit probe points (Angstrom)
$ runEDDB -i tms.fchk --table nmr                             # writes tms.nmrref
$ runEDDB -i benzene.fchk --table nmr --nmr-ref tms.nmrref     # + chemical-shift table, partitioned
$ runEDDB -i benzene.fchk --table nmr --nmr-log benzene_csgt.log   # CTOCD-DZ lane (Gaussian NMR=CSGT)

Every --nics* and --nmr-ref option implies --table nmr. For each ring the table lists NICS(0), NICS(1) and their zz components along the ring normal (the normal follows the right-hand rule of the atom order and is printed in the ring header); NICS = −σ (Schleyer’s sign convention). Under the caption every ring is listed by its atoms in the closed-cycle notation of -p and --ring-detect (R1 1-2-3-4-5-6-1), and the ring-current table lists the extra cycles of the --rcs / --bcs lists the same way, with the numbering continued (R4 5-4-13-...-5), so a line can be pasted back into a command. NICS(1) is the mean of the two probes at ±1 Å, which are always both computed (the loop-current fit uses both) and is the only NICS(1) printed, in every mode; for a planar ring the two faces are equal, and for a non-planar or one-sided situation --nics-scan or --nics x,y,z give any point above or below the ring. --nics-scan R,step scans from −R to +R along the normal (step 0 or omitted: only the two points ±R). The ring list of --nics-rings uses the same ring-graph cutoff as -r (--ring-thresh).

Reading the tables. By default every NMR table has four numeric columns with five decimals: the headline value (σ_iso, NICS, δ_iso or the shielding of an MO set), its EDLB_G and EDDB_G parts and other (the core + Rydberg share of the occupied space, zero with -b nao; the frontier-MO table adds rest). The terminal never shows more than these four columns. Under -v the same tables are written again to .NMR_DIAG.txt with eight columns at full precision: each EDLB_G / EDDB_G group then lists its total first, followed by the two parts that separate the density from the response; there the NICS table also lists the two faces NICS(+1) and NICS(−1) and the centroid and normal of every ring, and the file adds the probe-resolved sum rule, the loop-current fit and the bond-current edge terms (everything that is a check on the calculation rather than its result). The .nmrref, .NMR_NOBD.txt and .NMR_MO.txt files keep the full precision as before. The two parts are:

• GIAO: dia = Larmor circulation of the unperturbed density, evaluated with the gauge origin at the probe (the nucleus for σ, the probe point for NICS — the classic Lamb convention); para = the CPHF response of the density, referred to the same origin. The total is gauge-invariant by construction; with the probe as origin the dia / para split, and with it the EDLB_G / EDDB_G partition, is also independent of where the molecule sits in the coordinate frame (the raw GIAO integrals refer the split to the coordinate origin, which would make it change under a translation of the input geometry). This is the bookkeeping of natural chemical shielding (NCS) analysis and of GIAO-based current-density programs: every occupied orbital, deep ones included, carries a Larmor term proportional to its density.
• CTOCD-DZ: trans = translational response (momentum operator, Δλ = ±1, diatropic-type — it replaces the Larmor term); rot = rotational response (angular momentum, Δλ = 0, paratropic-type). The angular momentum is taken about the probe (the nucleus for σ, the probe point for NICS), so the trans / rot split does not depend on the common origin of the CPHF run; for NICS(0) and for every zz row on the ring axis this coincides with the ring-centre origin of Steiner and Fowler. The EDLB_G / EDDB_G partition itself is origin-independent in this lane whatever the origin. These are the selection rules of the ipsocentric model: only orbitals with symmetry-allowed low-lying excitations contribute, so in benzene the two e_1g π-NOBDs carry the whole diatropic ring current and the deep a_2u π-NOBD contributes nothing, whereas in the GIAO table the same a_2u NOBD shows a large dia entry. Both tables sum to the same NICS.

Which analysis space (-b) for which probe, and why the answer looks backwards. The recommendation below (NVB at the nuclei, -b nao for NICS) reads like the opposite of what one expects: the natural valence basis has no core and no Rydberg orbitals, so it seems unfit for anything close to a nucleus and natural for a probe far away. That reasoning rests on a wrong premise about what -b does. It does not choose the orbitals with which the shielding is computed. The shielding at every probe is always evaluated with the full density and the full response in the complete AO basis, all electrons and all functions. -b only chooses the analysis block: which NAOs are partitioned into EDLB_G and EDDB_G. The share of the occupied space that lies outside the block (core and Rydberg NAOs) is not discarded but collected in the other column, and the total σ is identical for every choice of -b. NVB is therefore not "unfit" near a nucleus; it merely labels the core contribution differently. At a nucleus the Larmor circulation of the 1s shell is large (benzene carbon: +57 ppm in the other column of the CTOCD-DZ table at def2-TZVP, about 200 ppm in the GIAO Larmor bookkeeping; hundreds of ppm for O or F) and it is a chemically distinct term that deserves its own column. With -b nao the core NAOs enter the block, they are perfectly localized, and their contribution dissolves into EDLB_G, which then no longer means "localized bonds and lone pairs" but "localized bonds, lone pairs and the core". Hence NVB at the nuclei: not because it describes the core better, but because it keeps it apart. Far from the nuclei the core contribution is negligible (benzene: 0.2 ppm at NICS(0), −0.3 ppm at NICS(1)); the issue there is a different one. The occupied space projected onto the valence block is not complete, a small part of the response of the occupied density sits on NAOs classified as Rydberg, that part depends on the basis set and on where the response to a distant point is anchored, and under NVB it lands in other or blurs the EDLB/other border. With -b nao the block contains every NAO, its projector equals the projector onto the occupied space, other ≡ 0 identically and the NICS is resolved exactly into EDLB_G + EDDB_G. The EDDB_G column is stable either way (benzene, RHF, NICS(0) of the π-NOBDs: −21.4, −22.5, −21.6 and −22.3 ppm for 6-31G, 6-31G*, def2-TZVP/NVB and def2-TZVP/NAO), so the recommendation is about the cleanliness of the remainder, not about the EDDB_G number. The rule is thus the reverse of the intuition: where the core is large we want to see it separately, so NVB; where it is small we want an exhaustive partition with no remainder, so NAO.

Equivalent loop current (I_loop) — a ring-size-independent measure. NICS values of rings of different size are not comparable with each other: the probe 1 Å above a five-membered ring sits closer to its current loop than above a six-membered one, and the induced field of a loop scales with its area, so a larger NICS does not mean a stronger current. The current itself does not depend on the ring area, so below every NICS table runEDDB fits to the NICS_zz(z) profile along the ring normal the field of a thin circular loop of radius R carrying the current I per unit external field,

NICS_zz(z) = −2π I R² / (R² + z²)^3/2 ⇒ I [nA/T] = −NICS_zz(z) [ppm] · (R² + z²)^3/2 / (2π R²) (R, z in Å)

with the sign chosen so that a diatropic (aromatic) ring gives a positive current. Since the Ampère integration below measures the ring current exactly, the loop fit is a diagnostic and no longer appears on the terminal: under -v the file .NMR_DIAG.txt lists one I_loop row per NICS ring in its RING CURRENT STRENGTHS (nA/T) table (the total, EDLB_G and EDDB_G currents, under -a/-s: the selected NOBDs, for the ring radius fixed to the mean centroid–atom distance) and the table LOOP-CURRENT FIT DETAILS adds, per column, I(R_ring) (a linear least-squares fit of I alone) and I(R_fit) with the radius fitted as well, each with the rms deviation of the fit in ppm. The points used are NICS(±1) and every --nics-scan point with |z| ≥ 1 Å (nearer the ring plane the thin-loop picture does not hold); without a scan only the two NICS(±1) points exist and just I(R_ring) is fitted, so use --nics-scan 4,0.25 (or finer) whenever R_fit matters.

• What to read. The EDDB_G column is the ring current of the delocalized electrons and is the number to compare between rings and with current-density integrations (GIMIC): benzene at B3LYP/def2-TZVPP gives I_loop^EDDB = 11.3 nA/T (rms 0.3 ppm, R_fit = 1.74 Å, i.e. just outside the nuclei, as expected for a π current) against the GIMIC ring current of 11.8 nA/T; the same fit orders benzene ≫ pyrrole > thiophene > furan (11.3, 7.0, 5.2, 4.4 nA/T), whereas their NICS(1)_zz values are all within a few ppm of each other.
• What not to read. The EDLB_G profile (local σ and lone-pair circulations) is not that of a loop: its R_fit usually runs to the edge of the search range (marked *) and its I is not a ring current. The total fit inherits this and is correspondingly less well defined (benzene: rms 3.9 ppm at R_ring); quote the EDDB_G current.
• Units. With NICS in ppm and lengths in Å the prefactor of the loop formula is exactly 1, so the number printed is directly the current per unit field in nA/T (positive = diatropic, negative = paratropic).

Ring and bond current strengths (--rcs [LIST], --bcs [LIST]) — the GIMIC number without a current density. The loop current above is a model fit; the current itself can be measured exactly from the same shielding machinery. GIMIC defines the bond current strength as the flux of the induced current density through a plane perpendicular to the bond, from the ring centre through the bond midpoint to far outside the molecule. By Stokes' theorem that flux is the line integral of the induced field around the boundary of the plane, μ₀ I = ∮ B_ind·dl (Ampère's law), and the induced field at any point is minus the shielding tensor there, B_ind/B₀ = −σ·n for the external field along the ring normal n. So --rcs / --bcs place shielding probes on the four edges of a rectangle through every ring bond (the ring axis, the bottom and top edges at ±h, and the outer edge L beyond the bond midpoint; for a bond shared by two rings the far edge is the other ring's axis), integrates them with Gauss–Legendre quadrature and prints, in ppm·Å turned into nA/T by the factor 1/4π,

I_bond [nA/T] = (1/4π) ∮ (−σ·n)·dl [ppm Å], positive = diatropic,

as I_bond rows of the same RING CURRENT STRENGTHS table (total, EDLB_G, EDDB_G; under -a/-s: the selected NOBDs), bond by bond, followed by the I_cyc row of the cycle (with --bcs), or as the I_cyc rows alone (with --rcs); the diagnostics file of -v adds the trans/rot (dia/para) rows under every bond, the mean and spread over the outward bonds, and the table BOND-CURRENT EDGE TERMS with the four line-integral terms. Every edge probe carries the full EDLB/EDDB partition, so the current of the delocalized electrons through a bond is obtained exactly, with no loop model and no current-density grid.

• What to read. For a divergence-free current every outer bond of a ring carries the same flux, which is the ring current of that ring (benzene at RHF/def2-TZVP, CTOCD-DZ: 13.04 nA/T total, 10.22 nA/T EDDB_G, the same for all six bonds; the GIAO lane of the same molecule gives 12.96 / 10.14). A bond shared by two rings carries the difference of the two ring currents (azulene: 6.6 nA/T around the five-membered ring, 5.0 around the seven-membered one, 1.6 through the shared bond, sign by ring). The spread over the outer bonds measures how far the finite-basis current fails to be divergence-free (0.00 for benzene, 0.05–0.07 nA/T on the GFN2-xTB lane). The result does not depend on the box: L = h = 3, 5 and 8 Å give 13.05, 13.04 and 13.04 nA/T for benzene, and 12 or 24 Gauss points per edge give the same digits.
• Why I_loop and I_bond are different numbers. They answer the same question with two different instruments, and for benzene they disagree by about 10 % in the EDDB_G column (11.6 vs 10.2 nA/T at RHF/def2-TZVP; total 15.9 vs 13.0). I_loop is a model: it assumes that the whole delocalized current is an infinitely thin circular loop of exactly the ring radius R_ring, and it reads the current off the NICS_zz values at |z| ≥ 1 Å through the loop formula above. I_bond is a measurement: Ampère's law makes no assumption about where the current flows, and in the limit of a large box it reduces to the integral of the NICS_zz(z) scan along the ring axis, I = (1/4π) ∫ −NICS_zz(z) dz, exactly and for any current distribution. The two agree only if the true profile has the loop shape (R² + z²)^−3/2, because that profile integrates to exactly 4πI. The real π current does not: it flows at R_fit ≈ 1.7 Å, outside the nuclei, with a finite width above and below the plane, so its NICS_zz falls off more slowly than the thin-loop formula with R_ring = 1.4 Å predicts, and the fixed-radius fit compensates with a larger current; the excluded |z| < 1 Å region adds to the bias. That the discrepancy is only 10 % is a point in favour of the simple model, and the free-radius fit closes most of the gap; but it costs nothing beyond the NICS(±1) probes only because it is a model, and that is why it lives in the diagnostics file. The number to quote or to compare with GIMIC is I_cyc (or I_bond, its EDDB_G column) from --rcs / --bcs; the two must never be mixed in one table.
• Two options, one integral: I_cyc and I_bond. --rcs reports ring current strengths: without a list every NICS ring, with a list only closed cycles in the -p syntax (first atom repeated at the end, e.g. --rcs 1-2-3-4-5-6-7-8-9-10-1; a cycle defines its own centroid and normal and need not be a NICS ring; a bond or an open chain is an error). For every cycle it prints one row, I_cyc, the minimum over its outward bonds: since any decomposition of the current into loop currents gives the loop around the whole cycle at most the current of its weakest bond, I_cyc is the largest current that a coherent circulation around the entire cycle can carry, the current analogue of the CDI (the minimum EDDB_P population along a cycle). --bcs is the same integral bond by bond: without a list every bond of every NICS ring; with a list, e.g. --bcs 1-2-3-4-5-6-7-8-9-10-1,4-9,1-2-3, every bond of a closed cycle (followed by its I_cyc row, so --bcs 1-2-3-4-5-6-1 prints what --rcs 1-2-3-4-5-6-1 prints plus the six bonds), the bonds of an open chain (1-2-3 = 1-2,2-3, no ring value) and single bonds. Every bond is integrated from the centre of the smallest cycle that contains it (a NICS ring first, then a listed cycle; a bond outside every cycle is skipped with a warning), so the perimeter of naphthalene is measured bond by bond through its two rings and a plane never starts on a bond the cycle straddles (a warning is printed when a cycle's centroid falls on an atom or bond). The two options may be combined; the cycles are numbered R1, R2, ... in one sequence, NICS rings first. In a monocycle I_cyc equals the ring current (Kirchhoff: the same flux through every bond, so the heterogeneity of conjugation is invisible in the net currents); in a fused system it separates global from local circulation: naphthalene's perimeter carries I_cyc = 7.0 nA/T with 0.0 through the central bond (GFN2-xTB lane), azulene's perimeter 4.9 nA/T while 6.6 circulate around the five-membered ring, 4.9 around the seven-membered one and 1.6 through the shared bond.
• Currents are not scalars: shared bonds, enclosed rings and the mesh picture. A bond current is a flux through an oriented plane, and its sign belongs to the plane, not to the bond: the plane starts in the centre of a cycle and the current is positive when it circulates diatropically around that cycle. A bond shared by two rings therefore appears twice with opposite signs (naphthalene 4–9: +0.00 under R1 and −0.00 under R2; the hub bonds of coronene: −3.34 nA/T under the hub, +3.34 under the outer rings), and the number itself is the net current of that bond, the difference of the two ring currents, never the current of one ring. The clean way to think about it is Kirchhoff's mesh picture: for a planar network the current decomposes uniquely into one loop (mesh) current per ring; an outward bond carries the mesh current of its ring, a shared bond the difference of two mesh currents. The mesh current of a ring is the flux through any plane from its centre, through one of its bonds, out to the vacuum. By Stokes' theorem it does not matter what the plane crosses on the way (another ring, its far bond, the hydrogens), and this is the quantity that I_cyc reports. For a ring with an outward bond the ordinary plane (L beyond the bond midpoint) already ends in the vacuum, so I_cyc is the minimum over those bonds. For a ring whose every bond is shared (the hub of coronene, the inner rings of larger nanographenes) no ordinary plane reaches the vacuum, so the program adds, per bond, a plane that continues through the neighbouring ring to L beyond the last atom in that direction; I_cyc is the minimum over these (they are the vac rows of the BOND-CURRENT EDGE TERMS table in the -v diagnostics file; the environment variable RUNEDDB_BCS_VACUUM=1 builds them for every ring as a consistency check: in azulene the plane through the shared bond, continued through the seven-membered ring, gives 6.57 nA/T against 6.55–6.60 through the outward bonds). Coronene on the GFN2-xTB lane shows why the distinction matters: the outer rings carry I_cyc = 9.74 nA/T, the hub 6.41 nA/T, and the hub bonds −3.34 nA/T seen from the hub (= 6.41 − 9.74; the same for L = h = 5 or 8 Å). The textbook “paratropic hub current” of coronene is the net current of the hub bonds, an I_bond; the mesh current of the hub is diatropic and is the I_cyc. Both are printed, under different names, and both are correct; when a number is quoted, say which one it is. The loop fit I_loop at the hub centre (2.0 nA/T) reads NICS_zz produced by all seven loops at once and is not comparable with either.
• Cost and options. Each bond adds 4n probes (the axis edge, 2n probes, is shared by all bonds of a cycle): benzene with n = 12 needs 312 probes. The contour probes are independent and run in parallel over probes (each with its own work arrays; the thread count follows -n and a 4 GB memory budget): benzene at def2-TZVP takes a minute on one core and about five seconds on sixteen. --bcs-box L,h,n sets the box and the quadrature; the -v diagnostics file adds the trans/rot (dia/para) split, the four edge terms and, on the GFN2-xTB lane, the EDDB current divided by the in-plane translational sum rule (EDDB/S_par), the same basis-completeness correction as for I_loop.

What a ring current is, and is not, evidence of. runEDDB reports I_loop, I_bond and I_cyc as what they are: magnetic properties of the molecule in a field, the same observables that determine NICS, proton shifts and the magnetizability exaltation. They are not an aromaticity index, and the partition into EDLB_G and EDDB_G does not make them one. The sign of a ring current is set by the symmetry of the frontier π–π* excitations (translationally allowed: diatropic; rotationally allowed: paratropic, the ipsocentric picture), its magnitude largely by the HOMO–LUMO gap. Neither is a measure of how far the bonds are delocalized, nor of the stabilization that IUPAC requires together with cyclic delocalization for a ring to be called aromatic. The familiar correspondence between the 4n+2 / 4n electron counts and diatropic / paratropic currents is a property of free-electron (Hückel-type) annulene models, where the two follow from the same orbital symmetry; outside that family it fails in both directions: bond-alternated (localized) annulenes still sustain a diatropic current, and the four phenalenyl systems of the tutorial share one delocalized 12-membered perimeter (the same EDDB maps and CDI) while the carbon-centred ions are diatropic and the isoelectronic B- and N-centred molecules paratropic. The π electrons of benzene are not more mobile than those of hexatriene; they are bound by the cyclic potential, and it is the degeneracy of the cyclic problem, the coherence it enforces around the ring and the wide gap it opens that stabilize the molecule, while the current stays moderate. A strong paratropic current signals a small gap, not destabilization. Use EDDB_G/EDDB_P and the CDI for delocalization, energies for stabilization, and the currents for the magnetic response; inferring bonding or aromaticity from a current alone has led to well-documented over-interpretations, down to the heaviest elements (the tri-thorium cluster: Angew. Chem. Int. Ed. 2022, e202204337; RSC Adv. 2023, 13, 34224).

Orbital resolution — files, not terminal tables. --table nmr prints only the two complete tables (nuclei and NICS). Two variants print the same tables and additionally write a sidecar file next to the input (honouring --out-prefix / --out-dir):

• --table nmr-nobd → <stem>.NMR_NOBD.txt: the EDDB_G column of every nucleus and NICS probe resolved by NOBD, one block per probe closed by the sum over the listed NOBDs (= the EDDB_G column of the terminal table). Without -a / -s the file lists all non-zero-occupied NOBDs for every probe and can become very large (hundreds of NOBDs × hundreds of NICS scan points); use it together with -a / -s, which restricts the file to the selected NOBDs.
• --table nmr-mo → <stem>.NMR_MO.txt: the full shielding of every probe resolved by occupied canonical MO (index, spin, occupation, orbital energy, the two parts and their sum), from the Gaussian .fchk or the ORCA .json MO coefficients. The partition assigns row and column of the trace to each occupied orbital and removes the double-counted occupied–occupied block, so the per-MO sums reproduce the totals exactly (the closure is printed under -v and in the file header). For UHF the file has alpha and beta sections; since Gaussian stores only the spin-summed response density, the exact partition over the natural orbitals of the total density is handed to the spin-orbitals by their weights |⟨NO|MO⟩|²/n_NO, which sums back exactly. ROHF NMR is not available in Gaussian. -a / -s select NOBDs, not MOs, so this file always resolves the full σ; under a dissection an extra w_sel column gives the weight of every MO inside the selected NOBD subspace (for a planar π system in the full NAO basis the σ MOs have exactly 0 and the π MOs exactly 1, i.e. there the NOBD dissection coincides with a π-MO selection).

Frontier orbitals (--table nmr-homo[N], nmr-somo, --nmr-mos LIST). In the ipsocentric (Steiner–Fowler) picture the ring currents are sums of occupied-orbital contributions dominated by the frontier orbitals. These options print one extra table with the shielding carried by a chosen set of occupied MOs: the set's density D_F and its share of the response density (row and column of the trace assigned to the set, minus the double-counted occupied–occupied block) are partitioned exactly like the total into EDLB_G, EDDB_G and other, and rest is everything carried by the other occupied MOs (σ_iso = σ(set) + rest exactly). The header lists every MO of the set with its energy, occupation and ED/EDLB/EDDB populations (the --table mo row), so the orbital's share of the delocalized density and of the magnetic response are read side by side; a NICS block follows when probes were requested. nmr-homo takes the HOMO, nmr-homoN the levels HOMO-N..HOMO (N ≤ 9), nmr-somo every singly occupied MO, --nmr-mos 5:12,20 an explicit set (numbered as in --table mo). For UHF a level pairs α_k with β_k, so for a triplet HOMO and HOMO-1 are the two singly occupied α orbitals and HOMO-2 is the first α/β pair; a5:12,b5:11 picks one spin. A degenerate HOMO (benzene's e_1g pair) triggers a warning suggesting the whole shell (nmr-homo1). Prefer the CTOCD-DZ lane for a ring-current interpretation: there the occupied-orbital contributions are exactly the ipsocentric ones, whereas GIAO carries an occupied–occupied response block that is shared symmetrically.

Chemical shifts (--nmr-ref) — step by step. Nothing has to be subtracted by hand:

1. Two NMR single points, one level of theory. Run the reference molecule (TMS, or the solvent molecule if you report shifts against a solvent signal) and the molecule of interest with the same method, basis set and gauge, and convert both checkpoints:

   #P B3LYP/def2-TZVP NMR=GIAO           ! tms.gjf and benzene.gjf, same route line
   $ formchk tms.chk tms.fchk
   $ formchk benzene.chk benzene.fchk

   For the CTOCD-DZ lane use NMR=CSGT IOp(10/33=2) in both jobs and keep the .log files.
2. Reference first. Run runEDDB on the reference: runEDDB -i tms.fchk --table nmr. Every --table nmr run writes <stem>.nmrref (honouring --out-prefix / --out-dir) with the element-averaged shieldings and their partition columns — compute it once, reuse it for every molecule.
3. Then the molecule of interest with the reference attached: runEDDB -i benzene.fchk --table nmr --nmr-ref tms.nmrref. Below the shielding table runEDDB prints a CHEMICAL SHIFTS table in which every column is already the difference δ = σ(ref) − σ for each nucleus whose element occurs in the reference: δ_iso and its EDLB (total|dia|para), EDDB (total|dia|para) and other parts. A positive entry means that component deshields the nucleus relative to the reference, a negative one shields it.

Variants and pitfalls.

• Explicit reference values. --nmr-ref H=31.7,C=185.4 (e.g. experimental or published shieldings) gives δ_iso only, no partition columns.
• Implied table. --nmr-ref switches --table nmr on by itself.
• Gauge must match. GIAO with GIAO, CTOCD-DZ (--nmr-log) with CTOCD-DZ; the .nmrref file records its gauge, and on a mismatch runEDDB warns and prints δ_iso only, because the partition columns of the two gauges are not comparable.
• Same analysis space. Use the same -b in both runs: the reference file stores the partition columns of the space it was computed in (NVB by default), and the δ columns are differences of those columns.
• No reference under dissection. With -a or -s no .nmrref is written and a --nmr-ref is ignored with a warning: a partition restricted to selected NOBDs is not a chemical-shift reference, and there is nothing to select in TMS. If you want both the shifts against TMS and the π dissection of the same molecule, run it twice — once with --nmr-ref and once with -a.
• Element averaging. The reference values are averaged over all atoms of an element in the reference molecule (the four methyl carbons and twelve protons of TMS); a nucleus whose element is absent from the reference is skipped in the δ table.
• Sign convention. δ = σ(ref) − σ column by column: a positive entry means that component deshields the nucleus relative to the reference, a negative one shields it.

Example (benzene vs. TMS, RHF/6-31G, GIAO). δ(H) = 7.85 ppm = 10.09 (EDLB: 4.70 Larmor + 5.39 response) − 2.27 (EDDB: −4.94 Larmor + 2.68 response) + 0.03 (other): the deshielding of the benzene protons comes from the C–H bonds and the σ framework, while the π density shields the protons in the GIAO bookkeeping; the π ring current is seen where it belongs — in NICS: NICS(0) = −8.6 ppm with −21.4 ppm from the three π-NOBDs, NICS(1)zz = −29.0 ppm with −24.9 ppm from π. In the CTOCD-DZ lane the π-NOBD contribution to the proton shielding is essentially zero and the e_1g pair alone carries the ring current.

Which lane for which question. The totals are gauge-independent; the partition is not. For the shielding of nuclei and for chemical shifts use the GIAO table (fast basis-set convergence, local “atomic” bookkeeping). For ring currents and their NICS — the π/σ question, the four-/two-electron rules of aromatic and antiaromatic rings — the CTOCD-DZ table is the physically transparent one. Reporting both side by side turns the gauge dependence of the “π contribution” into a result rather than an ambiguity. The default analysis basis is NVB; for NICS probes far from nuclei also check -b nao, because the response at a nucleus-free point spans the Rydberg NAOs and the EDLB/other split (not the EDDB column, which is stable to about 1 ppm) depends on the basis size.

Approximate magnetic partition from any SCF wavefunction (automatic; --nmr-kernel). When --table nmr, a NICS option or --bcs is requested on a wavefunction file that carries no NMR response data (a Gaussian .fchk from a plain SCF job, an ORCA .json or .molden without the pbfield densities, a Turbomole .tmol without XDENS, a GAMESS log), runEDDB does not stop: the same machinery that produces NICS on the GFN2-xTB lane below runs on the canonical orbitals of the file. The response densities D^(L) and D^(p) are built as the sum over states of the input MOs (all occupied and virtual MOs with their energies are required, SCF reference only; post-HF densities and NBO .49 archives still stop with an error), and the CTOCD-DZ tables, the EDLB_G/EDDB_G partition, the loop and bond currents and the frontier tables follow unchanged. The run says so everywhere: the stage line reads Computing APPROXIMATE NMR partition from the input MOs (CTOCD-DZ, B3LYP)... with the method it used, every table caption carries [APPROX., UNCOUPLED], a line under the first caption names the exchange kernel (exchange kernel: B3LYP a_x = 0.20 (.fchk header)), a counted WARNING states what is approximated and by how much the currents are expected to be off, the -v block has Lane: and Functional: lines, and the .nmrref written by such a run is marked so that --nmr-ref refuses to mix it with an exact one. When the file does carry NMR data the exact partition is always computed from them and a given --nmr-kernel is ignored with a WARNING. Use the approximate lane for signs, trends and the partition of any SCF file you already have; for numbers to be quoted run the NMR job. The worked examples below use the GFN2-xTB lane; the same tables appear for a wavefunction file without NMR data, with the [APPROX., UNCOUPLED] caption.

What is approximated, and the role of the method. For a purely imaginary perturbation the Coulomb and local exchange–correlation kernels vanish, so for pure functionals (LDA, GGA) the uncoupled sum over states is the coupled response and the tables agree with a full NMR job to every printed digit (PBE/def2-TZVP: benzene, naphthalene, pyrrole and planar cyclooctatetraene, ratio 1.000); for meta-GGAs a small current-dependent kernel survives (TPSS: 2 %). For hybrids and Hartree–Fock the exchange coupling is missing from the plain sum over states, and it is large: the ring currents of aromatic rings come out 1.25–1.4 times too small for B3LYP and PBE0, 2.2–2.5 times too small for HF, and much more for 4n rings (planar cyclooctatetraene: 2.7, 3.2 and 11.5), so no single factor per functional exists. runEDDB therefore adds an sTDA-type exchange kernel whose only input is the exact-exchange fraction a_x of the functional (0.20 B3LYP, 0.25 PBE0, 1.0 HF). In a Gaussian basis the Löwdin-charge kernel of the xTB lane is far too strong and its required scale changes with the basis (0.40, 0.125 and 0.06 for def2-SVP, TZVP and QZVPP), because the transition charges spread over polarisation and diffuse functions; runEDDB builds the kernel in the natural minimal basis of its own NAO construction (the core and valence NAOs, the basis-independent minimal basis of the molecule), where it needs no scaling at all and gives the same numbers in every basis. Validated on benzene, naphthalene, pyrrole and planar cyclooctatetraene at B3LYP, PBE0 and HF in def2-SVP, def2-TZVP and def2-QZVPP against NMR=CSGT, the ring currents of the hybrids are reproduced within 2–5 % for every ring type (the 4n ring included, where the uncoupled value alone is off by 2.7–11), whereas Hartree–Fock stays 10–13 % low for the hydrocarbons and 20 % low for pyrrole: the missing term is exact exchange, four to five times larger at HF than in a hybrid, so any imperfection of the monopole kernel scales with it, and the exchange involving the nitrogen lone pair is what the atomic kernel represents worst. This is not the π-donor deficit of the GFN2-xTB lane, which sits in the tight-binding Hamiltonian and appears with or without coupling. The results of the three basis sets agree with each other within 0.03. Treat the coupled numbers as a good estimate, not as a replacement for the NMR job.

Telling runEDDB the method (--nmr-kernel). The only thing the kernel needs is which method produced the orbitals, and runEDDB usually finds it by itself: the second line of a Gaussian .fchk records it (SP RB3LYP def2TZVP, also RHF/UHF/ROHF), failing that the route line of the .log next to the file; an ORCA .json records no functional at all (orca_2json writes HFTyp: RHF even for B3LYP), so for .json and .molden input the ORCA output <stem>.out (or .log) next to the file is read, whose Fraction HF Exchange line is the exact-exchange fraction itself (pure, hybrid, range-separated and double-hybrid runs and Hartree–Fock are told apart from the same block), so keep the output next to the export; a Turbomole .tmol carries its $dft section (functional b3-lyp; no $dft means HF); a GAMESS log echoes DFTTYP= (NONE means HF). The -v block reports the method and where it came from. When the method cannot be determined (no output next to a .json or .molden file, a functional missing from the table) the default DFT kernel a_x = 0.20 is used and the line under the first caption says so and how to be specific: exchange kernel: default DFT a_x = 0.20, method not identified (--nmr-kernel NAME|hf). The default is DFT-like and not Hartree–Fock on purpose: the HF kernel applied to DFT orbitals (small gap) is catastrophic (benzene B3LYP orbitals with a_x = 1: a loop current of −199 nA/T with the wrong sign, planar cyclooctatetraene unstable), whereas the DFT kernel applied to HF orbitals merely underestimates (benzene HF: 6.1 instead of 11.1 nA/T). If you know the method, say it: --nmr-kernel b3lyp (or pbe0, hf, tpssh, m06-2x, …) overrides whatever was detected, on your responsibility; --nmr-kernel @list prints the table; names are matched case-insensitively, the spellings of the supported programs are accepted (pbepbe, pbe1pbe, wb97xd, b3lyp/g, b3-lyp, b-p, s-vwn, b3lypv1r, camb3lyp, …) and dispersion suffixes (-d3, -d3bj, -d4) are ignored. For a pure functional the table gives a_x = 0 and the plain uncoupled response is used, which is exact. Range-separated functionals are accepted with their short-range fraction and double hybrids with their hybrid part, each with a WARNING that the long-range exchange or the MP2 part of the response is not represented by the monopole kernel. --nmr-kernel dft asks for the default explicitly; an explicit number (--nmr-kernel 0.30) sets a_x directly, and a_x,y also scales the kernel (research use; the validated scale is y = 1). Adding a functional is one line in the table (functional_ax in the response module; the developer notes docs/NMR_KERNEL_FUNCTIONALS.md list every functional with its a_x and its spelling in Gaussian, ORCA, Turbomole and GAMESS).

Method families at a glance. The table lists the families, the functionals currently in the built-in table with their exact-exchange fraction a_x, what the lane then computes and how well it was validated. The lane itself runs for every family whenever the input carries no NMR data; --nmr-kernel only tells it the method.

Approximate NICS from the GFN2-xTB lane (geometry input). With a geometry file as input (.xyz, .mol/.sdf, .mol2, .pdb or a built-in @model) any --nics* option, or --table nmr together with one of them, makes runEDDB compute NICS from its own GFN2-xTB wavefunction, with the full EDLB_G / EDDB_G partition, the trans | rot split, -a/-s, the NOBD- and MO-resolved files, the frontier-orbital tables, NICS scans and the equivalent loop current, at the cost of a few seconds:

$ runEDDB -i @arenes/benzene --nics-rings 6 --nics-scan 4,0.25 -v
$ runEDDB -i pyrrole.xyz --nics-rings 5 -a --table nmr-homo

The GFN2-xTB Hamiltonian has no exact exchange and its electrostatics depend on populations only, so an imaginary (antisymmetric) change of the density changes nothing in the Hamiltonian: the magnetic linear response of the model is uncoupled and is obtained exactly as a sum over states of the converged tight-binding orbitals for the angular-momentum and linear-momentum perturbations (spin-polarized runs: per spin, then summed). These two response densities feed the ipsocentric CTOCD-DZ lane unchanged, so every table looks exactly like the one produced from a Gaussian NMR=CSGT log. This is London’s tight-binding ring-current theory with all valence electrons, the ipsocentric gauge and the EDLB/EDDB partition on top.

What is and what is not computed. The GFN2 basis is a minimal valence basis without core electrons. It carries the translationally allowed valence π → π* excitations that make the ring current, but not the Larmor circulation around the nuclei that dominates nuclear shieldings. Nuclear shieldings and chemical shifts are therefore not computed: the nuclei are not probed, no nuclei table and no .nmrref are written, and --nmr-ref, --nmr-log and --nmr-xdens are rejected with an error, as is a run without any --nics* probe. A NOTE above the NICS table repeats this. Meaningful are the NICS values at ring centroids, NICS(±1), NICS_zz, scans and points away from the nuclei, their signs (4n+2 diatropic, 4n paratropic), their partition and their trends along a series.

Magnitude and the sum-rule scaling. The uncoupled minimal-basis response is too small in absolute terms: for benzene NICS(1)_zz is −17.3 ppm against −29.7 at B3LYP/def2-TZVPP (CTOCD-DZ), and paratropic 4n rings come out too paratropic (cyclobutadiene NICS(1)_zz +115 ppm). On a 35-molecule set with B3LYP-D3(BJ)/def2-TZVPP geometries and CSGT references the ratio between the reference and the tight-binding EDDB loop current factorises as k = k_coup × k_rest: k_coup = 1.23–1.37 is the response coupling that the hybrid functional has and GFN2 lacks by construction (recovered by --xtb-coupled, below), and k_rest = 1.1–1.2 for hydrocarbon, azine, borazine and thiophene rings is the residual Hamiltonian/basis deficit; rings with a π-donor N–H or O (pyrrole, furan, imidazole, oxazole, …) have k_rest = 1.4–1.6 and 4n rings k_rest = 0.3–0.6 (their GFN2 gaps are too small). runEDDB also measures the basis part directly: under -v the translational sum rule is printed (detail block and diagnostics file) as a tensor, S_ab = Tr[D^(p_a)∂_b]/N (1 = complete basis), as its principal values and, for every ring, as the in-plane average S_∥ and the normal component. The in-plane momentum responses carry the zz ring current, so dividing a ring-current quantity by S_∥ corrects for the basis:

• Validated so far (B3LYP/def2-TZVPP reference, same geometries): the EDDB_G equivalent loop current I_loop/S_∥ falls within 10–20 % of the reference for benzene, pyrrole, thiophene and furan (12.0 / 5.8 / 4.1 / 4.6 vs. 11.3 / 7.0 / 5.2 / 4.4 nA/T) and reproduces their order; the total NICS(1)_zz/S_∥ is within 2 % for benzene but 13–30 % low for the heterocycles, whose lone-pair (EDLB) contribution the minimal basis under-represents even after scaling. The scaled current is printed as the I/S_par column of the LOOP-CURRENT FIT DETAILS table and the EDDB/S_par column of the BOND-CURRENT EDGE TERMS table, both in the -v diagnostics file, so no arithmetic by hand is needed; the terminal shows the raw currents. On the 35-molecule set a single factor of 1.6 on the raw uncoupled EDDB loop current reproduces the reference within 5 % for neutral six-membered rings and 10 % over all in-sample rings (NICS(1)_zz: 1.4 and 4.5 ppm mean absolute deviation).
• Not to be scaled this way: NICS(0), NICS(0)_zz and isotropic values (their in-plane σ terms near the nuclei are not a uniform basis deficit; benzene NICS(0)_zz/S_∥ = −27.6 vs. −16.6 ppm), and any comparison between rings of different size, for which the loop current is the right quantity anyway.
• Probe-resolved sum rule (-v). In the CTOCD-DZ lane (any input) runEDDB also prints under -v, for every NICS probe, the common-origin Larmor term of the unperturbed density about that probe and s = trans / Larmor for the total, EDLB_G and EDDB_G columns: in a complete basis the translational term equals the Larmor term exactly, so s is the sum rule weighted with the NICS kernel at the very point of evaluation. Read the zz rows (in-plane fields mix σ and π in the partition, so the isotropic s is not a completeness measure). It is a rigorous diagnostic — B3LYP/def2-TZVPP gives s_zz(EDDB) = 0.98 for benzene, pyrrole, furan and thiophene, i.e. a converged π response, while GFN2 gives 0.41–0.85 with a π Larmor term within 1–10 % of the reference (the tight-binding π density is right, the deficit sits in the response). The loop-current block adds the fit of the locally corrected profile NICS_zz(z)/s_zz(z) (I/s(z)); on the four rings it is not more accurate than I/S_par (furan over-corrected, thiophene still low: a Hamiltonian rather than a basis deficit), so treat both scaled columns as heuristics.
• Unscaled tables. All NICS tables print the raw GFN2 values; the scaling is offered only where it has been tested (the loop-current block), and the sum-rule tensor is there for any other use. Unlike a GIAO calculation, the Gaussian-expanded vSTO-nG basis of GFN2 is the same for every molecule, so the deficit is systematic rather than random.

Coupled tight-binding response (--xtb-coupled [a_x]). For a purely imaginary perturbation only exact exchange couples the first-order density (Coulomb and local exchange–correlation kernels vanish on an antisymmetric density change), so the coupled response of a hybrid functional differs from its uncoupled sum over states by an exchange term alone. --xtb-coupled adds that term to the GFN2 lane in the zero-differential-overlap form of Grimme’s simplified Tamm–Dancoff approach (sTDA): the sum over states is iterated with X_eff = X + c_x K[D], K[D]_μν = γ(A_μ,B_ν) D_μν in the Löwdin basis, with the Ohno–Klopman kernel γ = (R² + η_AB⁻²)^−1/2 built from the GFN2 atomic hardnesses (no new atomic parameters) and scaled by one number a_x. The bare flag uses a_x = 0.20, the exact-exchange fraction of B3LYP, which is also what a least-squares fit on the 16 neutral six-membered rings of the set gives (0.22, the same value as a fit on all rings): with it the EDDB loop current of those rings agrees with B3LYP/def2-TZVPP within 7 % and NICS(1)_zz within 3 ppm without any scaling, and the coupling ratios of every aromatic ring (1.23–1.37) are reproduced. Two limits are printed by the program: rings with a π-donor N–H or O are tagged [pi-donor N-H/O] in the ring list of the NICS table (and listed in the -v detail block) (their current stays ~1.5× too small, coupled or not; this is the GFN2 Hamiltonian, which puts the π* orbitals too much on the heteroatom), and for 4n rings the coupling map becomes unstable at a_x = 0.2 (the tight-binding gaps are too small), in which case a WARNING is printed and the uncoupled response is used; a_x ≈ 0.1 reproduces their B3LYP amplification. Under -v the stage reports a_x, the kernel, the number of iterations, the spectral-radius estimate of the coupling map and the amplification of the p and L responses. --xtb-uncoupled restores the default uncoupled lane. The empirical factor of 1.6 remains the better single number for screening because it absorbs the Hamiltonian deficit uniformly; the coupled model is the physical one and shows where that deficit sits.

Limitations of xTB-NICS — read before using the numbers. Validated on 35 molecules against B3LYP-D3(BJ)/def2-TZVPP (CSGT) at the same geometries:

1. Where it works: rings without a π-donor heteroatom — hydrocarbons (benzenoid PAHs, annulenes, macrocycle paths), azines (pyridine, pyrazine, pyrimidine), borazine, thiophenes and thiazoles, tropylium. The EDDB_G loop current scaled by 1.6 is within ~5 % of B3LYP for neutral six-membered rings (10 % over all such rings); with --xtb-coupled the unscaled loop current is within 7 % and NICS(1)_zz within ~3 ppm. Ring-by-ring ordering inside fused systems (pyrene, triphenylene, phenanthrene, porphine) is reproduced.
2. π-donor N–H and O rings (pyrrole, furan, imidazole, pyrazole, oxazole, the pyrrole rings of porphyrins): the tight-binding current is ~1.5× too small, with or without coupling (GFN2 puts the π* orbitals too much on the heteroatom). The program tags such rings [pi-donor N-H/O] in the ring list of the NICS table and names them in the -v detail block. Use them for sign and trend only; the NH ring of porphine even flips sign.
3. 4n (antiaromatic) rings: the sign is always right, the magnitude is not (uncoupled: too paratropic by up to 2×; coupled at a_x = 0.2: unstable, the program falls back to the uncoupled response with a WARNING). Qualitative use only.
4. Anions (C₅H₅⁻, C₈H₈²⁻): ~1.5× too small, like the π-donor rings. Cations (tropylium) are fine.
5. Never read NICS(0), NICS(0)_zz, isotropic NICS or nuclear shieldings from this lane, and never compare rings of different size by NICS — use the loop current.
6. Non-aromatic references (cyclopentadiene, cycloheptatriene, tub COT) come out with the right small magnitude and sign, so the lane does not invent currents where there are none.
7. The lane is a screening tool: seconds instead of hours, systematic (the same minimal basis for every molecule), but a semi-empirical Hamiltonian without core electrons. Confirm any number that matters with a DFT reference through the Gaussian / ORCA / Turbomole lanes of the same partition.

Notes. The response is built in a stage of its own (Building the uncoupled GFN2-xTB magnetic response (L, p); -v reports the number of occupied–virtual pairs and the smallest gap). Fermi smearing and open shells (--xtb-spin, --xyz-mult) are handled through the occupation numbers of each spin block; near-degenerate open-shell pairs produce a strong paramagnetic response, as they should in an uncoupled model. The other column is empty in this lane (no core or Rydberg NAOs in the GFN2 basis). -b nao and -b nvb coincide here. The lane is selected by the input type; --xtb-coupled / --xtb-uncoupled choose the response (the stage line reads coupled or uncoupled accordingly).

ICSS shielding maps (--export-cube icss, icss-edlb, icss-eddb). ICSS (iso-chemical-shielding surfaces, Klod and Kleinpeter 2001) put the NICS probe on every point of a 3D grid: the cube holds the isotropic shielding σ_iso(r) in ppm that a ghost nucleus at r would see, visualised as isosurfaces (e.g. +5, +2, +1, +0.5, −0.5, −1 ppm) or as a colour-coded plane 1 Å above a ring. The sign is the shielding sign of the ICSS literature: positive = shielded (inside the cone of an aromatic ring), negative = deshielded (the in-plane region outside the ring); NICS at that point is −σ_iso. runEDDB adds what no other program has: the same partition as the NICS tables, point by point — icss-edlb is the map of the localized electrons (EDLB_G), icss-eddb that of the delocalized ones (EDDB_G), and at every point σ_iso = EDLB_G + EDDB_G + other (core and Rydberg outside the analysis block; zero in -b nao). Under -a / -s the EDDB_G map is that of the selected NOBDs (e.g. the π system alone). The three maps are computed in one pass, so requesting all of them costs the same as one.

Lanes, cost and grid. The maps come from the same response densities and integrals as the tables, so they exist for every lane: GIAO files, Gaussian CSGT logs (CTOCD-DZ), a wavefunction file without NMR data (the approximate --nmr-kernel lane) and a geometry file (GFN2-xTB). Approximate lanes say so in the stage line and in the first line of the cube; in the xTB lane, which has no core electrons, points within 1 Å of a nucleus are set to zero. The grid is the cube grid of --export-cube (-c, --cube-ref, scope), but with Kleinpeter’s lattice as the default instead of ~100³ points: an exact 0.5 Å step on every axis (an odd number of nodes, so the box centre — the ring centre of a symmetric molecule — is a node; -c N switches back to ~N³ points) and a box margin of 10 Å around the atoms instead of 3.2 Å (--cube-margin R changes it): the values at a grid point are the shielding integrals contracted once per primitive pair of the basis, so the cost grows with the square of the basis size and linearly with the number of points and is shared by all cores (-n). Benzene/6-31G takes seconds, a 24-atom def2-TZVP molecule about a minute on 16 cores at -c 30 and a few minutes at -c 50; the values do not depend on the number of threads. Under -v the stage lists the grid, the lane, the partition, the numbers of shell and primitive pairs after screening, a progress bar and the range of the map.

Tensor component (--icss-axis). The isotropic maps are Kleinpeter’s ICSS; ring currents are read more sharply from the component of the shielding tensor along the ring normal (the NICS_zz analogue, field and probe along n). --icss-axis z (or x, y) adds the maps of σ_zz for a frame axis, --icss-axis ring N takes the normal of NICS ring N (from --nics-rings / --nics-ring) and --icss-axis 0.3,0.4,0.87 any direction; the files are named <stem>.ICSS_zz.cube, <stem>.ICSS_zz_EDLB_G.cube, … (ICSS_nn for a ring normal or a vector) and the isotropic maps are written as well, from the same pass. At a ring probe of the tables the ICSS_nn value equals the NICS_zz row (sign reversed). For benzene the σ_zz cone is about three times deeper than the isotropic one (29.3 vs 11.1 ppm at 1 Å, B3LYP/6-311G(d,p)), its +0.5 ppm surface reaches 8.5 Å and its 0.1 ppm surfaces extend beyond the 10 Å box; the delocalized electrons carry 25.0 of the 29.3 ppm at 1 Å.

Long passes (--icss-checkpoint). A 3D map of a large molecule in a big basis can take hours. With --icss-checkpoint [MIN] the partial sums of the pass are saved every MIN minutes (default 10) to <stem>.ICSS*.tmp; a run started again with the same command finds the file, checks that it belongs to the same pass (grid, lane, components, memory blocks) and continues from the saved point — the stage line says so — and the result is identical to that of an uninterrupted pass, because every grid row is accumulated in a fixed order. The file is removed when the pass completes; a file that does not match the pass is ignored. The 2D maps of --export-map take seconds and have no checkpoint.

2D maps (--export-map). The classic ICSS picture is a colour-coded plane 1 Å above the ring; --export-map icss,icss-edlb,icss-eddb writes exactly that, on a 2D lattice (0.25 Å step by default, --map-step) in the plane given by --map-plane: ring N[,h] (h Å above NICS ring N along its normal), fit[,h] (the best-fit plane of the heavy atoms of the scope, the default with h = 1.0), atoms i,j,k[,h] (through three atoms) or x=V / y=V / z=V (a frame plane, h = 0). The density tokens ed, eddb, edlb, bdf put the corresponding density slice on the same lattice (e.g. the EDDB_G density above the ring next to the EDDB_G shielding it produces). Density slices are averaged over a slab by default. A density read at one height misrepresents heteroatoms: the 2p functions of N and O are more contracted than those of C, so their π density peaks closer to the plane and decays faster. In pyrrole, furan and pyridine (B3LYP/def2-TZVPP) the EDDB_G density above N or O exceeds that above C below 1 Å, matches it at 1.0–1.1 Å and falls below it from 1.25 Å on (furan O at 1.5 Å: 30 % less than C), although the EDDB populations of N and O are the largest in the ring. The mean of five planes between 0.5 and 1.5 Å restores the population ordering and ratios (pyridine N/C 1.25–1.36 vs. 1.25–1.37 from the populations), so this slab is the default of the density tokens; --map-slab h1,h2[,n] sets another slab, --map-slab none or an explicit height in --map-plane gives a single plane, and the ICSS maps stay single-plane (Kleinpeter’s convention) unless --map-slab is given. The file header records the slab. A 2D map costs seconds (10⁴ nodes) and is not restricted to planar molecules: the plane is defined locally and the map is a section through the 3D field. The files are plain text (.map.dat, see §5) for gnuplot or matplotlib; --icss-axis adds the nn maps from the same pass, and the value at the plane origin of a ring map equals the corresponding NICS row of the table with the sign reversed (benzene RHF/6-31G: 11.961 = −NICS(1), 28.972 = −NICS(1)_zz).

Reading the maps. The values near the nuclei are the nuclear shieldings (tens to hundreds of ppm) and dominate any automatic colour scale: set the isosurface levels or the colour range by hand (±0.5 to ±5 ppm is the usual window). At a NICS probe of the tables the cube value is identical to the table value (same integrals, same partition), so the tables remain the reference for numbers and the maps show where in space the shielding of each domain comes from. As with the NICS tables, the map is a magnetic response, not a measure of delocalization: a strong EDDB_G cone above a ring shows that the delocalized electrons carry the ring current there, and the EDLB_G map shows the local (Larmor and σ-bond) shielding that every molecule has.

Validation against the ICSS literature (benzene). Kleinpeter’s protocol places ghost atoms on a lattice from −10 to +10 Å with a 0.5 Å step (68,921 NICS points, GIAO at B3LYP/6-311G(d,p) on B3LYP/6-311G(d,p) geometries) and characterises a ring current by how far a given iso-chemical-shielding surface reaches from the ring centre. The table gives the reach (Å) of the shielding surfaces along the ring normal and of the −0.1 ppm deshielding surface in the ring plane along the C–H direction, read from runEDDB cubes on the same 41³ lattice (--cube-ref with a ±10 Å template grid) by trilinear interpolation of the icss map; the literature values are those quoted by Kleinpeter for benzene as the standard aromatic compound.

^a Kleinpeter reports that the +5 ppm surface of benzene “never exceeds 4.2 Å”, i.e. 2.1 Å above and below the ring; the +0.5 and +0.1 ppm values (5.1 and 8.9 Å) are from his 2024 papers (see §11).

Three checks are worth knowing. (i) At every lattice node the cube value is the NICS table value: the --nics-scan and --nics probes placed on the nodes agree with the map to 10⁻¹² ppm (same integrals, same partition, only the order of summation differs). (ii) The reach of a surface depends on the grid step only through the interpolation between nodes: 1.0, 0.5 and 0.25 Å lattices (21³, 41³, 81³ points) move the values above by at most 0.05 Å, so Kleinpeter’s 0.5 Å step is sufficient; the 81³ benzene map at B3LYP/6-311+G(d,p) took 3.7 minutes on 16 cores. (iii) The 0.1 ppm surfaces lie 7–9 Å from the ring centre, which is why the ICSS maps use a 10 Å box margin by default (the 3.2 Å margin of the density cubes would cut them off): the default grid (--export-cube icss, no -c, no --cube-margin) is exactly this lattice with the molecule’s bounding box plus 10 Å (benzene: 51×51×43 nodes, 70 s on 16 cores with the σ_zz maps included) and reproduces the table above (2.04 / 3.01 / 3.92 / 5.02 / 8.77 Å for the B3LYP/6-311G(d,p) map); --cube-margin R widens or narrows the box, and --cube-ref clones any lattice. The partition adds what the literature maps cannot show: on the ring axis the delocalized electrons carry essentially the whole shielding cone (EDDB_G 20.9 of 8.9 ppm at the centre and 10.5 of 11.1 ppm at 1 Å in the B3LYP/6-311G(d,p) map, the localized part being negative at the centre), while far from the ring the EDLB_G and EDDB_G contributions become comparable (0.74 and 0.15 ppm at 5 Å in the HF map) and the in-plane deshielding belt is a joint effect of both.

References: GIAO — Wolinski, Hinton, Pulay (1990); the NCS-type partition idea — Bohmann, Weinhold, Farrar (1997); NICS — Schleyer et al. (1996); ICSS maps — Klod and Kleinpeter (2001); CTOCD-DZ — Lazzeretti, Malagoli, Zanasi (1994); the ipsocentric selection rules — Steiner and Fowler (2001); CSGT response densities — Keith and Bader (1993); the exchange kernel of --xtb-coupled — Grimme (2013), Grimme and Bannwarth (2016). See §11; --cite adds them automatically to any run that uses --table nmr.

8. Understanding the Output

Analysis tables are printed only when requested with --table — no tables appear by default (the default view is the stage messages plus the EDDB summary block). All requested tables are collected into one captioned > Printing analysis details: section, printed just before the summary in the order mo → no → bo → nobd → ao → subshell → shell → npa: --table mo adds the per-MO ED / EDLB / EDDB table (when canonical MOs are present), --table no the NO table, --table bo the bond-order orbitals (aliases boo/nba), --table nobd the NOBD table, --table ao the per-orbital (atomic-orbital) population table, --table subshell the subshell-resolution population table, --table shell the shell-resolution population table, --table npa the condensed-atomic table (led by the natural charge), and --table all every table. The orbital-based tables (MO / NO / BOO) print before the NOBD-driven dissection; ao / subshell / shell / npa print after it. --verbose adds per-stage detail (and, under it, the per-stage Papers to be cited (DOI): lists) but not tables; -q (true-quiet) suppresses stage messages and timings too — only the banner, the EDDB summary block, and the two bookend stage lines around it are emitted.

NOBD-resolution table: appears with --table nobd (also shown by -s used without a list, so you can pick NOBD indices for a follow-up run). Each row is one NOBD (Natural Orbital for Bond Delocalization):

• Sel — star (*) marks selected NOBDs (all starred without dissection).
• Total — EDDB-density eigenvalue. Values near 2.0 mark strongly delocalized bonding orbitals.
• s, p, d, f + %s, %p, %d, %f — sub-shell contributions and percentages. High %p flags π-type orbitals; mixed s/p flags σ. Under --auto-pi a NOBD is selected as π when its (s+p) content is < 10% s and ≥ 50% of that content is p aligned with the local plane normal at each atom.

Per-MO decomposition (MO): appears with --table mo, when the input carries canonical molecular orbitals (HF / DFT references; rejected for densities with no MOs). One row per MO — its energy and occupation, and its contribution decomposed into ED (total), EDLB (localized) and EDDB (delocalized) — so you can read off which orbitals carry the delocalized density (e.g. the ring-current-active π set). Printed before the dissection, so the numbers reflect the full, undissected density.

Bond-order orbitals (BOO): appears with --table bo (aliases --table boo / --table nba; orbitals exportable with --export-fchk bo). Each row is one atom pair, whose (Wiberg) bond order is resolved into orthonormal bond-order orbitals — the natural orbitals of that single bond. The columns give the Index of the pair’s first BOO (matching the MO numbering in the exported FCHK; the pair’s last BOO is implicit — one below the next row’s Index), the atom-index and element pairs joined by ~, the Total bond order WBO (the full Wiberg index, summed over all BOOs), and the twelve largest BOO occupations (σ²) above the listing threshold, each tagged with a σ / π / δ / φ symmetry label (from the two-centre (n·L)² classifier about the bond axis). Only those are listed, so they need not sum exactly to the Total — the small remaining tail is omitted from the printout and from --export-fchk bo (orbitals with a weight below 10⁻⁴ are not exported; the file header says so). The Total itself is the complete trace of the bond-order operator over all components, so it never depends on the listing or export cutoffs; under -v the complete per-pair balance (every component with its coupling eigenvalue and label, the σ/π/δ/φ/other channel sums and the positive/negative/null branch sums) is written to <stem>.NBA_DIAG.csv and <stem>.NBA_MODES.csv. The resolution is spin- and correlation-consistent — identical for restricted / unrestricted and HF / DFT / post-HF densities — and is evaluated in the basis selected by -b. Only with -b nao does the Total equal the full (Mayer / layer-consistent) bond order of the molecule: in the default NVB (or NMB) block the Rydberg NAOs are excluded and, for a correlated density, the quasi-degenerate layers of the BOO construction are built from that block (unlike the BOP layer partition, which is always derived in full NAO — see Post-HF layer basis above; for the bond-order components the difference is a few thousandths), so the absolute Totals come out slightly smaller (benzene C–C: 1.440 in NAO vs 1.418 in NVB at RHF, 1.394 vs 1.371 at MP2). The σ/π resolution that the table is meant for is practically unaffected (the π share of the C–C bond is 0.316 vs 0.312), so the valence-block default is the sensible choice; reach for -b nao only when the exact textbook value of the total bond order matters. For a correlated density whose natural occupations leave the admissible range (relaxed MP2 densities, see the post-HF note in §2) the NBA stage prints a counted WARNING with the violation counts and the input → effective trace, applies no repair and continues; the CSV files carry the same diagnostics in their header.

Large systems — scope the bond-order analysis. Both --table bo and --export-fchk bo honour the analysis scope: they resolve BOOs only for the atom pairs in the active --fragment / --pathway / --ring-detect selection, using the same connectivity that drives the EDDB. On a large molecule, prefer one of these scopes over --global: a global run diagonalises a bond-order operator for every connected pair in the structure, which is rarely what you want when you are studying one bond, ring or fragment. Scoping also keeps the exported <stem>.BO.fchk small — it contains only the selected pairs’ BOOs (each still expanded in the full AO basis, so the file stays a valid “partial MO set” that most FCHK viewers accept), making it ideal for inspecting the bonding pattern of a specific fragment.

Orbital-resolution table (AO): appears with --table ao. One row per NAO in the active -b subset. The columns are: Index (the orbital number — it matches the MO index in the --export-fchk ao FCHK 1:1, so the table doubles as a navigation key in Avogadro2 / IQmol / GaussView), the NAO label (its angular type, e.g. 2s, 2p — the header carries the basis name NAO / NMB / NVB) and its Atom, then the orbital’s total electron population ED and its delocalized part EDDB (both tagged with the scope suffix _G / _F / _P / _R). Each of ED and EDDB has α and β sub-columns for open-shell inputs; for a closed-shell (restricted) reference these are not defined and print a centred dash (–).

Subshell-resolution table (SUBSHELL): appears with --table subshell — the intermediate resolution between ao and npa, new in Aug 2026 (until mid-Aug 2026 this table answered to --table shell, which now prints the coarser shell-resolution table below). One row per subshell (atom, n, L): the ED / EDDB populations (and their α / β sub-columns for open-shell inputs) are the sums over the 2L+1 magnetic components of the corresponding AO-table rows — e.g. the 2p row is 2p1 + 2p2 + 2p3. The summation is keyed by the NAO quantum labels, never by table position, so it is correct for any AO display ordering (the .49 lane and the different --nao-algorithm constructors order the AO rows differently). Rows print in canonical order — atom by atom, ascending n, then L (1s, 2s, 2p, 3s, 3p, 3d, …); there is no Index column (subshell rows do not map to AO / FCHK orbital indices — use the ao table for that). Works for every -b subset; with --nao-algorithm cc-layers the rows directly display the frozen shell-layer hierarchy.

Shell-resolution table (SHELL): appears with --table shell — the coarsest orbital-resolved view, new in Aug 2026. One row per principal shell (atom, n): the ED / EDDB populations (and α / β sub-columns for open-shell inputs) are the sums over all L of the subshell rows — e.g. shell 2 is 2s + 2p. Shells are labelled by the plain principal number 1, 2, 3, … (no spectroscopic K / L / M letters: those end at n = 8, while decontracted bases reach far higher n, and custom sets may lack the highest-L members). Because the NAO Rydberg populations fall roughly an order of magnitude per shell, this table exposes the Rydberg population ladder at a glance — a fast diagnostic of the shell truncation point when selecting compact, shell-structured virtual subspaces (pairs naturally with --nao-algorithm cc-layers, where each row is exactly one frozen layer). Same accumulation and ordering rules as the subshell table.

NO-resolution table: appears with --table no. One row per natural orbital with non-trivial occupation, sorted by descending NOON. Each row reports the NOON itself plus the breakdown of its occupation into Core (Cor), Valence (Val), and Rydberg (Ryd) NAO contributions, with percentages. Row index matches the MO index in the --export-fchk no FCHK 1:1.

Multireference / spin-contamination indices. These are reported in the NO stage’s verbose (-v) detail block — not inside the table — right after the occupation-spectrum summary, each with its source DOI in brackets (run --table no -v to see them):

• Yamaguchi y₀ biradical index — 0 for a closed-shell single reference, 1 for an ideal singlet diradical, intermediate for partial diradical character.
• Bofill-Pulay index — the per-electron idempotency deviation N_U^Mayer/N_e with N_U^Mayer = Σ_i n_i(2−n_i) (zero for a pure single Slater determinant; rises with correlation).
• Head-Gordon index — the smoothed effectively-unpaired-electron count Σ_i n_i²(2−n_i)² (less inflated by the small-correlation tail of large basis sets — the most-cited NOON-only multireference indicator).
• Spin contamination — the UHF ⟨S²⟩ − S_z(S_z+1); printed for open-shell inputs only (it is identically zero, and omitted, for closed-shell references).

Condensed-atomic table (NPA): appears with --table npa (the orbital-resolution table is now the separate --table ao token; when both are requested, AO prints just before NPA). Each row leads with the natural atomic charge q_N = Z_eff(A) − N_A, where N_A is the full natural electron population on the atom (core + valence + Rydberg) and Z_eff is the effective nuclear charge — so the charges are correct, and sum to the molecular charge, even with ECP-replaced cores or the valence-only GFN2-xTB core (computed from the undissected density, so independent of --auto-pi / --select-nobd). The remaining columns give the per-atom total electron population (ED) and electron delocalization (EDDB, the EDDB density integrated over the atom), each tagged with the scope suffix _G / _F / _P / _R. ED and EDDB each carry α and β sub-columns for open-shell inputs; for a closed-shell (restricted) reference these are undefined and print a centred dash (–).

Summary: three population numbers — ED (total electrons in the active subset), EDLB (localized; ED − EDDB), and EDDB (delocalized). For UHF / ROHF the α and β components are reported separately in the orbital and atomic tables.

Cyclic Delocalization Index (CDI): reported only for --pathway runs. CDI equals the smallest per-atom EDDB along the pathway — the maximum population of electrons that can delocalize uniformly around the cycle. CDI ≈ 0.9 per atom is strongly aromatic; CDI ≈ 0 is non- or anti-aromatic.

Closing line and the warning counter: the run terminates with a one-line time-stamp such as > Calculation completed on 2026-MM-DD. Total time: Xs. Whenever any non-fatal warning fired during the run, the closing line is annotated with a parenthetical count: > Calculation completed (with 1 warning!) on 2026-MM-DD. Total time: Xs. The annotation appears even in true-quiet (-q) mode where the warning text itself is suppressed — so a batch user always knows when something needs a closer look. To see the actual warning message, rerun the same command without -q. ECP-related inconsistencies (missing [Pseudo] in molden, Q-Chem GUI=2 quirk in fchk) are fatal errors — they stop the run rather than just counting; see §2 for the recovery steps.

9. Known Issues and Limitations

Tested QC programs and formats: validated against .fchk files from Gaussian 16 and Q-Chem 6 (HF / DFT and post-HF), .json and .molden files from ORCA 6 (HF / DFT and post-HF via orca_2json and orca_2mkl), .molden files from Multiwfn-converted wavefunctions, native .tmol bundles from Turbomole, and .gms output files from GAMESS-US (RHF / UHF / ROHF / MCSCF on RUNTYP = ENERGY / OPTIMIZE / SADPOINT). Cross-program agreement on benzene at cc-pVDZ is 4–5 significant figures for both MP2 and CCSD; the WF₆ ECP benchmark agrees to ~3 significant figures across Gaussian / ORCA fchk / json / molden at RHF / def2-TZVPP. The built-in GFN2-xTB engine for the .xyz lane is validated against the reference tblite implementation across 81 systems covering 26 elements (H…W) including first / second / third-row transition metals, lanthanides, and main-group hypervalent compounds — final per-atom Mulliken charges agree to ~10⁻⁶e at tight SCC convergence, and total energies match Grimme’s reference xtb binary to within ~10⁻⁴ E_h.

Supported basis functions. The native lane (.fchk / .json / .molden) supports angular momenta up to k-type (l = 7), covering every standard basis set in routine use including cc-pV6Z. Overlaps are evaluated by direct Gauss-Hermite quadrature (no Obara-Saika recursion), so each shell is treated explicitly. The legacy .49 visualisation lane supports up to h-type (l = 5) only; for cc-pV6Z and similar (i / k functions) please use the matching native input format.

Post-HF density requirements. The QC program must be told to write the relaxed (or linearized) correlated density:

• Gaussian: density=mp2 (MP2) or density=current (CCSD, CI, CASSCF) on the route line.
• ORCA: %mp2 density relaxed end for MP2 or %mdci density linearized end for the CCSD family.
• Q-Chem: cc_ref_prop = true in $rem (CC family only; for MP2 use Gaussian or ORCA, since Q-Chem does not emit the relaxed MP2 density to FCHK).

If the correlated block is missing, runEDDB silently falls back to the SCF reference density and the run is (correctly) classified as plain HF/DFT. The Reference type: line in the verbose-mode Stage 1 detail block reports the actual classification.

HLRL — NAO pictures at large basis sets, and a viewer warning. From cc-pVQZ onwards (worse at cc-pV5Z/6Z, def2-QZVPP) the rigorous NAO algorithm stores each NMB orbital as a large, near-cancelling coefficient expansion spread over high-L shells and neighbouring atoms — an inherent property of the NAO method, not a runEDDB artefact (NBO7 produces the same representation). The orbital function remains atom-localised; the fragile representation is the problem. Most molecular viewers evaluate high-L (f/g/h) pure harmonics approximately, which breaks the cancellation and displays the exported raw NAOs as heavily delocalised junk — verified (Jul 2026) in current Avogadro2 and IQMol, while Multiwfn (and exact grid evaluation) render the same files correctly. runEDDB’s fchk conventions are validated against Gaussian-written references. Mitigations, in order of preference: the default frozen-nmb algorithm has no L>3 content by construction — it renders correctly in every viewer and shows the physically polarised orbitals; orbital cubes (--export-cube ao[LIST]) are exact grid evaluations and immune to viewer basis handling for ANY algorithm (and are the only option beyond h functions, where no mainstream viewer reads the fchk at all); --nao-algorithm rigorous/nbo7legacy writes the raw standard coefficients for reproducing the viewer artefact itself. (The former Galerkin display projection of exported standard NAOs is opt-in via --export-lproj — its free-atom-like pictures over-correct by stripping the genuine polarisation tails, so it is no longer applied silently.) See the NAO construction algorithm paragraph in Section 4.

NBO archive (.49) — legacy lane. Population analyses (NOBD / orbital / atomic tables, summary) and all six visualisation outputs work via a parallel NBO-convention pipeline that mirrors the legacy runEDDB_nbo49 writers. The .49 path is provided for legacy compatibility only and will be removed once the native NAO driver completes its acceptance suite; for new projects, please use .fchk / .json / .molden.

Cube grid quality. The default grid (~100³ points) gives integration accuracy within ~1% for most systems. For publication images, use -c 150 or higher. The integrated density and its percentage of the analytical value are printed after each cube export for quality control.

10. Performance Tips

Parallel BOP analysis: use --ncores N to parallelize the BOP step. Speedup is nearly linear for large systems (100+ atoms).

$ runEDDB -i large_molecule.fchk --ncores 4

Using every core (-n 0): --ncores 0 (equivalently -n auto) asks runEDDB to detect how many cores it may use and to run on all of them, so the one command adapts itself to whatever machine it lands on — laptop, workstation or node. It behaves identically on Linux, Windows and macOS, and the detected count is echoed right below the banner (Auto-detected N available CPU cores), so every log records the thread count that was actually used. On Linux the number comes from the process’ CPU affinity mask rather than from the machine as a whole, so a batch allocation is honoured: a job given 8 cores of a 128-core node through SLURM / taskset resolves -n 0 to 8, not to 128.

NOTE — more cores is not always faster. -n 0 is deliberately literal: it takes every core it is allowed to use, which is what you want on a desktop, a laptop or a single-socket server, but is rarely what you want on a large many-core machine. The dense linear algebra that dominates the NAO/BOP stages stops scaling after a few tens of threads, and past that point the extra threads spend their time synchronising instead of computing — a cost that grows sharply on dual- or multi-socket systems, where threads on one socket must reach across the interconnect for memory attached to the other. The effect is not subtle: on a dual-socket 128-core EPYC, a large-basis (aug-cc-pV5Z) benchmark took 2.2 s with -n 8 and 41 s with -n 0 — roughly 20× slower while burning 75 CPU-minutes instead of 16 CPU-seconds. As a rule of thumb: use -n 0 on laptops, desktops and single-socket servers; on big multi-socket nodes set an explicit -n 8–32 (and, if the node is shared, confine the job with numactl / taskset so that -n 0 follows the affinity mask — but see the NUMA note below before confining it to a single socket for speed, which is counter-productive). Cube-grid evaluation is the exception — it parallelises almost perfectly, so grid-dominated jobs do benefit from high core counts. When in doubt, time a short run at a few settings rather than assuming.

Large NUMA machines — interleave memory only when the job spans several nodes. Linux places each memory page on the NUMA node of whichever thread first touches it, which for runEDDB is the input-parsing stage. That default is ideal as long as the worker threads stay on that same node, and harmful once they do not: the large matrices stay concentrated on one node while every thread elsewhere reaches for them across the interconnect. Which regime you are in depends on how -n compares with the cores per NUMA node (numactl --hardware).

Measured on a dual-socket AMD EPYC 9554 (128 cores, 8 NUMA nodes of 16 cores, SMT off):

• At -n 128, spanning all eight nodes, a large system took 15.9 s by default and 9.4 s with --interleave=all. Without it, more than 64 cores was slower than 32 — easy to mistake for a limit of the program.
• At -n 16, which fits inside a single 16-core node, interleaving was consistently about 4% slower (9.12 s vs 9.54 s). Here the default placement is already perfect, and interleaving only makes seven of every eight accesses remote.

So: leave the default alone when -n is at or below the cores per node, and add --interleave=all when you deliberately spread a big job across nodes. Note also that binding a large job to a single socket for speed is counter-productive — the workload is limited by memory bandwidth, and one socket has half the memory controllers.

GFN2-xTB on multiple cores: for the .xyz lane, --ncores N also parallelizes the GFN2-xTB SCC itself (the heaviest step for large molecules). The engine auto-sizes its thread pools to the cores you request, so requesting more cores than the job can use no longer wastes CPU. As a rule of thumb, the GFN2-xTB SCC keeps speeding up to roughly --ncores 16–32 on a few-hundred-atom system; beyond that the gains taper off (it is dominated by linear-algebra steps that do not scale indefinitely). For small molecules, --ncores 1–4 is plenty.

Analysis basis subset: the default -b nvb (Val-only) is the fastest option, since the BOP cycle scales with the number of NAOs in the working basis. -b nmb (Cor + Val) costs only marginally more (~10–30% depending on element). -b nao (Cor + Val + Ryd) is significantly slower at large basis sets — use it for completeness checks rather than routine analysis.

Cube grid sizing: the default --cube-size 100 (~1M points) is fine for routine work; for publication-quality figures and large systems we recommend at least 150. Grid evaluation can be demanding but parallelizes perfectly, so use as many cores as you can (--ncores / -n).

11. Citations

If you use runEDDB in your research, please consider to cite the following works:

Tip — let runEDDB pick the references for you. Add --cite both to any run to write the exact, de-duplicated set of papers for that run to <stem>.bib (BibTeX) and <stem>.refs.txt (plain list). The set is method- and option-dependent: the NAO / NPA (Reed et al.), BOP and core EDDB papers always; the NOBD paper with -s / -a; the EDDB_P (-p / -r) or EDDB_F (-f) paper; the AO-resolved EDDB paper with --table / --export-fchk ao; the Bond-Delocalization-Field paper for a bdf / dbdf cube; the frozen-NMB NAO paper (in preparation) whenever the default --nao-algorithm frozen-nmb is active; and, for the .xyz lane, GFN2-xTB + DFT-D4 (spGFN2-xTB on spin-polarised runs). Add -v to also see them on screen, grouped by analysis stage. The runEDDB-method catalogue follows (the NAO / NPA reference is Reed, Weinstock, Weinhold, J. Chem. Phys. 83, 735−746, 1985, DOI 10.1063/1.449486; GFN2-xTB / DFT-D4 are listed under it).

1. Electron delocalization index based on bond order orbitals.
   D.W. Szczepanik (), E.J. Zak, K. Dyduch, J. Mrozek
   Chemical Physics Letters  593 (2014) 154−159. DOI: 10.1016/j.cplett.2014.01.006.   URL 
   
   
2. A uniform approach to the description of multicenter bonding.
   D.W. Szczepanik (), M. Andrzejak, K. Dyduch, E.J. Zak, M. Makowski, G. Mazur, J. Mrozek,
   Physical Chemistry Chemical Physics  16 (2014) 20514−20523. DOI: 10.1039/c4cp02932a.   URL 
   
   
3. A new perspective on quantifying electron localization and delocalization in molecular systems.
   D.W. Szczepanik ()
   Computational and Theoretical Chemistry  1080 (2016) 33−37. DOI: 10.1016/j.comptc.2016.02.003.   URL 
   
   
4. On the three-center orbital projection formalism within the electron density of delocalized bonds method.
   D.W. Szczepanik ()
   Computational and Theoretical Chemistry  1100 (2017), 13−17. DOI: 10.1016/j.comptc.2016.12.003.   URL 
   
   
5. The electron density of delocalized bonds (EDDB) applied for quantifying aromaticity.
   D.W. Szczepanik (), M. Andrzejak, J. Dominikowska, B. Pawełek, T.M. Krygowski, H. Szatylowicz, M. Solà
   Physical Chemistry Chemical Physics  19 (2017) 28970−28981. DOI: 10.1039/c7cp06114e.   URL 
   
   
6. Electron delocalization in planar metallacycles: Hückel or Möbius aromatic?
   D.W. Szczepanik (), M. Solà ()
   ChemistryOpen  8 (2019) 219−227. DOI: 10.1002/open.201900014.   URL 
   
   
7. The electron density of delocalized bonds (EDDB) as a measure of local and global aromaticity
   D.W. Szczepanik (), M. Solà ()
   in
   Aromaticity: Modern Computational Methods and Applications
   Edited by I. Fernández
   Chapter 8 (pp. 259−283), Elsevier, 2021.
   ISBN: 978-0-12-822723-7
   DOI: 10.1016/b978-0-12-822723-7.00008-x.
   
   
8. Decoding chemical resonance in π-conjugated systems.
   D.W. Szczepanik (), P.A. Wieczorkiewicz
   ACS Physical Chemistry Au  (2026). DOI: 10.1021/acsphyschemau.6c00031.   URL 
   — the Bond-Delocalization-Field (BDF) paper, cited whenever a bdf / dbdf cube is exported.
   
   
9. Natural atomic orbitals at the complete-basis-set limit: a frozen-minimal-basis two-cycle construction.
   D.W. Szczepanik ()
   Journal of Computational Chemistry  (2026), in preparation.
   — the frozen-NMB NAO construction, cited whenever the default --nao-algorithm frozen-nmb is active.
   
   

When you analyse a plain .xyz geometry through the built-in GFN2-xTB lane, the underlying electronic-structure method is not ours — it is the work of Grimme and co-workers. Please also cite the relevant original method papers:

1. GFN2-xTB — An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions.
   C. Bannwarth, S. Ehlert, S. Grimme
   Journal of Chemical Theory and Computation  15 (2019) 1652−1671. DOI: 10.1021/acs.jctc.8b01176.   URL 
   — the GFN2-xTB Hamiltonian used for every .xyz calculation (minimal valence vSTO-nG basis, multipole / anisotropic electrostatics, density-dependent dispersion).
   
   
2. A generally applicable atomic-charge dependent London dispersion correction.
   E. Caldeweyher, S. Ehlert, A. Hansen, H. Neugebauer, S. Spicher, C. Bannwarth, S. Grimme
   The Journal of Chemical Physics  150 (2019) 154122. DOI: 10.1063/1.5090222.   URL 
   — the self-consistent DFT-D4 dispersion correction (charge-dependent C₆ coefficients + ATM three-body term) included in the GFN2-xTB energy.
   
   
3. High-throughput screening of spin states for transition metal complexes with spin-polarized extended tight-binding methods.
   H. Neugebauer, B. Bädorf, S. Ehlert, A. Hansen, S. Grimme
   Journal of Computational Chemistry  44 (2023) 2120−2129. DOI: 10.1002/jcc.27185.   URL 
   — the spin-polarized extension (spGFN2-xTB) used for open-shell systems (multiplicity > 1).
   
   

When you use --table nmr / --nics*, please also cite the methods behind the magnetic response and its partition:

1. Efficient implementation of the gauge-independent atomic orbital method for NMR chemical shift calculations.
   K. Wolinski, J.F. Hinton, P. Pulay
   Journal of the American Chemical Society  112 (1990) 8251−8260. DOI: 10.1021/ja00179a005.  [link]
2. Natural chemical shielding analysis of nuclear magnetic resonance shielding tensors from gauge-including atomic orbital calculations.
   J.A. Bohmann, F. Weinhold, T.C. Farrar
   The Journal of Chemical Physics  107 (1997) 1173−1184. DOI: 10.1063/1.474464.  [link]
3. Nucleus-independent chemical shifts: a simple and efficient aromaticity probe.
   P.v.R. Schleyer, C. Maerker, A. Dransfeld, H. Jiao, N.J.R. van Eikema Hommes
   Journal of the American Chemical Society  118 (1996) 6317−6318. DOI: 10.1021/ja960582d.  [link]
4. Ab initio calculation of the anisotropy effect of multiple bonds and the ring current effect of arenes — application in conformational and configurational analysis.
   S. Klod, E. Kleinpeter
   Journal of the Chemical Society, Perkin Transactions 2  (2001) 1893−1898. DOI: 10.1039/b009809o.  [link]
5. Identification of local benzenoid aromaticity and global aromaticity of polycyclic aromatic hydrocarbons (PAHs) via the corresponding ring current effects in ¹H NMR spectroscopy.
   E. Kleinpeter, A. Koch
   Arkivoc  2024 (7) 202412263. DOI: 10.24820/ark.5550190.p012.263 (open access; Kleinpeter’s ICSS protocol and the benzene +5 ppm reach).  [link]
6. ¹H and ¹³C NMR spectra of infinitene and the ring current effect of the aromatic molecule.
   E. Kleinpeter
   Magnetic Resonance in Chemistry  2024. DOI: 10.1002/mrc.5467 (benzene ICSS(0.5 ppm) = 5.1 Å, ICSS(0.1 ppm) = 8.9 Å).  [link]
7. The chemical structure of a molecule resolved by atomic force microscopy.
   L. Gross, F. Mohn, N. Moll, P. Liljeroth, G. Meyer
   Science  325 (2009) 1110−1114. DOI: 10.1126/science.1176210.  [link]
8. A simple model of molecular imaging with noncontact atomic force microscopy.
   N. Moll, L. Gross, F. Mohn, A. Curioni, G. Meyer
   New Journal of Physics  14 (2012) 083023. DOI: 10.1088/1367-2630/14/8/083023 (the Pauli-repulsion power law behind --export-map afm).  [link]
9. Theory of the scanning tunneling microscope.
   J. Tersoff, D. R. Hamann
   Physical Review B  31 (1985) 805–813. DOI: 10.1103/PhysRevB.31.805 (the constant-height STM image of --export-map stm as the local density of the states in the bias window)  [link]
10. Computational approach to molecular magnetic properties by continuous transformation of the origin of the current density.
    P. Lazzeretti, M. Malagoli, R. Zanasi
    Chemical Physics Letters  220 (1994) 299−304. DOI: 10.1016/0009-2614(94)00158-8.  [link]  (CTOCD-DZ lane)
11. Calculation of magnetic response properties using a continuous set of gauge transformations.
    T.A. Keith, R.F.W. Bader
    Chemical Physics Letters  210 (1993) 223−231. DOI: 10.1016/0009-2614(93)89127-4.  [link]  (CSGT response densities, CTOCD-DZ lane)
12. Patterns of ring currents in conjugated molecules: a few-electron model based on orbital contributions.
    E. Steiner, P.W. Fowler
    The Journal of Physical Chemistry A  105 (2001) 9553−9562. DOI: 10.1021/jp011955m.  [link]  (frontier-orbital NMR tables)
13. A simplified Tamm–Dancoff density functional approach for the electronic excitation spectra of very large molecules.
    S. Grimme
    The Journal of Chemical Physics  138 (2013) 244104. DOI: 10.1063/1.4811331.  [link]  (exchange kernel of --xtb-coupled)
14. Ultra-fast computation of electronic spectra for large systems by tight-binding based simplified Tamm–Dancoff approximation (sTDA-xTB).
    S. Grimme, C. Bannwarth
    The Journal of Chemical Physics  145 (2016) 054103. DOI: 10.1063/1.4959605.  [link]  (kernel parameters of --xtb-coupled)
15. Patterns of ring currents in conjugated molecules: a few-electron model based on orbital contributions.
    E. Steiner, P.W. Fowler
    The Journal of Physical Chemistry A  105 (2001) 9553−9562. DOI: 10.1021/jp011955m.  [link]  (ipsocentric selection rules)

runEDDB’s GFN2-xTB lane implements the GFN2-xTB and DFT-D4 methods described above and is validated against Grimme’s reference implementation. Portions of it are derived from the corresponding open-source reference code — dftd4 and tblite (LGPL-3.0-or-later) and mctc-lib (Apache-2.0) — and are used and redistributed under those licenses. See the THIRD-PARTY-LICENSES.txt file included with the program for full attributions and license texts.

12. Contact

Author and main contributor: dr hab. Dariusz W. Szczepanik, prof. UJ

Contributors: mgr Pawel A. Wieczorkiewicz, dr Ouissam El Bakouri.

Project e-mail: runeddb@gmail.com