Grid

WTP provides non-orthogonal 3D periodic grids that support

1. Three Indexing schemes and iteration.
2. Fourier transform.

The goals is to provide a uniform interface for dealing with

1. The reciprocal lattice.
2. The Brillouin zone.
3. The crystal lattice.
4. The homecell.

A grid is made of a set of basis vectors and the domain in units of the basis vectors. If you subtype this, the concrete type should have these two fields for things to work out of the box.

basis::NTuple{3, Vector3}
domain::NTuple{3, NTuple{2, Integer}}

make_grid(T, basis, domain)

Exmaple:

Make a grid of type T from basis and domain.

To make a home cell with basis vectors

$a_1 = \begin{pmatrix} 0\\ \sqrt{3}/2 \\ 0 \end{pmatrix} \quad a_2 = \begin{pmatrix} 1\\ 1/2 \\ \sqrt{3}/2 \end{pmatrix} \quad a_3 = \begin{pmatrix} 0\\ 0 \\ 1/2 \end{pmatrix}$

and domain $((-2, 1), (-3, 2), (-4, 3))$

julia> homecell_basis = (Vector3(0, sqrt(3)/2, 0), Vector3(1, 1/2, sqrt(3)/2), Vector3(0, 0, 1/2));

julia> homecell = make_grid(HomeCell3D, homecell_basis, ((-2, 1), (-3, 2), (-4, 3)))
type: HomeCell3D
domain: ((-2, 1), (-3, 2), (-4, 3))
basis:
    ket: 0.000, 0.866, 0.000
    ket: 1.000, 0.500, 0.866
    ket: 0.000, 0.000, 0.500

The home cell in 3D. The $u_{nk}$ orbitals in the real space is represented as function on this grid.

The 3D reciprocal lattice. The $u_{nk}$ orbitals in the frequency space are represented as functions on this grid.

This grid is the dual grid of HomeCell3D. Given a homecell, one can find its corresponding reciprocal lattice by

reciprocal_lattice = transform_grid(homecell)

The 3D Brillouin zone. The usage is the same of HomeCell3D.

The 3D crystal lattice. This is the dual grid of BrillouinZone3D.

basis(g)

Basis vectors for the grid. Depending on the type of grid, the basis is that of the grid such as real/reciprocal lattice vectors.

A tuple of vectors will be returned, and the vectors returned are always ket vectors.

Example:

julia> basis(homecell)

(ket: 0.000, 0.866, 0.000, ket: 1.000, 0.500, 0.866, ket: 0.000, 0.000, 0.500)

If it is a matrix that you want, use basis_matrix instead for better performance.

basis_matrix(g)

The basis vector as a matrix, where each column is a basis vector. The matrix is cached because this has an impact on performance.

Example:

julia> basis_matrix(homecell)

3×3 SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 0.0       1.0       0.0
 0.866025  0.5       0.0
 0.0       0.866025  0.5

domain(grid)

The domain of the grid are the extremal miller indices. The grid includes the end points of the domain.

The result is a tuple of n_dim tuples in the form ((x_min, x_max), (y_min, y_max), (z_min, z_max)).

Example:

julia> domain(homecell)

((-2, 1), (-3, 2), (-4, 3))

domain_matrix(grid)

Get the domain of the grid as a matrix, where each row corresponds to a dimension.

julia> domain_matrix(homecell)

3×2 SMatrix{3, 2, Int64, 6} with indices SOneTo(3)×SOneTo(2):
 -2  1
 -3  2
 -4  3

set_domain(grid, new_domain)

Create a copy of grid with its domain set to new_domain.

Example:

julia> odd_grid = set_domain(homecell, ((-0, 0), (-1, 1), (-2, 2)))

type: HomeCell3D
domain: ((0, 0), (-1, 1), (-2, 2))
basis:
    ket: 0.000, 0.866, 0.000
    ket: 1.000, 0.500, 0.866
    ket: 0.000, 0.000, 0.500

expand(grid, factors=[2, 2, 2])

Expand a grid by a factor of factor. The returned grid will have the save basis vector, but it will be larger.

julia> supercell = expand(homecell, [3, 3, 3])

type: HomeCell3D
domain: ((-6, 5), (-9, 8), (-12, 11))
basis:
    ket: 0.000, 0.866, 0.000
    ket: 1.000, 0.500, 0.866
    ket: 0.000, 0.000, 0.500

mins(grid::Grid)

The minimum indices of the grid along each direction as a tuple of integers.

Example:

julia> mins(homecell)

(-2, -3, -4)

maxes(grid)

Similar to mins(grid), but gives the maximum indices instead.

Example:

julia> maxes(homecell)

(1, 2, 3)

center(grid)

The center of a grid. This should be a tuple of zeros for a grid that complies to the centering convention. The result will be a tuple of n_dims integers. A tuple instead of a vector is used for performance reasons.

julia> center(homecell)

(0, 0, 0)

n_dims(T)

Gives the number of dimensions of the type of grid.

Example:

julia> n_dims(HomeCell3D)

3

size(g)

The size of the grid. The result is a tuple of n_dims(g) integers.

julia> size(homecell)

(4, 6, 8)

length(g)

The number of grid points in the grid.

julia> length(homecell)

192

snap(grid, point)

Snap a cartesian coordinate to a grid point.

julia> snap(homecell, [0, sqrt(3), 1.49])

GridVector{HomeCell3D}:
    coefficients: [2, 0, 3]

g[i, j, k]

Indexing a grid with a miller index gives the grid vector corresponding to that index.

julia> origin = homecell[0, 0, 0]

GridVector{HomeCell3D}:
    coefficients: [0, 0, 0]

g[x_1:x_2, y, z_1:z_2]

Ranges are supported, and the result will be a multi-dimensional array.

julia> homecell[-1:0, 0, -1:0]

2×2 Matrix{GridVector{HomeCell3D}}:
 GridVector{HomeCell3D}:
    coefficients: [-1, 0, -1]
  GridVector{HomeCell3D}:
    coefficients: [-1, 0, 0]

 GridVector{HomeCell3D}:
    coefficients: [0, 0, -1]
   GridVector{HomeCell3D}:
    coefficients: [0, 0, 0]

linear_index(g::Grid, i::Integer)

Alternative syntax: g(i).

This provides a mapping between 1D indices (which are to be used for matrix algorithms) and grid vectors.

Example:

julia> homecell(4)

GridVector{HomeCell3D}:
    coefficients: [-1, 0, 0]

julia> homecell(1:3)

3-element Vector{GridVector{HomeCell3D}}:
 GridVector{HomeCell3D}:
    coefficients: [0, 0, 0]

 GridVector{HomeCell3D}:
    coefficients: [1, 0, 0]

 GridVector{HomeCell3D}:
    coefficients: [-2, 0, 0]

The danger of using a one dimensional index is that it can be disassociated with the object that it's referring to. Therefore, a two way mapping must be provided.

julia> linear_index(homecell(4))

4

iterate(grid)

All grids are iterable. The iteration order is the same order as the one dimensional index.

julia> collect(homecell)[5] == homecell(5)

true

Example

For loops

julia> for r in homecell
           r == homecell(1) && println(r)
       end
    
GridVector{HomeCell3D}:
    coefficients: [0, 0, 0]

Higher order functions

julia> sum(r->r'*r, homecell)

1779.138438763306

Mapping over a grid gives you an OnGrid object, which can be indexed by grid vectors.

julia> r2 = map(r->norm(r)^2, homecell);
julia> typeof(r2)

WTP.SimpleFunctionOnGrid{HomeCell3D}

julia> r2[homecell[0, 0, 0]]

0

One can construct model orbitals this way.

julia> gaussian = map(r->exp(-r'*r), homecell);
julia> contourplot(gaussian[homecell[0, :, :]])
┌────────────────────────────────────────┐  ⠀⠀⠀⠀  
6 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│  ┌──┐ 1
  │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│  │▄▄│  
  │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│  │▄▄│  
  │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│  │▄▄│  
  │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│  │▄▄│  
  │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⣀⣀⣀⠤⠤⣀⡀⠀⠀⠀│  │▄▄│  
  │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⡠⠤⠤⠔⠒⠒⠒⠉⠉⠉⠉⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀│  │▄▄│  
  │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡠⠤⠒⠉⢀⣀⡠⠤⠔⠒⠒⠒⠒⠒⠒⠒⠒⢄⠀⠀⠀⠀⠀⠀⢱⠀⠀⠀│  │▄▄│  
  │⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠊⠁⠀⢀⠤⠒⠁⣀⡠⠤⠒⠒⠉⠑⠒⠢⠤⠤⣀⠀⠉⠢⡀⠀⠀⠀⢸⠀⠀⠀│  │▄▄│  
  │⠀⠀⠀⠀⠀⠀⠀⠀⠸⡀⠀⠀⠀⠉⠢⣀⠈⠑⠒⠒⠤⠤⣀⣀⠤⠤⠒⠊⠉⡠⠔⠊⠁⠀⣀⠔⠉⠀⠀⠀│  │▄▄│  
  │⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠑⠤⠤⠤⠤⠤⠤⠤⠤⠔⠒⠊⠉⢀⡠⠔⠒⠉⠀⠀⠀⠀⠀⠀│  │▄▄│  
  │⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⡠⠤⠤⠤⠔⠒⠒⠊⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│  │▄▄│  
  │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠑⠒⠒⠉⠉⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│  │▄▄│  
  │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│  │▄▄│  
1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│  └──┘ 0
  └────────────────────────────────────────┘  ⠀⠀⠀⠀  
  ⠀1⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀8⠀  ⠀⠀⠀⠀  

The result does not look isotropic even though our Gaussian is. This is because the grid is not orthogonal. Keep this in mind for all the pictures.

transform_grid(g)

If we apply a Fourier transform on a function defined on a grid, the result is a function defined on a different grid, which we refer to as the dual grid.

julia> guassian_reciprocal = fft(gaussian);
julia> reciprocal_lattice = grid(guassian_reciprocal)

type: ReciprocalLattice3D
domain: ((-2, 1), (-3, 2), (-4, 3))
basis:
    ket: -0.907, 1.814, -0.000
    ket: 1.047, 0.000, -0.000
    ket: -1.360, 0.000, 1.571

Example:

julia> reciprocal_lattice = transform_grid(homecell)

type: ReciprocalLattice3D
domain: ((-2, 1), (-3, 2), (-4, 3))
basis:
    ket: -0.907, 1.814, -0.000
    ket: 1.047, 0.000, -0.000
    ket: -1.360, 0.000, 1.571

The resulting basis vectors should satisfy aᵢᵀ bᵢ = 2π / sᵢ, where s is the size of the grid (number of grid points in each direction).

julia> basis_matrix(reciprocal_lattice)' * basis_matrix(homecell) * diagm([size(homecell)...])
3×3 Matrix{Float64}:
 6.28319   0.0          0.0
 0.0       6.28319      0.0
 0.0      -6.31568e-17  6.28319

HomeCell3D and ReciprocalLattice3D are dual grids and BrillouinZone3D and RealLattice3D are dual grids.