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     BLOCKSITER 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     44448710  1  0    -3003.9501726068   -83.0330587906  54.636432025   0.749401761        9790866     44449311  2  0    -2862.1853212282   141.7648513786 115.095938299   1.824897836        9796139     43952212  3  0    -2921.4368542174   -59.2515329892 115.495602218   0.722894426        9783664     44390313  4  0    -2848.2964423996    73.1404118178 115.538938253   2.394847457        9782038     44479514  5  0    -2905.4155965583   -57.1191541586 115.529755578   0.678236838        9781255     44535115  0  0    -2859.1068450197    46.3087515386  14.660262004   1.964358011        9783437     44464316  1  0    -3019.7713654862  -160.6645204665  14.032213984   0.801579670        9794434     44356317  2  0    -2934.9445251088    84.8268403773  99.891887745   1.748565491        9794346     44561018  3  0    -2905.6599083961    29.2846167127 100.392269592   1.167341691        9789953     44515719  4  0    -2864.9366576173    40.7232507788 104.001625891   1.797814244        9781805     44659120  5  0    -2855.1081148594     9.8285427580 103.862237046   1.447532031        9781260     44602621  6  0    -2865.3401581051   -10.2320432457  82.200569407   1.143844678        9780463     44605222  0  0    -2857.7261418857     7.6140162194  78.128990865   1.766124007        9781544     44544723  1  0    -3036.7926259046  -179.0664840190   5.185277571   0.905312779        9796083     44234424  2  0    -2977.1610659805    59.6315599241   5.419478057   1.193389134        9796601     44417225  3  0    -2859.5939849458   117.5670810347 114.799971576   1.724872824        9789000     44482126  4  0    -2913.7281172708   -54.1341323250 115.553208606   1.094842493        9785096     44366327  0  0    -2845.9860251842    67.7420920865  16.890315601   1.995942532        9783715     44417128  1  0    -2989.2091196251  -143.2230944409  16.260193167   1.244487455        9796844     43984629  2  0    -2911.6240721064    77.5850475187 114.870828357   1.711669748        9798193     43959030  3  0    -2897.8498052725    13.7742668339 115.443658127   0.704588469        9789911     442464SCF 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.

Noel O'Boyle 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 Jensen said...

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

Noel O'Boyle 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 Jensen 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 Jensen 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 A. Baldera-Aguayo 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 Jensen 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.

Unknown 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 Jensen 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 Jensen 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 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 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
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 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