Convergence in HartreeFock Calculations
The convergence of HartreeFock SCF calculations depends on a series of options. All of these have sensible defaults which lead to usually well behaved SCF calculations. The most obvious parameter is the number of SCF cycles which can be set with
SCF=(MaxCycle=N)
the default being N=64. This default is sufficient for most systems, but significantly higher numbers are required for systems with closely spaced orbitals such as some transition metal complexes. The convergence criterion for rms deviations in electron density matrix elements of the SCF algorithm can be set directly by
SCF=(Conver=n)
the default for n being 8 for geometry optimization and frequency jobs and n=4 for single point calculations. The latter default is rarely appropriate and should be adjusted to the same value as for geometry optimizations. Using SCF=tight is synonymous with SCF=(Conver=8) and SCF=SinglePoint is synonymous with SCF=(Conver=4). The following table illustrates how much the number of SCF cycles varies with the SCF convergence criterion for the HF/STO3G calculation on formaldehyde :
Conver=n  optimization cycles 
final energy (Hartree) 

4  6  112.354346245 
5  7  112.354347141 
6  8  112.354347141 
7  9  112.354347141 
8  10  112.354347141 
9  11  112.354347141 
One can clearly see that the total energy of the system is converged to within 10^{9} Hartree already with SCF=(Conver=5) and that each increase in the convergence criterion is accompanied by an increase in the number of SCF cycles. For larger and for less well behaved systems the consequences of changing the convergence criterion are usually much more pronounced. Calculations on large systems can be accelerated substantially by lowering various cutoffs and convergence parameters using SCF=sleazy. This is particularly helpful in geometry optimizations when still far away from the minimum.
The initial guess used to start the SCF calculation is of great importance for both the efficiency as well as the final result obtained in the SCF calculation. The latter point is of particular importance in symmetric species in which the symmetry of the initial guess decides over the symmetry of the converged SCF wavefunction. Whether the symmetry of the initial guess is retained during SCF cycles depends on the version of Gaussian in use. Newer versions do NOT retain the symmetry of the initial guess unless SCF=Symm is specified. Using the NH_{2} radical as an example, we will see how either the ^{2}A_{1} or the ^{2}B_{1} states can be obtained at the ROHF/STO3G level of theory depending on the initial guess being used. Starting the SCF calculation with the default (Harris) guess with:
#ROHF/STO3G scf=(symm,tight)
rohf/STO3G amide radical NH2, 2B1 state
0 2
N1
H2 1 r2
H3 1 r2 2 a3
r2=1.05726795
a3=100.0637224

will give the following result:
Harris functional with IExCor= 205 diagonalized for initial guess.
ExpMin= 1.69D01 ExpMax= 9.91D+01 ExpMxC= 9.91D+01 IAcc=1 IRadAn= 1 AccDes= 1.00D06
HarFok: IExCor= 205 AccDes= 1.00D06 IRadAn= 1 IDoV=1
ScaDFX= 1.000000 1.000000 1.000000 1.000000
Initial guess orbital symmetries:
Occupied (A1) (A1) (B2) (A1) (B1)
Virtual (A1) (B2)
The electronic state of the initial guess is 2B1.
Requested convergence on RMS density matrix=1.00D08 within 128 cycles.
Requested convergence on MAX density matrix=1.00D06.
Requested convergence on energy=1.00D06.
No special actions if energy rises.
Keep R1 and R2 integrals in memory in canonical form, NReq= 419815.
SCF Done: E(ROHF) = 54.8368134090 A.U. after 7 cycles
Convg = 0.5123D11 V/T = 2.0047
S**2 = 0.7500
Annihilation of the first spin contaminant:
S**2 before annihilation 0.7500, after 0.7500
**********************************************************************
Population analysis using the SCF density.
**********************************************************************
Orbital symmetries:
Occupied (A1) (A1) (B2) (A1) (B1)
Virtual (A1) (B2)
The electronic state is 2B1.
The Harris functional guess predicts a highest (singly) occupied orbital of "(B1)" symmetry leading to a doublet ^{2}B_{1} state in the converged wavefunction. The SCF algorithm retains the number of occupied orbitals in each irreducible representation (3*A1, 1*B2, 1*B1). The first unoccupied molecular orbital (LUMO) is of A1 symmetry. Interchange of the SOMO and LUMO will in this case lead to occupation of orbitals of different symmetry and thus to a different electronic state. Adjustment of the initial guess orbitals can be achieved with the guess=alter keyword and definition of the list of orbital exchanges after the geometry definition. In the following input orbitals 5 and 6 of the initial guess will be interchanged.
#ROHF/STO3G scf=(symm,tight) guess=alter
rohf/STO3G opt amide radical NH2, 2A1 state
0 2
N1
H2 1 r2
H3 1 r2 2 a3
r2=1.05726795
a3=100.0637224
5 6

In some cases it might be advantageous to inspect the initial guess before actually doing an SCF calculation with guess=only and also to check the intended alterations with guess=(only,alter). With the SCF calculation for the ^{2}B_{1} state in hand, we can compfortably skip this procedure. Here the SCF procedure gives the following final result:
Harris functional with IExCor= 205 diagonalized for initial guess.
ExpMin= 1.69D01 ExpMax= 9.91D+01 ExpMxC= 9.91D+01 IAcc=1 IRadAn= 1 AccDes= 1.00D06
HarFok: IExCor= 205 AccDes= 1.00D06 IRadAn= 1 IDoV=1
ScaDFX= 1.000000 1.000000 1.000000 1.000000
Pairs of Alpha orbitals switched:
5 6
Initial guess orbital symmetries:
Occupied (A1) (A1) (B2) (A1) (A1)
Virtual (B1) (B2)
The electronic state of the initial guess is 2A1.
Requested convergence on RMS density matrix=1.00D08 within 128 cycles.
Requested convergence on MAX density matrix=1.00D06.
Requested convergence on energy=1.00D06.
No special actions if energy rises.
Keep R1 and R2 integrals in memory in canonical form, NReq= 419815.
SCF Done: E(ROHF) = 54.3257900934 A.U. after 9 cycles
Convg = 0.7698D09 V/T = 1.9800
S**2 = 0.7500
Annihilation of the first spin contaminant:
S**2 before annihilation 0.7500, after 0.7500
**********************************************************************
Population analysis using the SCF density.
**********************************************************************
Orbital symmetries:
Occupied (A1) (A1) (B2) (A1) (A1)
Virtual (B1) (B2)
The electronic state is 2A1.
After performing the default Harris functional guess all orbital changes are listed (here only 5 and 6 are interchanged). This leads to an initial guess wavefunction with a SOMO of A1 symmetry. The relative order of orbital symmetries is retained during the SCF optimization and the final wavefunction therefore describes an ^{2}A_{1} state. This wavefunction is stored in the checkpoint file on completion of the calculation and can be read from there as the initial guess for yet another calculation with guess=read. The following input file reads, for example the geometry as well as the wavefunction obtained in the previous calculation from the checkpoint file (here named /scr1/testit.chk) and calculates the ROHF/STO3G energy within a single SCF cycle. Since the wavefunction stored in the checkpoint file is of ^{2}A_{1} symmetry, the SCF will converge to the ^{2}A_{1} state without any orbital interchanges.
%chk=/scr1/testit.chk
#ROHF/STO3G scf=(symm,tight) guess=read geom=check
rohf/STO3G opt amide radical NH2, 2A1 state
0 2
Comparison of the two states optimized here shows that the ^{2}B_{1} state based on the default Harris functional guess is energetically more favorable at 54.8368134090 au than the ^{2}A_{1} state at 54.3257900934 au. Thus, the default guess arrives at the electronically more favorable state. This is, however, not always the case.
Selection and manipulation of the initial guess for UHF wavefunctions is only slightly more challenging than that for the ROHF case. The added difficulty consists in deciding whether alpha or beta orbitals should be interchanged in order to arrive at the desired electronic state. For the current case of the NH_{2} radical, use of the default Harris functional guess yields again the energetically more favourable ^{2}B_{1} state at 54.8392961557 au. Please note the ROHF/UHF energy lowering of only 0.0025au (6.5 kJ/mol) in this case, suggesting that the wavefunction is well described already in the ROHF picture. The orbital listing obtained for the converged UHF wavefunction is:
SCF Done: E(UHF) = 54.8392961557 A.U. after 10 cycles
Convg = 0.4780D08 V/T = 2.0048
S**2 = 0.7585
Annihilation of the first spin contaminant:
S**2 before annihilation 0.7585, after 0.7500
**********************************************************************
Population analysis using the SCF density.
**********************************************************************
Orbital symmetries:
Alpha Orbitals:
Occupied (A1) (A1) (B2) (B1) (A1)
Virtual (A1) (B2)
Beta Orbitals:
Occupied (A1) (A1) (B2) (A1)
Virtual (B1) (A1) (B2)
The electronic state is 2B1.
Comparison of the occupied alpha and beta spin orbitals shows that the additional alpha electron is located in an orbital of B1 symmetry. The complementary beta spin orbital is not occupied. Excitation of the electron in the highest lying beta orbital (A1 symmetry) to the next higher beta orbital (of B1 symmetry) changes the situation such that in both the alpha and beta orbitals there is one B1 and one B2 occupied orbital while the number of alpha orbitals with A1 symmetry exceeds that of the beta orbitals by one:
In order to effect the interchange of beta orbitals 4 and 5 the following input file will be used. Please observe the two empty lines after the geometry input. The first empty line terminates the geomtry input section, the second empty line specifies no changes in the alpha spin orbitals, and the following line specifies the changes in beta spin orbitals:
#UHF/STO3G scf=(symm,tight) guess=alter
uhf/STO3G opt amide radical NH2, 2A1 state
0 2
N1
H2 1 r2
H3 1 r2 2 a3
r2=1.05808265
a3=100.14734623
one empty line!!
a second empty line!!
4 5
final empty line!!!
The final orbital ordering obtained after the SCF in this case is:
SCF Done: E(UHF) = 54.7318269247 A.U. after 10 cycles
Convg = 0.6550D08 V/T = 2.0029
S**2 = 0.7508
Annihilation of the first spin contaminant:
S**2 before annihilation 0.7508, after 0.7500
**********************************************************************
Population analysis using the SCF density.
**********************************************************************
Orbital symmetries:
Alpha Orbitals:
Occupied (A1) (A1) (B2) (A1) (B1)
Virtual (A1) (B2)
Beta Orbitals:
Occupied (A1) (A1) (B2) (B1)
Virtual (A1) (A1) (B2)
The electronic state is 2A1.
Please observe that the difference between the ROHF/STO3G (at 54.3257900934 au) and the UHF/STO3G (at 54.7318269247 au) energies of the ^{2}A_{1} state is much larger at 1066 kJ/mol than found for the ^{2}B_{1} state at the same theoretical levels.
There are two different SCF algorithms available in Gaussian. The default algorithm DIIS is quite fast and works well for most systems. For problematic cases the quadratically convergent (QC) algorithm is much more reliable but also much slower than DIIS. The latter option is used with the SCF=QC keyword and becomes quite usefull in the case of convergence problems with the default DIIS algorithm.
In the case of severe convergence problems, a few further strategies exist to approach the problem. As a first point it is important to check the geometry of the system under investigation. Errors in the ZMatrix definition, which position atomic centers too close together, often result in convergence problems during SCF calculations. Also, it is often found that self consistency is achieved with one method, but not another. It is therefore quite helpful to calculate a wavefunction first with one method and then use the converged wavefunction as the initial guess for another calculation. It is, for example frequently found that HartreeFock calculations converge more readily than DFT calculations and initial calculation of the HF wavefunction can therefore aid in getting the DFT calculation started. Ultimately, the convergence behaviour also improves with reduced electron numbers. If, for example, calculations for neutral radicals do not converge, a converged wavefunction can often be obtained for the corresponding cationic systems (lacking one electron). The cation wavefunction can then be used as the initial guess of the radical system. If still no convergence can be achieved: go look for another project!