SCF Convergence and Chaos Theory

What are Nonlinear Systems

		x = f(x)

The branch of mathematics called "chaos theory" which has been created and explored over the last decade is the study of equations and systems of equations of this type. To some extent the mathematicians are tacking new names onto behavior that computational chemists have been seeing for years, but they are also learning more about why these systems behave the way they do. It is hoped that, in time, the work being done by mathematicians will provide the insight necessary to create better SCF convergence methods.

The purpose of this document is two fold. First it attempts to make the connections between the field of chaos theory and SCF calculations. Second it gives suggestions on how to fix calculations with convergence problems.

What Nonlinear Systems Do

1. After a number of iterations, the value returned by the equation may be the value that was put in on that iteration. This is what chemists desire - a converged solution.
2. The values could be almost repeating but not quite. These are called Lorenz attractor systems in chaos theory.
3. The values produced from one iteration to the next may oscillate between 2 values, 4 values or any other power of 2. (I am not aware of any examples other than powers of 2.)
4. The values produced may be random within some fixed range. Random number generators use this property intentionally, but it is rather annoying when an SCF calculation does so.
5. The values produced may be random and not bounded within any upper or lower limits.

It stands to reason that the total energy may follow any of these patterns during SCF iterations. We have encountered oscillating and random behavior in the convergence of open shell transition metal compounds but have never tried to determine if the random values were bounded or not.

How to Control Nonlinear Behavior

Another way of changing the constants in an SCF calculation is to change the geometry a bit. Often pulling a bond length a bit shorter than expected is effective (say making the length 90% of the expected value). lengthening bond lengths a bit and avoiding eclipsed or gauche conformations are second and third best. Once you have a converged wave function move the geometry back where you want it and use the converged function as the initial guess.

Changing the initial value of variables is equivalent to using a different initial guess in an SCF calculation. The best initial guess is usually a converged SCF calculation for a different state of the same molecule or a slightly different geometry of the same molecule. This can be a very effective way to circumvent convergence problems. In the worst case it may be necessary to construct an initial guess by hand in order to ensure that the nodal properties of all of the orbitals are correct for the desired electronic state of the molecule. The construction of the virtual orbitals as well as the occupied orbitals can have a significant effect on convergence. Chaos theorists will try many starting points and color code a plot by which solution is obtained for each starting point, this is called fractal geometry (i.e. the famous Mandelbrot diagram).

There are quite a number of ways to effectively change the equation in an SCF calculation. These include switching computation methods, using level shifting and using forced convergence methods.

Switching between Hartree-Fock (HF), semiempirical, generalized valence bond (GVB), multi-configuration self consistent field (MCSCF), complete active space self consistent field (CASSCF) and Moller-Plesset calculations (MPn) will change the convergence properties. Configuration interaction (CI) and coupled cluster (CC) calculations usually start with an SCF calculation, thus will not circumvent problems with an SCF. In general higher levels of theory tend to be harder to converge. Ease of convergence as well as calculation speed are why lower level calculations are usually used to generate the initial guess for higher level calculations.

Oscillating convergence in an SCF calculation is usually an oscillation between wave functions that are close to different states or a mixing of states. Thus oscillating convergence can often be helped by level shifting, which artificially raises the energies of the virtual orbitals. Level shifting may or may not help in cases of random convergence.

Most programs will stop trying to converge a problem after a certain number of iterations. In a few rare cases the wave function will converge if given more than the default number of iterations.

Most SCF programs do not actually compute orbitals from the previous iteration orbitals in the way that is described in introductory descriptions of the SCF method. Most programs use a convergence excelleration method which is designed to reduce the number of iterations necessary to converge to a solution. The method of choice is usually Pulay's DIIS method (Direct Inversion of the Iterative Subspace). Some programs also give the user the capability to modify the DIIS method, such as putting in a dampening factor. These modifications can be useful for fixing convergence problems but they usually require a significant amount of experience to know how best to modify the procedure. Turning off the DIIS extrapolation can help a calculation converge, but usually requires more iterations.

Some programs contain alternative convergence methods that are designed to force even the most difficult problems to converge. These methods are often called direct covergence or quadratic convergence methods. While these methods almost always work, they often require a very large number of iterations and thus a very large amount of CPU time.

What to Try First

1. Try a different initial guess (using the "guess" keyword in Gaussian). For an open shell system, try converging the closed shell ion of the same molecule then using that as an initial guess for the open shell calculation. Adding electrons may give more reasonable virtual orbitals but as a general rule of thumb cations are easier to converge than anions.
2. Try level shifting ( "SCF=Vshift" in Gaussian).
3. Try changing the geometry. First slightly shortening a bond length then slightly lengthening a bond length then shifting the conformation a bit.
4. Consider trying a different basis set.
5. Consider doing the calculation at a different level of theory. This isn't always practical, but beyond this point the increased number of iterations may make the computation time as much as using a higher level of theory anyway.
6. Turn off the DIIS extrapolation ( "SCF=NoDIIS" in Gaussian ). You should probably give the calculation more iterations along with this.
7. Give the calculation more SCF iterations ( "SCF(MaxCyc=N)" where N is the number of iterations in Gaussian ). This seldom helps but the next option often uses so many iterations that its worth a try.
8. Use a forced convergence method. (In Gaussian "SCF=QC" is usually the best but on rare occasions "SCF=DM" will be faster). Don't forget to give the calculation an extra thousand iterations or so. The wave function obtained by these methods should be tested to make sure it is a minimum and not just a stable point ( see the "stable" keyword in Gaussian ).
9. See if the software documentation suggests any other ways to change the DIIS method. You may well have to run hundreds of calculations to get enough experience with this to know what works when and how much to change it by.

Further Information

There is a very small discussion of handling convergence problems in
T. Clark "A Handbook of Computational Chemistry" Wiley-Interscience (1985)

A good introduction to chaos theory is
J. Gleick "Chaos: Making a New Science" Viking (1987)

A more mathematical treatment of chaos theory is
S. H. Strogatz "Nonlinear Dynamics and Chaos With Applications to Physics, Biology, Chemistry and Engineering" Addison Wesley (1994)

If you have any suggestions for improvements or additions to this document, please contact me at