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.

22 comments:

baoilleach said...

This was one of the major stumbling blocks for me during my PhD, and the reason I couldn't use GAMESS. Does GAMESS still not allow projection from a smaller basis set onto a larger? If so, why not? It solves just about every problem in Gaussian.

Jan said...

GAMESS now allows projection of the MINI basis using
$guess guess=rdmini $end
And you're right, it often solves the problem.

baoilleach said...

For pasting in output that doesn't fit, you might want to checkout the pre tag with style overflow:auto. Check out my blog for a few examples here and there.

NUChem said...

Jan,

Again your posts are very informative and I look forward to each new post. I know you're always looking for ideas so I had one that I hope you could help me out on. Could you show how to model a TS of a bond migration with an overall charge of +1? For example a Mannich type reaction in which a bond migrates to an iminium ion which results in a carbocation? Is this possible in GAMESS?

Thanks again,
NUChem

Jan said...

Noel - That worked great. Thanks for the tip!

NUChem - Thanks! Do you mean a reaction like this?

NUChem said...

Jan,

Yes that is basically the idea. I have a system in which a C-C bond is migrating onto the iminium ion which then leaves a carbocation where the C-C bond use to be. I am having trouble getting the TS of the bond migration. Thanks for your help.

NUChem

Jan said...

NUChem - are you sure there is a transition state? i.e. have you optimized the carbocation-imine complex, and all other minima along the reaction path?

See the end of this post.

Gas phase chemistry is often very different from that in solution when ions are involved.

NUChem said...

Jan,

I'm in the process of minimizing everything. My problem is I have a somewhat large molecule (C22H29N) and a macbook pro, so the optimizations using higher theory levels takes a long time. My next task is to try the TS search with a smaller model system.

NUChem

Pedro Alexis said...

Hi Jan,
First of all thanks for your posts which are very informative. What I'm doing is trying to optimize clusters of transition metals like Ni,Pt,Pd. What i want to know if those algorithms for SCF convergence are needed:
$scf dirscf=.t. diis=.t. damp=.t. $end
And if they force I mean the method for 30 iterations in SCF to converge? Actually, what i wrote is:MAXIT=200 and I let the process of minimizing work.

Thanks

Alex,

Jan said...

Hi Pedro, happy to hear you find the posts informative.

With regard to your question: damp and diis could speed up convergence, but there's no way to know without trying.

mamin said...

Dear Jan,
I have a convergence problem with a manganese cluster formed by 4 manganese and one calcium ions connected through oxu-bridges. I tried many things but nothing works.
here is my input file.
$CONTRL SCFTYP=UHF RUNTYP=OPTIMIZE COORD=CART NZVAR=0
MULT=1 ICHARG=1 DFTTYP=B3LYP INTTYP=HONDO $END
$SYSTEM MEMORY=9000000 $END
$BASIS NGAUSS=6 GBASIS=N31 $END
$STATPT OPTTOL=1.0E-4 NSTEP=30 $END
$SCF DIRSCF=.TRUE. FDIFF=.FALSE. diis=.t. damp=.t. $END
$ELPOT IEPOT=1 WHERE=PDC $END
$PDC PTSEL=CONNOLLY CONSTR=CHARGE $END
$PCM SOLVNT=H2O $END
$DATA

C1
Ca 20. 29.09000 39.75700 66.82400
O 8.0 28.26900 38.72900 67.80000
O 8.0 29.86400 37.56000 68.57000
O 8.0 27.71000 36.97400 68.46500
O 8.0 29.17100 36.66200 65.76200
Mn 25. 28.49900 38.12000 69.54300
Mn 25. 29.02500 36.52000 67.43100
Mn 25. 26.31200 38.26900 68.06400
Mn 25. 29.38900 35.90400 64.26600
O 8.0 27.48915 35.79648 63.98412
H 1.0 26.93515 35.24248 63.43010
O 8.0 29.53406 34.30248 63.02473
H 1.0 28.98005 34.85648 62.47073
H 1.0 30.08806 33.74848 62.47073
O 8.0 29.02800 38.99052 64.26143
H 1.0 28.47400 38.43652 63.70743
H 1.0 29.58200 39.54453 63.70743
O 8.0 31.24017 34.76896 63.95595
H 1.0 30.68617 34.21496 63.40196
H 1.0 31.79417 34.21496 64.50993
O 8.0 30.96084 36.31232 66.49402
H 1.0 31.51484 35.75832 67.04802
$END
thanks

Jan said...

OK, that looks like a tough one. It will be a while before I have time to give it the attention it deserves. You may want to post your question to the GAMESS google list in the meantime.

Andre said...

Hi Jan,

I'm trying to study some adsorbed aminoacids, do you think that I can do this on GAMESS? I tried cysteine adsorbed on Pt, but I had some dificulties... How can I freeze the Pt atoms?

Thanks

Jan said...

You can use IFCART in $STATPT to freeze the Cartesian coordinates of atoms.

My main worry would be SCF convergence problems for large numbers of Pt atoms.

André said...

I've tried to use a lot of atoms but it was to much for the pc so I'm using now just a part of unit cell. I've tryed to use just one too but I also didn't get any advantage, GAMESS coudn't converge this input.
The best input I've tried was this with a part of the unit cell, but seems that was to much too... So I'm searching ways to impruve this one... Now i've used the UFF that I didn't know.
Do you think that this freezing could save calculation time and use less of the computer memory?

Andre said...

Hi Jan,

I've tried to use IFCART but GAMESS didn't understend this... Always the output says that are a problem in $STATPT... Could you teach me how to use this? And as you had said, the convergence is a problem... I've made a input w/ just 5 Pt (a face of the crystal) and the problem continious...Do you have any sugestion?

Thank You very much

Jan H. Jensen said...

I suggest you post your questions to the Google GAMESS mailing list: http://groups.google.com/group/gamess?hl=en

Include as much information as possible, such as input files and error messages from the output

Mridu said...

JAN,

I am very new to GAMESS.And I am trying hard to go with this. for last few days, I have performed to get the excitation energy of Heptacene in neutral, protonated and deuterated forms, I am done with these three, but when i tried to do the same with different protonation site, It's showing too many iterations and scf is not converged. My input is

!protonated heptacene
$CONTRL SCFTYP=RHF DFTTYP=BLYP TDDFT=EXCITE RUNTYP=GRADIENT COORD=CART $END
$CONTRL icharg=+1 $END
$CONTRL MOLPLT=.TRUE. PLTORB=.TRUE. NPRINT=9 $END
$SYSTEM MWORDS=250 $END
$BASIS GBASIS=N311 NGAUSS=6 NDFUNC=1 NPFUNC=1 $END
$GUESS GUESS=HUCKEL $END
$TDDFT NSTATE=20 $END
$SCF DIRSCF=.T. DIIS=.T. DAMP=.T. $END
$DATA
protonated heptacene TDDFT
cs

C 6.0 -7.3928601857 1.4098455534 0.0000000000
C 6.0 7.3928601857 1.4098455534 0.0000000000
C 6.0 -7.3928601857 -1.4098455534 0.0000000000
C 6.0 7.3928601857 -1.4098455534 0.0000000000
C 6.0 -8.5642863201 0.7172401590 0.0000000000
C 6.0 8.5642863201 0.7172401590 0.0000000000
C 6.0 -8.5642863201 -0.7172401590 0.0000000000
C 6.0 8.5642863201 -0.7172401590 0.0000000000
C 6.0 -6.1278524901 -0.7287038203 0.0000000000
C 6.0 6.1278524901 -0.7287038203 0.0000000000
C 6.0 -6.1278524901 0.7287038203 0.0000000000
C 6.0 6.1278524901 0.7287038203 0.0000000000
C 6.0 -4.9227845523 -1.4078474530 0.0000000000
C 6.0 4.9227845523 -1.4078474530 0.0000000000
C 6.0 -4.9227845523 1.4078474530 0.0000000000
C 6.0 4.9227845523 1.4078474530 0.0000000000
C 6.0 -3.6778874904 -0.7301125257 0.0000000000
C 6.0 3.6778874904 -0.7301125257 0.0000000000
C 6.0 -3.6778874904 0.7301125257 0.0000000000
C 6.0 3.6778874904 0.7301125257 0.0000000000
C 6.0 -2.4603204113 1.4090219507 0.0000000000
C 6.0 2.4603204113 1.4090219507 0.0000000000
C 6.0 -2.4603204113 -1.4090219507 0.0000000000
C 6.0 2.4603204113 -1.4090219507 0.0000000000
C 6.0 -1.2261544434 0.7310531305 0.0000000000
C 6.0 1.2261544434 0.7310531305 0.0000000000
C 6.0 -1.2261544434 -0.7310531305 0.0000000000
C 6.0 1.2261544434 -0.7310531305 0.0000000000
C 6.0 0.0000001454 1.4093991286 0.0000000000
C 6.0 -0.0000001454 -1.4093991286 0.0000000000
H 1.0 9.5111136776 1.2457649106 0.0000000000
H 1.0 -9.5111136776 1.2457649106 0.0000000000
H 1.0 9.5111136776 -1.2457649106 0.0000000000
H 1.0 -9.5111136776 -1.2457649106 0.0000000000
H 1.0 7.3930208595 -2.4951590275 0.0000000000
H 1.0 -7.3930208595 -2.4951590275 0.0000000000
H 1.0 7.3930208595 2.4951590275 0.0000000000
H 1.0 -7.3930208595 2.4951590275 0.0000000000
H 1.0 4.9236957234 -2.4938764573 0.0000000000
H 1.0 -4.9236957234 -2.4938764573 0.0000000000
H 1.0 4.9236957234 2.3466224202 0.8729421226
H 1.0 -4.9236957234 2.4938764573 0.0000000000
H 1.0 2.4607792579 -2.4949211058 0.0000000000
H 1.0 -2.4607792579 -2.4949211058 0.0000000000
H 1.0 2.4607792579 2.4949211058 0.0000000000
H 1.0 -2.4607792579 2.4949211058 0.0000000000
H 1.0 0.0000000776 -2.4952684669 0.0000000000
H 1.0 -0.0000000776 2.4952684669 0.0000000000
H 1.0 4.9236957234 2.3466224202 -0.8729421226
$END

Jan H. Jensen said...

your CH bond lengths to C16 are very long. Try shortening them to something more reasonable like 1.1 Å

Nagat said...

Hi Jan
I'm a new in Gamess. I've SCF convergence problem in dealing with heavy atoms by Gamess. For examples, for Pr free atom
The log.file read
ENERGY COMPONENTS
-----------------

WAVEFUNCTION NORMALIZATION = 1.0000000000

ONE ELECTRON ENERGY = -34.3239010426
TWO ELECTRON ENERGY = 34.3239010426
NUCLEAR REPULSION ENERGY = 0.0000000000
------------------
TOTAL ENERGY = 0.0000000000
My input file is
$CONTRL SCFTYP=ROHF RUNTYP=ENERGY MULT=6 QMTTOL=1.0E-07 MAXIT=199 ECP=READ $END
$SCF DIRSCF=.TRUE. DIIS=.TRUE. DAMP=.TRUE. SHIFT=.TRUE. $END
$GUESS GUESS=HCORE $END
$DATA

C1
Pm 61.0 -0.00000 0.00000 0.00000
S 1
1 2.3480000 1.0000000
S 1
1 0.3778000 1.0000000
S 1
1 0.3058000 1.0000000
S 1
1 0.0554400 1.0000000
S 1
1 0.0276800 1.0000000
S 1
1 0.0143500 1.0000000
P 1
1 2.4190000 1.0000000
P 1
1 0.3014000 1.0000000
P 1
1 0.1093000 1.0000000
P 1
1 0.0466300 1.0000000
P 1
1 0.0171300 1.0000000
P 1
1 0.0098150 1.0000000
D 1
1 4.2670000 1.0000000
D 1
1 0.6316000 1.0000000
D 1
1 0.2996000 1.0000000
D 1
1 0.1536000 1.0000000
D 1
1 0.0724700 1.0000000
D 1
1 0.0317500 1.0000000
F 1
1 70.8100000 1.0000000
F 1
1 23.1300000 1.0000000
F 1
1 9.1960000 1.0000000
F 1
1 3.6410000 1.0000000
F 1
1 1.4020000 1.0000000
F 1
1 0.4748000 1.0000000

$END
$ECP
PM-ECP GEN 54 3
6 ----- F POTENTIAL -----
-0.27581099 2 0.18629999
-3.67671704 2 0.72170001
-19.20204740 2 2.52950001
-80.51749420 2 10.04810050
-295.95886200 2 51.84069830
-32.22702030 1 154.48420700
8 ----- S-F POTENTIAL -----
0.02899800 2 0.12690000
59.00218200 2 0.63349998
-79.18837740 2 0.73229998
43.98767090 2 1.05410004
36.28269960 2 3.92600012
30.94334980 2 4.81040001
29.88560490 1 13.42710020
5.22673082 0 6.61800003
8 ----- P-F POTENTIAL -----
21.08691790 2 0.47749999
-63.05217740 2 0.74510002
87.14070130 2 1.02049994
-51.46301650 2 1.54830003
54.61180120 2 2.59549999
2.29255104 2 5.04430008
37.44199370 1 7.40789986
5.61524105 0 21.90730100
8 ----- D-F POTENTIAL -----
-0.40402901 2 0.43770000
30.99323080 2 1.01680005
-96.03415680 2 1.22340000
168.40460200 2 1.65579999
-190.01771500 2 2.43820000
150.01397700 2 3.39700008
36.84980780 1 9.34249973
7.18157816 0 21.38209920
$END

I don't know what is the error.
Please I need your help and many thanks in advance.
Best regards
Nagat

Anonymous said...

i am getting error = ENERGY=0 ...
i.e the iteration are too much and SCF is not converging .......i was trying for system of benzene and guanidinium...
my input file is
$CONTRL SCFTYP=RHF RUNTYP=ENERGY ICHARG=1 DFTTYP=B3LYP COORD=ZMT $END
$BASIS GBASIS=N31 NGAUSS=6 $END
$SCF DIRSCF=.TRUE. FDIFF=.FALSE. diis=.t. damp=.t. $END
$DATA
benzene-dimer-2.5A...1-A-1 state...RHF/STO-3G
C1
H
H 1 B1
H 2 B2 1 A1
H 3 B3 2 A2 1 D1 0
X 3 B4 2 A3 1 D2 0
N 5 B5 3 A4 2 D3 0
H 6 B6 5 A5 3 D4 0
H 6 B7 5 A6 3 D5 0
N 5 B8 3 A7 2 D6 0
N 5 B9 3 A8 2 D7 0
H 10 B10 5 A9 3 D8 0
H 10 B11 5 A10 3 D9 0
H 9 B12 5 A11 3 D10 0
H 9 B13 5 A12 3 D11 0
C 6 B14 5 A13 3 D12 0
C 15 B15 6 A14 5 D13 0
C 10 B16 5 A15 3 D14 0
C 9 B17 5 A16 3 D15 0
C 5 B18 3 A17 2 D16 0
C 5 B19 3 A18 2 D17 0
C 18 B20 9 A19 5 D18 0
H 17 B21 10 A20 5 D19 0
H 16 B22 15 A21 6 D20 0

B1 4.66999775
B2 4.66999775
B3 4.66999775
B4 4.66917310
B5 3.54260941
B6 1.90000000
B7 1.90000000
B8 3.54261091
B9 3.51679573
B10 1.90000000
B11 1.90000000
B12 1.89999999
B13 1.90000001
B14 2.53622578
B15 2.62999999
B16 2.46933576
B17 2.52588655
B18 2.62920747
B19 2.50000000
B20 2.63000000
B21 2.04000000
B22 2.04000000
A1 119.99999997
A2 119.99999981
A3 60.00539255
A4 69.75997982
A5 111.26957970
A6 111.71783511
A7 69.66833171
A8 136.14375513
A9 111.01703839
A10 111.01711949
A11 111.18819161
A12 111.79950189
A13 47.80528532
A14 88.28655094
A15 48.34549142
A16 47.83512626
A17 0.19365678
A18 90.61172695
A19 88.30508692
A20 92.44439034
A21 120.11347300
D1 -0.00000000
D2 0.28007305
D3 130.88589456
D4 0.76016322
D5 134.57579571
D6 49.17267056
D7 89.94039382
D8 -113.43090886
D9 113.62298405
D10 -0.72252791
D11 -134.53804670
D12 67.46098371
D13 60.05722290
D14 -179.90583492
D15 -67.40196804
D16 -103.48820490
D17 90.00900374
D18 -60.06516142
D19 -179.87639503
D20 92.46229480
$END

any help will be great...
thanx

Jan H. Jensen said...

I suggest you post your questions to the Google GAMESS mailing list: http://groups.google.com/group/gamess?hl=en

Include as much information as possible, such as input files and error messages from the output