Is there chaos in the SCF equations?

 

This article is a working preliminary draft, NOT yet submitted for peer review. Leave your comments on the discussion page (talk page) or contact the First Author, MartinY, at their talk page or by email.

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 c². 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

Wikipedia:Image:LogisticMap BifurcationDiagram.png

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 (CH₂O), 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 2ⁿ 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 CH₂O

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 sp² 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 q₁, q₂, q₃ and q₄ 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.