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

8

Norrby

by a lack of parameters. Metal systems are structurally more diverse than organic compounds (7). As an example, the C–O–C bond angle seldom deviates more than a few degrees from true tetrahedral, whereas observed P–Pd–P angles vary over a range of ca. 100° depending on coordination geometry and steric requirements. Thus, parameter transferability between different types of complexes is limited, and alternative functional forms may be required (vide infra) (8,9). Despite the apparent difficulties, several force fields exist that allow calculations to be performed for almost any type of complex (10). However, predictivity may well be low for complexes outside the set used in parameter generation (11).

An alternative approach, which will be pursued here, is to tailor a force field to one specific type of complex. For organometallic complexes, it is still possible to use existing parameters for the organic part of the system and to develop new parameters only for the coordination sphere. Many examples can be found in the literature (8,9,12), but the need to develop new parameters largely limits applications to force field experts, as opposed to the organic field, where practicing chemists can easily model the system with only basic computational experience. The goal of the current chapter is to simplify the process of producing a high-quality organometallic force field by providing a workable recipe for the procedure. Some examples from the literature are included, but the coverage is by no means complete.

1.1.Force Fields

A force field is essentially a relationship between the geometry of a molecule and the force on each atom. The force is a vector quantity, the derivative of the energy with respect to coordinates. To simplify the expressions, force fields are generally presented in the form of energy as a function of coordinates. The true zero of the energy is an unknown, different for each force field and molecule. Thus, the total energy calculated for any molecule cannot be interpreted in a physically meaningful way, and no special meaning should be attached to a calculated energy of zero (or a negative energy). However, when two energies are calculated from exactly the same functions (i.e., when the connectivities of two structures are identical), the unknown constants be considered identical, and the energies can be compared directly.*

One of the fundamental postulates of molecular mechanics is that the steric energy of a molecule can be separated into terms resulting from small, transferable moieties. For all bond lengths and angles, it is assumed that there exists an unstrained state with a steric energy of zero. All deviations from this ‘‘ideal’’

*Formal heats of formation can be calculated from steric energies by adding geometry-independent terms for several structural features; see Ref. 4. By this method, structural isomers with different connectivity can also be compared.

Recipe for an Organometallic Force Field

9

FIGURE 1 A simple force field.

value will give rise to an energy increase.* It is generally impossible for all interactions to achieve their unstrained state in the same geometry, and thus the ‘‘ideal values’’ will never be directly observed, but in organic molecules, the deviations from the unstrained state are usually small. Other contributions to the total energy of the molecule come from rotations around bonds as well as nonbonded interactions. In order to reproduce strained structures or vibrational data, it has also been found necessary to employ cross-terms in the force field. An example is the stretch–bend interaction, which can be described as the change in a bond angle function when the constituent bonds are distorted. For trigonal atoms, it is also common to employ a term that differentiates between planar and pyramidal form (an out-of-plane or inversion term).

The functional form of a simple example force field is shown in Figure 1. Most current force fields are substantially more complicated, but the additions take many different forms and will not be covered here. For more detailed accounts, see, for example, Refs. 1 and 4.

The basic unit of a force field is the atom type. In general, there is at least one atom type for each element, more if several chemical environments are to be considered. For example, all force fields differentiate between sp2 and sp3 hybridized carbons, assigning a distinct atom type to each. For organometallic modeling, it is frequently necessary to add new metal atom types to existing force fields. Even when the metal atom types exist in the force field, there is seldom any differentiation based on, for example, oxidation state.† Atom types are used to classify other interactions. Any unique pair of connected atom types identifies a bond type; an angle type is labeled by a unique set of three connected atoms, etc. Each unique interaction type needs its own set of parameters. Many atom-

*The l0 and θ0 parameters are also called reference values. But to avoid confusion with bond and angle ‘‘reference data,’’ the term ideal values will be used for these parameters throughout this chapter.

†One exception is the PCModel program, which allows at least a basic differentiation, see Ref. 10c.

10

Norrby

type combinations will not have existing parameters; in particular the torsions would require determination of millions of parameters for a complete set.

1.2.Parameters

A complete force field consists of a functional form, as exemplified in Figure 1, and a set of parameters. For example, for each type of bond in the example force field, two parameters are needed: an ideal length l0 (corresponding to the bond length in a hypothetical unstrained molecule), and a stretching force constant ks. The latter can be seen as the relative stiffness of the bond, and determines how much the energy increases upon a certain distortion. Some parameters, such as the ideal bond length, correspond closely to observables. However, the optimum set of parameters can rarely be identified by observation.

Take the torsional parameters for the central bond in butane as an example (Fig. 2). There are several observable energies that are closely related to the νn- parameter (Fig. 1), but each is also affected strongly by other parameters. The rotation barrier might be taken as the amplitude of a threefold cosine function (when n 3, ν3 ∆E‡rot/2 will give an energy difference of ∆E‡rot between the lowest and highest point on the torsional profile). But in reality part of the barrier is due to van der Waals (vdW) repulsion, so ν3 should be less than half the observed barrier. Likewise, the conformational difference between gauche and anti forms is largely determined by vdW interactions, but the remaining error might be reproduced using an added ν1 parameter in the force field (the ν3 parameter has no influence on the relative energy of gauche and anti forms, because the contribution from ν3 cos 3ω must always be equal at 60° and 180°).

From vibrational or microwave spectroscopy it is possible to obtain the curvature at the bottom of each well. Ignoring mixing with other structural elements, this corresponds to the second derivative of the energy with respect to the torsional angle, ∂2E/∂ω2. Fitting to this observable may require either sacri-

FIGURE 2 Butane torsional profile.

Recipe for an Organometallic Force Field

11

ficing some accuracy for other data points or adding more torsional terms (i.e., n 1, 2, and 3, Fig. 1). When a substantial amount of data is used, no term in the final force field corresponds to only one type of observable. Instead, the optimal value for each parameter is that which, together with all other terms in the force field, gives the best overall fit to all observables. This concept will be defined more rigorously in Section 3.1. However, it should be clear that changing any parameter might lead to a shift in the optimum value for several others.

Despite what was said in the preceding paragraph, most parameters depend closely on some specific type of data. Good starting values for further refinement can therefore be obtained by manual fitting of one parameter at a time to small subsets of the reference data. The most intuitive example is the ideal bond lengths and bond angles (l0 and θ0, Fig. 1). Averages of observed values (possibly after removal of outliers) are good initial estimates for these parameters. Other examples are given in subsequent sections.

1.3. Parameterization

Defining new force fields has long been as much an art as a science. In the literature, there are two major schools on how to derive force field parameters, manually (4,13) and automatically (14). The manual method has the advantage of creating a deep familiarity with the force field and data, but it requires great expertise. Moreover, when the parameter set grows large, it becomes slow and tedious to ensure that fitting to new data retains consistency with all previously optimized sets.

An automated parameterization may be difficult to set up. But when this has been accomplished, the process is substantially faster than the manual method, and much larger bodies of data can be fitted simultaneously. The main drawback of the automated scheme is that errors may remain undetected more easily than in manual parameterization. Automated parameterization therefore requires substantial validation to identify outliers in the data set and deficiencies in the force field. Statistical tools should be used to verify that each parameter is well defined by the chosen set of reference data, and any ill-fitting data points can be rationalized on sound physical grounds.

The necessary steps in executing an automated parameterization for an organometallic complex are outlined in Figure 3 (15). Each step will be detailed in later sections. Selecting the basic force field is possibly the most critical step. The functional form of the basic force field must be flexible enough to accommodate the variability in metal complexes (7). In addition, the existing parameters for organic moieties will usually not be modified and will therefore limit the accuracy that can be obtained for organic ligands.

The target for an automated parameterization sequence is to enable the force field to reproduce a set of reference data, such as structures and relative

12

Norrby

FIGURE 3 Parameterization flowchart.

energies. The quality of the reference data will therefore limit the attainable accuracy in the final force field. The accuracy can be improved by using large data sets, because random errors are expected to cancel to some extent. However, any systematic errors will be propagated into the final force field. It is also necessary to weight the reference data points, according to both quality and relevance to the intended use of the force field.

With the basic ingredients in hand, the next step is to set up a working force field. It is not necessary at this point to achieve a good fit, but the force field should allow calculations for all structures needed to reproduce the reference data set. This involves choosing functional forms for bonds and angles involving the metal and then guessing reasonable values for all previously undefined parameters. When all of this is accomplished, automatic procedures can vary the param-

Recipe for an Organometallic Force Field

13

eters and calculate all data points iteratively to obtain the best possible fit with the reference data. In the initial stages, it might be necessary to tether parameters and/or to divide them into subsets.

As a final step, the force field should be validated. In part, this is done by evaluating how well the reference data are reproduced and comparing that accuracy to the accuracy needed in the intended application of the force field. However, it is also advisable to apply the force field to an external test set, that is, data points that have not been used at any stage in the parameter refinement.

1.4. Force Fields for Catalysis

To predict reaction selectivities, a special type of force field is needed. Relative reactivities are determined in transition states, whereas most force fields are geared for calculating properties at energy minima. Only rarely have molecular mechanics methods been used for bond-breaking phenomena (16). However, an alternative method that has been successfully applied to selectivity predictions is to treat the transition state as an energy minimum and to develop a force field to reproduce the transition-state (TS) structure (17). This approach allows application of standard molecular mechanics tools such as conformational searching. In a recently developed method, the part of the potential energy surface (PES) perpendicular to the reaction path calculated by QM methods can be closely reproduced by force fields (18). The new method, dubbed ‘‘QM-guided molecular mechanics (Q2MM),’’ has been applied with good results to selectivity predictions in asymmetric synthesis (19) and catalysis (20).

1.5. Selecting a Force Field

There are many points to consider when selecting the program package and force field to be used as a basis for introduction of new organometallic moieties. The available modeling tools, the flexibility of the functional form, and the accuracy of the existing force field are all important. The intended use of the force field will dictate what tools must be present. Sometimes, all that is needed is the generation of good gas-phase structures. If so, the available tools need hardly be considered, since all existing MM packages allow energy minimization. In some situations, the graphical interface may be more important than performance, especially if the force field is to be used for visualization.

A very common use of force fields is to determine relative energies of isomeric forms, since most physical properties will depend on the relative energies of plausible isomers. In this case, it is very important that the underlying force field is already able to produce accurate energies (5). Prediction of thermodynamic properties, solvation, intermolecular interactions, etc. also requires that the basic force field already does well in calculating the particular property for organic molecules.

14

Norrby

Finally, the selected force field must include functionality for describing the coordination environment. Bonds and angles around a metal atom do not behave like organic structures and can only rarely be described by the same functions. Several models for describing coordination angles have been implemented

(8). In some types of complexes, the metal will exert only a weakly directing force. Such coordination can be implemented by replacing all metal-centered angles with nonbonded interactions while still retaining the metal–ligand bonds, as in the points-on-a-sphere (POS) model (21). Alternatively, even the metal– ligand bonds can be described by tailored nonbonded potentials (22), allowing also a variable coordination number.

For more rigid geometries, metal-centered angles are used. Depending on coordination geometry, it may be necessary to differentiate between, for example, cis and trans bond angles, with separate parameters for each. An alternative is to employ functional forms with multiple minima, such as trigonometric functions (23) or more complicated forms (24). More intricate problems are posed by π-ligands, in particular if rotation barriers around the metal–ligand axis are to be reproduced. This type of problem has frequently been addressed by bonding the ligand to the metal through a pseudoatom (8,25).

Coordination complexes frequently display trans-induction and Jahn– Teller distortions, which can been handled by specialized functional forms (9,26,27). In simple cases modified ligand–ligand interactions may suffice (28).

2. REFERENCE DATA

Molecular mechanics is essentially an interpolation method. Reliable predictions for a class of compounds usually require that the force field has been fitted to data of a similar type. When the functional form is physically sound and has been carefully parameterized, limited extrapolation can be successful, but generally only to new combinations of known structural moieties. Thus, an accurate and varied set of reference data is necessary for determination of a good force field. The exact selection depends on the intended use of the force field. Production of rough structures is easily accomplished, but selection of the most favored conformer requires accurate energetics. For prediction of vibrational frequencies, or strongly distorted structures, the shape of the local potential energy surface (PES) around minima must be well described. Solvation and docking requires a good set of nonbonded parameters. For each application, appropriate data must be included in the reference set.

2.1.Structures

The basis for all force field calculations is the generation of sound structures. Without consistent structures, no other properties can be reliably predicted. Thus,

Recipe for an Organometallic Force Field

15

FIGURE 4 Thermally induced oscillation resulting in offset nuclear positions.

the reference set must contain structural information. Depending on the intended application, it might also be necessary to consider the effect on bond lengths from differences in structural determination paradigms (4,29). Most computational methods will determine minima on the PES, that is, nuclear positions at zero Kelvin, without consideration of anharmonic vibrations.* All experimental determinations will take place at higher temperature and thus will include contributions from thermal vibrations, in effect lengthening most bonds slightly.

By far the most common source of structural data, particularly in organometallic chemistry, is X-ray crystallography. It must be noted here that atomic positions determined by X-ray are points of maximum electron density, not nuclear positions. For most atoms, it is a good approximation to consider that the electrons are centered on the nucleus. However, this is never true for hydrogens: the electron of a hydrogen atom always participates in bonding, and is thus offset from the nucleus. This is the major reason why X-ray structures should never be used for determining bond lengths to hydrogens. Other sources, such as neutron diffraction and QM structures, must be used for hydrogen positions.

A large majority of computational structures are determined in vacuo, corresponding most closely to experimental gas-phase structures. A fundamental difference between crystallographic and gas-phase structures is illustrated in Figure 4. Gas-phase methods generally determine bond lengths, whereas crystallographic methods find average atomic positions. Assuming that bonds are stiff and vary little in length, oscillations of rigid moieties in crystals can yield average positions that are closer together than any instantaneous bond length. This behavior is rather common, for example, in flexible chelate rings or in freely rotating phenyl groups, where the apparent Cipso–Cortho bond is shortened. In extreme cases, the reported structures may even be averages of several cocrystalizing conforma-

˚

tions, resulting in bond length errors exceeding 0.1 A (30).

If possible, crystal structure reference data should be compared to calculations in a crystal environment (31). However, many packages do not include

*Exceptions include MM2 and MM3; see Ref. 4.

16

Norrby

the necessary tools for solid-state calculations, necessitating the use of isolated structures in parameterization. When comparing in vacuo calculations with crystal structures, it should also be realized that the crystal structure need not be an energy minimum for the isolated molecule. Crystal packing can have a strong influence on torsions in particular, but also on any long interatomic distances. It is generally safe to compare lists of bond lengths and angles, for these interactions are strong compared to crystal packing. However, structural overlays or complete lists of interatomic distances should not be used as measures of force field quality.

As a validation tool, one may also measure the crystal distortion energy using the force field (11,28,30). This energy must be low, and certainly much lower than the total contribution from the packing forces. If the calculated energy of an error-free, nonionic crystal structure is high relative to the global minimum, the force field may be deficient.

2.2.Energies

Having obtained a good structure, the most important property to be calculated is the energy of the molecule. Except for completely rigid molecules, the structures are distributed among conformational forms where the population depends directly on the energy. Thus, to calculate any property, one must first know the relative energies of all conformers. It follows then that energies must be included in the reference set for any force fields that are not designed solely to yield crude geometries.

Comparing experimental and computational energies is not always straightforward. Molecular mechanics energies are ‘‘steric’’ or potential energies for a single fixed geometry. All experimental energies contain, at the very least, vibrational contributions and are therefore sensitive to the shape of the PES around the minimum. For relative energies it is frequently assumed that vibrational contributions cancel, allowing a direct comparison of MM potential energies with experimental enthalpy differences. A cruder but still common practice is to compare steric energies directly to experimental free energies, ignoring the effects of entropy and usually also of condensed-phase contributions.

Accurate comparisons to experimental enthalpies and free energies can be achieved in molecular mechanics by application of normal mode analysis, solvation models, solid-state calculations, and/or dynamic averaging over large ensembles. Such methods are time consuming and therefore are not easily implemented in a parameter refinement, where each data point is calculated multiple times with different trial force fields. However, the full calculation may be performed once, to derive a correction term allowing the use of the simple potential energy in further calculations. The correction term can be iteratively updated whenever the force field has changed substantially, allowing the use of rapid calculations in the parameter refinement.

Recipe for an Organometallic Force Field

17

2.3. PES data

The exact shape of the PES influences many properties, such as vibrational frequencies and the magnitude of distortions in strained structures. An exact representation of the PES will implicitly allow reproduction of structures and energies. It has been shown that a force field can be successfully derived from PES data alone (32).

Experimental information about the local PES around minima can be derived from vibrational spectroscopy. Employing a harmonic approximation, the vibrational modes and frequencies can be calculated by diagonalization of the mass-weighted Hessian (the matrix of Cartesian second derivatives of the energy) of a structure at an energy minimum. Unfortunately, it is by no means as easy to back-transform experimental frequencies to a Hessian. Experiments do not give any direct information about the vibrational modes. An exact assignment of all frequencies requires spectra of several isotopomers and extensive iterative fitting. Experimental frequencies are valuable in validation, but for parameterization of MM force fields, it is more efficient to find other types of PES data. A very attractive alternative is to use data from QM calculations. At correlated levels (e.g., MP2 or B3LYP), QM frequencies are close to the experimental results (33). Quantum mechanical methods also allow PES determinations at nonstationary points. Furthermore, both QM and MM methods determine structures as energy minima of nuclei on a PES, alleviating the need for conversion of bond-length types (4,29). Finally, the parameter refinement can be performed without time-consuming energy minimizations, because the energies and energy derivatives are calculated at fixed geometries (32).

2.4. Electrostatics

The largest difference between force fields is probably how they handle electrostatics. Each force field uses its own definition of what functions and data should be used. The well-known MM2 force field describes all electrostatic interactions by bond dipoles (4), but most other force fields utilize atomic point charges. The charges may in turn be obtained from fragment matching (34), from bond-type- dependent charge flux (35), or from more complex schemes that can also respond to the environment (36).

Neither atomic charges nor bond dipoles are observables. About the only experimental data for isolated molecules that can be used as parameterization reference are molecular dipoles and higher multipole moments. Substantial effort has also been expended to find electrostatic schemes that can rationalize the behavior of condensed phases (37). However, electrostatic data may be more conveniently obtained from QM calculations. Several schemes exist for partitioning the electron density into atomic charges (38). In general, methods that reproduce the QM-calculated electrostatic field outside the molecular surface are preferred,