Center and Spread

ApproximationScheme

The center and the spread are often approximated for performance. Concrete approximations should subtype this.

CosScheme

Approximations based on a single mode of cos.

CosScheme3D(ũ, n_shells = 1)

The function $r^2$ can be approximated with $w_{\mathbf{b}} \cos(\mathbf{b}^T \mathbf{r})$ functions, which are then used for approximating the convolution between $r^2$ and u (inverse fft of ũ). The CosScheme3D includes shells of \mathbf{b} vectors and their corresponding weights $w_{\mathbf{b}}$.

Other dimensions will be implemented as extensions.

Example:

julia> scheme = CosScheme3D(ũ);

julia> length(shells(scheme))
1

shells(scheme)

Shells of $\mathbf{b}$ vectors involved in the approximation scheme. Each shell is a vector of kpoints.

Example:

julia> shells(scheme)

1-element Vector{Vector{<:KPoint}}:
 KPoint[GridVector{BrillouinZone3D}:
    coefficients: [-1, -1, -1]
, GridVector{BrillouinZone3D}:
    coefficients: [-1, 0, 0]
, GridVector{BrillouinZone3D}:
    coefficients: [0, -1, 0]
, GridVector{BrillouinZone3D}:
    coefficients: [0, 0, -1]
, GridVector{BrillouinZone3D}:
    coefficients: [0, 0, 1]
, GridVector{BrillouinZone3D}:
    coefficients: [0, 1, 0]
, GridVector{BrillouinZone3D}:
    coefficients: [1, 0, 0]
, GridVector{BrillouinZone3D}:
    coefficients: [1, 1, 1]
]

weights(scheme)

The weights corresponding to each shell within a scheme. The weights are ordered from the inner-most to the outer-most shell.

julia> weights(scheme)
1-element Vector{Number}:
 5.336038037571918

The neighbor integrals indexed by two kpoints. For each pair of kpoint, the integrals are stored as a matrix. This is the same as the $M_{mn}^{k, b}$ matrix in MLWF. The matrix elements are accessed by M[k, k+b][m, n]

neighbor_basis_integral(scheme)

The integrals between neighboring k-points (The MMN matrix). The integral is amonst immediate neighbor because the $\cos$ approximation is truncated at the first mode.

julia> M = neighbor_basis_integral(scheme)
julia> M[brillouin_zone[0, 0, 0], brillouin_zone[0, 0, 1]]
4×4 Matrix{ComplexF64}:
     0.85246+0.512198im     0.0798774-0.00342381im  4.54441e-8-1.98313e-8im  8.32139e-9+4.06675e-8im
   0.0057819+0.021992im     -0.363241-0.324958im     -0.393272-0.187671im     -0.489628+0.127022im
   0.0112149-0.000481798im  -0.194858+0.141126im      0.665299+0.0837434im    -0.498015-0.0113525im
 -0.00992131-0.0216156im     0.432817+0.269315im     -0.319281+0.193723im     -0.524017-0.0193465im

getindex(M, k_1, k_2)

The integral matrix between the neighboring k-points k_1 and k_2. Can also write M[k_1, k_2]. Note that M[k_1, k_2] = M[k_2, k_1]'. So only one of the matrices is stored.

julia> M[brillouin_zone[0, 0, 0], brillouin_zone[0, 0, -1]] == M[brillouin_zone[0, 0, -1], brillouin_zone[0, 0, 0]]'
true

setindex!(M, value, g...)

Set the integral matrix between two k-points g[1] and g[2] to value. Can also write M[g...] = value.

haskey(M, k_1, k_2)

Check if the integral matrix between k_1 and k_2 has been computed and stored.

julia> haskey(M, brillouin_zone[0, 0, 1], brillouin_zone[0, 0, 0])
true

julia> haskey(M, brillouin_zone[0, 0, 1], brillouin_zone[0, 0, -1])
false

gauge_transform(M, gauge)

Perform a gauge transform on the neighbor integrals.

$U^{k \dagger} M^{k, k+b} U^{k+b}$

julia> M = gauge_transform(M, U);
julia> M[brillouin_zone[0, 0, 1], brillouin_zone[0, 0, 0]]
4×4 adjoint(::Matrix{ComplexF64}) with eltype ComplexF64:
    0.85246-0.512198im    0.0057819-0.021992im  0.0112149+0.000481798im  -0.00992131+0.0216156im
  0.0798774+0.00342381im  -0.363241+0.324958im  -0.194858-0.141126im        0.432817-0.269315im
 4.54441e-8+1.98313e-8im  -0.393272+0.187671im   0.665299-0.0837434im      -0.319281-0.193723im
 8.32139e-9-4.06675e-8im  -0.489628-0.127022im  -0.498015+0.0113525im      -0.524017+0.0193465im

Use the same algorithm as in the original Wannier90 (without guiding centers, fixed centers, etc.). Apply this algorithm if you like to cross validate with Wannier90.

Mostly the same algorithms as W90BranchCut, but make a potential more consistent choice of branch cut.

The truncated convolution algorithm. (Resta's first spread)

center(M, scheme, n, ::Type{W90BranchCut})

Compute the center of the nth Wannier orbital using the original Wannier90 approach. The replication is exact.

Example:

julia> center(M, scheme, 1, W90BranchCut)
3-element Vector{Float64}:
 -8.73495454038011
  3.9151488754936588
  4.021344063664258
julia> center(M, scheme, 2, W90BranchCut)
3-element Vector{Float64}:
  0.11705358388655285
 -1.4773864607361118
 -2.4684488265569624

center(M, scheme, n, ::Type{TruncatedConvolution})

Compute the center of the nth Wannier orbital using the turncated convolution algorihtm.

Example:

julia> center(M, scheme, 1, TruncatedConvolution)
3-element Vector{Float64}:
 -8.70234629129622
  4.016099860574682
  4.093053970246974
julia> center(M, scheme, 1, TruncatedConvolution)
3-element Vector{Float64}:
  1.3634042350006577
 -2.712474285697622
 -3.6252083684030336

One can compare this with the "exact" center

julia> wanniers = commit_gauge(ũ)(:);
julia> _, r̃2 = compute_r2(supercell(u));
julia> c_1, σ²_1 = center_spread(fft(abs2(ifft(wanniers[1])), false), r̃2)
([-8.831796622851812, 3.953063696051289, 3.940388714312329], 28.748708554079474)
julia> c_2, σ²_2 = center_spread(fft(abs2(ifft(wanniers[2])), false), r̃2)
([-9.192231000846695, 7.604578700842824, 6.189389913055544], 27.460489592320883)

spread(M, scheme, n, W90BranchCut)

Compute the spread of the nth Wannier orbital using the original Wannier90 approach.

julia> spread(M, scheme, 1, W90BranchCut)
15.742034477681969
julia> spread(M, scheme, 2, W90BranchCut) # failure.
133.09413338071354

spread(M, scheme, n, W90BranchCut)

Compute the spread of the nth Wannier orbital using the original Wannier90 approach.

julia> spread(M, scheme, 1, W90BranchCut)
15.742034477681969
julia> spread(M, scheme, 2, W90BranchCut) # failure.
133.09413338071354

spread(M, scheme, n, TruncatedConvolution)

Compute the spread of the nth Wannier orbital using the truncated convolution.

julia> spread(M, scheme, 1, TruncatedConvolution)
17.50438313709964
julia> spread(M, scheme, 2, TruncatedConvolution)
17.313972338201154

Compare this to the "exact" spread.

julia> σ²_1, σ²_2
(28.748708554079474, 27.460489592320883)

gauge_gradient(M, scheme, brillouin_zone, W90BranchCut)

Gauge gradient in the original Wannier90. The result is a OnGrid{<:BrillouinZone} with a matrix on each k-point.

julia> G_w = gauge_gradient(M, scheme, brillouin_zone, W90BranchCut);
julia> G_w[brillouin_zone[0, 0, 1]]
4×4 Matrix{ComplexF64}:
       0.0-0.030634im     -1.13757+0.493027im   0.145804-0.00940629im  -0.0627567-0.0302657im
   1.13757+0.493027im          0.0+1.76773im   -0.876742+0.121255im      0.885767+0.370595im
 -0.145804-0.00940629im   0.876742+0.121255im        0.0-0.158586im    -0.0557475+0.0957635im
 0.0627567-0.0302657im   -0.885767+0.370595im  0.0557475+0.0957635im          0.0-0.12397im

gauge_gradient(M, scheme, brillouin_zone, TruncatedConvolution)

Gauge gradient for the truncated convolution. The result is a OnGrid{<:BrillouinZone} with a matrix on each k-point.

julia> G_t = gauge_gradient(U, scheme, brillouin_zone, TruncatedConvolution);
julia> G_t[brillouin_zone[0, 0, 1]]
4×4 Matrix{ComplexF64}:
        0.0-0.0388009im  -0.00321511+0.158118im    0.160808-0.0451399im  -0.0646812+0.0334942im
 0.00321511+0.158118im           0.0-0.179264im   -0.158035-0.0368048im    0.205487+0.0206957im
  -0.160808-0.0451399im     0.158035-0.0368048im        0.0-0.108496im    -0.070136+0.157672im
  0.0646812+0.0334942im    -0.205487+0.0206957im   0.070136+0.157672im          0.0-0.0640979im