The trouble with transition metals II: SCF convergence

In a previous post I showed how to build two structural models containing transition metals: a small model of a zinc finger and a ferrocene. This post is about computing the energy using GAMESS.

The first challenge in computing the energy is choosing the correct charge and multiplicity. Most organic molecules are singlets (i.e. all orbitals are doubly occupied) and have an easily identifiable charge. This is often not the case with molecules containing transition metals. Many transition metal atoms have unpaired d electrons and can have several different charges (oxidation states). Also, the charge of the ligands is not always obvious, though they almost always are singlets.

The zinc finger model is relatively easy: zinc is almost always Zn²⁺ and has 5 doubly occupied d orbitals, so that's a pretty safe bet. The zinc finger model contains two neutral groups (the imidazole rings) and two negatively charged groups (CH3-S^-). So the charge of the entire system is zero and the multiplicity is 1.

The ferrocene model is a bit more complicated: iron is commonly found in both the +2 and +3 oxidation state. Before one starts the calculation it is important to find out which one is appropriate for ferrocene. Google says +2 (meaning 6 d electrons) and all doubly occupied orbitals (see the orbitals at the bottom of the page), so a singlet. Furthermore, it is also important to know that each ring has a -1 charge, so the total charge of the complex is 0 and the mutiplicity is 1.

Now, I was all set to talk about SCF convergence problems, i.e. that the zinc finger would converge fine but the ferrocene would give problems. But when I ran a RHF/6-31G(d) energy calculation

$contrl runtyp=energy $end
$basis gbasis=n31 ngauss=6 ndfunc=1 $end
$scf dirscf=.t. $end

both runs converged fine!

"Luckily" when using a smaller basis set (RHF/STO-3G) leads to convergence problems - for both. Here is the relevant part of the SCF output for zinc finger.

                                                                                                                                                                     NONZERO     BLOCKS
ITER EX DEM     TOTAL ENERGY        E CHANGE  DENSITY CHANGE     ORB. GRAD      INTEGRALS    SKIPPED
 1  0  0    -3064.6703808712 -3064.6703808712   1.733010527   0.000000000        9751402     457807
 2  1  0    -3059.0602240553     5.6101568159  16.724847696   0.794728909        9737106     471195
 3  2  0    -3053.0258026844     6.0344213709  16.846293078   0.656307634        9778274     459644
 4  3  0    -2935.9005712044   117.1252314800 115.104028436   1.822663152        9792060     446482
 5  4  0    -2923.8944526114    12.0061185930 115.740577390   0.620923864        9792989     441267
 6  5  0    -2768.2361993103   155.6582533011 115.747362137   5.442595435        9783876     443565
 7  6  0    -2921.1390148151  -152.9028155048 115.772970667   0.810505145        9782321     444765
 8  7  0    -2775.0326048138   146.1064100013 115.770983819   4.741230605        9783214     444582
 9  0  0    -2920.9171138162  -145.8845090024  60.570511833   0.796143103        9783475     444487
10  1  0    -3003.9501726068   -83.0330587906  54.636432025   0.749401761        9790866     444493
11  2  0    -2862.1853212282   141.7648513786 115.095938299   1.824897836        9796139     439522
12  3  0    -2921.4368542174   -59.2515329892 115.495602218   0.722894426        9783664     443903
13  4  0    -2848.2964423996    73.1404118178 115.538938253   2.394847457        9782038     444795
14  5  0    -2905.4155965583   -57.1191541586 115.529755578   0.678236838        9781255     445351
15  0  0    -2859.1068450197    46.3087515386  14.660262004   1.964358011        9783437     444643
16  1  0    -3019.7713654862  -160.6645204665  14.032213984   0.801579670        9794434     443563
17  2  0    -2934.9445251088    84.8268403773  99.891887745   1.748565491        9794346     445610
18  3  0    -2905.6599083961    29.2846167127 100.392269592   1.167341691        9789953     445157
19  4  0    -2864.9366576173    40.7232507788 104.001625891   1.797814244        9781805     446591
20  5  0    -2855.1081148594     9.8285427580 103.862237046   1.447532031        9781260     446026
21  6  0    -2865.3401581051   -10.2320432457  82.200569407   1.143844678        9780463     446052
22  0  0    -2857.7261418857     7.6140162194  78.128990865   1.766124007        9781544     445447
23  1  0    -3036.7926259046  -179.0664840190   5.185277571   0.905312779        9796083     442344
24  2  0    -2977.1610659805    59.6315599241   5.419478057   1.193389134        9796601     444172
25  3  0    -2859.5939849458   117.5670810347 114.799971576   1.724872824        9789000     444821
26  4  0    -2913.7281172708   -54.1341323250 115.553208606   1.094842493        9785096     443663
27  0  0    -2845.9860251842    67.7420920865  16.890315601   1.995942532        9783715     444171
28  1  0    -2989.2091196251  -143.2230944409  16.260193167   1.244487455        9796844     439846
29  2  0    -2911.6240721064    77.5850475187 114.870828357   1.711669748        9798193     439590
30  3  0    -2897.8498052725    13.7742668339 115.443658127   0.704588469        9789911     442464

SCF IS UNCONVERGED, TOO MANY ITERATIONS
  TIME TO FORM FOCK OPERATORS=      84.7 SECONDS (       2.8 SEC/ITER)
  FOCK TIME ON FIRST ITERATION=       2.6, LAST ITERATION=       3.0
  TIME TO SOLVE SCF EQUATIONS=       0.4 SECONDS (       0.0 SEC/ITER)

FINAL RHF ENERGY IS        0.0000000000 AFTER  30 ITERATIONS

Clearly, the this SCF is not converging and giving it more iterations is not going to make the problem go away. Instead, different algorithms for SCF convergence are needed and

$scf dirscf=.t. diis=.t. damp=.t. $end

leads to convergence in for both ferrocene and zinc finger.

DIIS stands to Direct Inversion of the Iterative Subspace, which is an extrapolation procedure, while DAMP "damps" any oscillations between in the SCF.

There are a few other tricks to SCF convergence but they are clearly not needed here. If DIIS and DAMP doesn't solve the problem for you, leave a message describing your molecule and I'll see what I can do.