Is there chaos in the SCF equations?Edit
Principal AuthorEdit
Martin Grayson MartinY 09:31, 10 September 2006 (UTC)
Preface Edit
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.
AbstractEdit
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.
IntroductionEdit
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^{[1]} ^{[2]} ^{[3]} 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^{2}. 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 ^{[4]} level shifters ^{[5]} and DIIS(diis). ^{[6]} 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 CascadeEdit
Wikipedia:Image:LogisticMap BifurcationDiagram.png
ResultsEdit
It was whilst using the SYSMO ^{[7]} with no convergence aids that the following behaviour was observed.
Without using level-shifting or Hartree Damping methanal or formaldehyde (CH_{2}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 ^{[8]}.) 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^{n} 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 ^{[9]} ^{[10]} . 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_{2}OEdit
Table 1 -Damping factor k and number of iterations to convergenceEdit
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.883Edit
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 upEdit
k | E^{1} | E^{2} |
---|---|---|
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 oscillationEdit
k | E^{1} | E^{2} | E^{3} | E^{4} |
---|---|---|---|---|
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 MethodEdit
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^{2} 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_{1}, q_{2}, q_{3} and q_{4} 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 $.
ConclusionEdit
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.
ReferencesEdit
- ↑ Peterson, I. (1988). "{{{title}}}". The Mathematical Tourist, p150, Freeman, New York.
- ↑ Devany, R. L., (1987). "Chaotic Bursts in Nonlinear Dynamical Systems". Science 235, 342-345.
- ↑ Geick, J. (1988). "{{{title}}}". Chaos, Heinemann.
- ↑ Firscht, D. and Pickup, B. T. (1977). "{{{title}}}". Int. J. Quantum Chem. 12, 765..
- ↑ Guest, M. F. and Saunders, V. R. (1974). "{{{title}}}". Mol. Phys. 28, 819..
- ↑ Pulay, P. (1982). "{{{title}}}". Chem. Phys. Lett. 74, 2384..
- ↑ Lazzerreti P. and Zanasi P.. "{{{title}}}". The Modena Molecular Properties System.
- ↑ W. J. Hehre, R. Ditchford and J. A. Pople (1972). "{{{title}}}". J. Chem. Phys., 56, 2257.
- ↑ Wheland, G. W.; Mann D. E. (1949). "{{{title}}}". J. Chem. Phys., 17, 264.
- ↑ Streitwieser, A. (1960). "{{{title}}}". J. Am. Chem. Soc., 82, 4123.