Density Fields, NCI Analysis, and Volumetric Data

A volumetric result is a scalar field sampled on a three-dimensional grid. Its scientific meaning comes from the field definition, units, source wavefunction, grid geometry, derivative method, and visualization threshold together. An isosurface is only the set of interpolated points where that sampled field equals a chosen value. Changing the isovalue changes the surface even though the underlying electronic calculation is unchanged.

Tako can generate electron density, a Gaussian-smoothed Mulliken charge field, non-covalent-interaction (NCI) fields, molecular orbitals and orbital densities, spin density, electrostatic potential, density derivatives, Fukui functions, a dual descriptor, and an electron-localization-function field from the GFN2 single-point path. It does not treat these fields as interaction energies or experimental observables.

Scalar fields and grids

A scalar field $f(\mathbf r)$ assigns a number to each spatial position. A regular grid stores samples

where

$\mathbf r_0$ is the origin; $\mathbf a_1$, $\mathbf a_2$, and $\mathbf a_3$ are voxel vectors; and $(N_1,N_2,N_3)$ is the shape. The parallelepiped represented by the samples extends from the origin to $(N_\alpha-1)\mathbf a_\alpha$ along each axis. The volume per grid point is

Resolution and extent are separate choices. A larger box at fixed shape has larger voxels and less accurate derivatives. A finer shape at fixed box costs more memory and evaluation time. A field may also decay beyond the finite box, so numerical integration can omit tails even when local surfaces appear plausible.

Electron density and related fields

For a molecular wavefunction expressed in basis functions $\chi_\mu$, a one-particle density can be written

with a density matrix $\mathbf P$. Electron density is nonnegative and has units of electrons per volume. Its integral over all space gives electron count for a properly normalized complete representation. A finite numerical grid approximates that integral and truncates the density outside its box.

Tako’s electron_density_cube is evaluated from the configured GFN2 wavefunction on the common grid. Charge, unpaired electrons, SCF convergence, implicit solvation, geometry, and method settings can therefore affect it. The cube does not contain an electron-count integral or discretization error estimate.

Tako’s charge_density_cube has a different definition. It is not the electronic density, not nuclear charge minus electron density, and not a charge-deformation density. Mulliken atomic net charges $q_A$ are spread as normalized isotropic Gaussians:

with hardcoded $\sigma=0.6$ bohr. Positive and negative lobes visualize the chosen population analysis after smoothing. They are not accumulations and depletions of the true electronic density unless an actual density difference was computed.

An orbital cube stores signed molecular-orbital amplitude $\psi_p(\mathbf r)$; an orbital-density cube stores $|\psi_p(\mathbf r)|^2$. An amplitude surface needs positive and negative lobes because phase sign carries nodal information. Squared density is nonnegative and discards phase. Tako can return a limited frontier window and optional squared cubes, together with orbital index, energy, and occupation metadata.

Spin density is $\rho_\alpha-\rho_\beta$. A positive/negative surface describes excess spin-channel density, not positive/negative electric charge. For a restricted closed-shell state it should vanish apart from numerical behavior; the runtime can create a companion spin-polarized calculation when necessary.

The electrostatic-potential cube is constructed from the current atomic multipole representation with short-range regularization. It is not a full continuous electrostatic potential evaluated directly from nuclear point charges and the exact electron density. Mapping this field onto a separate density isosurface can show relative electrostatic character, but colors remain model- and scale-dependent.

Reduced density gradient

The reduced density gradient (RDG) is the dimensionless quantity

It compares spatial density variation with the local density scale of a uniform electron gas. Near nuclei and covalent bonds, density is high and the interpretation differs from weak low-density contacts. In low-density regions between fragments, a low RDG can reveal slowly varying density associated with a contact surface.

Tako evaluates

using a hardcoded density floor of $10^{-8}$ in internal atomic units. The absolute value appears only in the denominator guard. The output RDG is dimensionless.

RDG alone identifies a slowly varying density region but does not classify its local curvature as attractive-like or repulsive-like. That classification uses the density Hessian

Let its ordered eigenvalues be $\lambda_1\le\lambda_2\le\lambda_3$. NCI analysis commonly examines

For nonnegative electron density, negative values correspond to $\lambda_2<0$ and are conventionally associated with attractive-like concentration; positive values correspond to $\lambda_2>0$ and repulsive-like depletion. Values near zero often accompany weak, diffuse contacts. These labels describe density topology, not a partitioned force or interaction energy.

Finite differences and boundary behavior

Tako evaluates the density on a regular Cartesian grid. Interior gradients use two-point central differences. At the outer grid planes, gradients use first-order one-sided differences. The derivative vectors are transformed with the inverse transpose of the voxel-axis matrix, so the calculation accounts for physical axis spacing.

Interior Hessian diagonal terms use centered second differences; mixed terms use four-corner stencils. Any pure or mixed second derivative touching the boundary is set to zero. The complete outer voxel shell therefore receives $\lambda_2=0$ and a zero signed-density value by construction. NCI interpretation must focus away from the box boundary.

Derivative error scales with voxel size. Since the grid shape is fixed while the box expands with molecular extent, a larger molecule receives coarser derivatives. Tako exposes no NCI grid-convergence control. Qualitative comparison between differently sized systems is therefore entangled with grid resolution.

Tako’s common volumetric grid

Every generated cube uses one hardcoded template:

The padding is 3 bohr, not 3 Å. Grid spacing is neither fixed nor user-configurable. For periodic inputs, the field box still follows Cartesian atom extrema rather than the lattice vectors, so it is not a periodic unit-cell density representation.

This policy is adequate for compact interactive visualization but weak for quantitative derivatives, large systems, diffuse density, cross-system voxel counts, or periodic topology. A surface disappearing at the box edge is a grid defect.

NCI summary counters

Tako returns four integer counters. They inspect only strict interior grid points, $26^3=17{,}576$ points for the current template:

The last three counters partition all interior points. They are not restricted to the low-RDG subset. They do not count connected interaction regions, and they are not volumes because voxel volume is not applied. They change with box extent and spacing. Use them as runtime summaries of one particular grid, not as cross-molecule measures of attraction or repulsion.

Gaussian cube representation

A Gaussian cube file begins with two text lines, then atom count and origin, three voxel-count/axis records, atom records, and scalar values. Tako-generated files use positive voxel counts and bohr geometry. The frontend converts origin, voxel vectors, and atom positions to ångström when parsing. Scalar values are left unchanged.

Tako writes tako_cube_role="..." in the comment line to identify electron density, charge density, NCI RDG, NCI signed density, molecular orbital, orbital density, spin density, ESP, density-gradient norm, density Laplacian, Fukui fields, dual descriptor, or ELF. The header does not record scalar units, method settings, charge/spin, SCF convergence, density floor, derivative stencil, grid-padding policy, or source structure identifier. Preserve result.json and the input/protocol with every cube.

The cube atom-record charge column is written as zero rather than Mulliken or nuclear charge. Do not read it as the charge used to make charge_density_cube.

Imported cube support is intentionally narrower than the complete cube standard. Negative voxel counts conventionally signal ångström geometry, but the parser takes their absolute values and still applies the bohr-to-ångström conversion, shrinking such files by 0.529177. Negative atom-count orbital-dataset extensions are not supported. Extra scalar tokens are ignored after the first shape product. Imported cubes are stored under a generic electron-density artifact kind even when header role inference recognizes another field. Verify imported geometry and field identity before interpretation.

Isosurfaces, sections, and the viewer

Building an isosurface, choosing its sign and default isovalue, mapping a sibling cube’s percentile-limited colors onto it, and slicing or clipping a section are display and interaction behavior, not additional physics. That behavior — including the standard NCI pairing of nci_rdg.cube and nci_sign_lambda2_rho.cube — is documented in Inspect Cube and Grid Data.

Other density-derived fields

Density-gradient norm and Laplacian cubes expose $|\nabla\rho|$ and $\nabla^2\rho$. They inherit the fixed grid and finite-difference boundary behavior. The Laplacian is the trace of the density Hessian and can identify concentration/depletion patterns, but is sensitive to resolution.

Fukui fields are finite differences between separately converged $N$, $N+1$, and $N-1$ states. Tako generates $f^+$, $f^-$, and a dual descriptor with a heuristic unpaired-electron choice for charged companions. They are conceptual-DFT screening fields, and failures in either companion SCF invalidate the stage.

ELF is computed only through the restricted-spin path from occupied-orbital kinetic density and electron density with an internal floor. It is dimensionless and bounded by its implemented formula, but an isovalue still needs stated justification. Spin-polarized systems lie outside this current public ELF implementation.

Interpretation and validation

An NCI surface shows where a chosen approximate density has low RDG and how the local Hessian curvature is signed. It does not supply interaction energy, binding enthalpy, free energy, population, bond order, or causal force. A broad attractive-colored surface can coexist with an unfavorable total association energy once deformation, entropy, solvation, and reference fragments are treated consistently.

For a defensible result:

1. Verify structure, charge, spin, GFN2 settings, SCF convergence, and any solvent model.
2. Record the fixed $28^3$ grid, 3-bohr padding, actual voxel vectors, and finite box.
3. Inspect both RDG and signed-density raw ranges and confirm the contact lies away from boundaries.
4. Open RDG as the surface and explicitly color it by the aligned signed-density sibling.
5. State the RDG isovalue and color scale; do not rely on the field-agnostic defaults.
6. Repeat the visualization across reasonable isovalues and confirm the interpretation is stable.
7. Use geometrical distances and a consistently defined energy cycle when interaction strength matters.
8. Avoid comparing raw summary voxel counts between systems with different extents.

The current grid is not configurable, so Tako cannot provide an internal grid-convergence study. For quantitative density topology, export the wavefunction/density to software with controllable high-resolution grids and established analysis procedures.

The interface sequence is documented in Run a Single-Point Calculation (request the NCI cubes) and Inspect Cube and Grid Data (build the surface). Programmatic parsing, derivatives, and file writing are kept in NCI and Volumetric Data with Tako Script.

References and credit

• E. R. Johnson, S. Keinan, P. Mori-Sánchez, J. Contreras-García, A. J. Cohen, and W. Yang, “Revealing Noncovalent Interactions,” J. Am. Chem. Soc. 132, 6498–6506 (2010), DOI: 10.1021/ja100936w.
• J. Contreras-García, E. R. Johnson, S. Keinan, R. Chaudret, J.-P. Piquemal, D. N. Beratan, and W. Yang, “NCIPLOT: A Program for Plotting Noncovalent Interaction Regions,” J. Chem. Theory Comput. 7, 625–632 (2011), DOI: 10.1021/ct100641a.