Cundari Th.R. -- Computational Organometallic Chemistry-0824704789

168

Maseras

because of the impossibility of having dangling bonds. The position of these atoms can be very critical if electronic effects of the MM region are going to be introduced in the calculation through the introduction of electrostatic charges or related parameters at their positions. Fortunately, things are not so dramatic when methods like IMOMM are used.

The approach chosen by IMOMM and related methods consists of associating the position of these additional link atoms to those of the atoms they are replacing. The two bonds are forced to lie in the same direction, and this solves the problem of two of the three degrees of freedom associated with the coordinates of each link atom. There remains the problem concerning the bond distance between the atoms. It must have different values for the QM atom–link atom and for the QM atom–MM atom distances.

Different answers to this problem have been proposed by IMOMM and derived methods (34). The simplest consists of freezing this value in both the QM and MM calculations, at corresponding different values, while other, more elaborate approaches define either a constant factor between the two values or a constant difference between them. In the simplest case, which is the one used in the examples presented throughout this chapter, there is the problem of choosing the particular values for both frozen distances. This choice, which will affect the numerical outcome of the calculation, is to a certain extent arbitrary. It would be particularly troublesome if it had a large effect on the outcome of the calculations.

In order to analyze this topic, a set of IMOMM(RHF:MM3) calculations has been carried out on the relative stabilities of the cis and trans isomers of Pt(PtBu3)2(H)2. This is a square planar complex that has been used before for validation tests of the IMOMM method (23). The QM/MM partition applied here uses Pt(PH3)2(H)2 for the QM part, as in the previous tests, and the basis set is also the same (23). In this IMOMM calculation there are six P–H distances to be frozen in the QM part, with the associated six P–C distances to be frozen in the MM part. The values used previously for these parameters were 1.420 and

˚

1.843 A, respectively, taken from the equilibrium values used by the MM3 force field for the corresponding atom types. Additional calculations presented here

˚

consist of making a rather large displacement of 0.05 A in each direction for each of the two parameters. Thus, the P–H distance has been assigned values of

˚

1.370, 1.420, and 1.470 A; and the P–C distance has been set to 1.793, 1.843,

˚

and 1.893 A. The energy differences for the cis/trans pair in each of the resulting nine cases are collected in Table 1. The differences in relative stabilities of the two isomers do not deviate by more than 1.5 kcal/mol from those obtained with the standard values, even taking into account these abnormally large displacements from the equilibrium distances.

The conclusion therefore is that the choice of the particular values for the frozen distances related to the connecting atoms has a very minor effect on the outcome of the calculation, as far as the values applied are reasonable. A standard

Quantifying Steric Effects via QM/MM Methods

169

TABLE 1 Dependence of the IMOMM(RHF:MM3) Computed Energy Difference (kcal/mol) Between cis and trans Isomers of Pt(PtBu3)2(H)2 with Respect to the Frozen Values of Distances in the QM/MM Boundary Region

	˚	˚	˚
	P–H 1.370 A	P–H 1.420 A	P–H 1.470 A
˚	26.4	26.6	26.8
P–C 1.793 A
˚	25.2	25.4	25.6
P–C 1.843 A
˚	24.1	24.3	24.5
P–C 1.893 A

˚

The different values (in A) for the P–C distance in the MM calculation are presented

˚

in the rows, and the different values (in A) for the P–H distance in the QM calculation are presented in the columns.

equilibrium value taken from a force field can in principle be an acceptable option. Care must be taken, however, always to use the same value throughout the whole set of calculations that must be compared.

4. PERFORMANCE CONSIDERATIONS

Apart from the initial setup of the calculation discussed in the previous section, there are a number of additional technical features that, without being specific to the application of QM/MM methods to transition metal chemistry, find particular importance in the performance of this type of calculation. They are briefly reviewed in this section

4.1. Use of Microiterations

In a typical pure QM or pure MM geometry optimization, all geometrical degrees of freedom are treated equally. At a given step of the optimization cycle, the energy is computed at the corresponding geometry, and so is the gradient, the first derivative with respect to each geometrical parameter, and, eventually, the hessian or second derivative. Geometry convergence is checked; and if it is not achieved, a new step is defined by the displacement of each of the variables.

The use of microiterations supposes a substantial breaking of this algorithm. The flux diagram of the alternative optimization scheme is shown in Figure 4. The geometrical degrees of freedom are divided into two sets, which can be labeled as M macrovariables and m microvariables. While the typical algorithm just described would consider the eventual convergence of the M m total variables at each geometry optimization step and generate a new value for each of them, the microiteration scheme works in a different way. At a given step in the optimization cycle the energy, gradient, and, eventually, hessian are computed; and the stationary-point condition (zero gradient) is checked for all variables, in

170

Maseras

FIGURE 4 Flux diagram of an optimization scheme making use of microiterations.

much the same way as in a standard optimizer. The microiterations appear in the process of generation of the new geometry. Instead of modifying all the M m variables at once, only the m microvariables are modified (microiteration) in a frozen environment of the M macrovariables. The micro-optimization continues in a conventional way until the m microvariables are converged. Only after this process is finished are new values for the M macrovariables generated (macroiteration). The convergence of the total gradient is checked again, and the microiterations then have to be repeated for the new value of the macrovariables, the whole process continuing until the macrovariables are converged.

In principle, the optimizing scheme just described is independent of the QM/MM methodology and could also be applied to pure QM or pure MM calcu-

Quantifying Steric Effects via QM/MM Methods

171

lations. However, its main effect would be to increase the total number of iterations and, consequently, the required computer time, thus making little sense. Things are very different in QM/MM methods like IMOMM and related methods. The QM energy is not affected by the position of atoms in the MM region, and therefore there is no need to recompute the QM energy during the micro-optimi- zation process. As a result, a thousand microiterations will usually cost less than one single macroiteration. The use of a small number of macrovariables will reduce the number of macrosteps as compared to those required in a conventional optimization, with a consequent saving in the number of QM calculations and amount of computer effort required.

The efficiency of microiterations is well exemplified by test IMOMM calculations on the cis isomers of Pt(P(tBu)3)2(H)2, a system already used earlier and taken from an IMOMM calibration study (23). In this molecule, there are a total of 83 atoms. Therefore, the total number of geometrical degrees of freedom in the system is 243 (83 3 coordinates per atom, minus the 6 degrees of freedom corresponding to translation and rotation). If one further subtracts the 6 frozen distances between the connecting atoms, there are a total of 237 variables to be optimized in the calculation. If one uses a standard optimization scheme, with no microiterations, after 99 steps and 144 minutes of computer time the calculation is still far from convergence, with the maximum gradient 10 times above the threshold.

When one applies microiterations to this same problem, with the QM/MM partition described earlier, there are 21 macrovariables and 216 microvariables. The optimization is finished after 20 macrosteps, and it consumes only 41 minutes of computer time on the same workstation. The number of microiterations is of course much larger, with some macrosteps taking over 100 microsteps. This is, however, not a problem, because the microiterations, requiring only MM calculations, take much less computer time.

Although microiterations are a powerful tool that can save substantial computational effort in QM/MM calculations, it can also, in certain cases, introduce specific problems to the calculation. The user must take care to prevent them.

The first problem concerns the optimization of transition states. The use of microiterations implies that the microvariables are going to be fully optimized to a minimum energy in each step. Therefore, components of these variables in the transition vector will be neglected. This is not a problem if the atoms involved in the transition vector are in the QM region, and care should be taken so this is accomplished. At any rate, the problem affects only the computation of the second derivative, and the identification of geometries as stationary points (zero gradient) will still be rigorously accurate regardless of the use of microiterations.

A second problem with the use of microiterations appears when this method leads to different local minima for the same distribution of the macrovariables. This type of situation will usually make the calculation fail. The geometry opti-

172

Maseras

mizer assumes that each geometry (defined here exclusively by the macrovariables) has a unique energy and a unique gradient. If this is not the case at some point, the optimization procedure usually collapses in endless loops. There is, however, one easy way out of this problem: the addition of selected geometrical variables from the MM region into the macrovariables set, because there is certainly no methodological requirement to restrict the set of macrovariables to those in the QM and link regions.

4.2.Other Technical Aspects

The QM/MM programs that commonly have been applied to organometallic systems make use of geometry optimizers taken from pure QM programs. These optimizers are reliable in terms of locating the local minimum (or transition state) closer to the starting geometry, but offer no warranty in terms of whether this is the absolute minimum. This problem can be critical if there are several possible conformations available to the chemical system. Although this is true for any computational method, the problem is more likely to appear when one deals with very large interlocked ligands. And these are precisely the kind of systems that will be more commonly calculated with hybrid QM/MM methods.

There are at least two possible solutions to this problem. One is to carry out a conformational search. A number of possible minima are examined through a certain algorithm, and only the lowest of them is taken. This procedure is available in a number of MM programs, and in principle it should be required only for the MM region. Therefore, one could carry this search with the QM region frozen at a reasonable computational cost. The second option is to start from a known crystal structure. If the program starts with a geometry in the correct conformation, it will stay there. In any case, one must be aware that the option of defining a geometry with a graphical program from scratch can lead to an incorrect conformation and can seriously hinder the validity of the computation of relative energies.

A final comment on the practical use of QM/MM methods concerns the need to apply sophisticated features in the optimizer. In general, the QM/MM geometry optimization will be more challenging to the program than that of the corresponding QM optimization of the model system, for a variety of reasons. For example, one should expect a larger coupling between the geometrical variables and a less accurate initial estimation of the hessian. As a result, the user will probably have to exploit more uncommon optimization features than when interested in a pure QM calculation. In the particular case of the implementation of IMOMM in the mmabin program, using the gaussian92/dft code for the geometry optimizer, the features often required by QM/MM calculations included special care in the definition of the Z-matrix, numerical calculations of hessian elements for selected geometrical parameters, and application of additional features

Quantifying Steric Effects via QM/MM Methods

173

like restricting both the smallest and the largest acceptable absolute value of the eigenvectors of the hessian.

5. ANALYSIS OF RESULTS

Hybrid QM/MM methods have been developed essentially as a tool for the calculation of reliable geometries and energies at a reasonable computational cost. However, the results of their application can also be analyzed in a way that can provide further insight into the properties of chemical systems.

5.1. Identification of Steric Effects

Like several other concepts in chemistry, the definition of steric effects is clear from an intuitive point of view but difficult to put into a rigorous mathematical definition. The traditional division between electronic and steric effects loses precision when one realizes that steric effects are also included in the wavefunction and that they are also ultimately caused by the presence of electrons and nuclei. A possible definition of electronic effects as those going through bonds and steric effects as those going through space is appealing. But this has also problems in a number of cases. Should one label as steric electrostatic interactions between charged parts of different ligands in an organometallic complex? Are the agostic intramolecular M H–C interactions electronic or steric?

This problem in the definition of electronic and steric effects has no clear solution, because it is essentially a problem of semantics, and this is not the subject of computational chemistry. However, hybrid QM/MM approaches like IMOMM and derived methods, neglecting the EQM(QM/MM) term (Sec. 2), can be useful in providing an arbitrary definition, as good as many others, that can be quantified in mathematical terms. Our proposal consists of defining steric effects as the perturbations introduced by the MM region in the properties of the QM region of the system. The rationale behind this definition is that these interactions consist exclusively of geometrical strains, and that corresponds precisely to the usual understanding of steric effects. This definition is not absolute, because it depends on the particular force field applied in each case, but can lead to a consistent quantitative separation of qualitatively different effects on chemical systems.

A good example of the possibilities of hybrid QM/MM methods in the separation of electronic and steric effects was provided by a study on the relative stability of the cis and trans isomers of Ru(CO)2(PR3)3 complexes (11). The two isomers considered, presented in Figure 5, have a trigonal bipyramidal geometry, and they are labeled following the arrangement of the two carbonyl ligands. The two forms have been observed experimentally (11,49) by changing the nature of the phosphine ligand PR3. In particular, when PR3 PEt3, the cis isomer is the

174

Maseras

FIGURE 5 The two observed isomers in Ru(CO)2(PR3)3 complexes.

species present in the crystal. When PR3 P(iPr)2Me, two independent molecules are present in the crystal, one of them trans and the other cis. In this second case, IR studies in solution seem to indicate a larger proportion of the trans isomer. IMOMM(MP2:MM3) calculations were carried out on these complexes using Ru(CO)2(PH3)3 for the QM part. Both isomers, cis and trans, were found to be local minima for each complex, with geometries in agreement with available X- ray data.

The most relevant part of this study was, however, the comparison of the relative energies. For Ru(CO)2(PEt3)3, the cis isomer, the only one existing in the crystal, was computed to be more stable than the trans isomer by 3.0 kcal/ mol. The relationship between the two isomers was reversed for complex Ru (CO)2(P(iPr)2Me)3, with the trans species being more stable by 2.8 kcal/mol. Therefore, the QM/MM calculation reproduced the experimental observation accurately. This result would likely also have been obtained through the performance of much more expensive pure QM calculations. But the use of a QM/ MM method had the additional advantage of allowing the clarification of the origin of the difference between both complexes, by a simple decomposition of the total QM/MM energy in its two terms, QM and MM. The QM contribution, representing the electronic effects, always favors the cis isomer, by 3.1 kcal/mol in the case of Ru(CO)2(PEt3)3 and by 1.7 kcal/mol in the case of Ru(CO)2(P (iPr)2Me)3. The MM contribution, representing the steric effects, marks the difference between both complexes. This part always favors the trans isomer, but it does so by quite different magnitudes: 0.1 kcal/mol when the phosphine is PEt3, and 4.5 kcal/mol when the phosphine is P(iPr)2Me. The conclusion from these results is straightforward. There is always a small electronic preference for the placement of the π-acidic carbonyl ligands in equatorial positions (cis isomer). There is a steric preference toward the placement of the bulkier phosphine ligands in the equatorial positions (trans isomer), with the weight of this preference depending on the nature of the phosphine ligand. Only in the case of bulkier phosphines are the steric effects strong enough to overcome the electronic preference for the cis isomer.

Quantifying Steric Effects via QM/MM Methods

175

It can certainly be argued that the conclusions for Ru(CO)2(PR3)2 complexes could have been deduced from the experimental results alone without need of calculation. Apart from the fact that the performance of calculations always strengthens the validity of qualitative arguments, it is worth mentioning that this example proves precisely the validity of the assignment of the MM part of the calculation to the steric effects.

A second example of identification of steric effects by QM/MM calculations concerns a case when their presence was not obvious a priori. This is the case of a joint experimental and theoretical study (19) on [Ir(biph)X(QR3)2)] (biph biphenyl-1,2-diyl; X Cl, I; Q P, As). These five-coordinate complexes have a distorted trigonal bipyramidal geometry, as shown in Figure 6. The phosphine (or arsine) ligands occupy the axial sites, and the chelating ligand biph and the halide are in the equatorial sites. One of the most intriguing features of these compounds is the deviation that the halide presents from the symmetrical arrangement between the two coordinating carbons of biph, a deviation that is characterized by values of φ (Fig. 6) different from zero. Previously reported calculations at the extended Hu¨ckel level seemed to indicate an electronic origin for the deviation. However, pure Becke3LYP calculations on the [Ir(biph) Cl(PH3)2] model system produced a symmetrical structure (φ 0), making unlikely an electronic origin for the distortion. The steric origin of the distortion was proved by IMOMM(Becke3LYP:MM3) calculations on the real system [Ir(biph)Cl(PPh3)2], which yielded a distortion angle φ of 11.4° (as compared with the experimental value of 10.1°).

Again, this result would probably also have been obtained through the performance of very expensive pure QM calculations on the full system. A pure QM calculation on the real system would nevertheless only prove the decisive role of the phenyl substituents of phosphine in the distortion but could not decide on whether their effect was electronic or steric. The fact that this distortion appears in the IMOMM calculation, where the QM effects of phenyls are neglected, in itself constitutes a direct proof that the origin of the distortion is purely steric.

FIGURE 6 Two views of the geometry of [Ir(biph)X(QPh3)2] complexes, with indication of the definition of the distortion angle φ.

176

Maseras

Another analysis showed that the symmetrical arrangement of chlorine corresponds to a sterically hindered position and that its distortion had a very low electronic cost.

In this particular case, the steric origin of the distortion indicated by the hybrid QM/MM calculation was further confirmed by new experiments. In one of them, PPh3 was replaced by AsPh3, and in the other Cl was replaced by I. The logic behind these tests was that the modification of the size of the ligands should affect the magnitude of the steric effects and, therefore, the value of φ, and this was indeed observed in the crystal structures. In this way, this work constituted a nice example where the importance of steric effects was unearthed by hybrid QM/MM calculations and afterwards confirmed by new experiments.

5.2.Quantification of Steric Effects

Not only do QM/MM calculations lead to the identification of steric effects, they can also lead to their quantification and to its assignment to specific regions of the chemical systems.

An example of the characterization of intraligand steric effects in transition metal complexes was provided by the IMOMM(MP2:MM3) study of [Re- H5(PR3)2(SiR3)2] systems (10). Two different complexes of this type, [Re- H5(PPh(iPr)2)2(SiHPh)2] and [ReH5(PCyp3)2(SiH2Ph)2] (Cyp cyclopentyl), were analyzed, with the QM part being constituted by the [ReH5(PH3)2(SiH3)2] model system. The two complexes, which had been experimentally characterized by X-ray and, in one case, neutron diffraction (50), are nine-coordinate and have a capped square antiprism structure, with qualitatively similar structures. However, there are some quantitative differences between the crystal structures, especially concerning the bond angles. The largest difference between the experimental structures is in the Si–Re–Si angle, which was well reproduced by the IMOMM calculations.

More to the point, the same study (10) exemplifies the possibilities of a quantitative analysis of steric effects through hybrid QM/MM calculations. This was carried out through a simple computational scheme consisting of several steps: 1) separation of the geometrical variables of the model system into two sets (A, consisting of the geometrical variables to be analyzed, and B, consisting of the other geometrical variables); 2) full IMOMM geometry optimization on the real system; 3) full QM optimization on the model system; 4) restricted IMOMM optimization on the real system, with the geometrical variables of set A frozen at the values of the MO optimization in the model system; 5) comparison of the results of steps 2 (full IMOMM optimization) and 4 (restricted IMOMM optimization). The basic idea in this scheme is that the restricted IMOMM optimization yields the geometry the system would take in the absence of steric effects. The initial separation of geometrical variables in step 1 is necessary to put aside

Quantifying Steric Effects via QM/MM Methods

177

geometrical variables that are deemed unimportant by chemical common sense but that could heavily influence the analysis. A typical example would be ligand rotation around M–P bonds.

This scheme was applied to both the [ReH5(PPh(iPr)2)2(SiHPh)2] and [Re- H5(PCyp3)2(SiH2Ph)2] complexes (10). The comparison between the total energies of the restricted and full IMOMM optimized geometries showed energy differences of 4.9 kcal/mol for the first complex and 7.3 kcal/mol in the case of the second. This total energy was decomposed into its QM and MM contributions. In the case of [ReH5(PPh(iPr)2)2(SiHPh)2], the energy difference of 4.9 kcal/mol was obtained from an MM stabilization of 5.3 kcal/mol and an MO destabilization of 0.4 kcal/mol. The corresponding numbers for [ReH5(PCyp3)2(SiH2Ph)2] were 10.8 and 3.5 kcal/mol, respectively. Steric effects are therefore clearly more important in the second complex, a result that was not obvious before the calculations.

The calculations can nevertheless go further into the analysis of the steric effects. This is possible because of the mathematical structure of the MM energy, which is just a summation of terms. In the case of complex [ReH5(PCyp3)2 (SiH2Ph)2], presented in Figure 7, the term-by-term decomposition of the 10.8 kcal/mol of MM difference between the restricted and the full optimization shows that there are two clearly dominant contributions, one of 7.9 kcal/mol associated with the ‘‘van der Waals’’ (VdW) nonbonding interaction, and another of 2.9 kcal/mol from the ‘‘bending’’ interaction. This is not surprising, because the traditional view of steric repulsions associates them precisely to nonbonding interactions through space, which are represented in the MM3 force field by the ‘‘van der Waals’’ term. This VdW term can be further decomposed into each of

FIGURE 7 IMOMM(MP2:MM3) optimized structure of complex [ReH5(PCyp3)2- (SiH2Ph)2]. Hydrogen atoms not directly attached to the metal are omitted.