# Grid

`WTP.Grid`

β TypeWTP provides *non-orthogonal 3D periodic* grids that support

- Three Indexing schemes and iteration.
- Fourier transform.

The goals is to provide a uniform interface for dealing with

- The reciprocal lattice.
- The Brillouin zone.
- The crystal lattice.
- 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}}
```

## Centering convention

When dealing with multiple periodic grids, two questions on the matter of convention give people a lot of headache. Consider the case of 1D, the questions are

- What are the left-most and the right-most grid points?
- Which grid point is considered the center?

We will deal with these problems by imposing a simple convention for all grids in our package, and we refer to it as the centering convention. Again, in the case of 1D, suppose that the step size of the grid is $a$.

- For a grid with $2 N$ grid points, the leftmost grid point is at $-N a$, and the rightmost grid point is at $(N-1) a$.
- For a grid with $2 N + 1$ grid points, the leftmost grid point is at $-N a$, and the rightmost grid point is at $N a$.
- The center is at the origin.

A grid that does not comply to the centering convention is said to have an offset, which is the same as the center of the grid (offset from the origin).

## Creating a grid

`WTP.make_grid`

β Function`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
```

We provide four different grids by default

`WTP.HomeCell3D`

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

`WTP.ReciprocalLattice3D`

β TypeThe 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)`

`WTP.BrillouinZone3D`

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

.

`WTP.RealLattice3D`

β TypeThe 3D crystal lattice. This is the dual grid of `BrillouinZone3D`

.

## Basic usage

### Basis

`WTP.basis`

β Method`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.

`WTP.basis_matrix`

β Function`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

`WTP.domain`

β Function`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))
```

`WTP.domain_matrix`

β Function`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
```

`WTP.set_domain`

β Function`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
```

`WTP.expand`

β Method`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
```

### Boundaries

`WTP.mins`

β Function`mins(grid::Grid)`

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

Example:

```
julia> mins(homecell)
(-2, -3, -4)
```

`WTP.maxes`

β Function`maxes(grid)`

Similar to `mins(grid)`

, but gives the maximum indices instead.

Example:

```
julia> maxes(homecell)
(1, 2, 3)
```

### Center

`WTP.center`

β Method`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)
```

### Dimension

`WTP.n_dims`

β Function`n_dims(T)`

Gives the number of dimensions of the *type* of grid.

Example:

```
julia> n_dims(HomeCell3D)
3
```

`Base.size`

β Method`size(g)`

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

integers.

```
julia> size(homecell)
(4, 6, 8)
```

`Base.length`

β Method`length(g)`

The number of grid points in the grid.

```
julia> length(homecell)
192
```

### Snap

`WTP.snap`

β Function`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]
```

## Indexing a grid

### Single index

`Base.getindex`

β Method`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]
```

### Range index

`Base.getindex`

β Method`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

`WTP.linear_index`

β Method`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
```

The correspondance between the one dimensional index and the grid point is an implementation detail, and one must not rely on it. In particular, one should not use it as a function argument. The two way mapping will be lost since indices by themselves have no context (the grid). one should pass a `AbstractGridVector`

instead most of the time.

## Iterations and maps over grids

`Base.iterate`

β Method`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 a grid

`WTP.transform_grid`

β Function`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.