Skip to content

F-SAPT

Functional-group symmetry-adapted perturbation theory (F-SAPT) decomposes the non-covalent interaction energy between two molecules ("monomers") into physically meaningful components — electrostatics, exchange (Pauli repulsion), induction, and dispersion — and then further attributes each component to pairs of functional groups, one on each monomer. This makes it possible to see exactly which chemical groups drive (or oppose) binding, which is especially valuable for ligand–protein and host–guest interactions.

TeraChem implements the GPU-oriented F-SAPT variant of Parrish, Thompson and Martínez,1 in which the analysis is written almost entirely in terms of Coulomb (\(J\)) and exchange (\(K\)) matrix builds, and two pragmatic choices make large systems tractable:

  • Ab-initio dispersion and exchange-dispersion are replaced by Grimme's empirical DFT-D3(BJ) dispersion correction.3
  • Augmented basis sets are avoided so that integral screening remains effective.

The method is SAPT0 (i.e. SAPT based on a Hartree–Fock description of the monomers).

Quick start

F-SAPT is requested with fsapt yes added to an ordinary closed-shell RHF input. The dimer is split into the two monomers by atom order: the first fsapt_na atoms are monomer A and the remainder are monomer B. The geometry must list all of monomer A's atoms first.

water_dimer.in
basis         cc-pvdz
method        hf
coordinates   coord.xyz   # 6 atoms: water A (1-3), then water B (4-6)
charge        0
spinmult      1
run           energy
gpus          1

fsapt         yes
fsapt_na      3           # first 3 atoms = monomer A; rest = monomer B
fsapt_chargea 0
fsapt_chargeb 0
end

This runs three RHF calculations — the dimer and each monomer in the dimer-centered basis set (the other monomer's nuclei ghosted) — then the SAPT0 terms, the IAO/IBO localization and functional-group partition, and finally the F-SAPT tables.

Order matters for performance: put the smaller monomer first

The result is independent of which monomer you call A — but the cost is not. The order-1/order-2 partition does per-fragment Coulomb and ESP builds over monomer A's fragments, so the analysis is cheapest when A is the smaller monomer (this is the algorithm's intended "A@B" regime1, ideal for e.g. a small ligand A against a large protein B). Concretely, for a 63-atom ligand interacting with a 3-atom water, listing the water first (A) ran ~18% faster than listing the ligand first, with identical energies; the gap grows with the size disparity between the monomers.

With no fragment files (as above) each monomer is treated as a single fragment, so the analysis reduces to ordinary SAPT0 with a functional-group table that has one entry per monomer. To resolve groups within a monomer, supply fragment files (next section).

Functional-group fragmentation: fA.dat / fB.dat

The "F" in F-SAPT is a two-level partition. The first level splits the dimer into monomers A and B (fsapt_na). The second level subdivides each monomer into functional groups, specified by two plain-text files referenced with fsapt_fraga and fsapt_fragb:

ethane_water.in (fragment files)
fsapt         yes
fsapt_na      8
fsapt_fraga   fA.dat
fsapt_fragb   fB.dat
end

Each file has one fragment per line: a label followed by the atom indices belonging to that fragment. Indices are 1-based global atom numbers (the same numbering used in the input geometry, so the same files can drive both TeraChem and the reference implementation). Whitespace-separated tokens may be single indices, inclusive ranges N-M, or any mix of the two; # comments and blank lines are ignored. The example below — for ethane (atoms 1–8) split into two methyl groups interacting with a water (atoms 9–11) — deliberately uses all three forms: a list of single indices (Me1), a single range (Me2), and a mix of the two (Wat):

fA.dat
Me1   1 2 3 4
Me2   5-8
fB.dat
Wat   9 10-11

These are just shorthand for the underlying atom lists, so Me1 1-4, Me2 5-8, Wat 9-11 would specify exactly the same fragments and give identical results.

fA.dat must list only monomer-A atoms and fB.dat only monomer-B atoms, and together they must cover every atom of their monomer exactly once (TeraChem checks this and aborts on a duplicate, missing, or out-of-range atom).

When a fragment boundary cuts a covalent \(\sigma\) bond inside a monomer (e.g. the C–C bond between Me1 and Me2 above), the shared localized bond orbital is detected automatically and split equally between the two fragments using the "50/50 link bond" rule,2 with a compensating \(\pm\tfrac12\) reassignment of nuclear charge across the bond so each fragment stays neutral. No special markup is needed in the fragment files.

Theory (what each term is)

After the three RHF solutions, the SAPT0 terms are evaluated as generalized Fock builds. Writing \(D_X\) for the half-density of monomer \(X\), \(V_X\) for its nuclear-attraction potential, and \(J[D]/K[D]\) for the Coulomb/exchange matrices:

\[ E^{(10)}_{\mathrm{elst}} = E^{\mathrm{nuc}}_{AB} + 2\langle D_A|V_B\rangle + 2\langle V_A|D_B\rangle + 4\langle D_A|J[D_B]\rangle , \]

the exchange repulsion \(E^{(10)}_{\mathrm{exch}}\) (computed both to second order in overlap, \(S^2\), and to infinite order, \(S^\infty\)), the induction \(E^{(20)}_{\mathrm{ind}}\) and exchange-induction \(E^{(20)}_{\mathrm{exch\text{-}ind}}\) (both the uncoupled estimate and the coupled/response value from the monomer CPHF equations), and the DFT-D3(BJ) dispersion

\[ E_{\mathrm{disp}} = -\sum_{A<B}\sum_{n=6,8} s_n\, \frac{C^{AB}_n}{r_{AB}^{\,n} + (a_1 R^{AB}_0 + a_2)^{n}} . \]

The occupied orbitals of each monomer are localized into intrinsic bond orbitals (IBOs) (see Localization), with core and valence orbitals localized separately. Each SAPT term is then partitioned over functional-group pairs \((\mathcal A,\mathcal B)\) using per-fragment weights \(w\) on the nuclei and localized orbitals; for electrostatics

\[ E^{(\mathcal A\mathcal B)}_{\mathrm{elst}} = w_A w_B (A|B) + 2 w_{\bar a} w_B (\bar a\bar a|B) + 2 w_A w_{\bar b} (A|\bar b\bar b) + 4 w_{\bar a} w_{\bar b} (\bar a\bar a|\bar b\bar b) . \]

By construction the fragment-pair matrix of every term sums back to the corresponding total. The exchange matrices are scaled from \(S^2\) to \(S^\infty\) by \(S_{\mathrm{exch}} = E^{(10)}_{\mathrm{exch}}(S^\infty)/E^{(10)}_{\mathrm{exch}}(S^2)\), and the (uncoupled) induction matrices are scaled to the full SAPT0 induction, accounting for coupling and the \(\delta\)HF correction \(E_{\mathrm{ind}} = E^{\mathrm{HF}}_{\mathrm{int}} - E^{(10)}_{\mathrm{elst}} - E^{(10)}_{\mathrm{exch}}\), where \(E^{\mathrm{HF}}_{\mathrm{int}}\) is the supermolecular Hartree–Fock interaction energy.1

Projection (MinAO) basis

The IAO construction projects the occupied orbitals onto a minimal "MinAO" basis. F-SAPT uses cc-pvdz-minao, shipped with TeraChem at $TeraChem/basis/cc-pvdz-minao (independent of the main basis you choose). It is a normal TeraChem basis file, so it may also be selected as a basis set in its own right.

Output

The decomposition is printed as labeled tables — an order-1 table per monomer (each functional group's total Elst/Exch/IndAB/IndBA/Disp and its sum) and an order-2 fragment-pair matrix for each term — to both standard output and a file fsapt.dat in the run directory, all in kcal mol⁻¹. The plain SAPT0 totals and a set of internal consistency markers (each FSAPT <name> = ...) are also printed; the fragment matrices sum to these totals.

For example, the order-1 monomer-A table for an ethane(2 methyls)···water calculation:

Order-1 F-SAPT, Monomer A fragments:
Frag            Elst         Exch        IndAB        IndBA         Disp        Total
Me1           1.6112       5.3404      -0.5272      -0.0314      -1.8029       4.5901
Me2          -2.9198       0.2541      -0.1050      -0.0010      -0.0866      -2.8583

Notes and limitations

  • Closed-shell RHF monomers only; the dimer and both monomers must be closed shell.
  • Monomer A is the first fsapt_na atoms — order your geometry accordingly. Set fsapt_chargea/fsapt_chargeb consistently with each monomer's electron count.
  • Dispersion is empirical (DFT-D3(BJ)); there is no ab-initio Disp20/Exch-Disp20 term, by design.1
  • The method targets large, asymmetric (e.g. ligand–protein) complexes; the cost is dominated by the dimer RHF and the \(J/K\) builds.

Summary of relevant keywords

Keyword Type Default Description
fsapt bool no Run the F-SAPT analysis
fsapt_na int — Number of leading atoms forming monomer A (the rest are monomer B)
fsapt_chargea int 0 Total charge of monomer A
fsapt_chargeb int 0 Total charge of monomer B
fsapt_fraga string none Path to the monomer-A fragment file (fA.dat); if unset, monomer A is one fragment
fsapt_fragb string none Path to the monomer-B fragment file (fB.dat); if unset, monomer B is one fragment

  1. R. M. Parrish, K. C. Thompson, T. J. Martínez, Large-Scale Functional Group Symmetry-Adapted Perturbation Theory on Graphical Processing Units, J. Chem. Theory Comput. 14, 1737 (2018). ↩↩↩↩

  2. R. M. Parrish, T. M. Parker, C. D. Sherrill, Chemical Assignment of Symmetry-Adapted Perturbation Theory Interaction Energy Components: The Functional-Group SAPT Partition, J. Chem. Theory Comput. 10, 4417 (2014). ↩

  3. S. Grimme, J. Antony, S. Ehrlich, H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys. 132, 154104 (2010). ↩