NCI and Volumetric Data with Tako Script

Tako Script has two distinct volumetric routes. tako.singlePoint(..., { properties: ["nci"] }) invokes the native fixed-grid GFN2 NCI implementation. tako.cube and tako.grid parse and transform cube data in JavaScript. The second route is useful for inspection and derived fields, but it does not reproduce native NCI automatically because coordinate units and boundary stencils differ. Theory belongs in Density Fields, NCI Analysis, and Volumetric Data; viewer actions belong in Inspect Cube and Grid Data.

Request native NCI output

const complex = await tako.input("input/complex.extxyz");
const calculator = tako.calculator.gfn2({
  dispersion: "d4",
  charge: 0,
  unpaired: 0,
});

const result = await tako.singlePoint(complex, {
  calculator,
  properties: ["energy", "forces", "nci"],
});

await tako.output("output/nci-result.json", result);
await tako.output("output/nci-rdg.cube", result.nci.rdg_cube);
await tako.output(
  "output/nci-sign-lambda2-rho.cube",
  result.nci.sign_lambda2_rho_cube,
);

This runs the configured GFN2 calculator, evaluates electron density on the hardcoded $28^3$ box, and calculates native RDG and signed-density fields. There is no top-level tako.analysis.nci API and no public NCI grid option.

Do not pass workflow-style densityFloor, gridSpacingAngstrom, paddingAngstrom, or isovalue expecting them to affect this call. The native density floor is $10^{-8}$, padding is 3 bohr, shape is fixed, and isovalue is viewer state rather than calculation input.

Read the native result

const { summary, rdg_cube, sign_lambda2_rho_cube } = result.nci;

console.log({
  lowRdgVoxels: summary.low_rdg_points,
  allInteriorNegativeSignedDensity: summary.attractive_points,
  allInteriorPositiveSignedDensity: summary.repulsive_points,
  allInteriorNearZeroSignedDensity: summary.near_zero_points,
});

low_rdg_points counts interior voxels with $s\le0.5$. The other three counters independently partition every interior voxel by a $\pm10^{-12}$ signed-density threshold. They are not restricted to low RDG. Do not normalize or regress them as interaction volumes without including voxel volume and a scientifically defined joint mask.

The two cube strings are complete Gaussian cube files. The result does not return parsed origin/axes/shape, scalar units, grid spacing, or source identifiers as separate fields.

Parse a cube

const rdg = tako.cube.parse(result.nci.rdg_cube);
const signed = tako.cube.parse(result.nci.sign_lambda2_rho_cube);

console.log(rdg.role);    // "nci_rdg"
console.log(rdg.shape);   // [28, 28, 28]
console.log(rdg.origin);  // Å
console.log(rdg.axes);    // voxel vectors in Å
console.log(rdg.values);  // plain number[] at Script boundary

tako.cube.parse converts cube origin, voxel vectors, and atom coordinates from bohr to ångström. Scalar values are unchanged. The core parser stores a typed array, but the Script boundary returns a plain number array, increasing memory use for large cubes.

Validate grid identity before combining fields:

function sameVector(left, right, tol = 1e-10) {
  return left.every((value, index) => Math.abs(value - right[index]) <= tol);
}

function sameGrid(left, right) {
  return left.shape.every((value, index) => value === right.shape[index])
    && sameVector(left.origin, right.origin)
    && left.axes.every((axis, index) => sameVector(axis, right.axes[index]));
}

if (!sameGrid(rdg, signed)) {
  throw new Error("NCI fields are not aligned");
}

Native NCI siblings should align. Imported or independently generated cubes may not.

Understand coordinate and derivative units

Native NCI differentiates density on a grid whose axes are in bohr. Script parsing converts axes to ångström, so tako.grid.gradient, hessian, laplacian, and derivatives operate per ångström or per ångström squared when given the parsed object.

If scalar density remains e bohr$^{-3}$ while coordinates are Å, inserting Script gradients directly into the atomic-unit RDG equation mixes units. Either convert axes back to bohr before differentiation or convert density and derivatives consistently.

For axis conversion without changing scalar order:

const ANGSTROM_PER_BOHR = 0.52917721092;
const atomicUnitGrid = {
  ...rdg,
  origin: rdg.origin.map((value) => value / ANGSTROM_PER_BOHR),
  axes: rdg.axes.map((axis) =>
    axis.map((value) => value / ANGSTROM_PER_BOHR)
  ),
};

This example uses an already dimensionless RDG field only to demonstrate geometry conversion. To reconstruct NCI, begin from electron density, preserve its e bohr$^{-3}$ values, and use bohr axes.

Use grid primitives

The grid API includes point/index access, scalar/vector/tensor access, first and second derivatives, gradient, Hessian, Laplacian, mapping, zipping, rectangular integration, and histograms. Representative calls are:

const cube = tako.cube.parse(cubeText);
const index = tako.grid.offset(cube.shape, 10, 12, 8);
const point = tako.grid.point(cube, 10, 12, 8);
const value = tako.grid.at(cube, 10, 12, 8);
const gradient = tako.grid.gradient(cube, 10, 12, 8);
const hessian = tako.grid.hessian(cube, 10, 12, 8);
const laplacian = tako.grid.laplacian(cube, 10, 12, 8);
const integral = tako.grid.integrate(cube);
const histogram = tako.grid.histogram(cube, 80);

integrate is a rectangular sum $\sum_i f_i|\det\mathbf A|$; it applies no endpoint weights and does not estimate truncation error. histogram is a distribution of sampled values, not a spatially connected-region analysis.

zip requires identical shape and origin/axes within $10^{-10}$. Grid operations accept value arrays at least as long as the shape product and ignore extras; cube writing requires exact length.

Native and Script boundary stencils differ

Native NCI uses one-sided first derivatives at boundaries and sets every second/mixed derivative touching the boundary to zero. Script Hessian evaluation instead clamps the boundary stencil centre to the nearest interior location. A Script reconstruction can match interior formulas but not native outer-shell values without implementing the native policy explicitly.

Exclude boundary voxels from comparisons:

function isInterior([nx, ny, nz], i, j, k) {
  return i > 0 && j > 0 && k > 0
    && i < nx - 1 && j < ny - 1 && k < nz - 1;
}

Even in the interior, confirm coordinate units and density floor before expecting parity.

Construct a joint NCI mask

A scientifically declared joint mask is more meaningful than the runtime’s independent counters:

const rdgGrid = tako.cube.parse(result.nci.rdg_cube);
const signedGrid = tako.cube.parse(result.nci.sign_lambda2_rho_cube);
if (!sameGrid(rdgGrid, signedGrid)) throw new Error("grid mismatch");

const selected = [];
for (let i = 1; i < rdgGrid.shape[0] - 1; i += 1) {
  for (let j = 1; j < rdgGrid.shape[1] - 1; j += 1) {
    for (let k = 1; k < rdgGrid.shape[2] - 1; k += 1) {
      const offset = tako.grid.offset(rdgGrid.shape, i, j, k);
      const s = rdgGrid.values[offset];
      const signedRho = signedGrid.values[offset];
      if (s <= 0.5 && signedRho < -1e-12) {
        selected.push(offset);
      }
    }
  }
}

This counts jointly selected voxels but still is not a connected-component count or interaction energy. A volume estimate would multiply by voxel volume and requires a grid-convergence policy that Tako’s native fixed grid cannot supply.

Write a derived cube

tako.cube.writeLike reuses template geometry/atoms and writes one scalar per voxel:

const maskValues = rdgGrid.values.map((s, index) =>
  s <= 0.5 && signedGrid.values[index] < -1e-12 ? 1 : 0
);

const maskCube = tako.cube.writeLike(rdgGrid, maskValues, {
  title: "Joint attractive-like low-RDG mask",
  comment: "RDG <= 0.5 and sign(lambda2)rho < -1e-12",
  role: "unknown",
});

await tako.output("output/nci-joint-mask.cube", maskCube);

The writer converts Å geometry back to bohr, preserves scalar values unchanged, and writes six values per line. It requires finite values and exact shape-product length. role adds tako_cube_role; it does not attach scalar units, method provenance, or arbitrary structured metadata. Put thresholds in the comment and preserve a separate report.

Handle imported cubes cautiously

The parser assumes bohr geometry even when voxel counts are negative. Convert standards-compliant Å-unit cubes externally before using Script values. Negative atom-count orbital-dataset variants are not supported. Extra datasets can be ignored or misread.

Role inference is not complete provenance. Confirm title/comment, units, field formula, atom coordinates, grid axes, and generation method.

Preserve workflow and failure boundaries

Native NCI is a post-processing stage of tako.singlePoint, not a standalone resumable calculation. Earlier cumulative partial results may remain if NCI or a later property fails. Cancellation hard-terminates the worker and provides no density-grid checkpoint.

The visual workflow’s NCI node exposes density floor, isovalue, spacing, and padding settings, but the current runner ignores all four and dispatches the same fixed native property request (Known Limitations). Those controls do not affect the calculation, so they do not belong in a record of what was computed.

The exhaustive single-point result schema is generated in tako.singlePoint Settings. The full cube/grid namespace list is generated in Math and Grid API.