Skip to content

Integration Grids

The exchange–correlation energy and potential cannot be evaluated analytically, so they are computed by numerical quadrature,

\[ E_{\mathrm{xc}} = \int f\big(\rho(\mathbf{r}), \nabla\rho(\mathbf{r}), \dots\big)\, d\mathbf{r} \;\approx\; \sum_{g} w_g\, f\big(\rho(\mathbf{r}_g), \nabla\rho(\mathbf{r}_g), \dots\big), \]

on an atom-centered molecular grid. TeraChem builds this grid from three standard ingredients:

  1. a Murray–Handy–Laming (MHL) radial transformation1 of a uniform quadrature in a mapped variable, scaled by a tabulated atomic radius;
  2. pruned Lebedev–Laikov angular shells2, with the angular order reduced in the core and in the far tail;
  3. Becke molecular partitioning3 of the atomic quadratures into a molecular one, with the atomic-size adjustment of Becke's Appendix A.

The grid is selected with a single keyword, dftgrid 0–5. These presets are TeraChem-specific — they are not the SG-n grids4 and they are not the "(Nrad, Nang)" grids of other codes, so grid levels do not transfer between programs. This page documents the construction exactly as it is implemented, so that a calculation can be reproduced or the parent grid described in a publication.

Note

Everything below concerns the molecular XC grid. Two other grids in TeraChem are documented elsewhere: the Lebedev surface grids used to discretize PCM cavities (see PCM Solvation) and the uniform real-space grid available for periodic DFT (see Periodic Calculations).

Choosing a grid

dftgrid Typical use Nominal points/atom (H–Ne)
0 Coarse; SCF pre-convergence, dynamic-grid first stage 1,328 / 1,368
1 Default. Energies and gradients for routine work 3,720 / 3,892
2 Recommended when tight gradients matter (optimizations, frequencies, meta-GGAs) 8,018 / 8,482
3 Tight 15,218 / 15,746
4 Very tight 39,752 / 40,748
5 Reference-quality; benchmarking grid errors 101,692 / 103,492

(H, He / Li–Ne; see Nominal grid sizes for the remaining rows.)

The number of points actually used is smaller, because shells beyond the range of the basis are dropped and negligible-weight points are discarded. How much smaller depends on the basis set and the element, so the only reliable count is the one the program reports: at printlevel 4 or higher TeraChem prints the grid level and the surviving point count. As a rough calibration, the figures quoted historically for light atoms — ≈800 points/atom at dftgrid 0, ≈3,000 at dftgrid 1 and ≈80,000 at dftgrid 5 — correspond to 60–80% of the nominal counts above.

DFT grid type: 1
Using dynamic DFT grids.
DFT grid points: 30124 (3012 points/atom)

dftgrid turns off dynamic grids

Specifying dftgrid explicitly also sets dynamicgrid no, even though the documented default for dynamicgrid is yes. If you want both a non-default grid and dynamic grids, request them together:

dftgrid     2
dynamicgrid yes

Construction

Atomic rows

Every atom is assigned a row index 0–4 that selects its radial count and angular order. Elements beyond the fifth row reuse row 4, since no grids are tabulated past that point.

Row Elements
0 H, He
1 Li–Ne
2 Na–Ar
3 K–Kr
4 Rb and heavier

Radial quadrature

For an atom with \(N_r\) radial shells and atomic radius \(R_A\), shell \(i\) sits at

\[ q_i = \frac{i}{N_r}, \qquad r_i = \left(\frac{q_i}{1-q_i}\right)^{2}, \qquad r^{\mathrm{phys}}_i = R_A\, r_i , \qquad i = 0,\dots,N_r-1 , \]

i.e. the MHL \(m=2\) map applied to a uniform grid in \(q\). The associated radial weight is

\[ w^{\mathrm{rad}}_i = \frac{1}{N_r}\, \left(\frac{dr^{\mathrm{phys}}}{dq}\right)_{\!i} \big(r^{\mathrm{phys}}_i\big)^2 , \qquad \frac{dr^{\mathrm{phys}}}{dq} = \frac{2 R_A\, q}{(1-q)^3} , \]

and the full pre-partitioning weight of a grid point is \(w_g = 4\pi\, w^{\mathrm{ang}}\, w^{\mathrm{rad}}_i\), with \(w^{\mathrm{ang}}\) the unit-sphere Lebedev weight. Note that \(i=0\) places a shell at the nucleus with zero weight; those points are removed by the screening step below.

Radii are taken from TeraChem's internal table (values in Å, converted to bohr internally):
H 0.53 He 0.31 Li 1.63 Be 1.09 B 0.81 C 0.65 N 0.54 O 0.47 F 0.41
Ne 0.36 Na 2.16 Mg 1.67 Al 1.36 Si 1.15 P 0.99 S 0.87 Cl 0.78 Ar 0.71

The same radii are reused for the Becke size adjustment.

Basis-dependent radial cutoff. Radial shells are generated outward and the loop stops at the first shell satisfying

\[ \big(r^{\mathrm{phys}}\big)^2 > \frac{50}{\alpha^A_{\min}} , \]

where \(\alpha^A_{\min}\) is the smallest Gaussian primitive exponent on atom \(A\). Beyond that radius every basis function on the atom has decayed below \(e^{-50}\approx 2\times10^{-22}\). Because the cutoff depends on the basis set, the number of surviving shells — and hence the total point count — is basis-dependent: the same dftgrid level yields different grid sizes for 6-31G and aug-cc-pVTZ.

The nominal shell counts \(N_r\) per row and grid level are:

Row dftgrid 0 1 2 3 4 5
0 (H, He) 30 50 75 85 100 140
1 (Li–Ne) 30 50 75 85 100 140
2 (Na–Ar) 45 75 95 110 125 175
3 (K–Kr) 75 95 110 130 160 210
4 (Rb+) 105 130 155 180 205 235

Angular quadrature and pruning

Each radial shell carries a Lebedev–Laikov sphere. The unpruned angular order for a given row and grid level is

\[ \ell_{\max} = \ell^{(0)}_{\mathrm{row}} + \Delta_{\texttt{dftgrid}} , \]

with \(\ell^{(0)} = \{17, 17, 21, 25, 31\}\) for rows 0–4 and \(\Delta = \{0, 6, 12, 18, 30, 42\}\) for dftgrid 0–5. Requesting order \(\ell\) returns the smallest tabulated Lebedev rule complete through \(\ell\) (e.g. \(\ell=17 \rightarrow 110\) points, \(\ell=23 \rightarrow 194\), \(\ell=31 \rightarrow 350\)).

Pruning divides the radial range into five regions, indexed \(j=0,\dots,4\), by comparing the reduced radius \(r_i\) (before multiplication by \(R_A\)) against row-dependent boundaries \(\alpha_j\):

Row \(r < \alpha_0\) \(\alpha_0 \le r < \alpha_1\) \(\alpha_1 \le r < \alpha_2\) \(\alpha_2 \le r < \alpha_3\) \(r \ge \alpha_3\)
0 (H, He) \(r<0.25\) \(<0.5\) \(<0.9\) \(<4.5\) rest
1 (Li–Ne) \(r<0.1667\) \(<0.5\) \(<0.8\) \(<4.5\) rest
2 (Na–Ar) \(r<0.1\) \(<0.4\) \(<0.7\) \(<2.5\) rest
3 (K–Kr) — — all \(r\) — —
4 (Rb+) — — all \(r\) — —

In region \(j\) the angular order is reduced to

\[ \ell_j = \max\!\big(\ell_{\max} - \delta_j,\ \mathrm{round}(f_j\,\ell_{\max})\big), \]
\[ \delta_j = \{26, 18, 12, 0, 12\}\ \text{(rows 0–1)}, \quad \delta_j = \{36, 24, 12, 0, 12\}\ \text{(rows 2–4)}, \quad f_j = \{0.1,\, 0.4,\, 0.6,\, 1.0,\, 0.6\} . \]

Region 3 is unpruned (\(\ell_3 = \ell_{\max}\)): it is the chemically important valence region. Regions 0–2 coarsen toward the nucleus, and region 4 coarsens again in the diffuse tail.

Rows 3 and 4 are effectively unpruned — and permanently coarsened

For K–Kr and Rb+ the boundary table degenerates (\(\alpha = \{0, 0, 10^{10}, 0\}\)), so every radial shell falls in region 2 and uses \(\ell \approx 0.6\,\ell_{\max}\). Heavy atoms therefore never see the full \(\ell_{\max}\) sphere at any radius — e.g. at dftgrid 1 a K atom uses 146 angular points on all 95 shells rather than up to 350. If you need converged XC integration on heavy elements, validate against dftgrid 3 or higher rather than assuming the nominal \(\ell_{\max}\) applies.

Resulting angular point counts per region:

dftgrid Row \(\ell_{\max}\) Region 0 1 2 3 4
0 0–1 17 6 26 50 110 50
0 2 21 6 38 74 170 74
0 3 25 6 50 86 230 86
0 4 31 6 74 146 350 146
1 0–1 23 6 38 86 194 86
1 2 27 6 50 110 266 110
1 3 31 6 74 146 350 146
1 4 37 14 86 230 590 230
2 0–1 29 6 74 110 302 110
2 2 33 6 74 170 434 170
2 3 37 14 86 230 590 230
2 4 43 26 146 350 770 350
3 0–1 35 38 110 194 434 194
3 2 39 14 110 266 590 266
3 3 43 26 146 350 770 350
3 4 49 74 230 590 974 590
4 0–1 47 170 302 434 770 434
4 2 51 86 266 590 974 590
4 3 55 146 350 770 1202 770
4 4 61 230 590 974 1454 974
5 0–1 59 434 590 770 1202 770
5 2 63 266 590 974 1454 974
5 3 67 350 770 1202 1730 1202
5 4 73 590 974 1454 2030 1454

(For rows 3–4 only the region-2 column is ever used.)

The 74-, 230- and 266-point rules have negative weights

Three of the Lebedev–Laikov rules appearing in the table above — 74, 230 and 266 points — are not all-positive quadratures. They occur in ordinary settings: dftgrid 1 uses the 266-point rule in the valence region of Na–Ar, and dftgrid 2 uses the 74-point rule for H–Ne.

For XC quadrature this is harmless. The weights enter linearly, and the weight screening compares \(|w_g|\) against the threshold, so negative-weight points are retained and integrate correctly.

It matters wherever a weight is raised to a fractional power. Two places in TeraChem do that, and both handle it:

  • PCM takes \(\sqrt{w}\) as a segment area, so the affected Lebedev orders must not be used for cavity discretization — see PCM Solvation.
  • THC forms its collocation factors as \(X^P_\mu = w_P^{1/4}\phi_\mu\) and therefore discards every grid point with \(w_g \le 0\) before taking the fourth root. A THC grid derived from a dftgrid preset is consequently a strict subset of the corresponding XC grid, with slightly fewer points.

Becke partitioning and weight screening

Atomic quadratures are combined into a molecular one with Becke's fuzzy-cell scheme. For a point \(\mathbf{r}_g\) belonging to atom \(A\),

\[ w_g \leftarrow w_g \, \frac{P_A(\mathbf{r}_g)}{\sum_B P_B(\mathbf{r}_g)}, \qquad P_A(\mathbf{r}) = \prod_{B\ne A} s(\nu_{AB}), \]

with the elliptical coordinate \(\mu_{AB} = (r_A - r_B)/R_{AB}\) and the atomic-size adjustment

\[ \nu_{AB} = \mu_{AB} + a_{AB}\,(1 - \mu_{AB}^2), \qquad a_{AB} = \frac{u_{AB}}{u_{AB}^2 - 1}, \quad u_{AB} = \frac{\chi_{AB}-1}{\chi_{AB}+1}, \quad \chi_{AB} = \frac{R_A}{R_B}, \]

with \(a_{AB}\) clamped to \([-0.5, 0.5]\). The switching function is Becke's iterated polynomial with \(k=3\):

\[ s(\nu) = \tfrac12\big(1 - p(p(p(\nu)))\big), \qquad p(x) = \tfrac32 x - \tfrac12 x^3 . \]

The cell-function product is truncated early once it drops below the weight threshold, and after partitioning all points with

\[ |w_g| < w_{\mathrm{thre}} = 10^{-11} \]

are discarded. threall tightens this: \(w_{\mathrm{thre}} = \min(10^{-11}, \texttt{threall})\). The surviving points are then binned into 4 bohr cubic boxes for GPU locality, which reorders the point list but does not change the grid.

Summary of the pipeline

for each atom A:
    N_r, l_max        <- row(A), dftgrid
    for i = 0 .. N_r-1:
        r  = (q/(1-q))^2 * R_A                    # MHL radial map, q = i/N_r
        if r^2 > 50/alpha_min(A): break           # basis-dependent cutoff
        pick pruned Lebedev sphere for region(r)  # 5 regions
        emit points with w = 4*pi*w_ang*w_rad
apply Becke partitioning (size-adjusted, k=3)
drop points with |w| < 1e-11
pack surviving points into 4 bohr boxes

Nominal grid sizes

Nominal points per atom, i.e. before the basis-dependent radial cutoff and the Becke weight screening. Actual counts are lower, typically by 20–40% for light atoms with a double-\(\zeta\) basis, and are reported at printlevel 4.

Row dftgrid 0 1 2 3 4 5
0 (H, He) 1,328 3,720 8,018 15,218 39,752 101,692
1 (Li–Ne) 1,368 3,892 8,482 15,746 40,748 103,492
2 (Na–Ar) 3,002 7,330 14,994 25,468 59,974 143,846
3 (K–Kr) 6,450 13,870 25,300 45,500 123,200 252,420
4 (Rb+) 15,330 29,900 54,250 106,200 199,670 341,690

Worked example: dftgrid 0 for H, C, N, O

\(N_r = 30\) and \(\ell_{\max} = 17\) for both rows, giving pruned spheres of 6, 26, 50, 110 and 50 points across the five regions. Because the region boundaries differ between row 0 and row 1, the shells distribute as

Atom Shells per region Nominal points
H 11, 2, 2, 6, 9 1,328
C, N, O 9, 4, 2, 6, 9 1,368

before the radial cutoff and weight screening.

Dynamic grids

With dynamicgrid yes (the default unless dftgrid is given explicitly), TeraChem builds two grids — level 0 and level dftgrid — and switches between them during the SCF. Iterations run on grid 0 while the SCF error measure (the DIIS error, or the gradient error when using SOSCF) exceeds gridthre (default 0.01); once it falls below, the calculation finishes on grid dftgrid. The switch is announced at printlevel 4:

|                 >>> SWITCHING TO GRID 0 <<<
|                 >>> SWITCHING TO GRID 1 <<<

The switch down to grid 0 happens at most once per SCF. Dynamic grids are not available for periodic systems.

Because the final iterations — and all energies, gradients and properties — use grid dftgrid, dynamic grids affect only the SCF path, not the converged result. They can, however, slightly change the number of iterations.

  • Periodic DFT can integrate the XC term on this molecular grid (default) or on a uniform real-space grid; see Integration grids (DFT).
  • THC builds its collocation factors on a quadrature grid that may be this same molecular grid — thcfitbox_grid dftgrid, with the level set by thcfitbox_dftgrid (defaults to dftgrid). Points with non-positive Becke weights are dropped first, so the THC grid is a subset of the XC grid at the same level. See Tensor Hypercontraction.
  • PCM cavity discretization uses independent Lebedev surface grids (pcmgrid_heavy, pcmgrid_h, …); see PCM Solvation.

Keyword summary

Keyword Type Default Description
dftgrid int 0–5 1 Molecular XC grid level. Setting it also disables dynamic grids unless dynamicgrid yes is given as well.
dynamicgrid yes/no yes (no if dftgrid is set) Converge the SCF on grid 0, then finish on grid dftgrid.
gridthre float 0.01 SCF error below which the dynamic grid switches to dftgrid.
threall float not set Global threshold; tightens the grid-weight cutoff to \(\min(10^{-11}, \texttt{threall})\).
thcfitbox_dftgrid int 0–5 dftgrid Grid level for the THC collocation grid when thcfitbox_grid dftgrid.
periodic_uniform_grid_size int,int,int not set Use a uniform real-space grid instead (periodic only).
periodic_uniform_grid_density int not set Uniform-grid point density per unit cell (periodic only).

References

  1. C. W. Murray, N. C. Handy, and G. J. Laming, "Quadrature schemes for integrals of density functional theory," Mol. Phys. 78, 997–1014 (1993). doi:10.1080/00268979300100651 ↩

  2. V. I. Lebedev and D. N. Laikov, "A quadrature formula for the sphere of the 131st algebraic order of accuracy," Dokl. Math. 59, 477–481 (1999). ↩

  3. A. D. Becke, "A multicenter numerical integration scheme for polyatomic molecules," J. Chem. Phys. 88, 2547–2553 (1988). doi:10.1063/1.454033 ↩

  4. P. M. W. Gill, B. G. Johnson, and J. A. Pople, "A standard grid for density functional calculations," Chem. Phys. Lett. 209, 506–512 (1993). doi:10.1016/0009-2614(93)80125-9 ↩