5. operators

5.1. O

Overlap operator.

Python

Julia

1def O(atoms, W):
2    return atoms.Omega * W

1function op_O(atoms::Atoms, W::Matrix{ComplexF64})
2    return atoms.Omega .* W
3end

5.2. L

Laplacian operator.

Python

Julia

1def L(atoms, W):
2    if len(W) == len(atoms.G2c):
3        G2 = atoms.G2c[:, None]
4    else:
5        G2 = atoms.G2[:, None]
6    return -atoms.Omega * G2 * W

1function op_L(atoms::Atoms, W::Matrix{ComplexF64})
2    if size(W)[1] == length(atoms.G2c)
3        G2 = atoms.G2c
4    else
5        G2 = atoms.G2
6    end
7    return -atoms.Omega .* G2 .* W
8end

5.3. Linv

Inverse Laplacian operator.

Python

Julia

 1def Linv(atoms, W):
 2    out = np.empty_like(W, dtype=complex)
 3    with np.errstate(divide="ignore", invalid="ignore"):
 4        if W.ndim == 1:
 5            out = W / atoms.G2 / -atoms.Omega
 6            out[0] = 0
 7        else:
 8            G2 = atoms.G2[:, None]
 9            out = W / G2 / -atoms.Omega
10            out[0, :] = 0
11    return out

 1function op_Linv(atoms::Atoms, W::Array{ComplexF64})
 2    out = W ./ atoms.G2 ./ -atoms.Omega
 3    out[1] = 0.0
 4    return out
 5end
 6function op_Linv(atoms::Atoms, W::Matrix{ComplexF64})
 7    out = W ./ atoms.G2 ./ -atoms.Omega
 8    out[1, :] .= 0.0
 9    return out
10end

5.4. I

Backwards transformation from reciprocal space to real-space.

Python

Julia

 1def I(atoms, W):
 2    n = np.prod(atoms.s)
 3    if len(W) == len(atoms.G2):
 4        Wfft = W
 5    else:
 6        if W.ndim == 1:
 7            Wfft = np.zeros(n, dtype=complex)
 8        else:
 9            Wfft = np.zeros((n, atoms.Nstate), dtype=complex)
10        Wfft[atoms.active] = W
11
12    if W.ndim == 1:
13        Wfft = Wfft.reshape(atoms.s)
14        Finv = np.fft.ifftn(Wfft).ravel()
15    else:
16        Wfft = Wfft.reshape(np.append(atoms.s, atoms.Nstate))
17        Finv = np.fft.ifftn(Wfft, axes=(0, 1, 2)).reshape((n, atoms.Nstate))
18    return Finv * n

 1function op_I(atoms::Atoms, W::Array{ComplexF64})
 2    n = prod(atoms.s)
 3    if size(W)[1] == length(atoms.G2)
 4        Wfft = W
 5    else
 6        Wfft = convert(Array{ComplexF64}, atoms.active)
 7        Wfft[Wfft.!=0] = W
 8    end
 9    Wfft = reshape(Wfft, atoms.s[1], atoms.s[2], atoms.s[3])
10    Finv = vec(ifft(Wfft))
11    return Finv .* n
12end
13function op_I(atoms::Atoms, W::Matrix{ComplexF64})
14    Finv = zeros(ComplexF64, length(atoms.G2), size(W)[2])
15    for i = 1:size(W)[2]
16        Finv[:, i] = op_I(atoms, W[:, i])
17    end
18    return Finv
19end

5.5. J

Forward transformation from real-space to reciprocal space.

Python

Julia

1def J(atoms, W):
2    n = np.prod(atoms.s)
3    if W.ndim == 1:
4        Wfft = W.reshape(atoms.s)
5        F = np.fft.fftn(Wfft).ravel()
6    else:
7        Wfft = W.reshape(np.append(atoms.s, atoms.Nstate))
8        F = np.fft.fftn(Wfft, axes=(0, 1, 2)).reshape((n, atoms.Nstate))
9    return F / n

 1function op_J(atoms::Atoms, W::Array{ComplexF64})
 2    n = prod(atoms.s)
 3    Wfft = reshape(W, atoms.s[1], atoms.s[2], atoms.s[3])
 4    Finv = vec(fft(Wfft))
 5    return Finv ./ n
 6end
 7function op_J(atoms::Atoms, W::Matrix{ComplexF64})
 8    Finv = zeros(ComplexF64, length(atoms.G2), size(W)[2])
 9    for i = 1:size(W)[2]
10        Finv[:, i] = op_J(atoms, W[:, i])
11    end
12    return Finv
13end

5.6. Idag

Conjugated backwards transformation from real-space to reciprocal space.

Python

Julia

1def Idag(atoms, W):
2    n = np.prod(atoms.s)
3    F = J(atoms, W)
4    F = F[atoms.active]
5    return F * n

 1function op_Idag(atoms::Atoms, W::Array{ComplexF64})
 2    n = prod(atoms.s)
 3    F = op_J(atoms, W)
 4    F = F[atoms.active[:, 1]]
 5    return F .* n
 6end
 7function op_Idag(atoms::Atoms, W::Matrix{ComplexF64})
 8    F = zeros(ComplexF64, length(atoms.G2c), size(W)[2])
 9    for i = 1:size(W)[2]
10        F[:, i] = op_Idag(atoms, W[:, i])
11    end
12    return F
13end

5.7. Jdag

Conjugated forward transformation from reciprocal space to real-space.

Python

Julia

1def Jdag(atoms, W):
2    n = np.prod(atoms.s)
3    Finv = I(atoms, W)
4    return Finv / n

 1function op_Jdag(atoms::Atoms, W::Array{ComplexF64})
 2    n = prod(atoms.s)
 3    Finv = op_I(atoms, W)
 4    return Finv ./ n
 5end
 6function op_Jdag(atoms::Atoms, W::Matrix{ComplexF64})
 7    n = prod(atoms.s)
 8    Finv = op_I(atoms, W)
 9    return Finv ./ n
10end