GitHub

N-body matrix elements

The matrix element of an $N$-body operator between two Slater determinants may be expanded according to the Löwdin rules (which reduce to the Slater–Condon rules if all single-particle orbitals are orthogonal):

\[\begin{equation} \label{eqn:matrix-element-expansion} \matrixel{\Phi_A}{\Omega_n}{\Phi_B} = \frac{1}{n!}\sum_p (-)^p \matrixel{k_1k_2...k_n}{\Omega_n}{l_1l_2...l_n} D^{AB}({k_1k_2...k_n}|{l_1l_2...l_n}) \end{equation}\]

where $D^{AB}({k_1k_2...k_n}|{l_1l_2...l_n})$ is the determinant minor of the orbital overlap determinant $D^{AB}$ with the rows ${k_1k_2...k_n}$ and columns ${l_1l_2...l_n}$ stricken out, and $p$ runs over all permutations.

In general, a term in the expansion is thus of the form

\[\begin{equation} \alpha\matrixel{k_1k_2...k_n}{\Omega_n}{l_1l_2...l_n}\braket{a}{b}\braket{c}{d}\dots\braket{y}{z}, \end{equation}\]

where $\alpha$ is a scalar. This is represented by NBodyTerm type.

EnergyExpressions.OrbitalOverlap — Type
OrbitalOverlap(a,b)

Represents the overlap between the orbitals a and b in a N-body matrix element expansion.

Examples

julia> EnergyExpressions.OrbitalOverlap(:a,:b)
⟨a|b⟩
source
EnergyExpressions.OrbitalMatrixElement — Type
OrbitalMatrixElement(a,o,b)

Represents the N-body matrix element between the sets of orbitals a and b.

Examples

julia> struct MyTwoBodyOperator <: TwoBodyOperator end

julia> EnergyExpressions.OrbitalMatrixElement((:a,:b), MyTwoBodyOperator(), (:c,:d))
⟨a b|MyTwoBodyOperator()|c d⟩
source
EnergyExpressions.numbodies — Function
numbodies(::NBodyOperator{N})

Returns the number of bodies coupled by the N-body operator, i.e. N.

source
numbodies(lco::LinearCombinationOperator)

Returns the maximum number of bodies coupled by any of the N-body operators in the LinearCombinationOperator.

source
numbodies(::OrbitalOverlap)

Returns the number of bodies coupled by the zero-body operator in the orbital overlap, i.e. 0.

source
numbodies(::OrbitalMatrixElement{N})

Returns the number of bodies coupled by the operator, i.e. N.

source
EnergyExpressions.isdependent — Function
isdependent(o::OrbitalOverlap, orbital)

Returns true if the OrbitalOverlap o depends on orbital.

Examples

julia> isdependent(OrbitalOverlap(:a,:b), :a)
false

julia> isdependent(OrbitalOverlap(:a,:b), Conjugate(:a))
true

julia> isdependent(OrbitalOverlap(:a,:b), :b)
true
source
isdependent(o::OrbitalMatrixElement, orbital)

Returns true if the OrbitalMatrixElement o depends on orbital.

Examples

julia> isdependent(EnergyExpressions.OrbitalMatrixElement((:a,), OneBodyHamiltonian(), (:b,)), :a)
false

julia> isdependent(EnergyExpressions.OrbitalMatrixElement((:a,), OneBodyHamiltonian(), (:b,)), Conjugate(:a))
true

julia> isdependent(EnergyExpressions.OrbitalMatrixElement((:a,), OneBodyHamiltonian(), (:b,)), :b)
true

julia> isdependent(EnergyExpressions.OrbitalMatrixElement((:a,:b,), CoulombInteraction(), (:c,:d)), :c)
true
source
isdependent(nbt::NBodyTerm, o)

Returns true if any of the factors comprising nbt is dependent on the orbital o. Not that the result is dependent on whether o is conjugated or not.

source
EnergyExpressions.transform — Function
transform(f::Function, nbt::NBodyTerm)

Transform integrals of the the N-body matrix element expansion term nbt according to the function f, which should accept a single NBodyTermFactor as its argument.

source
transform(f::Function, nbme::NBodyMatrixElement)

Transform integrals of the the N-body matrix element nbme according to the function f, which should accept a single NBodyTermFactor as its argument, and return a NBodyMatrixElement. This is useful for adapting energy expressions to specific symmetries of the system under consideration.

source
EnergyExpressions.overlap_matrix — Function
overlap_matrix(a::SlaterDeterminant, b::SlaterDeterminant[, overlaps=[]])

Generate the single-particle orbital overlap matrix, between the orbitals in the Slater determinants a and b. All orbitals are assumed to be orthogonal, except for those which are given in overlaps.

Examples

First we define two Slater determinants that have some orbitals in common:

julia> sa = SlaterDeterminant([:i, :j, :l,:k̃])
i(1)j(2)l(3)k̃(4) - i(1)j(2)l(4)k̃(3) - i(1)j(3)l(2)k̃(4) + i(1)j(3)l(4)k̃(2) + …  + i(4)j(1)l(3)k̃(2) + i(4)j(2)l(1)k̃(3) - i(4)j(2)l(3)k̃(1) - i(4)j(3)l(1)k̃(2) + i(4)j(3)l(2)k̃(1)

julia> sb = SlaterDeterminant([:i, :j, :k, :l̃])
i(1)j(2)k(3)l̃(4) - i(1)j(2)k(4)l̃(3) - i(1)j(3)k(2)l̃(4) + i(1)j(3)k(4)l̃(2) + …  + i(4)j(1)k(3)l̃(2) + i(4)j(2)k(1)l̃(3) - i(4)j(2)k(3)l̃(1) - i(4)j(3)k(1)l̃(2) + i(4)j(3)k(2)l̃(1)

The orbital overlap matrix by default is

julia> overlap_matrix(sa, sb)
4×4 SparseArrays.SparseMatrixCSC{EnergyExpressions.NBodyTerm,Int64} with 2 stored entries:
  [1, 1]  =  1
  [2, 2]  =  1

which has only two non-zero entries, since only two of the orbitals are common between the Slater determinants sa and sb.

We can then define that the orbitals k̃ and l̃ are non-orthogonal:

julia> overlap_matrix(sa, sb, [OrbitalOverlap(:k̃,:l̃)])
4×4 SparseArrays.SparseMatrixCSC{EnergyExpressions.NBodyTerm,Int64} with 3 stored entries:
  [1, 1]  =  1
  [2, 2]  =  1
  [4, 4]  =  ⟨k̃|l̃⟩

We can even specify that the orbital k̃ is non-orthogonal to itself (this can be useful when the k̃ is a linear combination of orthogonal orbitals):

julia> overlap_matrix(sa, sa, [OrbitalOverlap(:k̃,:k̃)])
4×4 SparseArrays.SparseMatrixCSC{EnergyExpressions.NBodyTerm,Int64} with 4 stored entries:
  [1, 1]  =  1
  [2, 2]  =  1
  [3, 3]  =  1
  [4, 4]  =  ⟨k̃|k̃⟩

Notice that this overlap matrix was calculated between the Slater determinant sa and itself.

source
EnergyExpressions.EnergyExpression — Type
EnergyExpression

An energy expression is given by an energy matrix, or interaction matrix, sandwiched between a vector of mixing coefficients: E = c'H*c, where c are the mixing coefficients and H the energy matrix.

source
Base.Matrix — Type
Matrix(a, op::QuantumOperator, b[, overlaps])

Generate the matrix corresponding to the quantum operator op, between the SlaterDeterminants a and b, i.e ⟨a|op|b⟩. It is possible to specify non-orthogonalities between single-particle orbitals in overlaps.

source
Matrix(op::QuantumOperator, slater_determinants[, overlaps])

Generate the matrix corresponding to the quantum operator op, between the different slater_determinants. It is possible to specify non-orthogonalities between single-particle orbitals in overlaps.

source

Calculation of determinants

Actually computing the matrix element expansion $\eqref{eqn:matrix-element-expansion}$ is a combinatorial problem, that grows factorially with the amount of non-orthogonal orbital pairs. Furthermore, of the $(n!)^2$ terms generated from the expansion, only $n!$ are distinct, due to the integrals being symmetric with respect to interchange of the coordinates [hence the normalization factor $(n!)^{-1}$]. Thankfully, there are few symmetries that can be employed, to generate only the distinct permutations, as well as the fact that the overlap matrix is very sparse.

Finding non-zero minors of the overlap determinant

The algorithm to find which minor determinants $\Gamma^{(N)}(k_1k_2...k_N|l_1l_2...l_N)$ do not vanish, and hence which $N$ orbitals $k_1k_2...k_N,l_1l_2...l_N$ the $N$-body operator should be contracted over, is described briefly below. It is devised to be optimal for orthogonal orbitals (i.e. linear complexity $\mathcal{O}(Nn)$ where $n$ is the number of orbitals), and near-optimal for a small amount of non-orthogonal orbitals.

Given:

  • An $N$-body operator (implying $N$ rows and $N$ need to stricken out), and
  • an $n\times n$ matrix, with coordinates of non-zero matrix elements: $I,J$ (from these vectors, vanishing rows/columns can easily be deduced $\implies$ $N_r$ row/$N_c$ column "rank", i.e. yet to be stricken out),

do

  • find all $N_r$-combinations of the remaining rows,
  • for each such combination, find all columns which would be affected if striking out that particular combination of rows,
    • if the "support" (i.e. the only non-zero elements) of any of those columns vanishes when striking out the rows, that column must be stricken out, too. Total number of these columns is named $N_{cm}$,
    • if more than $N_c$ columns must be stricken out ($N_{cm}>N_c$), that row combination is unviable,
    • find all $N_c - N_{cm}$-combinations of the remaining columns,
    • for each such combination of columns, find all rows which would be affected,
      • if the support of any of those rows vanishes when striking out the candidate columns, and the row is not in the candidate set of rows to be stricken out, the column combination is unviable.
EnergyExpressions.nonzero_minors — Function
nonzero_minors(N, overlap) -> (ks,ls)

Find all (distinct) minor determinants of order N of the orbital overlap matrix that do not vanish, i.e. all non-vanishing minors are guaranteed to be present, but not all of the returned minors are guaranteed to be non-zero. Vanishing minors returned arise when the overlap matrix is rank deficient, which is unlikely to happen when computing energy expressions, but must still be guarded against. This is most easily checked by actually calculating the cofactor, which is most likely desired anyway.

source

Contracting over the stricken out orbitals

We use Julia's built-in Base.Cartesian.@nloops iterators to span the space of all possible choices of orbitals for the contraction of orbitals. If two or more orbitals are the same, the matrix element is trivially zero (the “Fermi hole”). To avoid double-counting, we also only consider those indices that are above the hyper-diagonal.

EnergyExpressions.detaxis — Function
detaxis(i::CartesianIndex{N})

Generate the axis index vector for the determinant minor, whose rows or columns represented by the CartesianIndex i should be omitted. Implemented via complement.

source
EnergyExpressions.cofactor — Function
cofactor(k, l, A)

Calculate the cofactor of A, where the rows k and the columns l have been stricken out. The cofactor is calculated recursively, by expanding the minor determinants in cofactors, so this function should only be used in case it is known that the cofactor is non-zero.

source
LinearAlgebra.det — Function
det(A)

Calculate the determinant of the matrix A whose elements are of the NBodyTerm type, by expanding the determinant along the first column. This is an expensive operation, and should only be done with relatively sparse matrices.

source
EnergyExpressions.@above_diagonal_loop — Macro
above_diagonal_loop(N, itersym, imax, args...)

Generate N Cartesian loops for the iteration variables itersym_{1:N}, where itersym_N ∈ 1:imax, itersym_{N-1} ∈ itersym_N+1:imax, etc, i.e. above the hyper-diagonal of the N-dimensional hypercube with the side imax. args... is passed on to Base.Cartesian._nloops. above_diagonal_loop is nestable.

Examples

julia> @above_diagonal_loop 2 i 3 begin
           println("==================================")
           println("i = ", Base.Cartesian.@ntuple 2 i)
           @above_diagonal_loop 2 j 3 begin
               println("j = ", Base.Cartesian.@ntuple 2 j)
           end
       end
==================================
i = (2, 1)
j = (2, 1)
j = (3, 1)
j = (3, 2)
==================================
i = (3, 1)
j = (2, 1)
j = (3, 1)
j = (3, 2)
==================================
i = (3, 2)
j = (2, 1)
j = (3, 1)
j = (3, 2)
source
EnergyExpressions.@anti_diagonal_loop — Macro
anti_diagonal_loop(N, itersym, imax, args...)

Generate N Cartesian loops for the iteration variables itersym_{1:N}, where itersym_N ∈ 1:imax, itersym_{N-1} ∈ 1:imax\itersym_N, etc, i.e. no two iteration variables have the same values simultaneously. args... is passed on to Base.Cartesian._nloops; however, preexpr is already used to skip the diagonal elements. anti_diagonal_loop is nestable.

Examples

julia> @anti_diagonal_loop 3 i 3 begin
           println("-----------------------------")
           t = (Base.Cartesian.@ntuple 3 i)
           println("$t: ", allunique(t))
           @anti_diagonal_loop 2 j 2 begin
               u = (Base.Cartesian.@ntuple 2 j)
               println("$u: ", allunique(u))
           end
       end
-----------------------------
(3, 2, 1): true
(2, 1): true
(1, 2): true
-----------------------------
(2, 3, 1): true
(2, 1): true
(1, 2): true
-----------------------------
(3, 1, 2): true
(2, 1): true
(1, 2): true
-----------------------------
(1, 3, 2): true
(2, 1): true
(1, 2): true
-----------------------------
(2, 1, 3): true
(2, 1): true
(1, 2): true
-----------------------------
(1, 2, 3): true
(2, 1): true
(1, 2): true
source