Is there chaos in the SCF equations?

===Is there chaos in the SCF equations?===
===Principal Author===
Martin Grayson MartinY 09:31, 10 September 2006 (UTC)
===Preface ===
This is an in progress article to kick off
Interpretations in the Physical and Computational Sciences
and for the first author to develop an initial referencing style.
It will contain thoughts about some data which has been awaiting a full
interpretation for several years now.
There is one weak piece of graphics and the Wheland-Mann equation section at the end is still being worked out though I think these old equations might provide an interesting student computer project even if their practical importance is small in the age of fast computers.
===Abstract===
An example of an SCF cycle which might exhibit chaos has
been observed. This is clearly equivalent in behaviour to the
logistic equation, however no Feigenbaum numbers have been found.
An analytically tractable and simple example of chaos is being searched for in the equations of the \omega-technique of Wheland and Mann
which is a kind of hybrid of Hückel Theory with a self consistent iteration of the secular equations.
===Introduction===
The iterative solution of the SCF equations for a molecular system can contain
regions of mathematical chaos.
An example series of iterations
have been found which have similarities to the Feigenbaum Cascade in
the logistic equation.
This rarely impinges on the convergence of the SCF procedure
with modern convergence aids but it may be of interest to explore this
chaotic region for theoretical interest. Strangely the chaotic
region is adjacent to the region of fastest convergence
It is now well known that
very simple feedback systems such as the logistic equation
exhibit complex chaotic behaviour when iterated.
There are regions of stable convergence, regions of
monotonic divergence to infinity but also regions of bounded chaos,
or cycling repetitions of groups of density matrices.
The logistic equation itself gives the striking Feigenbaum Cascade where
continuous bifurcation gives the pattern in figure 1.
Whilst running SCF
computer programs on carbonyl containing groups
behaviour quite similar to this was observed in the SCF iterations.
This is not entirely surprising as the SCF equations are an example of a non-linear feedback system.
The Hamiltonian is quartic in the coefficients of the expansion
orbitals as the coefficients occur in the Fock operator as well
as in the wavefunction, whereas at each iteration the energy
evaluation is quadratic for the fixed density matrix.
{\bf \rho} = \sum_{k} {\bf c}_{k}{\bf c}_{k}^{\dagger}
E = tr {\bf \rho h} + \frac {1}{2} tr {\bf \rho G(\rho)}
This is analogous to iterating a quadratic equation in c2.
The nature of the interation process might not be seen with a modern electronic
structure package as they automatically change the conditions of
iteration according to algorithms invoking
damping factors
level shifters
and DIIS(diis).
This is done to ensure and accellerate convergence but it sometimes fails.
Normally changing some of the iteration parameters and a repeat run
gives convergence but
the various convergence environments switched in by
a modern MO package might obscure the reasons for non-convergence.
====Figure 1 - The Feigenbaum Cascade====
===Results===
It was whilst using the SYSMO
with no convergence
aids that the following behaviour was observed.
Without using level-shifting or Hartree Damping methanal
or formaldehyde (CH2O), with a split valence basis will
not converge. (The geometry used was the experimental one obtained by Bowen.
The basis set was the Pople 6-31G set
.)
If one plots the convergence against damping
parameter it is seen that quickest convergence requires a \Delta
of about 0.7 / 0.8 (Table 1) but above 0.8 a numerical catastrophe occurs.
As \Delta increases to 0.8825 convergence becomes slower and
slower until at 0.883 a burst into apparent chaos occurs.
a whole sequence of densities with different energies occurs (Table 2).
However above this an oscillation between two density matrices occurs (Table 3).
Above 1.0 a four density stable iteration occurs (Table 4).
These densities do not seem to be associated with physically realistic states as
they are not quantised and smoothly deform as a result of varying
\Delta.
oscillation between two density matrices is common in the SCF method.
Whether this is also associated systematically with 2n fold oscillations is open to question.
Later versions of the SYSMO SCF program give the same phenomena but at different values
of \Delta. The behaviour is clearly quite sensitively dependant upon
the nature of the SCF process.
A clearer and simpler example may be extractable from the
\omega-technique of Wheland and Mann
.
This is a kind of proto-SCF where one iterates Hückel Theory to optimise a parameter\omega.
For a collection of \pi-atoms, (e.g. a double bonded hydrocarbon), this will be a characteristic polynomial of order n,
where n is the number of \pi-atoms.
This can be a good source of iteratable equations of order greater than 4 which relate to physical phenomena.
===Tables - the data is at the Bowen molecular geometry of CH2O===
====Table 1 -Damping factor k and number of iterations to convergence====
k 	Number of iterations
0.1 	114
0.2 	58
0.3 	39
0.4 	29
0.5 	23
0.6 	19
0.7 	16
0.8 	16
0.883	bifurcation
 	to chaos
====Table 2 - Chaotic Iterations at k=0.883====
Iteration	Energy 	Density matrix
Number 	 	difference
50 	-107.0477669643	2.6736259862
51 	-109.7634079064	2.7156409421
52 	-107.0909997342	2.6724081722
53 	-109.8102354135	2.7192356793
54 	-107.1408308300	2.6694045835
55 	-109.8623866752	2.7215558452
56 	-107.1985847400	2.6638019352
57 	-109.9206495208	2.7220647808
58 	-107.2662657388	2.6543837819
59 	-109.9862982700	2.7200325311
60 	-107.3470953732	2.6392028968
61 	-110.0614625099	2.7143671367
62 	-107.4465949066	2.6148676033
63 	-110.1498366421	2.7032417355
64 	-107.5750089843	2.5748276579
65 	-110.2582062757	2.6831972914
66 	-107.7535317438	2.5046745319
67 	-110.4002025792	2.6466708354
68 	-108.0334179312	2.3667846479
69 	-110.6071786583	2.5737607270
70 	-108.5680930677	2.0390855905
====Table 3 - at 0.884 a two fold oscillation becomes set up====
k 	E1 	E2
0.884	-105.969714	-104.838991
0.9 	-105.411554	-104.276957
1.0 	-99.310476 	-98.341802
1.02 	-98.317586 	-97.434605
====Table 4 - above 1.02 a fourfold oscillation====
k 	E1 	E2 	E3 	E4
1.1	-93.495382	-93.480168	-93.444864	-93.145210
1.2	-85.894695	-85.448206	-84.489825	-84.438016
1.3	-73.685920	-73.170117	-66.948456	-66.846642
1.4	-53.576210	-52.945134	-43.893513	-43.734596
1.5	-17.154841	-16.382111	-4.239830 	-4.015124
1.5	positive 	energies 	 	
===The Wheland-Mann Method===
It may be possible that the
\omega-technique of Wheland and Mann
may be applied with values of \omega which exibit chaotic behaviour.
The Hückel problem in this formulation is a feedback solution of a
characteristic polynomial of order n, where n is the number of atoms.
The Wheland-Mann method modifies the Hückel theory \alpha
parameter thus:
 \alpha_\mu = \alpha_0 + ( 1 - q_\mu) \omega \beta 
Where q is the value of 2c^2 at each atom, i.e. for an sp2 carbon atom in an alternant hydrocarbon it is 1 and standard Hückel theory is recovered. \omega is conventionally 1.4 to give agreement with experiment.
The Hückel problem is then solved and the system is iterated until the valus of q_\mu stop changing, so here is the analogy with the SCF procedure. \omega is an optimisable electron repulsion parameter which damps down the tendency of HT to generate too much charge separation.
This model however has some of the the analytical simplicity of HT and the iterative procedure of SCF and so might give a simple model of a divergence / convergence / chaos situation.
Take the example of furan which has 5 \pi atoms. This leads to a characteristic polynomial which is a quintic with one root and two conjugate pairs by symmetry. The equations are well behaved near the conventional value of \omega which has a sort of physical meaning as an average electron-electron repulsion parameter. What values of \omega will misbehave and does this misbehaviour have any conceivable physical meaning? Furan however generates a quintic equation and three unique values of q, subject to the constraints that two are equal to others and the total adds up to zero for a neutral molecule. This is not likely to be analytically tractable.
An attempt to was found to find a smaller example which could be investigated by pen and paper rather than an algebra system. Though furan has quite a bit of symmetry it leads to a quintic equation. As any of these systems are iterated nonlinear equations they might generate an object like the Wikipedia:Mandelbrot_set but in these examples the chemistry constrains us not to use complex numbers as possible values of observables as all experimental observables, (except perhaps lifetimes), must be real. However there might be pieces of an object called a Wikipedia:Cantor_set along the real axis. 
C=C-C=O was investigated to see if it was analytically tractable, (butadiene is no good because the charges are all zero, being an alternant). It is an example in Streitweiser with the secular determinant:
\begin{vmatrix}
x & 1 & 0 & 0 \\
1 & x & 1 & 0 \\
0 & 1 & x & \sqrt 2 \\
0 & 0 & \sqrt 2 & x + 2 \\
\end{vmatrix}
= 0
This expands via
\begin{vmatrix}
x & 1 & 0 \\
 1 & x & \sqrt 2 \\
 0 & \sqrt 2 & x + 2 \\
\end{vmatrix}
-
\begin{vmatrix}
1 & 0 & 0 \\
0 & x & \sqrt 2 \\
0 & \sqrt 2 & x + 2 \\
\end{vmatrix}
to
 x^4 + 2 x^3 -4x^2 -4x +2 = 0 
with roots -2.826 \beta, -1.152 \beta, 0.386 \beta ~~~ {\rm and} ~~~~
1.593 \beta . Remembering that \beta is a negative energy 0.386 and 1.593 are the occupied orbitals.
Now unfortunately when we convert this to the Wheland-Mann equations the quartic has lots of products of q1, q2, q3 and q4 as they are independant other than being subject to the constraint that they all add up to the number of \pi electrons. 
\begin{vmatrix}
x + (1-q_1)\omega & 1 & 0 & 0 \\
1 & x + (1-q_2)\omega & 1 & 0 \\
0 & 1 & x + (1-q_3)\omega & \sqrt 2 \\
0 & 0 & \sqrt 2 & x + 2 + (1-q_4)\omega \\
\end{vmatrix}
= 0
Expanding this out, though easy enough with an algebra system does not lead to a very tractable equation. 
It looks as though the best way to tackle this problem is to program the Wheland-Mann equations and run the program for many values of \omega.
===Conclusion===
What has been presented is a mathematical curiosity. The
Feigenbaum numbers have been searched for in the data without success.
However we are clearly seeing the same chaos as exhibited
by the logistic equation in a more complex form. The observation that chaos
ensues in a region immediately adjacent to the region of swiftest
convergence may be of interest.
The obsolete Wheland-Mann method has been examined to see if it is a good source of meaningful iterable quartic equations. This work has caused the first author to re-examine the usefulness of the WM method for teaching iteration as it has links back to familiar organic chemistry but it seems on balance it a complex digression 
too far and performing a small SCF on helium and practicing the Wikipedia:Newton-Raphson method would be a better use of student time.
===References===
Category:Interpretations in the Physical and Computational Sciences submitted article