# 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}$:

Field-free, one-body Hamiltonian, i.e. $\operator{A}=\hamiltonian_0$. This is the contribution from the orbital $\ket{\chi}$ to itself. Rank 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.

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”.

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.

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()
ĥ + ĝ
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̃⟩ + ĥ|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̃⟩ + ĥ|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̃⟩ + ĥ|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̃⟩ + ĥ|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[1].terms
6-element Array{Pair{Int64,Array{EnergyExpressions.MCTerm,1}},1}:
0 => [MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, 1, 𝐈₁, :k̃, [1]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, 1, 𝐈₁, :k̃, [2]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, 1, 𝐈₁, :k̃, [3]), MCTerm{Int64,OneBodyHamiltonian,Symbol}(1, 1, 1, ĥ, :k̃, Int64[]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, -1, 𝐈₁, :k̃, [4]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, 1, 𝐈₁, :k̃, [5]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, -1, 𝐈₁, :k̃, [6]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, 1, 𝐈₁, :k̃, [7]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, -1, 𝐈₁, :k̃, [8]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 1, 1, 𝐈₁, :k̃, [9]), 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̃, [13]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 2, -1, 𝐈₁, :l̃, [14]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 2, 1, 𝐈₁, :l̃, [15]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 2, -1, 𝐈₁, :l̃, [16]), MCTerm{Int64,IdentityOperator{1},Symbol}(1, 2, 1, 𝐈₁, :l̃, [17])]
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.MCTerm`

— Type`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.

`EnergyExpressions.OrbitalEquation`

— Type```
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.

`EnergyExpressions.orbital_equation`

— Function`orbital_equation(E::EnergyExpression, orbital, integrals::Vector)`

Generate the `OrbitalEquation`

governing `orbital`

by varying the `EnergyExpression`

`E`

, and storing common expressions in `integrals`

.

`EnergyExpressions.MCEquationSystem`

— Type`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.

`EnergyExpressions.pushifmissing!`

— Function`pushifmissing!(vector, element)`

Push `element`

to the end of `vector`

, if not already present. Returns the index of `element`

in `vector`

.