INTRODUCTION

TOP ABSTRACT INTRODUCTION TESTS OF FUNCTIONAL FORM SMALL CLUSTER MINIMUM-ENERGY... CONCLUSION REFERENCES

Atomic-level simulations of proteins interacting with transition-metal ions have the promise of elucidating a wide variety of biophysical phenomena (Kaim and Schwederski, 1994). Such simulations can help establish both the native structure and reaction paths of metalloproteins and explain the roles of metal ions in protein folding and stability. Similar applications can be seen for simulations of metal ions interacting with DNA and RNA. However, the utility of such simulations depends critically on the availability of accurate but computationally tractable force fields for metals interacting with proteins. There are no quantitatively accurate force fields for these purposes, and it is probably not possible to construct quantitatively accurate force fields of a simple enough form to be useful. However, one can hope to generate force fields that include the basic physical effects at a level of accuracy sufficient to ascertain chemical trends as metals or ligands vary. The situation in this regard is best for "simple" metal ions that have no partly occupied shells. Here, an ionic picture based on electrostatics, supplemented by empirical repulsive forces, may hope to treat the most important physical effects. Even here, though, care should be taken, because the bond is not completely ionic but has some partial covalent character. For transition metals, the situation is much more difficult. The ligand-field splitting of the d-shell leads to important electronic contributions that cannot be ignored. These are manifested, for example, in the typical square or tetragonal coordination of Cu²⁺ complexes. Such effects cannot be described by radial interactions, but rather require the introduction of angular forces describing the energetics of the transition-metal d-shell. Most existing methods for including angular forces in simulations of metal ions have used assumed functional forms (Comba et al., 1995; Allured et al., 1991; Timofeeva et al., 1995; Wiesemann et al., 1994; Sayle et al., 1995, 1997) based, for example, on the observed structures of small complexes. A recently developed method (Landis et al., 1998) treats the environment-dependent energetics of transition metals via a valence-bond approach. This method appears promising for covalently bonded metals. However, the bonding state of transition metals in proteins may be rather different from the sdⁿ hybridized configuration assumed in that work. In the present analysis, we seek to develop a method suitable for a broader range of chemical environments ranging from ionic to covalent.

The angular overlap model (AOM) (Gerloch et al., 1981; Jørgensen et al., 1963) provides a conceptual basis for treating ionically bound transition metals. This model provides a straightforward way of determining the parameters of an effective quantum-mechanical d-d Hamiltonian h in terms of the positions of the transition metal's ligand neighbors. The Hamiltonian parameters are given as simple angular functions of these positions. However, the total energy in the AOM does not have a simple analytic form. The ligand-field part of the energy is given as

(1)

where the are the eigenvalues of h (and only the filled levels are included). The are not given as explicit functions of the h, but are instead obtained using numerical matrix methods, except in certain highly symmetric cases. A total-energy method based on diagonalization of a d-d Hamiltonian obtained from the AOM has been combined with a classical molecular-mechanics method to study the structures of small molecules (Burton et al., 1995).

Our goal is to find an expression for the total energy that is simpler than that obtained by straightforward application of the AOM or more complete quantum-mechanical methods. We feel that this simplification is useful for four reasons:

1.	Interpretability. The results of molecular-dynamics and energy-minimization calculations are much more useful if they can be easily interpreted in terms of the features of the underlying force fields. In standard methods, such as those based on the AOM, the matrix-diagonalization step provides a "black box" that obscures the energetic origins of the structural features that are found. By expressing the entire energy in a simple form, one can often associate characteristic atomic configurations with important minima in the real-space potential-energy function.
2.	Ease of programming. One often wants to extend an existing force field to include transition-metal ions. The energy functions that we derive here are of a type not too different from those already present in standard force fields. For this reason, inclusion of these energy functions in such force fields is simpler than for straightforward applications of the AOM. We note that the present approach also yields some saving in CPU time, but, given that most biomolecules will contain transition metals only as a minority component, this savings is probably not a major factor.
3.	Predicting structures of new molecules. A simple real-space picture of the force field allows one to have an a priori idea of what structures are likely to form around a transition-metal ion surrounded by a given combination of neighbors. In cases of a simple environment, such as a copper ion surrounded by four -bonded neighbors, one knows the likely structures even without examining the force field. However, in more complex environments, the constraints of the force field can provide very useful information.
4.	Justification of classical descriptions of quantum-mechanical energetics. Many papers have been published that use classical energy descriptions to describe open-shell transition-metal systems. The type of analysis given in the present paper allows one to assess the accuracy of such calculations depending on the angular form of the interactions that are used.

We have recently shown (Carlsson, 1998) how the ligand-field splitting effects of -bonded transition-metal ions in ionic configurations can be described by a force field of a fairly simple form. One starts with an explicitly quantum-mechanical form for the electronic energy, and then solves the quantum-mechanical problem with systematic approximation methods to extract the real-space description of the electronic bonding energy. A second-order treatment of the hybridization terms between the ligand and metal orbitals is used to generate a ligand-field Hamiltonian for the d-electrons. The electronic bonding energy of this Hamiltonian is analyzed in terms of its moments. The moment analysis of the energy gives a "semiclassical" energy function with a simple trigonometric angular form. It is not precisely an additive function of ligand-ligand interactions, but is nearly as simple computationally. This energy function gives quite accurate energy results for -bonded ligands, with less than 10% error in the ligand-field stabilization energy (defined below).

This approach is extended in the present paper in three ways. First, we develop energy functions that treat -bonded ligands as well. This results in new and distinct angular forms. Second, we extend the analysis to treat covalent environments. Finally, we present results for the lowest-energy structures of small clusters, which are a guide to understanding the bonding preferences of transition metals in proteins. The organization of the remainder of the paper is as follows. The next section develops the formalism underlying the angular forces. The following section gives tests of the functional form of the angular forces by comparing with results from diagonalization of simple cluster Hamiltonians. The subsequent section gives the small-cluster results for ideal examples and for an analog of the copper environment in blue-copper proteins.

Model

We model the electronic structure of the environment of a single transition-metal ion in terms of orbitals |L, µ localized on the ligands and |M,  localized on the single metal ion. The orbitals are taken to be orthogonal for simplicity of calculation. Here, L denotes a particular ligand atom and the |L, µ are distinct orbitals on that atom. The Hamiltonian takes the form,

(2)

Here, the  terms are on-site energies for the orbitals, and the h terms are hybridization energies between the ligand orbitals and the metal orbitals. We ignore electron-interaction terms and explicit electrostatic effects. The energy associated with this Hamiltonian can then be obtained by diagonalizing its matrix as given in the |L, µ  |M,  basis, and taking the sum of the energies of the occupied eigenstates.

To obtain a real-space description of the energetics of this Hamiltonian, we treat the ionic and covalent cases separately. In the ionic case, we make the simplifying approximation that the dimensionless ratios |h|/(_M,  L,µ) are small and can be used as expansion parameters. This allows the application of ligand-field theory, as in the AOM, which gives an effective Hamiltonian for the d-shell,

(3)

where

(4)

Note that we have assumed all of the d-orbitals on the transition-metal ion to have the same energies before ligand-field effects are "turned on." The angular overlap model provides formulas for the h^µ in terms of the angular and radial coordinates of the ligands.

In a fully quantum-mechanical ligand-field-theory calculation, one would numerically diagonalize the matrix _d, as in Burton et al. (1995). To obtain a force field with a simpler form, we instead work with the moments of _d, defined in terms of the traces of its powers, as follows. The first moment is the average energy of the d-complex, or

(5)

where

(6)

Thus is given as a sum of independent contributions from the ligands. Because the d-shell by itself has spherical symmetry, the contribution from each ligand is a radial function (no angular dependence) of the metal-ligand distance. Our major interest is in the angular terms resulting from the ligand-field splitting, so we do not consider the term further.

The width W of the d-complex corresponding to the ligand-field splitting is, in the simplest picture, described by the second moment or variance ² of the eigenvalues of _d. We expect that W  . We can obtain a simple real-space form for ² as follows. First, we note that ² =  _n(_n  )² = Tr(_d  Î)², where the _n are the eigenvalues of _d. From Eqs. 3, 4, and 6, we see that

(7)

Then the variance is

(8)

We will focus on the ligand-field stabilization energy, which is defined as the sum of the energies of the occupied orbitals relative to , or E_LFSE = '_n(_n  ), where the sum is over only the occupied eigenfunctions. We assume that E_LFSE is proportional to the ligand-field splitting W, and also identify W with the standard deviation . Thus, our expression for the ligand-field stabilization energy has the form

(9)

where A is a dimensionless constant that will later be used as a fitting parameter.

The accuracy of this form will be demonstrated in the next section through specific numerical experiments. At present, we will show how this form for E_LFSE results in an expression for the energy in terms of simple angular interactions between the ligands. We thus need to develop analytic forms for the g_'. First, we note that we can write

(10)

where _d = _µ|M, µM, µ| is the projection onto the d-subspace of the transition-metal ion. Here, we define | = _dH|L, , and we have used the relation  = _d, which holds because _d is a projection operator.

Consider first the case (Carlsson, 1998) where both the |L,  and |L', ' orbitals have -character with respect to the metal ion. Then they couple only to the d-orbitals |M,  and |M, ' that have -character with respect to the bond axes, so that | = h|M,  and |' = h'|M, ', where h and h' are the appropriate coupling strengths. Thus g_' = hh'M, |M, ', and the angular dependence is contained in the last inner product. However, this is simply a matrix element of a rotation about M, which carries L into L'. We note that |M, is equivalent to |M, m = 0, where m = 0 denotes the angular dependence of the spherical harmonic Y_2m. Choosing our coordinate system so that L is along the z-axis and L' is in the z-x plane, we find that

(11)

where the D-term is a matrix element of the l = 2 representation of the rotation group, and the second equality follows from the explicit formulas of the D-terms given in Wigner (1959). In summary, for -bonded ligands,

(12)

For the case of a -bonded orbital L and a -bonded orbital L', we consider only the case in which the axis of the orbital lies along the circle connecting L and L'; if it is perpendicular to this circle, L and L' have different inversion symmetries so their coupling vanishes. We write | = h|M,  where |M,  = (1/)(|M, m = 1 + |M, m = 1) and m is the usual azimuthal angular momentum index for the spherical harmonics. (The minus sign comes from the definition of the spherical harmonics). Then, again choosing a coordinate system in which L is along the z-axis and L' is in the z-x plane, and following reasoning parallel to the -bonded case, we see that

(13)

and

(14)

Finally, we turn to the case of two -bonded orbitals. To obtain the subsequent results, it is sufficient to consider the case in which the orbitals are parallel to the arc connecting L and L', and the case in which they are perpendicular to it. By reasoning similar to that above, one sees that, in the first case,

(15)

and

(16)

In the second case,

(17)

and

(18)

Thus, combining Eqs. 9, 12, 14, 16, and 18, we find the following explicit form for E_LFSE as a sum of ligand-ligand interactions:

(19)

where the ligand-ligand interaction U is given as:

For - interactions,

(20)

For - interactions,

(21)

For - interactions,

(22)

where

(23)

(24)

(25)

(26)

(27)

In each case, U is given as a product of radial terms involving the two ligands and a simple angular function. (Note that, although denoted an interaction here, U does not have units of energy because of the square root in Eq. 19). In the calculation of the terms involving - interactions, each term involves a sum over two -orbitals, parallel and perpendicular to the arc connecting the two ligands. In the latter case, g_' in Eq. 9 vanishes, but g and g_'' do not. In the calculation of the terms involving - interactions, one has a similar scenario except that one sums over two pairs of -orbitals.

These forms are plotted out in Fig. 1 A. Note that the - interactions are fairly similar in form to the - interactions, both having pronounced minima at 180° (as well as the physically irrelevant one at 0°) and a shallower minimum at 90°. The - interaction is complementary to these, having minima at 45° and 135°. We shall see later that these differences lead to large differences in ground-state structures of small clusters. For comparison, we show, in Fig. 1 B, an empirical angular interaction curve (Comba et al., 1995), assumed in some previous calculations of small transition-metal structures. The angular dependence is based on the observed square structure of small complexes of the transition metals of interest, and is proportional to sin²2. Because the metals that have square coordination generally have predominantly -bonds to their neighbors, the most relevant comparison is to the - curve in Fig. 1 A. We see that the behavior is quite different. The empirical curve has equivalent minima at 180° and 90°, whereas, in the theoretical curve, the 180° minimum is much deeper.

View larger version (12K): [in this window] [in a new window]   FIGURE 1   (A) Interaction between ligands of central transition metal ion. The dimensionless quantity u (cf. Eqs. 19-21) is either u, (solid line), u, (dashed line), or u, (dotted line). The quantity u is plotted so that negative values will correspond to lower energies. (B) Empirical force field from Comba et al. (1995).

The above analysis provides theoretical support for the use of the form of Eq. 9 for the ligand-field stabilization energy in the ionic limit where the electronic structure can be described accurately by an effective Hamiltonian for the d-complex. We now derive an energy function of equal simplicity for covalently bonded systems. We begin with a Hamiltonian of the form of Eq. 2, in the covalent limit in which _L,µ = _M, = 0. (The choice of the common energy to be zero affects no physical results and is taken for algebraic convenience.) We derive the energy function using a moment analysis, which has been used previously to develop energy functions for elemental metals (Carlsson, 1990). The moments of the electronic density of states, projected on the d-shell, are given as

(28)

where the | are the eigenfunctions of and the factor of simplifies the algebra. In the completely covalent case considered here, all of the eigenfunctions have equal squared amplitudes on the ligand and d-states, so that

(29)

For the special case of no ligand-field splitting, in which all of the bonding states have the same energy , and all of the antibonding states have the energy , there are altogether 10 eigenfunctions (plus potentially other ligand orbital combinations that are not coupled to the metal ion). Using Eq. 29, one sees that, in this case,

(30)

It has been established (Cyrot-Lackmann, 1968) that the moments can be written as a real-space sum over neighbor positions as

(31)

(32)

The moments of odd order vanish when both on-site energies vanish. These expressions involve functions quite similar in form to those that appear in the ionic limit. Specifically,

(33)

We now derive an approximate energy function from these quantities. In the covalent case, the ligand-field stabilization energy is not as well defined as in the ionic case, because the d-complex is not well separated from the ligand orbitals. For this reason, we focus on the total interaction energy E_int between the metal ion and the ligands. With our choice of energy zero, this is simply the sum of the eigenvalues of all of the occupied states. To derive the energy function, we note that the assumption of only one bonding and one antibonding eigenvalue, described by Eq. 30, corresponds to the case in which ⁴_d = ⁴ = ². Thus a natural form for the energy would be one that begins with Eq. 30 and adds a correction based on ⁴_d  ². The simplest form for the total energy based on dimensional analysis is

(34)

where A and B are coefficients to be determined. The second term is taken to have a square-root dependence because, in systems with narrow bonding and antibonding complexes, the width of these complexes is proportional to . (One readily shows this by taking a model density of states consisting of two rectangular subbands, each of width W, centered at energies +₀ and ₀. In this case, straightforward calculation shows that ⁴_d  ² = 4W² + 4W⁴/45, which reduces to 4W² for narrow bands.) We also find that this provides an excellent fit to the quantum-mechanical results for small clusters given in the following section.

Using Eqs. 33, we have for the total energy

(35)

This form is quite similar to Eqs. 9 and 19 but differs in two ways. First, the term based on ²_d makes a contribution involving a sum of single-ligand contributions. This term is not sensitive to the angular character of the local environment, because it is determined entirely by the distances and types of the ligands. Second, the ⁴_d term is divided by a factor involving the ²_d term. Computationally, these differences incur no significant extra costs. We note that, in distinguishing potential structures with similar distances and types of neighbors, the energy functions, Eqs. 9 or 19 and 35, are equivalent, because (_L, g) is determined by the neighbor distances and types. Thus, the dependence of the structure on the angular aspects of the environment is given by a function of the form Eq. 19 in either case, and this is used in the next section.

A physical system will be somewhere between the limits of complete ionicity and complete covalency. A suitable function for interpolating between these two limits is obtained by replacing the quantity (_L, g_,) in the first term and in the denominator of the second term of Eq. 35 by (_L, g_,) + (_d  L,)². With this form, one sees that, in the ionic limit (large _d  L,), the second term of Eq. 35 becomes equivalent to Eq. 9, and the first term simply describes the average energy of the d-complex. One can also show (Carlsson, unpublished) that, for the case in which the widths of the bonding and antibonding complexes go to zero (no ligand-field splitting), the first term of Eq. 35 thus modified gives the correct interaction energy for all values of the ionicity.

	TESTS OF FUNCTIONAL FORM

TOP ABSTRACT INTRODUCTION TESTS OF FUNCTIONAL FORM SMALL CLUSTER MINIMUM-ENERGY... CONCLUSION REFERENCES

To evaluate the accuracy of the semiclassical form of Eq. 19 for the ligand-field stabilization energy, we have performed explicit tests for small clusters. These clusters consist of a central transition-metal ion with four neighbors placed at random orientations at random distances relative to the central ion. We consider only the minority-spin orbitals, because these determine E_LFSE for high-spin late-transition metals. Four electrons are placed in these states, corresponding to Cu²⁺; other band-filling values give similar results. The random distances are taken into account by varying the couplings h and h uniformly over a finite interval ranging from zero to h or h. We take h = 0.5h. The energies are obtained by explicit diagonalization of a tight-binding Hamiltonian for this cluster. All of the ligand orbitals are taken to have the same value of _L,. The only dimensionless variable that enters the results is then  = |h_max/(_d  L,)|. For small values of , the bonding is primarily ionic, and, for larger values, it acquires more covalent character.

In our tests of the energy function for ionic systems, we use  = 0.1. We fit the energies of the clusters to the semiclassical form for the energy,

(36)

involving two parameters A and B (where U contains A). We use a database of 10,000 clusters to determine the parameters, and then test them on a set of 1000 clusters not included in the "training" set. Typical results are shown in Fig. 2, which shows the semiclassical energies versus exact energies for 1000 clusters. Figure 2 A corresponds to three -bonded ligands and one -bonded one, whereas Fig. 2 B corresponds to four -bonded ligands. The rms errors in E_LFSE for these two cases (and for the cases of the two and three -bonded ligands) are between 9% and 10%. This is to be compared with rms errors of 25-30% that are obtained with empirical force fields (Carlsson, 1998).

View larger version (22K): [in this window] [in a new window]   FIGURE 2   Energies obtained by semiclassical force field versus exact energies for small model clusters. Energies given in units of hmax, maximal coupling strength between ligands and transition metal. (A) One -bonded ligand and three -bonded ligands, ionic system. (B) Four -bonded ligands, ionic system. (C) One -bonded and three -bonded ligands, covalent system. (D) Four -bonded ligands, covalent system.

In the tests for covalent systems, we use the same clusters but with  =  corresponding to vanishing energy differences. In this case, we use an energy function of the form given by Eq. 35, and vary A and B to obtain the best description of E_int. For clusters with zero, one, two, three, and four -bonded ligands, the rms errors in E_int are 2.6%, 3.3%, 3.2%, 3.8%, and 3.6%, respectively. Figure 2 C and D, show the energy data for the cases of three -bonded ligands and one -bonded one, and four -bonded ones. In comparing these results with those for the ionic limit, one should note that, although the percent errors are smaller, the values of E_int are, in general, larger in magnitude than those of E_LFSE, so that the absolute errors may not necessarily be smaller. Nevertheless, an accuracy of a few percent in E_int is quite encouraging. To compare the errors in E_LFSE and E_int, consider a case in which the eigenvalues coming from the ligand and metal orbitals are uniformly spread out from the top to the bottom of the complex. If one considers the d-complex as being the top half of this complex, and the holes are concentrated in the top quarter of the complex, then the magnitude of E_LFSE will be roughly a third of E_int. The absolute errors in the covalent and ionic limits cases are then roughly comparable.

	SMALL CLUSTER MINIMUM-ENERGY STRUCTURES

TOP ABSTRACT INTRODUCTION TESTS OF FUNCTIONAL FORM SMALL CLUSTER MINIMUM-ENERGY... CONCLUSION REFERENCES

In this section, we describe some of the implications of the angular forms developed above for the structure of small model clusters consisting of a metal atom and four ligands. Such cluster calculations cannot treat a protein environment accurately; this would await parameterization and incorporation of the force field into protein codes, which is in progress in our group. However, from the small-cluster calculations it is possible to see the structural preferences of the angular forces by themselves. These can have a large impact on the choice of binding sites by particular ligands, and on the corresponding affinities. We emphasize that the goal of these calculations is to obtain the simplest description of the energetics that is both reasonably accurate and physically justifiable. Therefore, we expect that more elaborate methods, such as a complete application of the AOM, would give results similar to the present ones.

To keep the calculations as simple as possible, we place the ligands at frozen bond lengths from the central metal ion. Here "frozen" means that they are fixed at a given set of values, but these values are not necessarily the same for all of the ligands. The energy terms include the ligand-field stabilization energy as described by Eq. 19, and a repulsive radial interaction between ligands forbidding close approaches. (We note that, because the bond lengths are fixed, the ² terms in Eq. 34 do not vary during the simulation, and, as pointed out in the second section, Eq. 19 can then be used to describe both covalent and ionic systems.) The repulsive ligand-ligand term is needed because otherwise spurious structures involving 45° bond angles can appear when -bonding ligands are present. To evaluate the magnitudes of the couplings, we assume a value of 100 kJ/mol = 1.04 eV for the ligand-field stabilization energy of a cluster with four -ligands in square coordination; this number would correspond to experimental values for the more strongly stabilized complexes (Cotton and Wilkinson, 1972). For clusters with unequal bond lengths, we assume an exponential decay, so that e(r) = e(r₀)exp[(r  r₀)], where r₀ is the reference distance and  is a decay parameter. Because A multiplies the couplings e and e, any changes in A can be compensated for by changing the overall magnitudes e and e. For simplicity, we therefore take A = 1. Using Eq. 19 and the above value of E_LFSE, we then find that e(r₀) = 0.79 eV. In the blue-copper protein fragment that we treat, the -bonding effects are quite strong, and, in the absence of reliable estimates for e, we take e = e. The effects of varying the couplings are discussed below. Because the most prominent case that we consider of a distant ligand is sulfur, we identify  with the spatial decay rate coming from the measured first ionization energy E₁ of sulfur, 10.4 eV (= 1000 kJ/mol), using the formula ²²/2m = E₁. This yields  = 1.65 Å¹. The ligand-ligand terms contain an exponential term taken as the repulsive part of the van der Waals interactions as given in the MM2 force field (Sprague et al., 1987). In our model clusters, we use parameters and bond lengths typical for copper or nickel interacting with nitrogen or sulfur ligands. This is because copper and nickel have the largest ligand-field energies among the 3d transition metals, and typical ligands for these metals are nitrogen and sulfur. We have examined three simple geometries; the corresponding minimum-energy clusters are shown in Fig. 3 A-C.

1.	A cluster with all four ligands -bonded and placed at equal distances of 2.0 Å. We take repulsion parameters appropriate for nitrogen. The result is a square-planar cluster (cf. Fig. 3 A), consistent with the observed structures of many copper and nickel complexes. With this parameterization, the square structure is quite strongly favored over the tetrahedral one, by 0.31 eV or 7.0 kcal/mole. The square structure remains stable under increases of the ligand-ligand repulsion strengths of up to a factor of over three. We note, however, that the stability of the square structure is strongly sensitive to the choice of repulsive terms. For example, we find that, when the van der Waals term from the OPLS-II force field (Jorgensen and Tirado-Rives, 1988) is used instead of the MM2 form, an increase of the repulsion strengths of only 20% is enough to stabilize the tetrahedral structure. We expect that input from fully quantum-mechanical total-energy calculations will be necessary to precisely pin down the competition between the electronic energy and the repulsion terms in determining the structural energy differences.
2.	A cluster similar to that of 1, but with one -bonded ligand. This leads to a substantial deviation from planarity, as shown in Fig. 3 B, where the -bonded ligand is atom 4. As expected from the maximum of the - potential in Fig. 1, the angle between atom 4 and atoms 1-3 is greater than 90°, about 123°. The angle among the -bonded atoms 1, 2, and 3, is slightly greater than the ideal value of 90°, because of the repulsive energy term.
3.	A cluster with two -bonded nitrogen ligands 1 and 2, one -bonded sulfur ligand 3, and a -bonded sulfur ligand 4. The distances are 2.04 Å for the nitrogen ligands, 2.18 Å for the -bonded sulfur ligand, and 2.64 Å for the -bonded sulfur ligand. This geometry is motivated by the observed geometry of copper sites in blue-copper proteins. In these proteins, the nitrogen ligands and a cysteine sulfur ligand are close to the copper, and it is believed that the cysteine sulfur ligand is primarily -bonded. An additional -bonded methionine sulfur ligand is farther from the copper. We note that experiments and calculations indicate that the copper binding in this system has a high degree of covalency (Solomon et al., 2000). However, because Eqs. 19 and 35 are equivalent for energy minimization with fixed bond lengths, this is not a problem. The values of the ligand distances here are taken from cluster calculations for blue-protein models (Ryde et al., 1996). As above, the van der Waals parameters are taken from the OPLS-II parameter set, but we do not have a reliable procedure for determining the differences between the ligand-ion interactions of the nitrogen and the sulfur ligands. Because the greater distance of the sulfur at 2.18 Å from the copper will be compensated to some extent by the greater size of sulfur relative to nitrogen, we simply assume that the three close ligands have the same coupling strength to the central atom. We have varied the ratio of the sulfur-to-nitrogen coupling strengths by up to 30% in both directions and found changes of only a few degrees in the bond angles of the cluster. For the far ligand, we use the scaling procedure described above, which leads to a coupling that has 50% of the strength of the close ligands. We have varied this coupling from 0% to 75% of the close-ligand value, and again found bond-angle changes of only a few degrees. The lowest-energy structure for this cluster is indicated in Fig. 3 C. The geometry is trigonal, with bond angles of 125° between the cysteine sulfur and the nitrogens; the nitrogen-nitrogen bond angle is 91°. By comparison, the optimal values for protein models obtained by quantum-chemical calculations (Ryde et al., 1996) are 125° for the cysteine sulfur-nitrogen bond angle and 103° for the nitrogen-nitrogen bond angle. The results obtained by the present method, of course, do not have quantitative accuracy, but the overall structure of the coordination shell is quite similar to that obtained by Ryde et al. (1996), which is shown in Fig. 3 D. The formation of this structure is not due to the bond lengths that we used as input. To demonstrate this, we have considered a cluster with the same bond lengths as the blue-protein fragment, with all ligands -bonded. In this case, the structure is planar, as shown in Fig. 3 E. Thus the formation of the trigonal structure is directly related to the special character of the angular interactions associated with -bonding.

View larger version (36K): [in this window] [in a new window]   FIGURE 3   Lowest-energy structures obtained in relaxations of five-atom cluster using angular forces. (A) Four equivalent -bonded neighbors. (B) Three -bonded neighbors and one -bonded neighbor 4, all at the same distance. (C) Two -bonded nitrogen neighbors 1 and 2, a near -bonded sulfur neighbor 3, and a far -bonded neighbor 4. (D) Structure for blue-protein model obtained by quantum-chemical calculations (Ryde et al., 1996). (E) Three near -bonded neighbors and a far -bonded neighbor.

	CONCLUSION

TOP ABSTRACT INTRODUCTION TESTS OF FUNCTIONAL FORM SMALL CLUSTER MINIMUM-ENERGY... CONCLUSION REFERENCES

We have seen that it is possible to extract the functional form of angular forces around transition metals in biomolecules by using an approximate treatment of the ligand-field Hamiltonian in ionic systems or the cluster Hamiltonian for covalent systems. The resulting angular forms have simple trigonometric forms, which are very different for -bonding as compared to -bonding. Ligand-field energies in ionically bonded systems are represented well by these angular forces, and much better than by force fields with assumed functional forms. The interaction energies of covalently bonded systems are also obtained accurately. Future work should aim to include the functional forms derived here in widely used biomolecular simulation packages.