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:
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:
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:
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:
The positions of all relevant centers are given by:
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:
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)[source]
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 attributelocalgrid.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 usefnperiodic[localgrid.indices] += fnlocal
, because this will give wrong results whenlocalgrid.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.htmlThe 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)[source]
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)[source]
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
- 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
- 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 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)[source]
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,)