qstack.fields.dori

Density Overlap Regions Indicator (DORI) computation.

qstack.fields.dori.compute_dori(rho, drho_dr, d2rho_dr2, eps=0.0001)[source]

Compute Density Overlap Regions Indicator (DORI) analytically.

Parameters:

  • rho (numpy ndarray) – 1D array (ngrids,) of electron density.

  • drho_dr (numpy ndarray) – 2D array (3, ngrids) of density first derivatives.

  • d2rho_dr2 (numpy ndarray) – 3D array (3, 3, ngrids) of density second derivatives.

  • eps (float) – Density threshold below which DORI is set to zero. Defaults to 1e-4.

Returns:

1D array (ngrids,) of DORI values ranging from 0 to 1.

Return type:

numpy ndarray

qstack.fields.dori.compute_dori_num(mol, coords, dm=None, c=None, eps=0.0001, dx=0.0001)[source]

Compute DORI using numerical differentiation (semi-numerical approach).

Alternative to analytical DORI calculation using finite differences for derivatives of k², where k=dρ/dr. Useful for validation or when analytical gradients are problematic.

Parameters:

  • mol (pyscf.gto.Mole) – pyscf Mole object.

  • coords (numpy ndarray) – 2D array (ngrids, 3) of grid coordinates in Bohr.

  • dm (numpy ndarray, optional) – 2D density matrix in AO basis. Conflicts with c.

  • c (numpy ndarray, optional) – 1D density fitting coefficients. Conflicts with dm.

  • eps (float) – Density threshold below which DORI is zero. Defaults to 1e-4.

  • dx (float) – Finite difference step size in Bohr. Defaults to 1e-4.

Returns:

(dori, rho) where both are 1D arrays (ngrids,) containing

DORI values and electron density respectively.

Return type:

tuple

qstack.fields.dori.compute_rho(mol, coords, dm=None, c=None, deriv=2, eps=0.0001)[source]

Calculate electron density and derivatives efficiently.

Computes density and its spatial derivatives on a grid from either a density matrix or fitting coefficients, with optimizations for numerical stability.

Parameters:

  • mol (pyscf.gto.Mole) – pyscf Mole object.

  • coords (numpy ndarray) – 2D array (ngrids, 3) of grid coordinates in Bohr.

  • dm (numpy ndarray, optional) – 2D density matrix in AO basis. Conflicts with c.

  • c (numpy ndarray, optional) – 1D density fitting coefficients. Conflicts with dm.

  • deriv (int) – Maximum derivative order (0, 1, or 2). Defaults to 2.

  • eps (float) – Minimum density threshold below which derivatives are set to zero. Defaults to 1e-4.

Returns:

Depending on deriv value:

  • deriv=0: rho as (ngrids,) numpy ndarray,

  • deriv=1: (rho, drho_dr) where drho_dr is (3, ngrids) numpy ndarray,

  • deriv=2: (rho, drho_dr, d2rho_dr2) where d2rho_dr2 is (3, 3, ngrids) numpy ndarray.

Return type:

tuple

Raises:

RuntimeError – If both or neither of dm and c are provided.

qstack.fields.dori.compute_s2rho(rho, d2rho_dr2, eps=0.0001)[source]

Compute signed density based on second eigenvalue of the density Hessian.

Useful for distinguishing bonding vs. non-bonding regions. The sign of the second eigenvalue of the Hessian indicates local density topology.

Parameters:

  • rho (numpy ndarray) – 1D array (ngrids,) of electron density values.

  • d2rho_dr2 (numpy ndarray) – 3D array (3, 3, ngrids) of density second derivatives (Hessian).

  • eps (float) – Density threshold below which values are set to zero. Defaults to 1e-4.

Returns:

1D array (ngrids,) containing rho * sign(λ₂) where λ₂ is the

second eigenvalue of the Hessian, or 0 where rho < eps.

Return type:

numpy ndarray

qstack.fields.dori.dori(mol, dm=None, c=None, eps=0.0001, alg='analytical', grid_type='dft', grid_level=1, nx=80, ny=80, nz=80, resolution=None, margin=3.0, cubename=None, dx=0.0001, mem=1, progress=False)[source]

High-level interface to compute DORI with automatic grid generation and file output.

Reference:

P. de Silva, C. Corminboeuf, “Simultaneous visualization of covalent and noncovalent interactions using regions of density overlap”, J. Chem. Theory Comput. 10, 3745–3756 (2014), doi:10.1021/ct500490b

Computes the Density Overlap Regions Indicator (DORI). Automatically generates appropriate grids and optionally saves results to cube files for visualization.

DORI is a density-based descriptor for identifying covalent bonding regions, with values close to 1 indicating strong electron sharing (covalent bonds).

DORI(r) = γ(r) = θ(r) / (1 + θ(r)), where: θ = |∇(k²)|² / |k|⁶, and k(r) = ∇ρ(r) / ρ(r)

Parameters:

  • mol (pyscf.gto.Mole) – pyscf Mole object.

  • dm (numpy ndarray, optional) – 2D density matrix in AO basis. Conflicts with c.

  • c (numpy ndarray, optional) – 1D density fitting coefficients. Conflicts with dm.

  • eps (float) – Density threshold for DORI. Defaults to 1e-4.

  • alg (str) – Algorithm: ‘analytical’ or ‘numerical’. Defaults to ‘analytical’.

  • grid_type (str) – Grid type: ‘dft’ for DFT quadrature grid or ‘cube’ for uniform grid. Defaults to ‘dft’.

  • grid_level (int) – For DFT grid, the grid level (higher = more points). Defaults to 1.

  • nx (int) – For cube grid, number of points in x direction. Defaults to 80.

  • ny (int) – For cube grid, number of points in y direction. Defaults to 80.

  • nz (int) – For cube grid, number of points in z direction. Defaults to 80.

  • resolution (float) – For cube grid, grid spacing in Bohr. Conflicts with nx/ny/nz.

  • margin (float) – For cube grid, extra space around molecule in Bohr. Defaults to BOX_MARGIN.

  • cubename (str, optional) – For cube grid, base filename for output cube files. If None, no files saved.

  • dx (float) – For numerical algorithm, finite difference step in Bohr. Defaults to 1e-4.

  • mem (float) – Maximum memory in GiB for AO evaluation. Defaults to 1.

  • progress (bool) – If True, displays progress bar. Defaults to False.

Returns:

(dori, rho, s2rho, coords, weights) containing: - dori (numpy ndarray): 1D array of DORI values - rho (numpy ndarray): 1D array of electron density - s2rho (numpy ndarray): 1D array of signed density (None if numerical) - coords (numpy ndarray): 2D array (ngrids, 3) of grid coordinates - weights (numpy ndarray): 1D array of grid weights

Return type:

tuple

Note

When cubename is provided with cube grid, creates three files: - <cubename>.dori.cube: DORI values

  • <cubename>.rho.cube: electron density

  • <cubename>.sgnL2rho.cube: signed density (analytical only)

qstack.fields.dori.dori_on_grid(mol, coords, dm=None, c=None, eps=0.0001, alg='analytical', mem=1, dx=0.0001, progress=False)[source]

Compute DORI on a user-specified grid with memory management.

Main computational function for DORI evaluation. Handles large grids by chunking based on available memory.

Parameters:

  • mol (pyscf.gto.Mole) – pyscf Mole object.

  • coords (numpy ndarray) – 2D array (ngrids, 3) of grid coordinates in Bohr.

  • dm (numpy ndarray, optional) – 2D density matrix in AO basis. Conflicts with c.

  • c (numpy ndarray, optional) – 1D density fitting coefficients. Conflicts with dm.

  • eps (float) – Density threshold for DORI calculation. Defaults to 1e-4.

  • alg (str) – Algorithm choice: ‘analytical’ or ‘numerical’. Defaults to ‘analytical’.

  • mem (float) – Maximum memory in GiB for AO evaluation. Defaults to 1 GiB.

  • dx (float) – Step size in Bohr for numerical derivatives. Defaults to 1e-4.

  • progress (bool) – If True, displays progress bar. Defaults to False.

Returns:

(dori, rho, s2rho) where: - dori (numpy ndarray): 1D array (ngrids,) of DORI values - rho (numpy ndarray): 1D array (ngrids,) of electron density - s2rho (numpy ndarray): 1D array (ngrids,) of signed density (None if numerical)

Return type:

tuple

qstack.fields.dori.eval_rho_df(ao, c, deriv=2)[source]

Compute electron density and its derivatives from density-fitting coefficients.

Parameters:

  • ao (numpy ndarray) – 3D array of shape (*, ngrids, nao) where: - ao[0]: atomic orbital values on the grid, - ao[1:4]: first derivatives (if deriv>=1), - ao[4:10]: second derivatives in upper triangular form (if deriv==2).

  • c (numpy ndarray) – 1D array of density fitting/expansion coefficients.

  • deriv (int) – Maximum derivative order to compute (0, 1, or 2). Defaults to 2.

Returns:

Depending on deriv value:

  • deriv=0: rho as (ngrids,) numpy ndarray,

  • deriv=1: (rho, drho_dr) where drho_dr is (3, ngrids) numpy ndarray,

  • deriv=2: (rho, drho_dr, d2rho_dr2) where d2rho_dr2 is (3, 3, ngrids) numpy ndarray.

Return type:

tuple

qstack.fields.dori.eval_rho_dm(mol, ao, dm, deriv=2)[source]

Compute electron density and its derivatives from a density matrix.

Modified from pyscf/dft/numint.py to return full second derivative matrices needed for DORI calculations.

Parameters:

  • mol (pyscf.gto.Mole) – pyscf Mole object.

  • ao (numpy ndarray) – 3D array of shape (*, ngrids, nao) where - ao[0]: atomic orbital values on the grid, - ao[1:4]: first derivatives (if deriv>=1), - ao[4:10]: second derivatives in upper triangular form (if deriv==2).

  • dm (numpy ndarray) – Density matrix in AO basis.

  • deriv (int) – Maximum derivative order to compute (0, 1, or 2). Defaults to 2.

Returns:

Depending on deriv value:

  • deriv=0: rho as (ngrids,) numpy ndarray,

  • deriv=1: (rho, drho_dr) where drho_dr is (3, ngrids) numpy ndarray,

  • deriv=2: (rho, drho_dr, d2rho_dr2) where d2rho_dr2 is (3, 3, ngrids) numpy ndarray.

Return type:

tuple