Monday, November 9, 2009

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 Zn2+ 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.
Post a Comment