grid.periodicgrid module

Grids with periodic boundary conditions.

Algorithm

To understand how the algorithm in get_localgrid works, it may be useful review a few basic crystallographic concepts:

• https://en.wikipedia.org/wiki/Bravais_lattice

• https://en.wikipedia.org/wiki/Reciprocal_lattice

• https://en.wikipedia.org/wiki/Fractional_coordinates

The goal of the algorithm is to find all grid points of a periodically extended grid, which lie in a cutoff sphere. The following ASCII art is an example of such a problem in the 1D case:

lattice:       |          |          |          |          |          |
periodic grid: |          |          |  1   2 3 |          |          |
cutoff sphere:                           (             c             )

The first line represents the crystal planes of the lattice. The second line are grid points in one primitive cell, with the lattice repeated for clarity. The last line is the cutoff “sphere” with center c and the extent of the sphere is shown by the parenthesis.

The conventional approach to build the local grid considers all the relevant periodic images of the grid points and retains only those that lie within the cutoff sphere:

lattice:       |          |          |          |          |          |
periodic grid: |          |          |  1   2 3 |          |          |
image -2:      |  1   2 3 |          |          |          |          |
image -1:      |          |  1   2 3 |          |          |          |
image  0:      |          |          |  1   2 3 |          |          |
image  1:      |          |          |          |  1   2 3 |          |
image  2:      |          |          |          |          |  1   2 3 |
cutoff sphere:                           (             c             )

For each periodic image, one can apply a nearest-neighbor lookup within the sphere to find the relevant grid points. In this module, the cKDTree from SciPy is used for this purpose. See https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.KDTree.html

Obviously, it is not necessary to consider an infinite number of periodic images. In this example, only images 0, 1 and 2 need to be constructed, because only these contain grid points that lie in the cutoff sphere.

After retaining only the grid points of the periodic images within the cutoff sphere, one obtains a local grid suitable for integration or evaluation of a function whose domain is limited (approximately) to the cutoff sphere:

lattice:       |          |          |          |          |          |
periodic grid: |          |          |  1   2 3 |          |          |
local grid:                                 2 3    1   2 3    1   2 3
cutoff sphere:                           (             c             )

The algorithm in this module is slightly different. In order to build the cKDTree instance only once for a fixed set of points, periodic images of the cutoff sphere are constructed:

lattice:       |          |          |          |          |          |
periodic grid: |          |          |  1   2 3 |          |          |
c.s. image  0:                           (             c             )
c.s. image -1:                (             c             )
c.s. image -2:     (             c             )

In this example, only images 0, -1 and -2 of the cutoff sphere result in non-zero “untranslated” grid points that lie within the translated sphere. For each image, the found grid points in the periodic grid need to be translated back, by the negative displacement of the image of the cutoff sphere, to find the grid points in the final local grid. This translation is carried out after the neighbor search and requires no change to the cKDTree object.

To implement the second algorithm, one needs to find all possibly relevant displacements of the cutoff sphere. The left-most position of the center of the cutoff sphere still containing a grid point, without worrying about periodic boundary conditions, is:

\[c_\text{left} = \min_i x_i - r\]

where \(x_i\) is the position of grid point \(i\) and \(r\) is the radius of the sphere. We only need to consider translations by multiples of the lattice vector:

\[c_\text{left} = c + a \left\lceil \frac{\min_i x_i - c - r}{a} \right\rceil\]

where \(a\) is the length of the 1D lattice vector (and also the spacing between adjacent crystal planes). Similarly, for the right-most center, we have:

\[c_\text{right} = c + a \left\lfloor \frac{\max_i x_i - c + r}{a} \right\rfloor\]

The positions of all relevant centers are given by:

\[c_j = c + a j \quad \forall \, j \in \mathbb{Z} \quad \text{with} \quad \left\lceil \frac{\min_i x_i - c - r}{a} \right\rceil \le j \le \left\lfloor \frac{\max_i x_i - c + r}{a} \right\rfloor\]

Note that this algorithm also works when the points in the periodic grid are not confined to one primitive cell (as long as no duplicates are present). This situation will just result in more images of centers to be considered. While this will just work, it implies more iterations at the Python level, which could be potentially slow. To avoid such efficiency issues, the wrap option can be set to True when creating a PeriodicGrid instance.

The generalization to higher-dimensional periodic grids is implemented with a combination (a superposition) of all possible displacements of the center along each lattice vector:

\[\vec{c}_{j_1, \ldots, j_K} = \vec{c} + \sum_{k=1}^K \vec{a}_k j \quad \forall \, j \in \mathbb{Z} \quad \text{with} \quad \left\lceil \min_i \tilde{x}_i - \tild{c} - \frac{r}{s_k} \right\rceil \le j \le \left\lfloor \max_i \tilde{x}_i - \tild{c} + \frac{r}{s_k} \right\rfloor\]

where \(\vec{a}_k\) is lattice vector \(k\), \(\tilde{x}_i\) is the fractional coordinate of grid point \(i\), \(\tilde{c}\) is the fractional coordinate of the center and \(s_k\) is the spacing between adjacent crystal planes along the :math:`k`th lattice vector.

class PeriodicGrid(points, weights, realvecs=None, wrap=False)

Bases: Grid

Grid with support for periodic boundary conditions.

The purpose of this class is to support certain operations on integration grids with periodic boundary conditions. It does not construct such grids and it assumes the grid points and weights (and also the lattice vectors) of a single unit cell are provided when creating an instance of the class.

The dimensionality of the grid points is not fixed but the most common use cases are 1D, 2D, 3D grids. The number of cell vectors can be anything from zero up to the dimensionality of the grid points. The cell vectors cannot form a singular matrix.

The whole-grid integration works in exactly the same way as in the base class. Through the get_localgrid method, the following two operations can be carried out easily:

1) The integration of a local function on a part of the periodic grid. The integrand may cover several unit cells and may contain contributions from a periodic function, as long as it decays sufficiently fast to zero near the cutoff radius.

2) The addition of the periodic repetition of a local function to the whole grid.

In both cases, one uses get_localgrid to construct an aperiodic grid, whose grid points lie within a given (hyper)sphere enclosing the local function of interest. This grid is a periodic repetition of the single-unit-cell grid of which only the points within the (hyper)sphere are retained. On this aperiodic grid, one can carry out operations as usual. The attribute localgrid.indices can be used for two purposes:

1) To transfer a periodic function on the parent grid to the localgrid, one simply uses NumPy slicing: fnlocal = fnperiodic[localgrid.indices]

2) To add a periodic repetition of a local function to the parent grid, one uses np.add.att(fnperiodic, localgrid.indices, fnlocal). One should not use fnperiodic[localgrid.indices] += fnlocal, because this will give wrong results when localgrid.indices contains the same index multiple times, which happens when the cutoff (hyper)sphere covers multiple primitive cells. See https://docs.scipy.org/doc/numpy/reference/generated/numpy.ufunc.at.html

The efficiency of the get_localgrid method will be optimal when all grid points fall within one primitive cell, e.g. the fractional coordinates should lie in [0, 1[ or [-0.5, 0.5[. The algorithm will also work when the grid points of the primitive cell accidentally overflow into neighboring images, but in this case it becomes less efficient.

__init__(points, weights, realvecs=None, wrap=False)

Construct a PeriodicGrid instance.

Parameters

• points (np.ndarray(N,) or np.ndarray(N, M)) – An array with N positions of M-dimensional grid points.

• weights (np.ndarray(N,)) – An array of weights associated with each point on the grid.

• realvecs (np.ndarray(1,) or np.ndarray(Nlv, M,), optional) – Nlv Lattice vectors in real space defining the periodic boundary conditions. Each row contains one vector. These are sometimes also called period vectors or cell vectors. When not given, this class behaves identically to the Grid base class.

• wrap (boolean) – When True, the points are wrapped back into the primitive cell, by adding the appropriate integer linear combination of real-space lattice vectors to each point. (This correction is stored in a copy of the points array provided by the user.) As a result, all points will have fractional coordinates in the interval [0, 1[. When False, a warning is raised when the fractional coordinates span an interval wider than 1.1, because this implies a degradation of efficiency of get_localgrid, which can be easily avoided.

property frac_intvls

Intervals containing the grid points in fractional coordinates.

Type

np.ndarray(Ncv, 2)

get_localgrid(center, radius)

Create a non-peroidic local grid within the given distance from the 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.

Returns

Instance of LocalGrid.

Return type

LocalGrid

integrate(*value_arrays)

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

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

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 m`th 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 points

Positions of the grid points.

Type

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

property realvecs

Real-space lattice vectors.

Type

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

property recivecs

Reciprocal-space lattice vectors.

Type

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

save(filename)

Save the points and weights as a npz file.

Parameters

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

property size

the total number of points on the grid.

Type

int

property spacings

Spacings between adjacent primitive crystal planes.

Type

np.ndarray(Ncv, 2)

property weights

the weights of each grid point.

Type

np.ndarray(N,)

exception PeriodicGridWarning

Bases: Warning

Raised when the fractional coordinates span an interval wider than 1.1.

__init__(*args, **kwargs)

args

with_traceback()

Exception.with_traceback(tb) – set self.__traceback__ to tb and return self.