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

the ‘‘erroneous’’ H–C–H ideal angle would be noticeable as a decrease in the C–C–H angle, but an automated procedure is just as likely to ‘‘correct’’ by increasing the C–C–H ideal angle as well.

What then is the solution to this problem? Nothing inherently says that the

H–C–H ideal angle should be close to tetrahedral. It might very well be possible to set all ideal angles in a force field to 180° and to reproduce the entire data set

by fitting the force constants. However, it has been seen that the predictive ability is enhanced if ideal values are set close to perceived ‘‘unstrained’’ states, and it is definitely more pleasing to the chemist. Therefore, the penalty function should be modified to favor the unstrained state, possibly by tethering.

3.5.Parameter Tethering

It is possible to bias the parameterization to specific parameter values by tethering. In essence, tethering is a way of telling the refinement ‘‘I know what this parameter should be, don’t deviate too much from it.’’ A ‘‘preferred’’ value for the parameter is set before refinement. The squared deviation of the parameter from the preferred value is then added to the penalty function. Another way to look at it (and the simplest way to implement tethering) is simply to see the ‘‘preferred’’ values of the tethered parameters as reference data and the actual values as the corresponding calculated data. The ‘‘weight factor’’ is set, according to the rules outlined previously, to the inverse of the acceptable deviation from the ‘‘preferred’’ values and then increased somewhat to compensate for the fact that only a few ‘‘data points’’ are included.

A weak tethering is generally beneficial for most parameters. If the parameter is well determined by the reference data, the effect of tethering will be negligible, as it should be. On the other hand, if the parameter is very badly determined, even a weak tethering potential will suffice to keep it close to a pleasing value. Ideal bond angles in particular should be tethered to perceived ‘‘unstrained’’ values (see Sec. 3.4). Torsional parameters and cross-terms may also be tethered, usually to a value of zero. Inclusion of QM charges as reference data may in a way be seen as a tethering of the electrostatic parameters, because many other data errors might otherwise have been ‘‘corrected’’ automatically by introduction of physically unrealistic charges. Tethering is especially valuable for achieving a balance in the initial stages of the refinement.

3.6.Strategies for Initial Refinement

Automated adjustments leading to ‘‘unnatural’’ parameter values are most frequently observed in the initial stages of the refinement, when errors are still large. Several techniques are available for minimizing unwanted deviations, including tethering, subset refinement, and analysis of outliers.

Tethering is a logical correction procedure for the ‘‘unnatural’’ parameter deviations observed in the initial refinement stages. The only problem is then to find ‘‘preferred’’ values for all parameters. Hopefully, any preferred values have been set in the initial force field, so the entire parameter set can be tethered to these values, with weight factors set from the confidence one has in each parameter value. As the refinement progresses, tethering weights can be lowered, to avoid biasing the final force field. However, for reasons already discussed, bond angle tethering could be retained throughout.

Initial refinement can be made more efficient by dividing parameters and data into subsets. For example, electrostatic parameters could be adjusted as a group to QM electrostatic data only (and possibly excluded from further refinement altogether). If ideal bond lengths and angles have been set to reliable values initially, force constants could be refined in a group using a penalty function based solely on QM Hessian data. Torsional parameters could be grouped with the force constants, but then the penalty function should be extended to include also conformational energies. When the force-related parameters have been balanced, ideal bond lengths and angles could be refined using structural data only. Unless a very large body of PES data is available, cross-terms should be given low values and left out of the refinement until a late stage.

Error detection should be attempted after a few refinement cycles. Erroneous reference data can frequently be detected by a complete failure of a partially refined force field to fit the data point. On the other hand, such failures are even more frequently a result of deficiencies in the functional form. In any case, the error must be corrected before refinement is finalized, either by removal of the offending reference data or by changes in the functional form of the force field. One can easily identify the data points that are most likely to affect the refinement adversely as the largest weighted contributions to the penalty function. All large deviations must be manually scrutinized to elucidate whether they result from data errors or force field deficiencies. If the latter, it should also be estimated whether a failure to reproduce the data point will adversely affect the intended use of the force field (this is not always the case). If so, the force field setup must be modified before continuing; otherwise the data point might simply be removed.

4. VALIDATION

New parameter sets should always be validated before use. The simple fact that no further improvement can be found is not a sufficient condition for accepting the force field. To verify that the new force field is accurate and predictive, several tests should be performed. First, it should be verified that all the reference data are indeed reproduced with an acceptable accuracy. Second, the precision of each parameter determination should be checked. If a large range of parameter values

can give the same results, the parameter is not well determined by the data set. Finally, the predictivity of the force field should be tested against an external data set.

4.1.Internal Validation

The first step in validation is simply to verify that the remaining errors in the reproduction of the reference data are acceptably small. If the weight factors have been set, as suggested earlier, to the inverse of the acceptable error for each data type, the test is particularly simple. If the final penalty function is lower than the number of data points, the root mean square (rms) error will automatically fall within the acceptable range. The data should also be divided by type and retested, to make sure that the proper balance has been obtained. As before, outliers should be carefully scrutinized. Any errors in the reference data or deficiencies in the functional form are most easily detected at this stage. Plots of calculated vs. reference data can also give valuable information on trends in remaining deviations and possible systematic errors (20).

The penalty function derivatives calculated in Figure 9 give information about how well a parameter is determined by the available data. First of all, it should be verified that the second derivative of the penalty function with respect to each parameter is positive. If not, a lowering of the penalty function can always be obtained by a slight change in the parameter—indeed, one of the test points must have been better for the numerical derivative to be negative. The expected response of the penalty function to a small parameter change can now be calculated from a truncated Taylor expansion (Fig. 10). The parameter change that would be needed to effect a given penalty function change is available from solving the second-order equation.

At this point, it is necessary to decide a maximum ‘‘allowed’’ change in the penalty function, ∆χ2. This choice is necessarily arbitrary, but should reflect a change in the data that corresponds to either the expected input error or the largest deviation that could be accepted in the intended use of the force field (28). Using this value with the equation in Figure 10, it is now possible to calculate the

FIGURE 10 Relationship between small changes in a parameter and the penalty function.

maximum change in the parameter that could be accommodated without significant deterioration of the fit. This range could be reported together with the final parameter value as an indication of the quality in parameter determination. Notice that the equations in Figure 10 do not take account of possible errors in the data (except possibly in the choice of ∆χ2). Thus, the calculated range is a measure not of the accuracy of the parameter but of the precision in the determination (11). In particular, if the reference data contain systematic errors, the ‘‘real’’ value of the parameter may well fall outside the calculated range.

4.2. External Validation

The final test for each force field is how well it reproduces data that have not been included in the reference set. The test set should reflect the intended future use of the force field. In most cases, the goal is to predict experimental properties. Thus, QM data are not needed in the test set, even if they are used extensively in parameter refinement. On the other hand, it is now possible to include data points that cannot practically be included in the parameterization. Examples are experimental IR spectra, where peak assignment is a problem with preliminary force fields, or equilibria, which must be calculated with time-consuming dynamic methods.

In a general treatise, it is not possible to state what discrepancies can be accepted for the final force field. However, if the deviations are substantially worse than the experimentally achievable accuracy, the usefulness of the force field will be very limited. As always, large deviations should be identified, and if possible the underlying causes should be rationalized.

It will sometimes happen that as a result of the final test a redefinition of the functional form with subsequent reparameterization becomes necessary. If the test data that pinpointed the failure are of a type that can be used in refinement, they should be transferred to the reference set. In any case, the test set should be discarded and a new test set selected. The reason for this is that the force field is no longer independent of the initial test set, since it has been used to influence a design decision. One can therefore argue that a reapplication of the same set tests for only internal, not external, predictivity.

5. EXAMPLES

The procedures described herein have been implemented as a package of small C programs and Unix scripts (15).* This implementation has principally been

*Updated versions of the programs and scripts are available from the author on request. Several force fields and example structures can be found at http://compchem.dfh.dk/PeO/.

FIGURE 11 The three types of structures that have been parameterized in the examples.

developed in conjunction with the MacroModel package (35), but all procedures that are specific to one force field have been collected in a few scripts that are easily modified to accommodate other formats. Only minor modifications have been necessary for the procedures to work with other packages. Following are a few recently published examples, reflecting various development stages of the parameterization methods and several types of complexes. Typical structures of each type are depicted in Figure 11.

5.1.Rigid Octahedral Geometry: Ru(II)(bipy)3

Octahedral geometry is among the easiest coordination modes to handle with force field methods, as evidenced by the large numbers of force fields that have been published (8,12). The functional forms of standard organic force fields must be extended to allow differentiation between cis and trans bond angles about the metal center, but the structures are generally rigid enough that standard bond and angle functions can be used. A point worth noticing is that the ideal angle parameter for the trans angle should be exactly 180°, to avoid creating a cusp. Observed deviations from this value should, as far as possible, be reproduced by a lowering of the force constant. The ideal cis angle should be close to 90°, but because there is no discontinuity in the function derivatives here, small variations are allowed. Alternative functional forms that could reproduce the octahedral geometry include POS (since the octahedron minimizes steric repulsion for six ligands) and trigonometric functions with minima at both 90° and 180° (23).

The current implementation used the standard bond and angle functions available in MacroModel MM3*, with cis and trans parameters assigned by a

test of the actual angle in the input geometry (in rare cases with strongly distorted structures, the test may erroneously assign one ligand as being trans to either two or no other ligands). A further refinement included a fourfold torsional term to describe rotation about the metal–ligand bond. The reference data consisted of several X-ray structures, together with QM normal modes and CHelpG charges (48) determined for one small model system. It could be shown that all parameters were well determined by the included data (28).

Several observed complexes with terpyridine display a trans induction. The bond from Ru to the central N is very short, resulting in an elongation of the trans Ru–N bond. This effect could possibly be described as a stretch–stretch cross-term in force field modeling, but no such function is available in MacroModel. It was found that the observed distortions could be reproduced by a direct interaction between two trans ligands. This solution should not be given any physical significance, but should be considered only a working model for reproducing observed geometries (28).

5.2. Flexible Points-on-a-Sphere Model: ( 2–alkene)Pd(0)L2

Alkene coordination to Pd(0) is rather loose, with the preferred in-plane geometry easily distorted by modest steric interactions. As a consequence, the barrier to rotation about the Pd–alkene axis is also low. The C–Pd–L angle is thus highly variable and cannot be well represented by a standard angle-bending function. The two most frequently used force field models for this type of coordination are the dummy atom approach and the POS model. The latter was used in the current model. The QM-calculated barrier to rotation could be reproduced by addition of a nonphysical dihedral angle parameter including all ligands but not Pd.

Very few relevant structures could be found: the reference set of four X-ray structures was therefore augmented with QM structures in addition to the QM normal modes and CHelpG charges. The final force field, although fully converged, could not fully reproduce the reference structures. A few bond angles deviated by several degrees, due mainly to the inadequacies of the POS model. The model could probably have been improved by treating the L–Pd–L angle by a standard function, but this type of differentiation was not easily achieved. On the positive side, the external predictivity was similar to the internal. The final model was judged adequate for the intended application, determination of relative energies of intermediates in the palladium-assisted alkylation reaction. Comparison to other available methods showed that the accuracy in structure determination was similar to the semiempirical method PM3(tm), albeit with different systematic errors, and substantially better than two force fields based on general metal parameters (11).

5.3.Transition-State Model: Os-Catalyzed Asymmetric Dihydroxylation (AD)

The Sharpless AD reaction is an almost ideal test case for a Q2MM study. The selectivity is determined in one well-defined step, which has been well characterized by a combination of high-level QM methods and isotope studies (49). Experimentally, the reaction is not overly sensitive to reaction conditions, tolerating a wide range of solvents, from toluene to water/alcohol mixtures (50). Selectivity data are available in the literature for a wide range of ligand–substrate combinations. Transition-state structures were obtained at the B3LYP level for 59 small model structures. Hessians and CHelpG data for three of the structures and relative energies of several distorted structures were also included in the parameterization, but no experimental data were employed. The final force field was tested on a range of substrates by extensive conformational searches for all low-energy reaction paths. The selectivities could then be calculated from the Boltzmann populations of diastereomeric structures and compared to experimental enantioselectivities. The results were very good, with most deviations below 2 kJ/mol (20).

6. SUMMARY

Force fields for organometallic complexes can be derived using a fast and consistent parameterization method. Structures, conformational energies, QM-derived vibrational modes, and charges should be used as basic reference data. Force fields derived from such data are useful in the prediction of structures and relative energies. For other types of property predictions, specific reference data may be added. The final parameter set should be tested for precision and internal and external predictivity. The method is complementary to the increasingly popular QM/MM methodology (see Chap. 6), in that extensive conformational searches can be performed rapidly, but the response to drastic electronic changes cannot be reliably predicted.

ABBREVIATIONS

REFERENCES

1.F Jensen. Introduction to Computational Chemistry. Chichester, UK: Wiley, 1999.

2.A Hinchliffe. Chemical Modeling from Atoms to Liquid. Chichester, UK: Wiley, 1999.

3.J Goodman. Chemical Applications of Molecular Modelling. London: Royal Society of Chemistry, 1998.

4.U Burkert, NL Allinger. Molecular Mechanics, ACS Monograph 177. Washington, DC: ACS, 1982.

5.(a) K Gundertofte, T Liljefors, PO Norrby, I Pettersson. J Comput Chem 17:429– 449, 1996. (b) I Pettersson, T Liljefors. In: KB Lipkowitz, DB Boyd, eds. Reviews in Computational Chemistry. Vol. 9. New York: VCH, 1996, pp 167–189.

6.(a) JM Goodman, WC Still. J Comput Chem 12:1110–1117, 1991. (b) I Kolossvary, WC Guida. J Am Chem Soc 118:5011–5019, 1996.

7.(a) TR Cundari. J Chem Soc, Dalton Trans 2771–2776, 1998. (b) P Comba, M Zimmer. J Chem Educ 73:108–110, 1996.

8.CR Landis, DM Root, T Cleveland. In: KB Lipkowitz, DB Boyd, eds. Reviews in Computational Chemistry. Vol. 6. New York: VCH, 1995, pp 73–148.

9.P Comba. Coord Chem Rev 182:343–371, 1999.

10.(a) AK Rappe´, CJ Casewit. Molecular Mechanics Across Chemistry. Sausalito, CA: University Science Books, 1997. (b) AK Rappe´, CJ Casewit, KS Colwell, WA Goddard III, WM Skiff. J Am Chem Soc 114:10024–10035, 1992. (c) JJ Gajewski, KE

Gilbert, TW Kreek. J Comput Chem 19 1167–1178, 1998.

˚

11. H Hagelin, B Akermark, M Svensson, PO Norrby. Organometallics 18:4574–4583, 1999.

12. (a) BP Hay. Coord Chem Rev 126:177–236, 1993. (b) P Comba. Comments Inorg Chem 16:133–151, 1994. (c) M Zimmer. Chem Rev 95:2629–2649, 1995.

13. JP Bowen, NL Allinger. In: KB Lipkowitz, DB Boyd, eds. Reviews in Computational Chemistry. Vol. 2. New York: VCH, 1991, pp 81–97.

14. U Dinur, AT Hagler. In: KB Lipkowitz, DB Boyd, eds. Reviews in Computational Chemistry. Vol. 2. New York: VCH, 1991, pp 99–164.

15. PO Norrby, T Liljefors. J Comput Chem 19:1146–1166, 1998.