# System of equations

When dealing with large energy expressions resulting from many configurations, when deriving the orbital equations (see N-body equations and Calculus of Variations – Implementation), many of the integrals involved will be shared between the different terms of the equations. It is therefore of interest to gather a list of integrals that can be calculated once at every iteration (in a self-consistent procedure for finding eigenstates, or in time-propagation), and whose values can then be reused when solving the individual equations. The routines described below, although primitive, aid in this effort.

The idea is the following: With an energy expression on the form

$$$\begin{equation} E(\vec{P},\vec{c}) = \frac{\vec{c}^H \mat{H}\vec{c}}{\vec{c}^H\vec{c}} \end{equation}$$$

where $\vec{P}$ is a set of orbitals, $\vec{c}$ is a vector of mixing coefficients, and

$$$\begin{equation} \mat{H}_{ij} \defd \matrixel{\chi_i}{\Hamiltonian}{\chi_j}, \end{equation}$$$

all terms of the equations can be written on the form

$$$\begin{equation} \label{eqn:equation-term} \operator{A}\ket{\chi} \underbrace{\left[ \sum_n \alpha_n \conj{c}_{i_n} c_{j_n} \prod_k \int_{n_k} \right]}_{\textrm{Rank 0}}, \end{equation}$$$

(or on its dual form with conjugated orbitals) where $\operator{A}$ is some one-body operator, $\alpha_n$ a coefficient due to the energy expression, $\conj{c}_{i_n} c_{j_n}$ the coefficient due to the multi-configurational expansion, and finally $\int_{n_k}$ all the extra integrals (which are necessarily zero-body operators) that may appear due to non-orthogonal orbitals (and which may be shared between many equations). The operator $\operator{A}$ can have ranks 0–2; formally, it can only have rank 0 (multiplication by a scalar) or 2 (multiplication by a matrix). What we somewhat sloppily to refer by “rank 1” is an operator of rank 2, but which is diagonal in the underlying coordinate, such as a local potential.

Thus, to efficiently perform an iteration, one would first compute all common integrals $\int_k$, and then for every equation term of the form $\eqref{eqn:equation-term}$, form the coefficient in the brakets, calculate the action of the operator $\operator{A}$ on the orbital $\chi$, multiplied the coefficient.

There are a few important special cases of $\eqref{eqn:equation-term}$:

1. Field-free, one-body Hamiltonian, i.e. $\operator{A}=\hamiltonian_0$. This is the contribution from the orbital $\ket{\chi}$ to itself. Rank 2.

2. One-body Hamiltonian, including field, i.e. $\operator{A}=\hamiltonian$. This is the contribution from another orbital $\ket{\chi'}$ to $\ket{\chi}$ via some off-diagonal coupling, such as an external field (e.g. an electro–magnetic field). Rank 2.

3. Direct interaction, i.e. $\operator{A}=\direct{kl}$, where two other orbitals $\ket{\chi_k}$ and $\ket{\chi_l}$ together form a potential acting on $\ket{\chi}$. “Rank 1”.

4. Exchange interaction, i.e. $\operator{A}=\exchange{kl}$, where another orbital $\ket{\chi_k}$ and $\ket{\chi}$ together form a potential acting on a third orbital $\ket{\chi_l}$. Rank 2.

5. Source term, i.e. a contribution that does not involve $\ket{\chi}$ in any way. This term arises from other configurations in the multi-configurational expansion. Case 2. is also formulated in this way. Rank 0–2. If for some reason the source orbital and/or the operator acting on it is fixed, it may be possible to precompute the effect of the operator on the source orbital and reduce the computational complexity to a rank 0-equivalent operation.

## Example

We start by defining a set of configurations and specifying that a few of the constituent orbitals are non-orthogonal:

julia> cfgs = [[:i, :j, :l, :k̃], [:i, :j, :k, :l̃]]
2-element Array{Array{Symbol,1},1}:
[:i, :j, :l, :k̃]
[:i, :j, :k, :l̃]

julia> continua = [:k̃, :l̃]
2-element Array{Symbol,1}:
:k̃
:l̃

julia> overlaps = [OrbitalOverlap(i,j) for i in continua for j in continua]
4-element Array{OrbitalOverlap{Symbol,Symbol},1}:
⟨k̃|k̃⟩
⟨k̃|l̃⟩
⟨l̃|k̃⟩
⟨l̃|l̃⟩

We then set up the energy expression as before:

julia> H = OneBodyHamiltonian() + CoulombInteraction()
ĥ + ĝ

julia> E = Matrix(H, SlaterDeterminant.(cfgs), overlaps)
2×2 Array{EnergyExpressions.NBodyMatrixElement,2}:
(i|i)⟨k̃|k̃⟩ + (j|j)⟨k̃|k̃⟩ + (l|l)⟨k̃|k̃⟩ + (k̃|k̃) - G(i,j)⟨k̃|k̃⟩ + F(i,j)⟨k̃|k̃⟩ + … - G(i,k̃) + F(i,k̃) - G(j,k̃) + F(j,k̃) - G(l,k̃) + F(l,k̃)  …  (l|k)⟨k̃|l̃⟩ - [i l|k i]⟨k̃|l̃⟩ + [i l|i k]⟨k̃|l̃⟩ - [j l|k j]⟨k̃|l̃⟩ + [j l|j k]⟨k̃|l̃⟩ - [l k̃|l̃ k] + [l k̃|k l̃]
(k|l)⟨l̃|k̃⟩ - [i k|l i]⟨l̃|k̃⟩ + [i k|i l]⟨l̃|k̃⟩ - [j k|l j]⟨l̃|k̃⟩ + [j k|j l]⟨l̃|k̃⟩ - [k l̃|k̃ l] + [k l̃|l k̃]                                     (i|i)⟨l̃|l̃⟩ + (j|j)⟨l̃|l̃⟩ + (k|k)⟨l̃|l̃⟩ + (l̃|l̃) - G(i,j)⟨l̃|l̃⟩ + F(i,j)⟨l̃|l̃⟩ + … - G(i,l̃) + F(i,l̃) - G(j,l̃) + F(j,l̃) - G(k,l̃) + F(k,l̃)

Finally, we derive the coupled integro-differential equation system for the continuum orbitals k̃, l̃:

julia> eqs = diff(E, Conjugate.(continua))
EnergyExpressions.MCEquationSystem(EnergyExpressions.OrbitalEquation{Symbol,SparseArrays.SparseMatrixCSC{LinearCombinationEquation,Int64}}[OrbitalEquation(k̃):
[1, 1]  =  + (i|i)𝐈₁|k̃⟩ + (j|j)𝐈₁|k̃⟩ + (l|l)𝐈₁|k̃⟩ + ĥ|k̃⟩ - G(i,j)𝐈₁|k̃⟩ + F(i,j)𝐈₁|k̃⟩ - G(i,l)𝐈₁|k̃⟩ + F(i,l)𝐈₁|k̃⟩ - G(j,l)𝐈₁|k̃⟩ + F(j,l)𝐈₁|k̃⟩ - [i|k̃]|i⟩ + [i|i]|k̃⟩ - [j|k̃]|j⟩ + [j|j]|k̃⟩ - [l|k̃]|l⟩ + [l|l]|k̃⟩
[1, 2]  =  + (l|k)𝐈₁|l̃⟩ - [i l|k i]𝐈₁|l̃⟩ + [i l|i k]𝐈₁|l̃⟩ - [j l|k j]𝐈₁|l̃⟩ + [j l|j k]𝐈₁|l̃⟩ - [l|l̃]|k⟩ + [l|k]|l̃⟩
, OrbitalEquation(l̃):
[2, 1]  =  + (k|l)𝐈₁|k̃⟩ - [i k|l i]𝐈₁|k̃⟩ + [i k|i l]𝐈₁|k̃⟩ - [j k|l j]𝐈₁|k̃⟩ + [j k|j l]𝐈₁|k̃⟩ - [k|k̃]|l⟩ + [k|l]|k̃⟩
[2, 2]  =  + (i|i)𝐈₁|l̃⟩ + (j|j)𝐈₁|l̃⟩ + (k|k)𝐈₁|l̃⟩ + ĥ|l̃⟩ - G(i,j)𝐈₁|l̃⟩ + F(i,j)𝐈₁|l̃⟩ - G(i,k)𝐈₁|l̃⟩ + F(i,k)𝐈₁|l̃⟩ - G(j,k)𝐈₁|l̃⟩ + F(j,k)𝐈₁|l̃⟩ - [i|l̃]|i⟩ + [i|i]|l̃⟩ - [j|l̃]|j⟩ + [j|j]|l̃⟩ - [k|l̃]|k⟩ + [k|k]|l̃⟩
], Any[(i|i), (j|j), (l|l), G(i,j), F(i,j), G(i,l), F(i,l), G(j,l), F(j,l), [i|i]  …  [j k|l j], [j k|j l], [k|k̃], [k|l], (k|k), G(i,k), F(i,k), G(j,k), F(j,k), [k|k]])

We can investigate the MCEquationSystem object eqs a bit. It consists of two coupled equations:

julia> eqs.equations
2-element Array{EnergyExpressions.OrbitalEquation{Symbol,SparseArrays.SparseMatrixCSC{LinearCombinationEquation,Int64}},1}:
OrbitalEquation(k̃):
[1, 1]  =  + (i|i)𝐈₁|k̃⟩ + (j|j)𝐈₁|k̃⟩ + (l|l)𝐈₁|k̃⟩ + ĥ|k̃⟩ - G(i,j)𝐈₁|k̃⟩ + F(i,j)𝐈₁|k̃⟩ - G(i,l)𝐈₁|k̃⟩ + F(i,l)𝐈₁|k̃⟩ - G(j,l)𝐈₁|k̃⟩ + F(j,l)𝐈₁|k̃⟩ - [i|k̃]|i⟩ + [i|i]|k̃⟩ - [j|k̃]|j⟩ + [j|j]|k̃⟩ - [l|k̃]|l⟩ + [l|l]|k̃⟩
[1, 2]  =  + (l|k)𝐈₁|l̃⟩ - [i l|k i]𝐈₁|l̃⟩ + [i l|i k]𝐈₁|l̃⟩ - [j l|k j]𝐈₁|l̃⟩ + [j l|j k]𝐈₁|l̃⟩ - [l|l̃]|k⟩ + [l|k]|l̃⟩

OrbitalEquation(l̃):
[2, 1]  =  + (k|l)𝐈₁|k̃⟩ - [i k|l i]𝐈₁|k̃⟩ + [i k|i l]𝐈₁|k̃⟩ - [j k|l j]𝐈₁|k̃⟩ + [j k|j l]𝐈₁|k̃⟩ - [k|k̃]|l⟩ + [k|l]|k̃⟩
[2, 2]  =  + (i|i)𝐈₁|l̃⟩ + (j|j)𝐈₁|l̃⟩ + (k|k)𝐈₁|l̃⟩ + ĥ|l̃⟩ - G(i,j)𝐈₁|l̃⟩ + F(i,j)𝐈₁|l̃⟩ - G(i,k)𝐈₁|l̃⟩ + F(i,k)𝐈₁|l̃⟩ - G(j,k)𝐈₁|l̃⟩ + F(j,k)𝐈₁|l̃⟩ - [i|l̃]|i⟩ + [i|i]|l̃⟩ - [j|l̃]|j⟩ + [j|j]|l̃⟩ - [k|l̃]|k⟩ + [k|k]|l̃⟩

The first equation consists of the following terms:

julia> eqs.equations.terms
6-element Array{Pair{Int64,Array{EnergyExpressions.MCTerm,1}},1}:
0 => [MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, 1, 𝐈₁, :k̃, ), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, 1, 𝐈₁, :k̃, ), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, 1, 𝐈₁, :k̃, ), MCTerm{Int64,OneBodyHamiltonian,Symbol}(1, 1, 1, ĥ, :k̃, Int64[]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, -1, 𝐈₁, :k̃, ), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, 1, 𝐈₁, :k̃, ), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, -1, 𝐈₁, :k̃, ), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, 1, 𝐈₁, :k̃, ), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, -1, 𝐈₁, :k̃, ), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, 1, 𝐈₁, :k̃, ), MCTerm{Int64,ContractedOperator{1,2,1,Symbol,CoulombInteraction,Symbol},Symbol}(1, 1, -1, [i|k̃], :i, Int64[]), MCTerm{Int64,ContractedOperator{1,2,1,Symbol,CoulombInteraction,Symbol},Symbol}(1, 1, -1, [j|k̃], :j, Int64[]), MCTerm{Int64,ContractedOperator{1,2,1,Symbol,CoulombInteraction,Symbol},Symbol}(1, 1, -1, [l|k̃], :l, Int64[]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 2, 1, 𝐈₁, :l̃, ), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 2, -1, 𝐈₁, :l̃, ), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 2, 1, 𝐈₁, :l̃, ), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 2, -1, 𝐈₁, :l̃, ), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 2, 1, 𝐈₁, :l̃, )]
10 => [MCTerm{Int64,ContractedOperator{1,2,1,Symbol,CoulombInteraction,Symbol},Symbol}(1, 1, 1, [i|i], :k̃, Int64[])]
19 => [MCTerm{Int64,ContractedOperator{1,2,1,Symbol,CoulombInteraction,Symbol},Symbol}(1, 2, 1, [l|k], :l̃, Int64[])]
11 => [MCTerm{Int64,ContractedOperator{1,2,1,Symbol,CoulombInteraction,Symbol},Symbol}(1, 1, 1, [j|j], :k̃, Int64[])]
12 => [MCTerm{Int64,ContractedOperator{1,2,1,Symbol,CoulombInteraction,Symbol},Symbol}(1, 1, 1, [l|l], :k̃, Int64[])]
18 => [MCTerm{Int64,ContractedOperator{1,2,1,Symbol,CoulombInteraction,Symbol},Symbol}(1, 2, -1, [l|l̃], :k, Int64[])]

where the MCTerm objects indicate which components of the mixing coefficient vector $\vec{c}$ need to be multiplied, and all the integers are pointers to the list of common integrals:

julia> eqs.integrals
32-element Array{Any,1}:
(i|i)
(j|j)
(l|l)
G(i,j)
F(i,j)
G(i,l)
F(i,l)
G(j,l)
F(j,l)
[i|i]
[j|j]
[l|l]
(l|k)
[i l|k i]
[i l|i k]
[j l|k j]
[j l|j k]
[l|l̃]
[l|k]
(k|l)
[i k|l i]
[i k|i l]
[j k|l j]
[j k|j l]
[k|k̃]
[k|l]
(k|k)
G(i,k)
F(i,k)
G(j,k)
F(j,k)
[k|k]

From this we see that the 𝐈₁ (one-body identity operator) contribution to |k̃⟩ can be written as a linear combination of |k̃⟩ and |l̃⟩, weighted by different components of the mixing coefficient vector $\vec{c}$ and various other integrals. This is all the information necessary to set up an efficient equation solver.

## Implementation

EnergyExpressions.MCTermType
MCTerm(i, j, coeff, operator, source_orbital, integrals=[])

Represents one term in the multi-configurational expansion. i and j are indices in the mixing-coefficient vector c (which is subject to optimization, and thus has to be referred to), coeff is an additional coefficient, and integrals is a list of indices into the vector of common integrals, the values of which should be multiplied to form the overall coefficient.

source
EnergyExpressions.OrbitalEquationType
OrbitalEquation(orbital, equation,
one_body, direct_terms, exchange_terms, source_terms)

Represents the integro-differential equation for orbital, expressed as a linear combination of the different terms, with pointers to the list of common integrals that is stored by the encompassing MCEquationSystem object.

source
EnergyExpressions.MCEquationSystemType
MCEquationSystem(equations, integrals)

Represents a coupled system of integro-differential equations, resulting from the variation of a multi-configurational EnergyExpression, with respect to all constituent orbitals. All integrals that are in common between the equations need only be computed once per iteration, for efficiency.

source
EnergyExpressions.pushifmissing!Function
pushifmissing!(vector, element)

Push element to the end of vector, if not already present. Returns the index of element in vector.

source