J . Phys. Chem. 1990, 94, 6619-6626
6619
Parametrization of Molecular Mechanics Calculations for the Accurate Description of Hydrogen-Bonding Interactions James R. Damewood, Jr.,*,+Robert A. Kumpf, Wolfgang C. F. Muhlbauer, Joseph J. Urban,+ and John E. Eksterowicz Department of Chemistry and Biochemistry, University of Delaware, Newark, Delaware 1971 6 (Received: December I I , 1989)
A new method for the parametrization of molecular mechanics force fields for the accurate calculation of hydrogen-bonded structures and energies is presented. The method is general, and applicable to a wide variety of force fields. Applications using the programs AMBER(~.O) and MM2(85) are presented and compared to both high-level ab initio calculations and results from M M ~ ( s ~ ) .
Introduction The accurate description of hydrogen-bondingl interactions is of central importance in molecular mechanics2 modeling. Since many chemical and most biological systems contain these types of bonding interactions, their special significance is easily appreciated. Despite this obvious importance, significant debate still surrounds the best approach to modeling hydrogen bonding within the molecular mechanics framework. In this work, we present an approach to the derivation of parameters necessary for the accurate description of hydrogen-bonding interactions that is not dependent upon the addition of special hydrogen-bonding functions to the force field. Our method is applicable to a wide variety of molecular mechanics force fields ( A M B E R ( ~ . M o)M , ~ ~ ( s ~etc.); ),~ it is straightforwardly implemented and differs from other previously reported approaches in that it requires only accurate (experimental’ or ab initio) binding energies and ground-state intermolecular separations in order to derive useful parameters. A wide variety of approaches to the description of intermolecular hydrogen bonding within the molecular mechanics framework have been reported p r e v i o ~ s l y . These ~ ~ ~ approaches may be categorized into those that employ special functions6 to describe the hydrogen-bonding interaction and those that rely only on the force field equation employed for intramolecular interaction^.^ The successful application of any special function approach requires a priori knowledge of which atoms in the system under study are involved in hydrogen-bonding interactions, something that cannot be unambiguously determined in every case. The types of functions that have been added to the force fields in the first category are varied, and include Lennard-Jones 10-12 potentials,6c~h~l~o exponential or Morse-like potentials,6a~b-g~i-nJ and trigonometriC6a.b,e,i.l.n,p.q.s-w and/or switching functions6dfikto describe the angular or distance dependence of the hydrogen-bonding interaction. For a number of reasons, force fields that do not employ special functions, i.e., those in the second category,’ offer considerable advantages. All trigonometric and/or switching function approaches to special functional forms require a decision (which is sometimes made arbitrarily) regarding the cutoff or angular dependence of the supplementary function. In addition, the use of computationally time consuming special functions for the description of the hydrogen-bonding interactions makes application of the resulting force field to bulk simulation8 difficult given the current limitations on computational resources. The accurate description of hydrogen-bonding interactions without the use of these additional functions, if possible, avoids these difficulties. We have therefore decided to pursue this direction in the current research. The position that special functions are not required for an accurate description of hydrogen bonding in molecular mechanics ‘Current address: IC1 Pharmaceuticals Group Structural Chemistry Section, Medicinal Chemistry Department, IC1 Americas, Inc., Wilmington, DE 19897.
0022-3654/90/2094-6619$02.50/0
calculations has been presented previo~sly.~However, simple parametrization methods for general application of the approach
(1) We define a hydrogen bond for the purpose of this work as a groundstate, attractive interaction with a near linear (Le., X-H-Y) geometry between a polarized hydrogen atom (X-H) and a heteroatom (Y). (2) For reviews of the molecular mechanics method, see: Burkert, U.; Allinger, N. L. Molgcular Mechanics; American Chemical Society; Washington, DC, 1981. Osawa, E.; Musso, H. Top Stereochem. 1982, 13, 117; Angew. Chem., Int. Ed. Engl. 1983, 22, 1. (3) Singh, U. C.; Weiner, P. K.; Caldwell, J. W.; Kollman, P. -4.University of California, San Francisco, 1986. The form of van der Waals equation for this force field is VvDw= c[(R/Rij)’*- 2(R/R,.)6]. (4) Allinger, N. L. MMZ(SS); Quantum Chemistry Program Exchange, Department of Chemistry, Bloomington, IN 47405. The form of the van der Waals equation for this force field is
(5) For extensive reviews of the available experimental data on hydrogen-bonded complexes see: (a) Curt& L. A.; Blander, M. Chem. Reu. 1988, 88, 827. (b) Hobza, P.; Zahradnik, R. Chem. Rev. 1988, 88, 871. (c) Buckingham, A. D.; Fowler, P. W.; Hutson, J . M. Chem. Reo. 1988.88, 963. (d) There exists a great deal of experimental data for hydrogen-bonded complexes, and therefore decisions as to what values to use in parameter development had to be made. In making these decisions we attempted to follow the guidelines for the critical evaluation of experimental data set forth in ref 5a. It is, of course, possible that one could use the methodology presented here to develop molecular mechanics parameters based on other experimental data. (6) See, for example: (a) Scott, R. A.; Scheraga, H. A. J . Chem. Phys. 1966,45, 2091. (b) Ooi, T.; Scott, R. A.; Vanderkooi, G.; Scheraga, H. A. J . Chem. Phys. 1967, 46, 4410. (c) McGuire, R. F.; Momany, F. A.; Scheraga, H. A. J . Phys. Chem. 1972, 76,375. (d) Stillinger, F. H.; Rahman, A. J . Chem. Phys. 1974, 60, 1545. (e) Singh, R. D.; Ferro, D. R. J . Phys. Chem. 1974, 78, 970. (f) Caillet, J.; Claverie, P.; Pullman, B. Acfa Crystallogr. 1976, 832, 2740. (g) Kitiagordsky, A. I. Chem. Soc. Rev. 1978, 7, 133. (h) Gelin, B. R.; Karplus, M. Biochemistry 1979,18, 1256. (i) Taylor, R. J . Mol. Struct. 1981, 1 1 , 31 I . 6 ) Levitt, M. J . Mol. Biol. 1981, 145, 251. (k) Levitt, M. Annu. Rev. Biophys. Bioeng. 1982, 1 1 , 251. (I) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J . Comput. Chem. 1983, 4, 187. (m) Levitt, M. J . Mol. Biol. 1983, 168, 595. (n) Kroon-Batenburg, L. M. J.; Kanters, J. A. J . Mol. Sfruct. (THEOCHEM) 1983, 105,417. (0)Weiner, S. J.; Kollman, P. A,; Case, D. A,; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S., Jr.; Weiner. P. J . Am. Chem. Soc. 1984, 106,765. (p) Goodford, P. J. J . Med. Chem. 1985.28.849. (9) Vedani, A,; Dunitz, J. D. J . Am. Chem. SOC.1985, 107, 7653. (r) Bocquet, J. F.; Chhor, K.; Lucazeau, G.; Dianoux, A. J. Can. J . Chem. 1985,63, 1940. ( s ) Nilsson, L.; Karplus, M. J . Comput. Chem. 1986, 7, 591. (t) Vedani, A.; Dobler, M.; Dunitz, J . D. J . Comput. Chem. 1986, 7, 701. (u) Masek, B. Personal communication, 1988. (v) Vedani, A. J . Comput. Chem. 1988, 9, 269. (w) Boobbyer, D. N . A.; Goodford, P. J.; McWhinnie, P. M.; Wade, R. C. J . Med. Chem. 1989, 32, 1083. (x) Empirical methods not derived within the framework of a general molecular mechanics force field were not included in this description. For example, see: Schroeder, R.; Lippincott, E. R. J . Phys. Chem. 1957, 61, 921.
0 1990 American Chemical Society
6620
The Journal of Physical Chemistry. Vol. 94, No. 17, 1990
TABLE I: Experimental and ab Initio Binding Energies (kcal/mol) and Heteroatom Separations (A) for Use in Parametrization complex AHcrp(ref) PE,,; Rb (ref) AE,,,' -3.59 f 0.5 (20a) -5.6 2.980' (2Ob,c) -5.5 (H20)2 (CH,OH)> -3.5 f 0.2 (21) -5.4 2.954' -5.4 -15.2 (HCO,H), -1 1.7 f 0.1 (22a) -13.7 2.696' (22b) (HCONH2)2 c 2.995' -13.4 -6.8 3.023' -6.6 CHJNH2. -5.59 f 0.2 (24) HOCH, NH,*H,O 2.983' (25) -6.4 CH,O*H20 3.042 -4.6 oCalculated from AH,,, using 6-31GL*//6-31G** frequencies. bHydrogen-bonded heteroatom separation. c6-31G**//6-31G** calculated values. dExperimental value. eSee ref 23 in text.
to readily available force fields have not been presented.I0 This study provides a simple parametrization scheme and allows for an assessment of the accuracy that can be attained by using this type of computational approach.
Molecular Mechanics Model for the Hydrogen Bond Many possible ways exist to describe hydrogen-bonding interactions even within the limitations of not employing special functions. The number and location of intermolecular interaction sites, their charges, and van der Waals parameters remain possible options within the framework of such an approach. Our parametrizations employ atom-centered charges without lone pairs and use no van der Waals radius on hydrogens attached to heteroatoms." Our atomic charges are determined by electrostatic potential surface fit to 6-31G**//6-31G** ab initio calculations.I2 This places our atomic charges on a firm quantum mechanical base and leaves the van der Waals R and 6 parameters for the heteroatoms of the system as the only parameters in need of derivation. While it may be possible to improve on the accuracy of computational results for hydrogen-bonded systems by including non-atom-centered charges or special functions, we have elected to explore how well the simple molecular mechanics model employed in this work can perform when the nonbonded parameters are optimized. This not only allows us to avoid possibly arbitrary decisions regarding the placement of additional interaction sites, but also preserves the capability of applying our parameters to economical molecular dynamics simulations. (7) See, for example: (a) De Santis. P.: Giglio, E.; Liquori, A. M.; Ripamonti, A. Nature 1965, 206,456. (b) Hagler, A. T.; Huler, E.; Lifson, S . J . Am. Chem. Soc. 1974, 96, 5319. (c) Hagler, A. T.: Lifson, S. J . Am. Chem. SOC.1974, 96, 5327. (d) Shipman, L. L.; Burgess, A. W.; Scheraga, H. A. froc. Natl. Acad. Sei. USA 1975, 72, 543. (e) Clementi, E.; Cavallone, F.; Scordamaglia. R. J . Am. Chem. Soc. 1977, 99, 5531. (0 Snir, J.; Nemenoff. R. A.; Scheraga, H. A. J . Phys. Chem. 1978, 82, 2497. (g) Robson, B.; Osguthorpe, D. J. J . Mol. Biol. 1979, 132, 19. (h) Hagler, A. T.; Lifson, S.; Dauber, P. J . Am. Chem. Soc. 1979, 101, 5122. (i) Allinger, N. L.: Chang, S. H.-M.; Glaser, D. H.: Honig, H. Isr. J . Chem. 1980, 20, 51. G) Jorgensen, W. L. J . Am. Chem. Soc. 1980, 102, 543. (k) Allinger, N. L.; Kok, R. A,: Imam, M. R. J . Comput. Chem. 1988, 9, 591. ( I ) Kim, S.; Jhon, M. S.; Scheraga, H. A. J . Phys. Chem. 1988, 92, 7216. (8) For recent reviews and leading references, see: McCammon, J. A.; Harvey, S. C. Dynamics of Proteins and Nucleic Acids: Cambridge: New York, 1987. Computer Simulation of Biomolecular Systems. Theoretical and Experimental Applications; van Gunsteren, W. F., Weiner, P. K., Eds.: Escom: Leiden, 1989. Jorgensen, W. L. Ace. Chem. Res. 1989, 22, 184. (9) See. for example: Ferro, D. R.: Hermans, J., Jr. Biopolymers 1972, / I . 105. Froimowitz. M. J . Med. Chem. 1982, 25, 689. Wolf, R. M.; Suter, U. W. Macromolecules 1984, 17, 669. Islam, S . A,; Neidle, S.;Grandecha, B. M.; Partridge, M.; Patterson. L. H.; Brown, J. R. J . Med. Chem. 1985, 28, 857. References 7b,d,h,k. (10) More elaborate parametrization schemes have been presented; see: Maple, J. R.; Dinur, U.: Hagler, A. T. Proc. Natl. Acad. Sci. USA 1988, 85, 5350. Dinur, U.: Hagler, A. T. J . Am. Chem. SOC.1989, 1 1 1 , 5149. ( I I ) A similar approach is taken in the T I P ~ model P for bulk water; see: Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J . Chem. Phys. 1983, 79, 926 and references therein. ( I 2) For an introduction to ab initio methods, see: Hehre, W. J.; Radom, L.;Schleyer. P. v. R.; Pople. J. A. Ab. Initio Molecular Orbital Theory; Wiley: New York. 19x6.
Damewood et al.
Computational Methodology Molecular mechanics calculations were performed with the ~ previous programs AMBER(3.0),3 MM2(85)? and M M z ( S ~ ) . ' ~Our research has d e m ~ n s t r a t e d 'the ~ need to modify the original electrostatic potential in MM2(85) in order to obtain accurate descriptions of intermolecular interactions. All calculations with the program MM2(85) therefore employed modified electrostatic parameters. Point charges, rather than the default bond moments, were used and were derived from the electrostatic potential surface fitting procedure described below. All calculations with M M ~ ( s ~ ) employed standard bond moments and the default parameters described in the literature. Ab initio calculations were performed with complete optimization at the 6-31G**//6-31G** level of sophistication using the program GAUSSIAN 82 (c82).I5 All frequencies were calculated at the 6-31G**//6-31G** level under the harmonic approximation of G82. Point chargesI6 for AMBER(3.0) and MM2(85) were determined by electrostatic potential surface (EPS) fit to the 631G**//6-31G** wave function by using a modified version of the program C H E L P . ~ ~ All geometries were completely optimized (without constraint) starting from either experimental or ab initio geometries (Table I ) . Parameters employed for complete geometry optimization for water were as follows: for H-0 stretching 553.0 kcal/(mol A*) (AMBER(3.0)) and 4.6 mdyn/A (MM2(85)) and for H-0-H bending 47.0 kcal/(mol rad2) (AMBER(~.O))and 0.35 mdyn A/rad2 (MM2(85)). For formaldehyde, the MM2 parameters employed were 0.370 mdyn A/rad2 for H-C-H bending. Other force field parameters derived for the systems under consideration are reported in the tables that follow. Cubic splines were calculated by using the ESSL Subroutine Library of Programs." All plots were generated with the program DI-3000.18
Experimental and ab Initio Data for Use in Parametrization Using a combination of the best available experimental dataSd and high-level a b initio results on the structure and energy of association for molecular dimers, we established a set of complexes for use in our initial parametrizations. The choice of these systems was determined by both the availability (or accessibility) of relevant data and their inherent theoretical interest. Previously described complexes for which experimental structures are not well-defined (e.g., ammonia dimert9)were not considered in this (13) (a) Reference 7k. (b) Lii, J.-H.; Gallion, S.;Bender, C.; Wikstrom, H.; Allinger, N. L.; Flurchick, K. M.: Teeter, M. M. J . Compur. Chem. 1989, IO, 503. (14) Damewood, J. R., Jr.: Anderson, W. P.; Urban, J. J. J . Comput. Chem. 1988, 9, 1 1 1 , Damewood, J. R., Jr.; Urban, J. J.; Williamson, T.C.; Rheingold, A. L. J . Org. Chem. 1988, 53, 167. (15) GAUSSIAN 82 (G82). Binkley, J. S.; Frisch, M. J.; DeFrees, D. J.; Raghavachari, K.; Whiteside, R. A.; Schlegel, H. B.; Fluder, M. J.; Pople, J. A. A copy of the program may be obtained from Carnegie-Mellon University. (16) Chirlian, L. E.: Francl, M. M. J . Comput. Chem. 1987,8, 894. The program CHELP as described in this reference does not, in some cases, obtain atomic charges that accurately reflect molecular symmetry. Two modifications were made in the program in order to correct this deficiency. A uniform distribution of points used for sampling the electrostatic potential was created by placing the molecule in a box of sample points (with specified, uniform density distribution) such that the edge of the box was no less than 2.8 A from any atom. All points within the van der Waals radius of an atom were then removed prior to derivation of the electrostatic potential. In addition, the Lagrange multiplier subroutine that fits the atomic charges to the electrostatic potential was modified so that symmetry-related atoms were constrained to have the same atomic charge. The charges so derived are as follows: for water, H, 0.3967; 0, -0.7934; for methanol, C, 0.3033; 0, -0.6929; 0-H, 0.4259; C H (anti to OH) 0.0289; C-H (staggering OH) -0.0326; for ammonia, N, -1.0506; H, 0.3502: for formaldehyde, C, 0.4084; 0,-0.4513; H, 0.0215: for formic acid, C, 0.6342; C=O, -0.5536, C-0, -0.6266; C-H, 0.0721; 0-H 0.4739; for methylamine, C, 0.5089; N, -1.0515; N-H, 0.3746; C-H (distal to N-H) -0.1 149; C-H (proximal to N-H), -0.0458; for formamide C, 0.5938: 0, -0.5577: N, -0.8651; C-H, 0.0394, N-H (syn to O), 0.4131: N-H (anti to 0),0.3764. ( I 7) IBM Engineering and Scientific Subroutine Library, Release 3, 1988. ( 1 8) Precision Visuals DI-3000 Contouring System Subroutines, Release 3, 1985.
Description of Hydrogen-Bonding Interaction
The Journal of Physical Chemistry, Vol. 94, No. 17, 1990 6621
C
Figure 1. Binding energy vs R vs c (left) and oxygen atom separation (calc dist) vs R vs c (right) grids for A M B E R ( ~ . O )calculations on the water dimer. R and calc dist in A and binding energy in kcal/mol.
work. The set of compounds employed in the parametrization include water dimer,20 methanol dimer,21formic acid dimer,22 formamide dimer,23methylamine-methanol,24 and ammoniawater.2s Relevant data on these complexes are shown in Table I. Complexation energies ( A E ) needed for the parametrization were calculated from experimental binding enthalpies ( A H ) by using standard methods26 and 6-31G**//6-31G** ab initio harmonic frequency data. Parameters for formaldehyde were also derived in this study after our initial parametrization efforts for the molecules listed above proved to be successful. These parameters were derived by using the ab initio 6-31G**//6-31G** structure and energy for formaldehyde-water complex la27(Figure 4) and allow for a comparison of molecular mechanics calculations on formaldehyde-water complexes with the extensive ab initio data that has been reported p r e v i ~ u s l y . ~ ~
General Method for Parameter Derivation The approach that we employ in the derivation of the van der Waals R and c parameters is dependent upon the complexity (dimensionality) of the parametrization problem. While our (19) See: Nelson, D. D., Jr.; Fraser, 6.T.; Klemperer, W. Science 1987, 238, 1670. (20) (a) Experimental binding enthalpy: Curtiss, L. A.; Frurip, D. J.; Blander, M. J. Chem. Phys. 1979, 71, 2703. (b) Experimental structure: Fredin, L.; Nelander, B.; Ribbegard, G. Chem. Phys. Lett. 1975,36, 375. (c) Dyke, T. R.; Muenter, J. S.J. Chem. Phys. 1974, 60, 2929. (21) Inskeep, R. G.; Dickson, F. E.; Olson, H. M. J. Mol. Spectrosc. 1960, 5 , 284. (22) (a) Experimental binding enthalpy: Henderson, G. J. J . Chem. Educ. 1987,64,88. (b) Experimental structure: Almenningen, A,; Bastiansen, 0.; Motzfeldt, T. Acta Chem. Scand. 1969, 23, 2848; 1970,24,747. We use the value of 2.696 A as calculated from the above experimental data by: Chang, Y.-T.; Yamaguchi, Y.; Miller, W. H.; Schaefer, H. F., 111 J . Am. Chem. SOC. 1987, 109, 7245. (23) In ref 13b, a value of 14.0 kcal/mol was reported (without reference)
as the experimental dimerization energy for the formamide dimer. Previo ~ s l y ,a~binding ~* energy of 7-8 kcal/mol was suggested as an average value for a single amide hydrogen bond. Based on this, some have used 14 kcal/mol as the dimerization energy of formamide See: (a) Pimentel, G. C.;McClellan. A. L. The Hydrogen Bond: W. H. Freeman: San Francisco, 1960. (b) Poland, D.; Scheraga. H. A. Biochemistry 1967, 6, 3791. (c) Dreyfus, M.; Maigret, 8.;Pullman, A. Theor. Chim. Acta 1970, 17, 109. For the experimental structure, see: Stevens, E. D. Acta Crystallogr. 1978, B34, 544. This experimental data was not used in our parametrization since it corresponds to an infinite hydrogen-bonded chain. (24) Millen, D. J.; Mines, G. W. J . Chem. SOC.,Faraday Trans. 2 1974, 70, 693. (25) Herbine, P.; Dyke, T. R. J. Chem. Phys. 1985, 83, 3768. (26) For example, see: (a) Reference Sa. (b) Del Bene, J. E. In Molecular
Structure and Energetics; Liebman, J. F., Greenberg, A,. Eds.; VCH: Deerfield Beach, FL. 1986; Vol. 1. (27) Kumpf, R. A.; Damewood, J. R., Jr. J . Phys. Chem. 1989, 93,4478.
1.80
J I 4
0.45
I
0.80
I
0.75
I
0.80
1
1.05
EPSILON Figure 2. Graphical representation of the location of the point of intersection of the appropriate splines for the water dimer from the grids in Figure 1. Binding ener y contours of -5.580 kcal/mol and oxygen atom separation of 2.980 are shown.
x
parametrization method is presented within the confines of our selected hydrogen bond model, the approach could, in principle, be extended to other molecular mechanics force fields.
Three-Dimensional Systems The parametrization of homodimers (water dimer, methanol dimer, etc.) for R and c is a three-dimensional problem in both binding energy and heteroatom separation. For the parametrization of the 3-D systems, a grid of 10 possible R and t values equally spaced between 1.5 and 1.95 8, (1.55 and 2.00 8, for nitrogen) for R and 0.2 and 1.55 (0.15 and 1.50 for the oxygen of water) for c is generated for the system under investigation. For each possible parameter combination, a complete geometry optimization calculation using the force field to be parametrized is performed. The data obtained from these full optimization calculations allow the generation of two three-dimensional surfaces, binding energy vs R vs c and heteroatom separation vs R vs t. The two graphs obtained from the AMBER(~.o) calculation for the water dimer are shown in Figure 1. Cubic spline functions are used to interpolate between the data points in the two respective graphs. The dimensionality of the problem is then reduced by finding the
6622
The Journal of Physical Chemistry, Vol. 94, No. 17, I990
spline contour on the respective graph that corresponds to the experimental (or ab initio calculated for other systems, Table I) binding energy and heteroatom separation (Figure 2). The analytically determined intersection of these two spline functions (Figure 2) determines R and t. In this work, all intersections were located for AE of complex formation. For heterodimers in which one of the partners has parameters derived from a homodimer calculation, it is possible to derive the necessary parameters for the second component using this three-dimensional procedure. For this type of calculation, the R and t for only the new partner are varied while the R and t for the other, previously determined monomer, remain at their homodimer values. We refer to this approach as the parameter buildup technique.
Damewood et al. TABLE 11: Derived R (A) and 6 (kcal/mol) Parameters for Hydrogen-Bonded Systems
AMBER(3.0) R
t
1.689d l.794d 1.550d 1.827' 1.633' 1.89p/ 1.726 1.6781 1.7119
0.71 1 0.269 0.212 0.218 0.924 0.223 1.509 1.238 0.606
complex (H20h4 (CH3OH)z" (HC02HLb (HCONH2)2b CH3NHZ.HOCH3' N H3*H20C C H 20.HzO'
MM2(85)
R 1.699d 1.878d 1.512d 1.868' 1.660' 1.834' 1.8 16 1.771' 1.7328
I
0.784 0.226 0.325 0.327 1.059 0.456 0.896 0.758 0.570
Parameters obtained from three-dimensional derivation. Parameters obtained from five-dimensional derivation. cParameters obtained by buildup technique. dParameters for hydroxyl oxygen. Parameters for carbonyl oxygen. /Parameters for nitrogen. CParameters for carbonyl oxygen. @
Five-Dimensional Systems The direct (i,e., nonbuildup) parametrization of heterodimers, or of homodimers with two different heteroatoms involved in hydrogen-bond formation (e.g., formamide dimer) for R, c and R f , e' is a five-dimensional problem in both dimerization energy and intermolecular separation. I n order to determine parameters for these five-dimensional systems, surfaces were generated for energy vs R, t and R', e', (tl)and heteroatom separation vs R, c and R', 8 (si). Five equally spaced values were employed for R, t and R', t' with R values ranging from 1.5 to I .9 A and t values ranging from 0.2 to 1.4. All possible combinations of parameter values were submitted for complete geometry optimization to generate the data for the two surfaces. These two surfaces were then fit to a cubic function of the form
< = (c,R + c2t + c3R' + c4t' + c
~ ) ~
with the corresponding 35 coefficients determined by least-squares techniques. As in the 3-D systems, contours were then selected and { = heteroatom that correspond to [ = binding energy separation (12). The values for the binding energy and the heteroatom separation are then subtracted from both sides of the corresponding r s so that the value for both functions is equal to zero. Solutions for R, t and R', t', can then be determined by finding the minimum in the function 6 (below):
(rl)
This approach is necessary since a standard root finding technique is precluded by our system of two equations with four unknowns. Conjugate gradient minimization of 6 results in a surface ( t ) in which = f2 and 6 = 0; Le., there exists an infinite number of nontrivial solutions for R, t, R', and E' that satisfy this condition. To sample t , we then select 10000 points from a regular grid in order to perform minimizations to locate the regions of closest approach between
