# Self-Consistent Iteration

`SCF.scf!`

— Function`scf!([fun!, ]fock[; ω=0, max_iter=200, tol=1e-8, verbosity=1])`

Perform the SCF iterations for the `fock`

operator. One iteration is performed by `scf_iteration!`

. The optional `fun!(P̃,c̃)`

argument allows for extra operations to be performed every SCF cycle (such as plotting, etc).

`ω`

is a relaxation parameter; the orbitals and mixing coefficients are updated as `wᵢ₊₁ = (1-ω)w̃ + ωwᵢ`

where `w̃`

is the solution in the current iteration and `wᵢ`

the previous solution. The default (`ω=0`

) is to only use this solution. A larger value of `ω`

makes it easier to achieve convergence in oscillatory problems, but a too large value decreases the convergence rate. It is therefore desireable to reduce the value of `ω`

if the convergence is deemed to be monotonous, the criterion for which is whether the change over the `monotonous_window`

last iterations is steadily decreasing. If this is the case, `ω`

will be reduced by multiplying by `ωfactor`

. Conversely, if the change is non-monotonous, `ω`

will be increased, but not beyond `ωmax`

.

The SCF procedure continues until either the amount of iterations equals `max_iter`

or the change in the coefficients is below `tol`

.

`SCF.scf_iteration!`

— Function```
scf_iteration!(fock, P, H, c, okws,
[; verbosity=0, method=:arnoldi, σ=nothing, tol=1e-10,
update_mixing_coefficients=true])
```

Perform one step of the SCF iteration. This will

Improve each of the orbitals in turn, using the values of the orbitals

`P`

and mixing coefficients`c`

from the previous step as input.Solve the secular problem to find new values for the mixing coefficients,

`c`

, unless`update_mixing_coefficients==false`

.

`okws`

is a vector of `OrthogonalKrylovWrapper`

s, one for each orbital.

`SCF.solve_orbital_equation!`

— Function`solve_orbital_equation!(Pj, okw, method, tol)`

Solve the orbital equation for the `j`

th orbital, by solving the eigenproblem of the `OrthogonalKrylovWrapper`

operator `okw`

, and store the result in the radial orbital vector of expansion coefficients `Pj`

. The solution is computed using `method`

; valid choices are

`:arnoldi`

, which uses ArnoldiMethod.jl. Really small grid spacings/large grid extents can be problematic with this method, due to high condition numbers.`:lobpcg`

, which uses IterativeSolvers.jl. Too coarse grids can cause never-ending loops or errors from the Cholesky decomposition, complaining about non-positive definite subproblems.

Both methods are controlled by the stopping tolerance `tol`

.

`SCF.solve_secular_problem!`

— Function`solve_secular_problem!(H, c, fock)`

Form the energy matrix, store it in `H`

, and then solve the secular problem `Hc = Ec`

for the lowest eigenvalue.

## KrylovWrapper

The orbital improvements are performed via diagonalization of the integro-differential equations for each orbital. The diagonalization procedure is the Arnoldi method, which requires the repeated action of the Hamiltonian on a vector. The `KrylovWrapper`

object is used to wrap Hamiltonians which are not simple `<:AbstractMatrix`

objects, but which require a specialized implementation of `mul!`

to act on a `<:AbstractVector`

of coefficients. An example is the kind of Hamiltonians implemented in the Atoms.jl library, which are based on the function space algebra from ContinuumArrays.jl.

`SCF.KrylovWrapper`

— Type`KrylovWrapper(hamiltonian)`

Proxy object used in the Krylov iterations, during orbital improvement. This is useful, since `hamiltonian`

may be defined to act on objects such as vectors living in function spaces (as e.g. implemented using ContinuumArrays.jl), whereas the SCF iterations act on the coefficients directly.

`Base.size`

— Method`size(::KrylovWrapper)`

Returns the dimension of the `KrylovWrapper`

. For Hamiltonians which are not `<:AbstractMatrix`

, this needs to be overloaded.

Missing docstring for `LinearAlgebra.mul!(y::V₁, A::KrylovWrapper{T,Hamiltonian}, x::V₂) where {V₁,V₂,T,Hamiltonian<:AbstractMatrix}`

. Check Documenter's build log for details.

## Gershgorin's circle theorem

Gershgorin's circle theorem can be used to estimate the bounds on the eigenspectrum of a matrix. It is always correct, but only a useful estimate when the matrix is diagonally dominant. Furthermore, `gershgorin_bounds`

is limited to the case of Hermitian matrices, whose eigenvalues lie on the real line. If an experimental energy is known, that is almost always a better lower bound to use for the Shift-and-invert strategy.

`SCF.gershgorin_disc`

— Function`gershgorin_disc(i, row)`

Return the Gershgorin disc for the `i`

th `row`

of a matrix.

`SCF.gershgorin_discs`

— Function`gershgorin_discs(A)`

Compute all Gershgorin discs of a matrix `A`

.

`SCF.gershgorin_bounds`

— Function`gershgorin_bounds(A)`

For a Hermitian matrix `A`

, compute the bounds on the real axis of the eigenspectrum. NB that this bound will only be useful for diagonally dominant `A`

s.