# grid.cubic module

Hyper Rectangular Grid In Either Two or Three Dimensions.

class Tensor1DGrids(oned_x, oned_y, oned_z=None)[source]

Bases: _HyperRectangleGrid

Tensor product of two/three one-dimensional grids.

__init__(oned_x, oned_y, oned_z=None)[source]

Construct Tensor1DGrids by tensor product of two (or three) one-dimensional grids.

Parameters
• oned_x (OneDGrid) – One-dimensional grid representing the grids along x-axis.

• oned_y (OneDGrid) – One-dimensional grid representing the grids along y-axis.

• oned_z (OneDGrid, optional) – One-dimensional grid representing the grids along z-axis.

coordinates_to_index(indices)[source]

Convert (i, j) or (i, j, k) integer coordinates to the grid point index.

Assumes the grid is ordered moving in the last-coordinate (with other coordinates fixed, e.g. k) followed by the next coordinate to the left (e.g. j), continuing to the first coordinate (e.g. i) (i.e., lexicographical ordering).

Parameters

indices ((int, int) or (int, int, int)) – The $$i-th$$, $$j-th$$, (or $$k-th$$) positions of the grid point.

Returns

index – Index of the grid point.

Return type

int

Create a grid contain points within the given radius of center.

Parameters
• center (float or np.array(M,)) – Cartesian coordinates of the center of the local grid.

• radius (float) – Radius of sphere around the center. When equal to np.inf, the local grid coincides with the whole grid, which can be useful for debugging.

Returns

Instance of LocalGrid.

Return type

LocalGrid

get_points_along_axes()[source]

Return the points along each axes.

Returns

The points in the x, y, if in three dimensions, z axes respectively, where $$M_i$$ is the number of points in the ith-direction.

Return type

np.ndarray(M_x,), np.ndarray(M_y,) np.ndarray(M_z)

index_to_coordinates(index)[source]

Convert grid point index to its (i, j) or (i, j, k) integer coordinates in the grid.

For 3D grid it has a shape of $$(N_x, N_y, N_z)$$ denoting the number of points in $$x$$, $$y$$, and $$z$$ directions. So, each grid point has a $$(i, j, k)$$ integer coordinate where $$0 <= i <= N_x - 1$$, $$0 <= j <= N_y - 1$$, and $$0 <= k <= N_z - 1$$. Two-dimensional case similarly follows. Assumes the grid enumerates in the last coordinate first (with others fixed), following the next coordinate to the left (i.e., lexicographical ordering).

Parameters

index (int) – Index of the grid point.

Returns

indices – The corresponding $$(i, j)$$ or $$(i, j, k)$$ integer coordinates in the grid.

Return type

(int, int) or (int, int, int)

integrate(*value_arrays)[source]

Integrate over the whole grid for given multiple value arrays.

Product of all value_arrays will be computed element-wise then integrated on the grid with its weights:

$\int w(x) \prod_i f_i(x) dx.$
Parameters

*value_arrays (np.ndarray(N, )) – One or multiple value array to integrate.

Returns

The calculated integral over given integrand or function

Return type

float

interpolate(points, values, use_log=False, nu_x=0, nu_y=0, nu_z=0, method='cubic')[source]

Interpolate function and its derivatives on cubic grid.

Only implemented in three-dimensions.

Parameters
• points (np.ndarray, shape (M, 3)) – The 3D Cartesian coordinates of $$M$$ points in $$\mathbb{R}^3$$ for which the interpolant (i.e., the interpolated function) is evaluated. If method=”cubic”, then M must be equal to one.

• values (np.ndarray, shape (N,)) – Function values at each of the $$N$$ grid points.

• use_log (bool, optional) – If True, the logarithm of the function values are interpolated. Can only be used for interpolating derivatives when the derivative is not a mixed derivative.

• nu_x (int, optional) – If zero, then the function in x-direction is interpolated. If greater than zero, then the “nu_x”th-order derivative in the x-direction is interpolated.

• nu_y (int, optional) – If zero, then the function in y-direction is interpolated. If greater than zero, then the “nu_y”th-order derivative in the y-direction is interpolated.

• nu_z (int, optional) – If zero, then the function in z-direction is interpolated. If greater than zero, then the “nu_z”th-order derivative in the z-direction is interpolated.

• method (str, optional) – The method of interpolation to perform. Supported are “cubic” (most accurate but computationally expensive), “linear”, or “nearest” (least accurate but cheap computationally). The last two methods use SciPy’s RegularGridInterpolator function.

Returns

The interpolation of a function (or of it’s derivatives) at a $$M$$ point.

Return type

float

moments(orders, centers, func_vals, type_mom='cartesian', return_orders=False)[source]

Compute the multipole moment integral of a function over centers.

The Cartesian type moments are:

$m_{n_x, n_y, n_z} = \int (x - X_c)^{n_x} (y - Y_c)^{n_y} (z - Z_c)^{n_z} f(r) dr,$

where $$\textbf{R}_c = (X_c, Y_c, Z_c)$$ is the center of the moment, $$\f(r)$$ is the density, and $$(n_x, n_y, n_z)$$ are the Cartesian orders.

The spherical/pure moments with $$(l, m)$$ parameter are:

$m_{lm} = \int | \textbf{r} - \textbf{R}_c|^l S_l^m(\theta, \phi) f(\textbf{r}) d\textbf{r},$

where $$S_l^m$$ is a regular, real solid harmonic.

The radial moments with $$n$$ parameter are:

$m_n = \int | \textbf{r} - \textbf{R}_c|^{n} f(\textbf{r}) d\textbf{r}$

The radial combined with spherical/pure moments $$(n, l, m)$$ are:

$m_{nlm} = \int | \textbf{r} - \textbf{R}_c|^{n+1} S_l^m(\theta, \phi) f(\textbf{r}) d\textbf{r}$
Parameters
• orders (int) – Generates all orders with Horton order depending on the type of the multipole moment type_mom.

• centers (ndarray(M, 3)) – The centers $$\textbf{R}_c$$ of the moments to compute from.

• func_vals (ndarray(N,)) – The function $$f$$ values evaluated on all $$N$$ points on the integration grid.

• type_mom (str) – The type of multipole moments: “cartesian”, “pure”, “radial” and “pure-radial”.

• return_orders (bool) – If true, it will also return a list of size $$L$$ of the orders corresponding to each integral/row of the output.

Returns

Computes the moment integral of the function on the mth center for all orders. If return_orders is true, then this also returns a list that describes what each row/order is, e.g. for Cartesian, [(0, 0, 0), (1, 0, 0) ,…].

Return type

ndarray(L, M), or (ndarray(L, M), list)

property ndim

Return the dimension of the grid.

property origin

Cartesian coordinates of the grid origin.

property points

Positions of the grid points.

Type

np.ndarray(N,) or np.ndarray(N, M)

save(filename)[source]

Save tensor product of three one-dimensional grids attributes as a npz file.

Parameters

filename (str) – The path/name of the .npz file.

property shape

Return number of grid points in the $$(x, y)$$ or $$(x, y, z)$$ directions.

property size

the total number of points on the grid.

Type

int

property weights

the weights of each grid point.

Type

np.ndarray(N,)

class UniformGrid(origin, axes, shape, weight='Trapezoid')[source]

Bases: _HyperRectangleGrid

Uniform grid (a.k.a. rectilinear grid) with evenly-spaced points in each axes.

__init__(origin, axes, shape, weight='Trapezoid')[source]

Construct the UniformGrid object in either two or three dimensions.

Three dimensions is presented below, but similarly will work in two-dimensions. Grid whose points in each (x, y, z) direction has a constant step-size/evenly spaced. Given a origin $$\mathbf{o} = (o_x, o_y, o_z)$$ and three directions forming the axes $$\mathbf{a_1}, \mathbf{a_2}, \mathbf{a_3}$$ with shape $$(M_x, M_y, M_z)$$, then the $$(i, j, k)-\text{th}$$ point of the grid are:

$\begin{split}x_i &= o_x + i \mathbf{a_1} \quad 0 \leq i \leq M_x \\ y_i &= o_y + j \mathbf{a_2} \quad 0 \leq j \leq M_y \\ z_i &= o_z + k \mathbf{a_3} \quad 0 \leq k \leq M_z\end{split}$

The grid enumerates through the z-axis first, then y-axis then x-axis (i.e., lexicographical ordering).

Parameters
• origin (np.ndarray, shape (D,)) – Cartesian coordinates of the grid origin in either two or three dimensions $$D$$.

• axes (np.ndarray, shape (D, D)) – The $$D$$ vectors, stored as rows of axes array, whose rows define the Cartesian coordinate system used to build the cubic grid, i.e. the directions of the “(x,y)(x,y,z)”-axis whose norm tells us the distance between points in that direction.

• shape (np.ndarray, shape (D,)) – Number of grid points along each axis.

• weight (str) –

String indicating weighting function. This can be:

Rectangle :

The weights are the standard Riemannian weights,

\begin{split}\begin{align*} w_{ij} &= \frac{V}{M_x \cdot M_y} \tag{Two-Dimensions} \\ w_{ijk} &= \frac{V}{M_x\cdot M_y \cdot M_z} \tag{Three-Dimensions} \end{align*}\end{split}

where $$V$$ is the volume or area of the uniform grid.

Trapezoid :

Equivalent to rectangle rule with the assumption function is zero on the boundaries.

\begin{split}\begin{align*} w_{ij} &= \frac{V}{(M_x + 1) \cdot (M_y + 1)} \tag{Two-Dimensions} \\ w_{ijk} &= \frac{V}{(M_x + 1) \cdot (M_y + 1) \cdot (M_z + 1)} \tag{Three-Dimensions} \end{align*}\end{split}

where $$V$$ is the volume or area of the uniform grid.

Fourier1 :

Assumes function can be expanded in a Fourier series, and then use Gaussian quadrature. Assumes the function is zero at the boundary of the cube. In three-dimensions it is

$w_{ijk} = \frac{2^3}{(M_x + 1) \cdot (M_y + 1) \cdot (M_z + 1)} \bigg[ \bigg(\sum_{p=1}^{M_x} \frac{\sin(ip \pi/(M_x + 1)) (1 - \cos(p\pi)}{p\pi} \bigg) \bigg(\sum_{p=1}^{M_y} \frac{\sin(jp \pi/(M_y + 1)) (1 - \cos(p\pi)}{p\pi} \bigg) \bigg(\sum_{p=1}^{M_z} \frac{\sin(kp \pi/(M_z + 1)) (1 - \cos(p\pi)}{p\pi} \bigg) \bigg]$

whereas in two-dimensions it is

$w_{ij} = \frac{2^3}{(M_x + 1) \cdot (M_y + 1)} \bigg[ \bigg(\sum_{p=1}^{M_x} \frac{\sin(ip \pi/(M_x + 1)) (1 - \cos(p\pi)}{p\pi} \bigg) \bigg(\sum_{p=1}^{M_y} \frac{\sin(jp \pi/(M_y + 1)) (1 - \cos(p\pi)}{p\pi} \bigg) \bigg]$
Fourier2 :

Alternative weights based on Fourier series. Assumes the function is zero at the boundary of the cube.

$\begin{split}w_{ijk} &= V^\prime \cdot w_i w_j w_k, \\ w_i &= \bigg(\frac{2\sin((j - 0.5)\pi) \sin^2(M_x\pi/2)}{M_x^2 \pi} + \frac{4}{M_x \pi} \sum_{p=1}^{M_x - 1} \frac{\sin((2j-1)p\pi /n_x) sin^2(p \pi)}{\pi} \bigg)\end{split}$
Alternative :

This does not assume function is zero at the boundary.

\begin{split}\begin{align*} w_{ij} &= V \cdot \frac{M_x - 1}{M_x} \frac{M_y - 1}{M_y} \tag{Two-Dimensions}\\ w_{ijk} &= V \cdot \frac{M_x - 1}{M_x} \frac{M_y - 1}{M_y} \frac{M_z - 1}{M_z} \tag{Three-Dimensions} \end{align*}\end{split}

property axes

Return the axes of the uniform grid.

closest_point(point, which='closest')[source]

Identify the index of the closest grid point to a given point.

Imagine a point inside a small sub-cube. If closest is selected, it will pick the corner in the sub-cube that is closest to that point. if origin is selected, it will pick the corner that is the bottom, left-most, down-most in the sub-cube.

Parameters
• point (np.ndarray, shape (3,)) – Point in $$[-1,1]^3$$.

• which (str) – If “closest”, returns the closest index of the grid point. If “origin”, return the left-most, down-most closest index of the grid point.

Returns

index – Index of the point in points closest to the grid point.

Return type

int

coordinates_to_index(indices)[source]

Convert (i, j) or (i, j, k) integer coordinates to the grid point index.

Assumes the grid is ordered moving in the last-coordinate (with other coordinates fixed, e.g. k) followed by the next coordinate to the left (e.g. j), continuing to the first coordinate (e.g. i) (i.e., lexicographical ordering).

Parameters

indices ((int, int) or (int, int, int)) – The $$i-th$$, $$j-th$$, (or $$k-th$$) positions of the grid point.

Returns

index – Index of the grid point.

Return type

int

classmethod from_cube(fname, weight='Trapezoid', return_data=False)[source]

Initialize UniformGrid class based on the grid specifications of a cube file.

Parameters
• fname (str) – Cube file name with *.cube extension.

• weight (str) – Scheme for computing the weights of the grid. See the acceptable values in the __init__() method.

• return_data (bool) –

If False, only the grid is returned. If True a tuple with the grid and the cube data is returned. The cube data is a dictionary with the following keys:

• atnums: atomic numbers of the atoms in the molecule.

• atcorenums: Pseudo-number of $$M$$ atoms in the molecule.

• atcoords: Cartesian coordinates of $$M$$ atoms in the molecule.

classmethod from_molecule(atcorenums, atcoords, spacing=0.2, extension=5.0, rotate=True, weight='Trapezoid')[source]

Construct a uniform grid given the molecular pseudo-numbers and coordinates.

Parameters
• atcorenums (np.ndarray, shape (M,)) – Pseudo-number of $$M$$ atoms in the molecule.

• atcoords (np.ndarray, shape (M, 3)) – Cartesian coordinates of $$M$$ atoms in the molecule.

• spacing (float, optional) – Increment between grid points along $$x$$, $$y$$, and $$z$$ direction.

• extension (float, optional) – The extension of the length of the cube on each side of the molecule.

• rotate (bool, optional) – When True, the molecule is rotated so the axes of the cube file are aligned with the principle axes of rotation of the molecule. If False, generates axes based on the x,y,z-axis and the spacing parameter, and the origin is defined by the maximum/minimum of the atomic coordinates.

• weight (str, optional) –

String indicating weighting function. Denoting the volume/area of the uniform grid by $$V$$, the weighting function can be:

Rectangle :

The weights are the standard Riemannian weights,

$w_{ijk} = \frac{V}{M_x\cdot M_y \cdot M_z}$
Trapezoid :

Equivalent to rectangle rule with the assumption function is zero on the boundaries.

$w_{ijk} = \frac{V}{(M_x + 1) \cdot (M_y + 1) \cdot (M_z + 1)}$
Fourier1 :

Assumes function can be expanded in a Fourier series, and then use Gaussian quadrature. Assumes the function is zero at the boundary of the cube.

$w_{ijk} = \frac{8}{(M_x + 1) \cdot (M_y + 1) \cdot (M_z + 1)} \bigg[ \bigg(\sum_{p=1}^{M_x} \frac{\sin(ip \pi/(M_x + 1)) (1 - \cos(p\pi)}{p\pi} \bigg) \bigg(\sum_{p=1}^{M_y} \frac{\sin(jp \pi/(M_y + 1)) (1 - \cos(p\pi)}{p\pi} \bigg) \bigg(\sum_{p=1}^{M_z} \frac{\sin(kp \pi/(M_z + 1)) (1 - \cos(p\pi)}{p\pi} \bigg) \bigg]$
Fourier2 :

Alternative weights based on Fourier series. Assumes the function is zero at the boundary of the cube.

$\begin{split}w_{ijk} &= V^\prime \cdot w_i w_j w_k,\\ w_i &= \bigg(\frac{2\sin((j - 0.5)\pi) \sin^2(M_x\pi/2)}{M_x^2 \pi} + \frac{4}{M_x \pi} \sum_{p=1}^{M_x - 1} \frac{\sin((2j-1)p\pi /n_x) sin^2(p \pi)}{\pi} \bigg)\end{split}$
Alternative :

This does not assume function is zero at the boundary.

$w_{ijk} = V \cdot \frac{M_x - 1}{M_x} \frac{M_y - 1}{M_y} \frac{M_z - 1}{M_z}$

generate_cube(fname, data, atcoords, atnums, pseudo_numbers=None)[source]

Write the data evaluated on grid points into a cube file.

Parameters
• fname (str) – Cube file name with *.cube extension.

• data (np.ndarray, shape=(npoints,)) – An array containing the evaluated scalar property on the grid points.

• atcoords (np.ndarray, shape (M, 3)) – Cartesian coordinates of $$M$$ atoms in the molecule.

• atnums (np.ndarray, shape (M,)) – Atomic numbers of $$M$$ atoms in the molecule.

• pseudo_numbers (np.ndarray, shape (M,), optional) – Pseudo-numbers (core charges) of $$M$$ atoms in the molecule.

Create a grid contain points within the given radius of center.

Parameters
• center (float or np.array(M,)) – Cartesian coordinates of the center of the local grid.

• radius (float) – Radius of sphere around the center. When equal to np.inf, the local grid coincides with the whole grid, which can be useful for debugging.

Returns

Instance of LocalGrid.

Return type

LocalGrid

get_points_along_axes()[source]

Return the points along each axes.

Returns

The points in the x, y, if in three dimensions, z axes respectively, where $$M_i$$ is the number of points in the ith-direction.

Return type

np.ndarray(M_x,), np.ndarray(M_y,) np.ndarray(M_z)

index_to_coordinates(index)[source]

Convert grid point index to its (i, j) or (i, j, k) integer coordinates in the grid.

For 3D grid it has a shape of $$(N_x, N_y, N_z)$$ denoting the number of points in $$x$$, $$y$$, and $$z$$ directions. So, each grid point has a $$(i, j, k)$$ integer coordinate where $$0 <= i <= N_x - 1$$, $$0 <= j <= N_y - 1$$, and $$0 <= k <= N_z - 1$$. Two-dimensional case similarly follows. Assumes the grid enumerates in the last coordinate first (with others fixed), following the next coordinate to the left (i.e., lexicographical ordering).

Parameters

index (int) – Index of the grid point.

Returns

indices – The corresponding $$(i, j)$$ or $$(i, j, k)$$ integer coordinates in the grid.

Return type

(int, int) or (int, int, int)

integrate(*value_arrays)[source]

Integrate over the whole grid for given multiple value arrays.

Product of all value_arrays will be computed element-wise then integrated on the grid with its weights:

$\int w(x) \prod_i f_i(x) dx.$
Parameters

*value_arrays (np.ndarray(N, )) – One or multiple value array to integrate.

Returns

The calculated integral over given integrand or function

Return type

float

interpolate(points, values, use_log=False, nu_x=0, nu_y=0, nu_z=0, method='cubic')[source]

Interpolate function and its derivatives on cubic grid.

Only implemented in three-dimensions.

Parameters
• points (np.ndarray, shape (M, 3)) – The 3D Cartesian coordinates of $$M$$ points in $$\mathbb{R}^3$$ for which the interpolant (i.e., the interpolated function) is evaluated. If method=”cubic”, then M must be equal to one.

• values (np.ndarray, shape (N,)) – Function values at each of the $$N$$ grid points.

• use_log (bool, optional) – If True, the logarithm of the function values are interpolated. Can only be used for interpolating derivatives when the derivative is not a mixed derivative.

• nu_x (int, optional) – If zero, then the function in x-direction is interpolated. If greater than zero, then the “nu_x”th-order derivative in the x-direction is interpolated.

• nu_y (int, optional) – If zero, then the function in y-direction is interpolated. If greater than zero, then the “nu_y”th-order derivative in the y-direction is interpolated.

• nu_z (int, optional) – If zero, then the function in z-direction is interpolated. If greater than zero, then the “nu_z”th-order derivative in the z-direction is interpolated.

• method (str, optional) – The method of interpolation to perform. Supported are “cubic” (most accurate but computationally expensive), “linear”, or “nearest” (least accurate but cheap computationally). The last two methods use SciPy’s RegularGridInterpolator function.

Returns

The interpolation of a function (or of it’s derivatives) at a $$M$$ point.

Return type

float

moments(orders, centers, func_vals, type_mom='cartesian', return_orders=False)[source]

Compute the multipole moment integral of a function over centers.

The Cartesian type moments are:

$m_{n_x, n_y, n_z} = \int (x - X_c)^{n_x} (y - Y_c)^{n_y} (z - Z_c)^{n_z} f(r) dr,$

where $$\textbf{R}_c = (X_c, Y_c, Z_c)$$ is the center of the moment, $$\f(r)$$ is the density, and $$(n_x, n_y, n_z)$$ are the Cartesian orders.

The spherical/pure moments with $$(l, m)$$ parameter are:

$m_{lm} = \int | \textbf{r} - \textbf{R}_c|^l S_l^m(\theta, \phi) f(\textbf{r}) d\textbf{r},$

where $$S_l^m$$ is a regular, real solid harmonic.

The radial moments with $$n$$ parameter are:

$m_n = \int | \textbf{r} - \textbf{R}_c|^{n} f(\textbf{r}) d\textbf{r}$

The radial combined with spherical/pure moments $$(n, l, m)$$ are:

$m_{nlm} = \int | \textbf{r} - \textbf{R}_c|^{n+1} S_l^m(\theta, \phi) f(\textbf{r}) d\textbf{r}$
Parameters
• orders (int) – Generates all orders with Horton order depending on the type of the multipole moment type_mom.

• centers (ndarray(M, 3)) – The centers $$\textbf{R}_c$$ of the moments to compute from.

• func_vals (ndarray(N,)) – The function $$f$$ values evaluated on all $$N$$ points on the integration grid.

• type_mom (str) – The type of multipole moments: “cartesian”, “pure”, “radial” and “pure-radial”.

• return_orders (bool) – If true, it will also return a list of size $$L$$ of the orders corresponding to each integral/row of the output.

Returns

Computes the moment integral of the function on the mth center for all orders. If return_orders is true, then this also returns a list that describes what each row/order is, e.g. for Cartesian, [(0, 0, 0), (1, 0, 0) ,…].

Return type

ndarray(L, M), or (ndarray(L, M), list)

property ndim

Return the dimension of the grid.

property origin

Return the Cartesian coordinates of the uniform grid origin.

property points

Positions of the grid points.

Type

np.ndarray(N,) or np.ndarray(N, M)

save(filename)[source]

Save uniform cubic grid attributes as a npz file.

Parameters

filename (str) – The path/name of the .npz file.

property shape

Return number of grid points in the $$(x, y)$$ or $$(x, y, z)$$ directions.

property size

the total number of points on the grid.

Type

int

property weights

the weights of each grid point.

Type

np.ndarray(N,)