J. Phys. Chem. 1993,97, 3122-3121
3722
Geometry Optimization of Organometallic Compounds Using a Modified Extended-Hiickel Formalism Franqois Savary and Jacques Weber' Department of Physical Chemistry. University of Geneva, 30 quai Ernest-Ansermet, CH-I 21 1 Geneva 4, Switzerland
Cion Calzaferri Institute for Inorganic and Physical Chemistry. University of Bern, Freiestrasse 3, CH-3000 Bern 9, Switzerland Received: December 2, 1992
The modified extended-Huckel method suggested recently by Calzaferri and co-workers has been parametrized for organometallic compounds in order to have at hand a semiempirical tool allowing for rapid optimization of their geometries. The model is based on (i) the introduction of an approximate two-body repulsiveelectrostatic energy term and (ii) the use of a distance-dependent Wolfsberg-Helmholz formula, K = 1 Ked(RAH-do),for 1-2 and 1-3 interactions. The optimum values for the K , 6 parameters determined for given classes of metalligand interactions have been found to be transposable to systems exhibiting several types of ligands such as aromatic rings, carbonyls, phosphines, etc. The overall performance of the model is satisfactory, as it leads for a series of more than 30 organometallic compounds to average error bars of 0.058 and 0.039 A on metal-ring and metal-carbonyl bond distances, respectively, and of 2.6' on bond angles.
+
I. Introduction Molecular modeling has recently evolved from the stage of an esoteric tool reserved to computational chemists into a mature field enabling practically every experimentalist to perform his own applications. Indeed, several user-friendly packages can be used today to study, at different levels of sophistication, molecular structure, function, and interaction.' However, whereas such investigations can be performed on a routine basis for organic compounds, the situation is different for transition metal complexes, such as organometallics,due to the well-known difficulties associated with an adequate ab initio quantum chemical description of the structure and properties of these compound~.2-~ Density functional methods represent undoubtedly in this field a very interesting alternative to the traditional ab initioapproach, specially from the point of view of the quality of the results,@ but a full geometry optimization of organometallics using these formalisms cannot yet be qualified as a routine task. Empirical force field methods, which are powerful tools for building the structure of organics? cannot be easily generalized for d-electron compounds, although some recent results in this area are promising.lsl4 There is therefore still a need for simple computational models, based, for example, on semiempirical quantum chemical methodologies, allowing for rapid and efficient optimization of the geometry of organometallic compounds. Interesting solutions in this direction have recently been suggested: Zerner et al. have parametrized the so-called ZINDO model for organometallic species leading to satisfactory results for the geometries of various complexes;l5 Blyholder et al. have derived MIND013 parameters for iron, which allowed them to rationalize the structure of several iron complexes;'6 finally, Li and Jug have developed an extended version of the SINDO1 method to optimize the geometry of transition metal compounds with an accuracy comparable to that of 3-21G ab initio SCF calculations. On the other hand, the extended-Hiickel method (EHM0)18 in its improved ASED (atom superposition and electron delocalization) form has been shown to describe reasonably well the structural propertiesoforganometalli~compounds.~~-2~ Recently, Calzaferri et a1.22and Tupper et a1.23 havediscussedand improved
this approach, which led to better structural predictionsfor small molecules. These authors, however, do not report any parameter development for organometallics. We therefore present results obtained for the geometry optimization of such compounds using an adequatelyparametrized EHMO-ASED procedure. This model has been chosen so as to be coherent with our reaction potential method calculating a reactivity index for electrophilic or nucleophilic attack on organometallic compounds and based on extended-Hiickelwave functions.24-2' 11. ASED Formalism The EHMO model is the simplest quantum chemical method consideringall the valence electrons of a system. It is well known that the EHMO model in its original form is, in principle, not able to optimize geometries correctly as it lacks repulsive electrostatic interactions. Andersonand HoffmannI9have shown that by introducing a two-body electrostatic correction term, this deficiency can be overcome. The total energy is calculated asz2
Et,, = UEHMO + E~cp
(1)
where AEEHMO is the extended-Hackel binding energy calculated as
E E H M is O the 'classical" extended-Hiickel total energy; bi and I$ are occupation number and valence-state ionization energy
(VSIE) of atomic orbital p, respectively, the summation over p running on all the atomic orbitals of the complex. The repulsive energy ERep is expressed as
(3) with EABoriginally calculated by Anderson and Hoffmanni9as (4)
where EABis the repulsive electrostatic energy for the interaction
0022-3654/93/2091-3122504.~0~0 @ 1993 American Chemical Society
The Journal of Physical Chemistry, Vol. 97, No. 15, 1993 3723
Geometry Optimization of Organometallic Compounds of atoms A and B, separated by distance RAB,with ZA,2, and p A being the nuclear charge, position, and electron density of the
atom A, respectively, B being more electronegative than A. Calzaferri et al. have suggested reformulating E A B as follows:22
where the second right-hand term of eq 5 is the arithmetic mean of A-B and B-A attractive interaction energies. On the other hand, Anderson has suggested improving the model by modifying the off-diagonal EHMO matrix elements in the following way:20 1
H,," = ~ K ( H , ,+ , H ~ ~ ) s ~ ~ ~ - ~ ~ A B(6) where K = 2.25, 6 is an adjustable parameter, whose optimal value is 0.13 A-I, and RABis the distance between the atoms bearing the orbitals p and Y , respectively. The deficiencies of this approach have been discussed and eliminated by including a distance-dependent exponential factor into the K constant as follows:22
K~~ = 1 + Ke-d(RAB-do)
(7)
where KABis the EHMO K parameter used in the nondiagonal with orbitals p and Y belonging Hamiltonian matrix elementsH,,, to atoms A and B, respectively; K and 6 are positive empirical parameters such as 0.75 IK 5 1.25 and 0.13 A-1 5 6 50.35 A-1, with do being equal to the sum of atomic orbital radii rn(A) r,(B) which are defined by
TABLE I: Values of A-B Interactions'
I[,
6 Parameters for Classes of Atomic
atom A
atom B
metal metal metal metal metal metal carbon (Carbonyl) carbon carbon
1-2 Interactions carbon (cyclobutadienyl) carbon (cyclopentadienyl) carbon (benzene) carbon (naphthalene) carbon (carbonyl) phosphorus oxygen (carbonyl) hydrogen phosphorus
carbon (benzene) carbon (carbonyl)
1-3 Interactions carbon (carbonyl) carbon (carbonyl)
K
6A-'
0.40 0.50 0.40 0.50 0.10
0.35 0.40 0.50 0.60 1.90
0.20 1.20
2.00 0.00 0.90 0.00
0.40 0.60
0.00 1.60
1.20
0.60
1-2 and 1-3 interactions refer to bond-stretching and angle-bending interactions, respectively.
parameters in order to obtain the best agreement with experimental geometries. To this end, we have modified the EHMO program30 by implementing (i) a quasi-Newton-Raphson gradient-based geometry optimizer3' and (ii) a scanning procedure over a range of K , 6 values. Then for test compounds whose experimental structural properties are known, the relative root mean square (RRMS) of the error on calculated geometricparameters Pi(calc) relative to the observed ones Pi(exp) has been minimized:
+
For single-t Slater-type orbitals, one has
whereas for double-f Slater-type orbitals, with exponents f l , J; and coefficients C I ,c2
n
rn =
2'"(f
c:rI
+ c:f2 +
f )n+1/2
(9b)
1 2
(ZI + f2P
where a0 is the Bohr radius. This reformulation of the EHMO method leads to a significant improvement of the potential energy curvesof diatomic and small organic molecules, especially after optimization of the K , 6 parameters for given classes of compounds.22 An investigation on the geometry of organic molecules containing C, H, 0,N, and S has led to the result that bond lengths and bond angles can well be described with 0.85 < K < 1.215 and 6 = 0.35 A-I. As an example the error on the calculated bond lengths is on the order of 1% only if ~(C,sp)= 0.85, K(C,SP~) = 0.975, K(C,SP~) = 1.215, and 6 = 0.35 A-I is applied.28 Our purpose here is to extend this parametrization to organometallic compounds. As a first step, we have decided to simplify the calculation of do by approximating it as the sum of the covalent radii R c o vand ~ R c o v(taken ~ from ref 29), which transforms eq 7 as
KAB 1 + , ~ - ~ ~ R A B - ( R , ~ A + R C ~ B ) ]
(10) In a second step, geometry optimizations have been carried out for given classes of organometallicsby modifying the constants K and 6 of eq 10, and it has been found that the calculated structural properties were relatively sensible on the basis of the values of these parameters. It seemed thus possible to optimize these
Our strategy in the search of optimum K and 6 values is the following: in order to be efficient and of a general use, the EHMO model builder is allowed tocontaindifferent setsof K , 6 parameters for the description of different atom-atom interactions and in particular metal-ligand interactions. In other words, a single set of K , 6 values should be able to lead to good predictions for a given type of metal-ligand interactions (i.e., carbonyls, aromatic rings, phosphinegroups,etc.). In addition, wedecided to limit ourselves to the first-row transition metal atoms and to use for all of them the same set of K , 6 parameters when coordinated to the same ligand. This procedure limits the number of K , 6 parameters to a minimum while being flexibleenoughtodifferentiate thevarious types of metal-ligand interactions encountered in organometallic chemistry; it is well known for instance that carbonyl ligands are better r acceptors than cyclopentadienyl rings and this should somehow be reflected by different choices of K , 6 parameters. However, in view of the deliberate strategy of using the same parameters for all the first-rowtransition metal atoms, the model could be more approximate when examining the variations of metal-ligand distances within series of compoundswith the same ligands such as carbonyls or metallocenes. It is clear, however, that the K , 6 parameters of KAB(q 10) correspondingto a long-range interaction which is not optimized should not be taken as variables. Together with the case of K values of both diagonal H,,and nondiagonal H,"elements on the same center, they are chosen as standard values: K = 0.75, 6 = 0.30. The values of the optimized K , 6 parameters are reported in Table I. In all the cases, we assumed that the internal geometry of the ligands (Le., the C-O bond distance in metal carbonyls, the C-C and C-H bond lengths in metallocenes, etc.) may be kept frozen to standard values. Indeed, we have noticed that these structural parameters exhibit small variationswithin a given
Savary et al.
3724 The Journal of Physical Chemistry, Vol. 97, No. 15, 1993
Q
_-
Figure 1. Three-dimensional mapof RRMSvalues calculated asa function of K and 6 in the case of the optimizationof Fe-C and C-O bond distances in Fe(C0)S assuming a Djh structureand inequivalentaxial and equatorial geometrical parameters. Experimental values are reported in Table 11. The valley on the left correspondsto small RRMS values, Le., to optimum K. 6 values.
series of organometallics and their optimization is much less important than that of the parameters involving the metal atom. The choice of K, 6 parameters for intraligand atom-atom interactions is therefore of little influence on the results and their values are reported in Table I. However, the values of these parameters may differ for interligand atom-atom nondiagonal EHMO matrix elements (see Table I) as they may reflect 1-3 interactions that are responsible for angle bending at the metal center.
F i e 2. Three-dimensionalmapof RRMS values calculated as a function of K and 6 in the case of the optimization of the ring-metal-carbonyl dihedral angle (Scheme 11) in (rl6-C6H,j)Cr(C0),. The exptrimental value is reported in Table IV. The valley at the center corresponds to the optimum K, 6 values. Kra 2.6 0.41 Z'
-6
6-0.30 i;'
.....-..6-0.13
22
A"
III. Results and Discussion Figure 1 presents a three-dimensional(3D)map corresponding to RRMS values calculated as a function of K, 6 parameters for Fe(C0)s while keeping the D3h symmetry and optimizing separately the axial and equatorial distances. It is seen that a very narrow valley is obtained between the region corresponding roughly to = KO.1,2.0> to = X0.5, O.O>. A detailed examination of the data correspondingto Figure 1 shows that the absolute minimum is found at = .We assume that this map can be considered as the one obtained if we had optimized the Fe-C distances (axial and equatorial) only as the C-O distances do not vary much. It is interesting to see that this valley extendsalong a region where K, 6 are approximately inversely proportional and this point will be treated below. A similar relation has been obtained for several complexes in which bond distances have been optimized. A similar surface of RRMS values calculated as a function of K, 6 parameters is presented in Figure 2 in the case of (C6H6)Cr(CO)3 when optimizing the dihedral angle (Scheme I). A valley is again observed but, contrasting with Figure 1, it is very wide without exhibiting a particular absolute minimum. The relation between K, 6 along the valley is different from that found in Figure 1: increasing K corresponds to increasing 6 values along the bottom of the valley. The different behavior exhibited by the RRMS values in Figure 1 and 2 can be easily rationalized when looking at Figure 3, which displays KABas a function of RABfor the Hz molecule by using various 6 values. It is seen in Figure 3 that while K A B is exponentiallydecreasing as a function of RABfor any value of 6 # 0, all the curves cross at the point . Indeed, when RAB= do, the exponential of eq 10 is equal to 1 and K = 1 + K, Le., 1.75. However the dependence of KABas a function of 6 is different for RAB< do with respect to case for R A B > do: in the first case increasing 6 from 0 to 1 amounts to an increase in K A B and therefore covalency; for RAB> do, Le., for A-B atoms exhibiting long-range
0.0
OSu*
1.0
15
2.0
25
3.0
35
4.0
Figure 3. KABcurves (eq 10) calculated as a function of RABfor various values of the 6 arameter in the case of the H2 molecule with K = 0.75 and do = 0.64
1.
SCHEME I
interactions, the opposite trend is observed and increasing 6 amounts to a decrease in KABand therefore covalency. In the case of bonded atoms with RAB do, Le., for nonbonded A-B atoms, any attempt to increase the covalent interaction between A-B atoms leads to a simultaneous increase of 6 as in this case the model tries again to reduce an unreasonable increase of covalency. The fact that the interpretation of the behavior of optimum RRMS values as a function of K, 6 is rather easy is due to the simplicity of the EHMO approach and to the clear dependence of the results upon the parameters used. The narrow valley displayed in Figure 1 suggests that for bonded atoms a relation may be found between K, 6 that would lead to a single parameter in eq 10. However this relation has been found to be typically dependent on the complex investigated and
Geometry Optimization of Organometallic Compounds
The Journal of Physical Chemistry, Vol. 97, No. 15, 1993 3725
TABLE 11: Geometry Optimization of Metal Carbonyl Complexes with K = 0.10 and 6 = 1.90 A-l distance (A).
V(CO)bh
Cr(C0)6'
Fe(CO)sd
Ni(C0)4c
M-C,,(calc) M-C,,(exp) M-C,,(calc) M-CJexp) C-0,, (calc) C-O,,(exp) C-0, (ca IC) C-Odexp)
1.909 2.003
1.856 1.910
1.305 1.838
1.160 1.139
1.146 1.141
1.834 1.827 I .806 1.806 1.150 1.153 1.140 1.153
RRMS
0.053
0.024
0.01 1
0.017
SCHEME I1
/ n \ .
/ I
\ \
.
'$$' C L
1.154 1.140
i
metal-ring distance
most of these compounds, the K , 6 parameters used for 1-2 interactions are those reported in Tables I1 and 111, a few cases appliesonly to Fe(C0)s. Reference 32. Reference 33. Reference 34. excepted (see footnotes of Table IV). In addition, the parameters Reference 35. for 1-3 interactions, Le., C-ring-M-CO (Scheme I) and ringwe have decided therefore to keep the two K , 6 parameters as M-CO (Scheme 11) angles, have been optimized in the case of independent variables, which leads in any case to more accurate (C6H6)Cr(C0)3. It is seen in Table I v that, as an average, the predictions. error bar on the M-ring distance (0.058 A) is much larger than Table I1 summarizes the results of the geometry optimization that on the M-CO distance (0.039 A), which may be ascribed of metal carbonyl complexes with a unique set of K , 6 values for to (i) the significant number of compounds with substituted rings, the metal-carbonyl interaction. It is seen that the performance for which the same K , 6 parameters as those of complexes with of the EHMO method in its improved form is reasonably good. unsubstituted ligands have been used, and (ii) the variety of For Fe(C0)S the model reproduces well the shortening of the n-membered ligand rings encountered in the series of systems axial metal-carbonyl distance with respect to the equatorial one. investigated. The M-L distances exhibit the largest deviation On the other hand we tried without success to optimize the K , 6 from experimental values, but it is of the same order as the M-ring values for a carbon-carbon interaction in order to reproduce the distances. This can be explained by the fact that, in some cases, 120' C-Fe-C equatorial angle of the bipyramidal complex. the same K , 6 parameters have been employed for different ligands. Apparently, this is due to the small value used for the K parameter However, taken as a whole, the results presented in Table IV are of the metal-carbonyl interaction. Indeed, when increasing this rather satisfactory taking into account the approximations parameter, we have noticed that proper K , 6 values for C-C inherent to the semiempirical model used. In addition, no electron interaction can be found so as to lead to a 120° C-Fe-C equatorial correlation nor multiplet structure description is explicitly angle. In this case, however, the Fe-C distance is significantly introduced in the present calculations. A comparison is provided underestimated. Finally, it is seen in Table I1 that the metalby the recent paper of Sosa et al.59which presents results obtained carbon distance is poorly predicted for V(CO)6. This could be for the geometry optimization of organometallic complexes, related to the fact that theground stateof thecompound is known generally carbonyls, using both local density functional (LDF) to exhibit a Jahn-Teller distortion. and S C F (3-21G basis set) methods. The average error bar The results obtained for the geometry optimization of metexhibited by these results is 0.038 8, for LDF and 0.068 A for allocenes and dibenzene-metal complexes are reported in Table SCF models, which is comparable to the errors reported in Table 111. As mentioned in section 11, the geometry optimizations have IV. been performed in this case using a single set of K , 6 for all On the other hand, examination of Table IV reveals that both metallocenes and another one for dibenzene complexes. It is bond and dihedral angles are particularly well reproduced by the seen in Table 111 that again the predictions of the EHMO model EHMO model: the average error is 2.6' for the former ones compare reasonably well with experiment, the cases of ferrocene (with a large discrepancy of about 9O observed for (CsH5)Coandvanadocene excepted. The qualitative trends in metal-ligand (C0)2 and (C5H5)V(CO),) and the experimental dihedral angles distance lengthening or shortening when going from vanadium are perfectly reproduced, although the energy differencebetween to nickel are similar to those observed experimentally. The case eclipsed and staggered conformations is generally of the order of ofmanganoceneisinterestingasthiscompound is known toexhibit 1 kcal/mol. Again the comparison with EHMO results obtained a high-spin configuration in the ground state whereas its for organics shows that the present predictions are of alike decamethyl analog is characterized by a low-spin c~nfiguration.~~ a c c ~ r a c y . ~ ~Similarly, * ~ ~ J 8 Sosa et al. report average error bars By using different populations for the predominantly metal 3d of 4.9O and 6.9' for bond angles calculated using the LDF and orbitals of manganocene so as to simulate a low-spin versus highS C F models, respectively, which allows us to consider the spin ground-state configuration, the modified EHMO model is performance of the EHMO predictions as satisfactory. able to reproduce adequately the lengthening of metal-ring Finally, we present in Table V the results of a comparison distance accompanying the low-spin to high-spin transition. The between EHMO and both LDF and SCF results obtained for lengthening is underestimated, however, probably because it is selected organometallic compounds. It is seen that again the not possible to carry out spin-unrestricted calculations in the EHMO model compares more favorably with experiment than EHMO framework. the SCFone, while providing results of roughly the same accuracy The results of the geometry optimization of organometallic as the LDF scheme. It has been shown recently, however, that complexes with mixed ligands are presented in Table IV. For the introduction of nonlocal corrections to the exchange and Eq corresponds to equatorial: axcorresponds toaxial. Thisdistinction
TABLE III: Ceometrv Optimization of Metallocenes and Dibenzene-Metal ComDlexes' (CrHs)2Vh
(C5HhCrh
distance (A) M-ring(calc) M-ring(exp)
1.841 1.928
1.795 1.798
RRMS
0.045 1
0.0017
(GHs)2Mnb
(CsHshFeh
(CsHs)d%'
(CsHs)2Nib
(C6Hb)2V1
(C6H6)2Crd
1.767/1.906' 1 .767/2.046r
1.745 1.661
1.756 1.739
1.78 1 1.828
1.682 1.660
I .630
0.0000/0.0684
0.0506
0.0098
0.0257
0.0133
0.0105
1.613
The values of K , d parameters are 0.5, 0.4 and 0.4,0.5 for metallocenes and dibenzene-metal complexes, respectively. Reference 36. Reference 37.
Reference 38. Metal-ring distance for low-spin/high-spin complexes.
3726 The Journal of Physical Chemistry, Vol. 97, No. IS, 1993
Savary et al.
TABLE IV: Geometry Optimization of Various Organometallic Complexes with Mixed Ligands distance (A) angle (degrees) compounds
ref
M-ring calc
M-ring exv
M-CO calc
M-CO exv
M-L calc
1.889 1.646 1.689 1.683 1.674 1.676 1.653 1.656 1.674 1.685 1.684 1.739 1.665 1.812 1.807 1.799 1.738 1.617 1.767
1.930 1.732 1.728 1.725 1.724 1.740 1.738 1.730 1.725 1.719 1.717 1.709 1.706 1.745 1.756 1.731 1.798 1.700 1.813
1.880 1.820 1.794 1.814 1.814 1.815 1.818 1.836 1.819 1.814 1.814 1.830 1.817 1.834 1.834 1.835 1.935 1.835 1.778 1.8441 1.83Sd 1.8401 1.848d 1.871 1 1.878d 1.837 1.799
1.900 1.784 1.813 1.824 1.825 1.827 1.830 1.842 1.840 1.842 1.845 1.838 1.861 1.821 1.824 1.843 1.800 1.835 1.679 1.8361 2.01(rJ 1.80od 1.7701 1.947e 2.00od 1.7631 1.87Se 1.81Id 1.770 1.819
1.874
1.817
M-L exv
ring-M- ring-MCOcalc C O e x v
dihedral (degrees) C-ring-MCOcalc
C-ring-MCOexD
~
39 40 40 40 40 40 40 41 40 40 40 42 40 43 44 45 46 47 48 49 50
51 52 1.864 53 1.661 54 1.7181 1.7938 54 1.704h 55 1.758' 56
1.740 1.725 1.681/ 1.6608 1.68Ih 1.764'
57 57
58 1.826
1.731
0.058
2.3101 2.319 2.479' 2.28Ik,'
1.867 1.838 1.8641 1.8931 1.868" 1.85od 1.872 1.817 2.323' 0.039
0.060
109.6
118.0
124.4"
125.W
123.5
125.0
126.1 126.1 126.3 126.0 124.4 128.7
125.0 125.0 125.8 123.8 124.5 119.1
123.4 124.0
120.0 121.3
127.1
126.7
30.0"
30.W
30.0 30.0 30.0
30.0 30.0 30.0
30.0
30.0
1.990'J 1.98W 1.968'
2.2821 2.33W 2.346k 2.366kJ 2.252'
2.6
0.0
1-3 interactions optimized for this compound ( K . 6 for carbonykarbonyl and benzene-carbonyl interactions). C14HIO corresponds to an anthracene ligand. ' C14H corresponds to a phenanthrene ligand. Corresponds to equatorial/axial metal-carbonyl distances, respectively. Compound chosen to optimize K , 6 parameters for metal-ethylene 1-2 interaction. /Corresponds to m e t a l - ( ( 2 4 ) distances. 8 Corresponds to metal-cyclobutadienyl/ metakyclopentadienyl distances, respectively. Compound chosen tooptimize K , 6 parameters for metakyclobutadienyl l-2 interaction, the experimental geometry has been taken equal to the metal-cyclobutadienyl distance of the ( C ~ H ~ ) C O ( C ~ complex, HS) see ref 53. Compound chosen to optimize K , 6 parameters for metal-naphthalene 1-2 interaction. Corresponds to the cisltrans metal-phosphorusdistances. Corresponds to the metal-phosphorus distances. Compound chosen to optimize K , 6 parameters for metal-phosphorus 1-2 interaction.
TABLE V: Comparison between Results Obtained with the Improved EHMO Model for Various Complexes and Those of More Elaborate Quantum Chemical Models compound
Cr(CO),, Fe(CO)
