Cubic Grids

Cubic grids are commonly found in many other packages and are useful for visualization, interpolation and integration purposes.

In Grid, it offers two classes for constructing these cubic (or Hyper-rectangular) grids in three-dimensions:

• Tensor1DGrids : Tensor product of three one-dimensional grids.

• UniformGrid : Evenly spaced grid in three axes.

See their API documentation for more information about the classes. This example illustrates how to use these grid classes for both visualization, interpolation (including differentiation), and integration, with a focus on usage in quantum chemistry. It is recommended for general usage to use the Tensor1DGrids and for quantum chemistry applications, to use the UniformGrid class.

Tensor 1D Grids

The first grid is a tensor combination of three grids with points \(\{p^1_i\}, \{p^2_i\}, \{p^3_i\}\) and weights \(\{w^1_i\}, \{w^2_i\}, \{w^3_i\}\) such that the new set of points are \(\{(p^1_i, p^2_j, p^3_k)\}\) with weights \(\{w^1_i \times w^2_j \times w^3_k\}\).

%matplotlib inline
from grid.onedgrid import UniformInteger
from grid.rtransform import LinearInfiniteRTransform
from grid.cubic import Tensor1DGrids
import numpy as np

limit = 10
# Construct a grid between -10 and 10 on the real line.
npoints = 50
oned_gridx = UniformInteger(npoints)
oned_gridx = LinearInfiniteRTransform(-limit, limit).transform_1d_grid(oned_gridx)

# Construct another grid between -0.5 and 0.5 on the real line
npoints = 45
oned_gridy = UniformInteger(npoints)
oned_gridy = LinearInfiniteRTransform(-limit, limit).transform_1d_grid(oned_gridy)

# Construct another grid between -10 and 10 on the real line
npoints = 55
oned_gridz = UniformInteger(npoints)
oned_gridz = LinearInfiniteRTransform(-limit, limit).transform_1d_grid(oned_gridz)

# Constructing the tensor product grid.
tensor_grid = Tensor1DGrids(oned_gridx, oned_gridy, oned_gridz)

Each of the individual grids axes are plotted.

import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.scatter(oned_gridx.points, [0] *  len(oned_gridx.points), [0] *  len(oned_gridx.points), color="k", label="Points on x-axis")
ax.scatter( [0] *  len(oned_gridy.points), oned_gridy.points, [0] *  len(oned_gridy.points), color="m", label="Points on y-axis")
ax.scatter( [0] *  len(oned_gridz.points), [0] *  len(oned_gridz.points), oned_gridz.points , color="b", label="Points on z-axis")
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")
plt.legend()
plt.show()

The Tensor1DGrids can now be easily constructed by passing in each of the three one-dimensional grids. These are then plotted

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.scatter(tensor_grid.points[:, 0], tensor_grid.points[:, 1], tensor_grid.points[:, 2], color="m", alpha=0.01)
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")
plt.title("Tensor Product of Three One-Dimensional Grids")
plt.show()

print(f"Points of the grid: {tensor_grid.points}")
print(f"Weights of the grid: {tensor_grid.weights}.")
print(f"The number of points: {tensor_grid.size}")
print(f"The number of dimensions: {tensor_grid.ndim}.")
print(f"The origin of the grid: {tensor_grid.origin}.")
print(f"The shape of the grid: {tensor_grid.shape}.")

# Conversion from indices i to coordinates (i, j, k) to points
coordinate = tensor_grid.index_to_coordinates(5)
print(f"Coordinate of the fifth point: {coordinate}.")
index = tensor_grid.coordinates_to_index((0, 0, 5))
print(f"Index of the (0, 0, 5) coordinate: {index}.")

Points of the grid: [[-10.         -10.         -10.        ]
 [-10.         -10.          -9.62962963]
 [-10.         -10.          -9.25925926]
 ...
 [ 10.          10.           9.25925926]
 [ 10.          10.           9.62962963]
 [ 10.          10.          10.        ]]
Weights of the grid: [0.06871435 0.06871435 0.06871435 ... 0.06871435 0.06871435 0.06871435].
The number of points: 123750
The number of dimensions: 3.
The origin of the grid: [-10. -10. -10.].
The shape of the grid: (50, 45, 55).
Coordinate of the fifth point: (0, 0, 5).
Index of the (0, 0, 5) coordinate: 5.

Example 1: Interpolation

This example will attempt to interpolate a Gaussian, first on a subset of the points used for interpolation then later on random set of points.

# Interpolation of a Gaussian
gaussian = lambda pts: np.exp(-np.linalg.norm(pts, axis=1)**2.0)
gaus_vals = gaussian(tensor_grid.points)

# Interpolation evaluated on subset of the same points used for interpolation
subset = tensor_grid.points[np.random.choice(range(tensor_grid.size), 10)]
true_vals = gaussian(subset)

print("Interpolate using linear method.")
interpol = tensor_grid.interpolate(subset, gaus_vals, use_log=False, method="linear")
print("The errors are: ", np.abs(interpol - true_vals))
print(f"Maximum error is {np.max(np.abs(interpol - true_vals))} \n")

print("Interpolate using linear method applying the logarithm.")
interpol = tensor_grid.interpolate(subset, gaus_vals, use_log=True, method="linear")
print("The errors are: ", np.abs(interpol - true_vals))
print(f"Maximum error is {np.max(np.abs(interpol - true_vals))} \n")

print("Interpolate using cubic splines")
interpol = tensor_grid.interpolate(subset, gaus_vals, use_log=False, method="cubic")
print("The errors are: ", np.abs(interpol - true_vals))
print(f"Maximum error is {np.max(np.abs(interpol - true_vals))} \n")

print("Inteprolate using cubic splines appling the logarithm")
interpol = tensor_grid.interpolate(subset, gaus_vals, use_log=True, method="cubic")
print("The errors are: ", np.abs(interpol - true_vals))
print(f"Maximum error is {np.max(np.abs(interpol - true_vals))} \n")

Interpolate using linear method.
The errors are:  [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
Maximum error is 0.0

Interpolate using linear method applying the logarithm.
The errors are:  [171.46716217  70.30215474  78.32720365  32.11924387 130.53013235
  28.19417151  35.37443308  85.77553134  85.69913673  26.65532941]
Maximum error is 171.4671621664289

Interpolate using cubic splines
The errors are:  [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
Maximum error is 0.0

Inteprolate using cubic splines appling the logarithm
The errors are:  [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
Maximum error is 0.0

# Interpolate at random set of points
random_pts = np.random.uniform(-1, 1, (25, 3))
true_vals = gaussian(random_pts)

print("Interpolate using linear method.")
interpol = tensor_grid.interpolate(random_pts, gaus_vals, use_log=False, method="linear")
print("The errors are: ", np.abs(interpol - true_vals))
print(f"Maximum error is {np.max(np.abs(interpol - true_vals))} \n")
# The logarithm is first applied to the Gaussian functions to be easier to interpolate at.
print("Interpolate using linear method applying the logarithm.")
interpol = tensor_grid.interpolate(random_pts, gaus_vals, use_log=True, method="linear")
print("The errors are: ", np.abs(interpol - true_vals))
print(f"Maximum error is {np.max(np.abs(interpol - true_vals))} \n")
# Rather then interpolate using linearly, this used the cubic
print("Interpolate using cubic splines")
interpol = tensor_grid.interpolate(random_pts, gaus_vals, use_log=False, method="cubic")
print("The errors are: ", np.abs(interpol - true_vals))
print(f"Maximum error is {np.max(np.abs(interpol - true_vals))} \n")
# Same as above but the logarithm is applied.
print("Inteprolate using cubic splines appling the logarithm")
interpol = tensor_grid.interpolate(random_pts, gaus_vals, use_log=True, method="cubic")
print("The errors are: ", np.abs(interpol - true_vals))
print(f"Maximum error is {np.max(np.abs(interpol - true_vals))} \n")

Interpolate using linear method.
The errors are:  [0.01233929 0.00172326 0.01414473 0.00302718 0.03765561 0.00854703
 0.00573863 0.0096147  0.06938564 0.00885821 0.03244877 0.00313253
 0.00879784 0.02016751 0.00573694 0.014182   0.00269118 0.01196634
 0.00416348 0.05229362 0.00944945 0.00436619 0.02369888 0.01221546
 0.03681424]
Maximum error is 0.0693856422879573

Interpolate using linear method applying the logarithm.
The errors are:  [1.5443592  1.71594526 1.34273699 1.80571319 1.18542439 1.62956061
 1.63051993 1.45611032 1.11876818 1.40452868 1.24498067 2.56010903
 1.54114751 1.13349723 1.88518539 1.45877593 1.86029613 1.40771558
 2.29762085 1.07745488 1.60149404 2.39484291 1.27974207 1.19860893
 1.19903943]
Maximum error is 2.560109034575345

Interpolate using cubic splines
The errors are:  [3.86592728e-04 4.22040488e-04 1.66738913e-04 4.13390441e-04
 1.01639001e-03 4.34750401e-05 5.16147392e-05 9.61299995e-05
 1.64110278e-03 5.10977886e-04 1.88314046e-04 1.11520856e-04
 5.16837466e-04 8.08227555e-04 1.77906599e-04 2.64324388e-05
 1.57194487e-04 3.27801023e-04 2.26379162e-04 1.25083661e-03
 2.33331268e-04 4.52028716e-04 3.51733200e-04 8.65339717e-04
 8.68365576e-04]
Maximum error is 0.0016411027780264265

Inteprolate using cubic splines appling the logarithm
The errors are:  [5.55111512e-17 5.55111512e-17 1.66533454e-16 1.11022302e-16
 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.11022302e-16
 1.11022302e-16 0.00000000e+00 0.00000000e+00 4.16333634e-17
 1.11022302e-16 0.00000000e+00 8.32667268e-17 1.66533454e-16
 5.55111512e-17 1.11022302e-16 5.55111512e-17 1.11022302e-16
 1.11022302e-16 1.11022302e-16 1.11022302e-16 1.11022302e-16
 0.00000000e+00]
Maximum error is 1.6653345369377348e-16

# Derivative can be interpolated
# X-derivative
deriv_x = tensor_grid.interpolate(random_pts, gaus_vals, use_log=False, nu_x=1, nu_y=0, nu_z=0, method="cubic")
print(f"Derivative in x-direction {deriv_x}.")
# Y-derivative
deriv_y = tensor_grid.interpolate(random_pts, gaus_vals, use_log=False, nu_x=0, nu_y=1, nu_z=0, method="cubic")
# Z-derivative
deriv_z = tensor_grid.interpolate(random_pts, gaus_vals, use_log=False, nu_x=0, nu_y=0, nu_z=1, method="cubic")
# XY- derivative
deriv_xy = tensor_grid.interpolate(random_pts, gaus_vals, use_log=False, nu_x=1, nu_y=1, nu_z=0, method="cubic")
# The logarithm can be used to interpolate derivatives on single points only (mixed derivative are not supported)
derivs = []
for pt in random_pts:
    deriv_x = tensor_grid.interpolate(np.array([pt]), gaus_vals, use_log=True, nu_x=1, nu_y=0, nu_z=0, method="cubic")
    derivs.append(deriv_x[0])
print(f"Derivative (using logarithm) in x-direction {np.array(derivs)}.")

Derivative in x-direction [ 0.30734258  0.24753392 -0.10164385 -0.35502454  0.62256378 -0.4806336
 -0.27490521  0.64421173  0.50708235 -0.62589627 -0.02271022 -0.18075862
 -0.46666906 -0.35434171 -0.07220761  0.37371807 -0.07209604  0.47252084
  0.18407512 -0.08490725  0.31419978 -0.18812879  0.11677048  0.27183554
  0.00227194].
Derivative (using logarithm) in x-direction [ 0.30673882  0.24715189 -0.1054587  -0.35469932  0.62461378 -0.48053782
 -0.27527774  0.64168648  0.50556283 -0.62694789 -0.02390025 -0.17964052
 -0.46912167 -0.35324203 -0.07340006  0.37428065 -0.07338858  0.4747938
  0.18456838 -0.08932495  0.31453276 -0.18754882  0.12127639  0.27184459
  0.00239627].

Example 2: Visualization of \(2p_{z}\) Orbital using Tensor 1D Grid

Given a function \(f(x, y, z)\), (e.g. an atomic orbital or electron density) we can evaluate it on a cubic grid and visualize the distribution of the function in space using isosurface values. Here we define a function that plots every point in the grid with a corresponding absolute value of the function greater than a threshold value. The function is latter used to visualize the \(2p_{z}\) orbital.

i. Utility Function For Plotting Isosurfaces

# Function to plot grid points bigger than a given threshold
def plot_isosurface(
    grid, vals, isovalue, *, both_signs=True, at_data=None, noticks=False, title=None
):
    """Plot grid points with values bigger than isovalue.

    The function plots the grid points with absolute values bigger than isovalue.

    Parameters
    ----------
    grid : Grid
        Cartesian grid (it can be uniform or TensorProductGrid) with 3 dimensions.
    vals : ndarray
        Values of the function on the grid points.
    isovalue : float
        Value of the isosurface to plot.
    both_signs : bool, optional
        If True, plot points with negative values too. The default is True.
    at_data : tuple, optional
        Tuple containing the atomic coordinates and atomic labels. If not None,
        it should be a tuple containing (atcoords, atlabels) where atcoords is
        a ndarray of shape (natoms, 3) containing the atomic coordinates in
        bohr and atlabels is a list of length natoms containing the atomic labels.
    title : str, optional
        Title of the plot. The default is None.
    """
    # indices of points with values close to isovalue
    idx_p_vals = np.where(vals > isovalue)
    idx_n_vals = np.where(vals < -isovalue)

    # points with values close to isovalue
    p_pts = grid.points[idx_p_vals]
    n_pts = grid.points[idx_n_vals]

    xmin, xmax = np.min(grid.points[:, 0]), np.max(grid.points[:, 0])
    ymin, ymax = np.min(grid.points[:, 1]), np.max(grid.points[:, 1])
    zmin, zmax = np.min(grid.points[:, 2]), np.max(grid.points[:, 2])

    fig = plt.figure()
    ax = plt.axes(projection="3d")
    if noticks:
        ax.set_xticks([])
        ax.set_yticks([])
        ax.set_zticks([])

    if at_data is not None:
        if len(at_data) != 2:
            raise ValueError("at_data should be a tuple containing (atcoords, atlabels)")

        atcoords, atlabels = at_data
        if atcoords.shape[0] != len(atlabels):
            raise ValueError("atcoords and atlabels should have the same length")

        # if there are more than 2 atoms set the view point perpendicular to the molecular plane
        if len(atlabels) > 2:
            # Calculate the centroid of atomic coordinates
            centroid = np.mean(atcoords, axis=0)
            # Center the atomic coordinates by subtracting the centroid from each point
            centered_coords = atcoords - centroid
            # Perform Singular Value Decomposition (SVD) on the centered coordinates
            _, _, V = np.linalg.svd(centered_coords)
            # The normal vector the molecular plane is the last column of the V matrix
            normal_vector = V[-1]
            # Set the view point perpendicular to the plane
            ax.view_init(
                elev=np.degrees(np.arcsin(normal_vector[2])),
                azim=np.degrees(np.arctan2(normal_vector[1], normal_vector[0])),
            )

        # plot the atomic centers
        for i in range(len(atlabels)):
            ax.text(
                atcoords[i, 0],
                atcoords[i, 1],
                atcoords[i, 2],
                atlabels[i],
                color="k",
                fontsize=12,
                horizontalalignment="center",
                verticalalignment="center",
            )
    # set axis limits
    ax.set_xlim((xmin, xmax))
    ax.set_ylim((ymin, ymax))
    ax.set_zlim((zmin, zmax))

    # plot points with values bigger than isovalue
    ax.scatter(
        p_pts[:, 0], p_pts[:, 1], p_pts[:, 2], alpha=0.01, color="r", label=f"Positive values"
    )
    if both_signs:
        ax.scatter(
            n_pts[:, 0], n_pts[:, 1], n_pts[:, 2], alpha=0.01, color="b", label=f"Negative values"
        )
    legend = plt.legend(loc="upper right")
    for lh in legend.legendHandles:
        lh.set_alpha(1)
    if title is not None:
        plt.title(title)
    plt.show()

ii. Plot \(2p_{z}\) orbital from isosurface value

import matplotlib.pyplot as plt
import numpy as np

# define 2pz orbital function
def psi_2pz(x, y, z, a_0):
    r = np.sqrt(x**2 + y**2 + z**2)
    return np.exp(-r / (2 * a_0)) * z


# Evaluate the function on the grid.
psi_vals = psi_2pz(*tensor_grid.points.T, a_0=1.0)

point_idx = target_x = 0.15
tolerance = 1e-2

plot_isosurface(tensor_grid, psi_vals, target_x, title="$2p_{z}$ Isosurface (0.15 a.u.)")

Uniform Grid

This is type of cubic grid (a.k.a. rectilinear grid) with evenly-spaced points in each axes. It supports the same properties and functions as Tensor1DGrids.

Example 1: Constructing Uniform Grid Around Formaldehyde Anion

This example illustates the use of Grid with the IOData package to construct a default grid around a molecule. The easiest method to do so is the “from_molecule” method. This example will showcase this for Formaldehyde anion. The .fchk files of which are read using the IOData package and the .fchk files are included in doc/notebooks/ch2o_q-1.fchk.

from iodata import load_one

# Load the wavefunction information about ch2o
mol_anion = load_one("ch2o_q-1.fchk")

# Get the atomic coordinates and numbers
atcoords = mol_anion.atcoords
atnums = mol_anion.atnums
print(f"Atomic Coordinates \n {atcoords}")
print(f"Atomic Numbers: \n {atnums}")

fig = plt.figure()
ax = plt.axes(projection='3d')
ax.scatter(atcoords[:, 0], atcoords[:, 1], atcoords[:, 2], color="k", s=50)
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")
plt.title("Atomic coordinates of Formaldehyde.")
plt.show()

Atomic Coordinates
 [[-1.23259516e-32  1.00530930e-32  1.27229733e+00]
 [ 0.00000000e+00 -1.32254311e-32 -9.95825669e-01]
 [ 0.00000000e+00  1.77311490e+00 -2.10171229e+00]
 [-2.17143949e-16 -1.77311490e+00 -2.10171229e+00]]
Atomic Numbers:
 [8 6 1 1]

from grid.cubic import UniformGrid

# Construct uniform grid
uniform_grid = UniformGrid.from_molecule(atnums, atcoords, spacing=0.1, extension=2.0, rotate=True)


print(f"The number of points: {uniform_grid.size}")
print(f"The number of dimensions: {uniform_grid.ndim}.")
print(f"The shape of the grid: {uniform_grid.shape}.")
print(f"The origin of the grid: {uniform_grid.origin}.")
print(f"The axes of the grid: {uniform_grid.axes}.")


fig = plt.figure()
ax = plt.axes(projection='3d')
ax.scatter(atcoords[:, 0], atcoords[:, 1], atcoords[:, 2], color="k", s=50)
ax.scatter(uniform_grid.points[:, 0], uniform_grid.points[:, 1], uniform_grid.points[:, 2], color="m", alpha=0.01)
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis")
ax.set_zlabel("Z-axis")
plt.title("Atomic coordinates of Formaldehyde.")
plt.show()

The number of points: 224960
The number of dimensions: 3.
The shape of the grid: [74 76 40].
The origin of the grid: [-2.   3.8 -3.7].
The axes of the grid: [[ 0.00000000e+00  0.00000000e+00  1.00000000e-01]
 [ 0.00000000e+00 -1.00000000e-01  0.00000000e+00]
 [ 1.00000000e-01  1.11022302e-17  0.00000000e+00]].

Example 2: Calculate and Visualize Dual Descriptor of Formaldehyde

This example will illusrate on how to use UniformGrid to compute and visualize the dual descriptor of formaldehyde. The dual descriptor is a function of electron density and its derivatives and is defined as:

\[f^{(2)}(r) = \Delta f(r) = (f^{+}(r) - f^{-}(r)\]

where \(f^{+}(r)\) and \(f^{-}(r)\) are the Fukui functions defined as:

\[f^{+}(r) = \left(\frac{\partial \rho(r)}{\partial \rho^2}\right)^{+}_{\upsilon(r)} = \rho_{N+1}(r) - \rho_{N}(r)\]

\[f^{-}(r) = \left(\frac{\partial \rho(r)}{\partial \rho^2}\right)^{-}_{\upsilon(r)} = \rho_{N}(r) - \rho_{N-1}(r)\]

then \(f^{(2)}(r)\) can be written as:

\[f^{(2)}(r) = \rho_{N+1}(r) - \rho_{N-1}(r)\]

i. Define utility function to calculate the dual descriptor

The dual descriptor is evaluated at each grid point as the difference between the electron density of the anion and cation of formaldehyde. The electron density can be evaluated in the grid points using the GBasis package with the wavefunction read, as above.

from gbasis.evals.density import evaluate_density

# define function to compute dual descriptor on a grid
def compute_dual_descriptor(grid, rdm_an, rdm_cat, basis_an, basis_cat):
    """Compute dual descriptor on a grid.

    Parameters
    ----------
    grid : Grid
        Grid object.
    rdm_an : np.ndarray
        1-RDM of the anion.
    rdm_cat : np.ndarray
        1-RDM of the cation.
    basis_an : list
        Basis set of the anion.
    basis_cat : list
        Basis set of the cation.
    """
    # evaluate electron density of cation and anion on the grid
    dens_vals_cat = evaluate_density(rdm_cat, basis_cat, grid.points)
    dens_vals_an = evaluate_density(rdm_an, basis_an, grid.points)

    # compute dual descriptor on the grid
    return dens_vals_an - dens_vals_cat

ii. Calculate and visualize dual descriptor

from iodata import load_one
from gbasis.wrappers import from_iodata

# load formaldehyde cation and anion fchk files
mol_cation = load_one("ch2o_q+1.fchk")

# Construct molecular basis from wave-function information read by IOData
basis_cat = from_iodata(mol_cation)
basis_an = from_iodata(mol_anion)

# get rdms for cation and anion
rdm_cat = mol_cation.one_rdms["scf"]
rdm_an = mol_anion.one_rdms["scf"]

# compute dual descriptor on the grid
dd = compute_dual_descriptor(uniform_grid, rdm_an, rdm_cat, basis_an, basis_cat)

# plot dual descriptor isosurface
mol_data = (atcoords, ["O", "C", "H", "H"])
plot_isosurface(
    uniform_grid, dd, 0.03, at_data=mol_data, title="Dual Descriptor Isosurface", noticks=True
)

print(f"Positive values (red) correspond to nucleophilic regions")
print(f"Negative values (blue) correspond to electrophilic regions")

Positive values (red) correspond to nucleophilic regions
Negative values (blue) correspond to electrophilic regions

Example 3: Integrate The Dual Descriptor Domains

All classes in Grid offers an easy method to integrate. Here we integrate the different lobes of the dual descriptor. This can be used to assess the relevance of electrophilic and nucleophilic regions of the molecule.

i. Define utility function to extract dual descriptor domains

# define utility function for grouping grid points in domains
from itertools import combinations
from scipy.sparse.csgraph import connected_components


# function to get grid points indices for same sign connected points with absolute
# value  greater than a threshold (select domain defined by isovalue)
def get_domains(grid, vals, isovalue):
    """Get grid points indices for domains with same sign connected points.

    The function returns the indices of the grid points that belong to the same domain. A domain
    is defined by a set of connected points with the same sign and absolute value greater than
    isovalue.

    Parameters
    ----------
    grid : Grid
        Cartesian grid (it can be uniform or TensorProductGrid) with 3 dimensions.
    vals : ndarray
        Values of the function on the grid points.
    isovalue : float
        Value of the isosurface used to define the domains.

    Returns
    -------
    tuple
        Tuple containing two lists. The first list contains domains with positive values and the
        second list contains domains with negative values. Each domain is a np.array with the
        indexes of the grid points that belong to the same domain.
    """
    # indices of points with value modulus greater than isovalue
    idx_p_vals = np.where(vals > isovalue)[0]  # positive values
    idx_n_vals = np.where(vals < -isovalue)[0]  # negative values

    # create adjacency matrix selected points with same sign
    p_vals_adj = np.zeros((len(idx_p_vals), len(idx_p_vals)))
    n_vals_adj = np.zeros((len(idx_n_vals), len(idx_n_vals)))

    # try all combinations of selected points with same sign (positive values)
    for i, j in combinations(range(len(idx_p_vals)), 2):
        # get coordinate indices of the points
        i_coord_idx = np.array(grid.index_to_coordinates(idx_p_vals[i]))
        j_coord_idx = np.array(grid.index_to_coordinates(idx_p_vals[j]))
        # two points are adjacent if they differ by at most 1 in each coordinate index
        if np.max(np.abs(i_coord_idx - j_coord_idx)) == 1:
            p_vals_adj[i, j] = 1

    # try all combinations of selected points with same sign (negative values)
    for i, j in combinations(range(len(idx_n_vals)), 2):
        i_coord_idx = np.array(grid.index_to_coordinates(idx_n_vals[i]))
        j_coord_idx = np.array(grid.index_to_coordinates(idx_n_vals[j]))
        # two points are adjacent if they differ by at most 1 in each coordinate
        if np.max(np.abs(i_coord_idx - j_coord_idx)) == 1:
            n_vals_adj[i, j] = 1

    # returns an array of integers, each integer represents domain and the indexes where
    # it appear are the points belonging to that domain
    # example: [0,1,0,0,1,1] -> points 0,2,3 belong to domain 0
    #                        -> points 1,4,5 belong to domain 1
    p_groups = connected_components(p_vals_adj, directed=False)[1]
    n_groups = connected_components(n_vals_adj, directed=False)[1]

    # creates a list with the domains. Each element of the list is a np.array with the indexes
    # of the selected points that belong to the same domain.
    # example: [0,1,0,0,1,1] -> [[0,2,3],[1,4,5]]
    p_domains = list([np.where(p_groups == i)[0] for i in range(max(p_groups) + 1)])
    n_domains = list([np.where(n_groups == i)[0] for i in range(max(n_groups) + 1)])

    # transform indexes of selected points to indexes of all grid points for each domain
    p_domains = [idx_p_vals[domain] for domain in p_domains]
    n_domains = [idx_n_vals[domain] for domain in n_domains]

    return p_domains, n_domains

ii. Select and integrate the positive (electrophilic) domains of the dual descriptor

# get dual descriptor domains, integrate and plot them
p_domains, n_domains = get_domains(uniform_grid, dd, 0.03)

for domain in p_domains:
    domain_vals = dd.copy()
    domain_mask = np.isin(np.arange(len(dd)), domain, invert=True)
    domain_vals[domain_mask] = 0.0

    plot_isosurface(
        uniform_grid,
        domain_vals,
        0.001,
        at_data=mol_data,
        title=f"Domain   Integral: {uniform_grid.integrate(domain_vals):.3f}",
        noticks=True,
    )

iii. Select and integrate the negative (nucleophilic) domains of the dual descriptor

for domain in n_domains:
    domain_vals = dd.copy()
    domain_mask = np.isin(np.arange(len(dd)), domain, invert=True)
    domain_vals[domain_mask] = 0.0

    plot_isosurface(
        uniform_grid,
        domain_vals,
        0.001,
        at_data=mol_data,
        title=f"Domain Integral: {uniform_grid.integrate(domain_vals):.3f}",
        noticks=True,
    )