J. Phys. Chem. 1992,96, 5138-5145
5138
for the vitreous and freezable water fraction with glass transition behavior reported in this work is related to the hydration value known to be necessary for full activity of the protein. Measurement of the water attached to a protein by hydrodynamic methods consistently gave higher values than those obtained by most other techniques. Kuntz and Kauzmann4 conclude that hydrodynamic hydration, by including trapped water, is expected to be greater than (or equal to) hydration measured by most other techniques. For hemoglobin hydration, values of 0.75.0.62, and 0.52(g of water)/& of protein) were determined by sedimentation, visaxsity, and diffusion measurements (from Table XVI in ref 4). As pointed out in ref 4, there is no simple explanation for this lack of agreement. Notwithstanding these difficulties, we propose that the two at fmt sight unrelated techniques, i.e. hydrodynamic methods and calorimetric measurements at subzero temperatures of vitreous water consisting both of unfreezable and of freezable water as described in this work, “see” the same type or types of water in the protein. In proteins containing up to =0.4 (g of water)/(g of protein), a dynamical transition is observed at =180 K’O3” which is absent in the dry protein.M It would be important to probe the additional influence of the vitreous, but freezable, type of water described in this work on the protein’s dynamics by making similar studies with protein samples containing higher amounts of (vitrified) water.
Acknowledgment. We are grateful to the Austrian Science Foundation for financial support and to Prof. G. P. Johari for suggesting this study and for many discussions. Registry No. MetHb, 12646-21-8; H20, 7732-18-5.
References and Notes (1) Rupley, A.; Gratton, E.; Careri, G. Trends Biol. Sci. 1983, 8, 18. (2) Finney, J. L. In Water & Aqueour Solurionr; Neilson, G. W., Enderby,
J. E., Eds.; Colston Papers No. 37; Adam Hilger: Bristol, Boston, 1985; p 227. (3) Frauenfelder, H.; Gratton, E. Methods Enzymol. 1986, 127, 207. (4) Kuntz, 1. D., Jr.; Kaunnann, W. Adu. Protein Chem. 1974, 28, 239. (5) Cooke. R.; Kuntz, I. D. Inr. Rev. Biophys. Bioeng. 1974, 3, 95. (6) Kuntz, I. D., Jr.; Brassfield, T. S.; Law, G. D.; Purcell, G. V. Science 1969, 163, 1329. (7) Kuntz, I. D., Jr.; Brassfield, T. S. Arch. Biochem. Biophys. 1971,112, 660. (8) Kuntz, I. D., Jr. J . Am. Chem. Soc. 1971,93,514; Ibid. 1971,93,516. (9) Hofer, K.; Mayer, E.; Johari, G. P. J. Phys. Chem. 1990, 94, 2689. (10) Zwart, A.; Buursma, A.; van Kampen, E. J.; -burg, B.; van der Ploeg, P. H. W.; Zijlstra,W. G. J. Clin. Chem. Clin. Biochem. 1981,19, 457. (11) Poole, P. L.; Finney, J. L. Methods Enzymol. 1986, 127, 284. (12) Hallbrucker, A.; Mayer, E. J . Phys. Chem. 1987, 91, 503. (13) Enke, C. G.; Nieman, T. A. Anal. Chem. 1976,48, 705A. (14) Stephens, R. B. J . Non-Cryst. Solids 1976, 20, 75. (15) Feltz, A. Amorphe und glasartige anorganische Festkgrper; Akademiaverlag: Berlin, 1983; p 51. (16) Franks, F. In Water, a Comprehensiue Treatise; Franks, F.,Ed.; Plenum Press: New York, 1982; Vol. 7, p 265. (17) Johari, G. P.; Hallbrucker, A.; Mayer, E. Nature 1987, 330, 522. (18) Hallbrucker, A,; Mayer, E.; Johari, G. P. Philos. Mag. B 1989,60, 179. (19) Hodge, I. M.; Berens, A. R. Macromolecules 1982, 15, 762. (20) Hofer, K.; Mayer, E.; Hodge, I. M. J. Non-Cryst. Solids 1992,139, 78. (21) Angell, C. A.; Tucker, J. C. J . Phys. Chem. 1980,84, 268. (22) Wunderlich, B. J . Phys. Chem. 1960, 64, 1052. (23) Brandts, J. I. J . Am. Chem. Soc. 1964.86,4291. (24) Franks,F.;Hatley, R. H. M.; Friedman, H. L. Biophys. Chem. 1988, 31, 307. (25) Franks, F.; Hatley, R. H. M. Cryoletters 1985, 6, 171. (26) Careri, G.; Gratton, E.; Yang, P.-H.; Rupley, J. A. Nature 1980,284, 572. (27) Rupley, J. A.; Gratton, E.; Careri, G. Trends Biochem. Sci. 1983, 18. (28) Dater, W.; Bachleitner,A.; Dunan, R.; Hiebl, M.; Ltischer, E. Biophys. J . 1986,50, 213. (29) Hofer, K.; Hallbrucker, A.; Mayer, E.; Johari, G. P. J . Phys. Chem. 1989, 93, 4674. (30) Parak, F. Methods Enzymol. 1986, 127, 196, and references cited
therein.
(31) Dater, W.; Cusack, S.; Petry, W. Nature 1989, 337, 754.
Application of the Diffusion Equation Method of Global Optimization to Water Cluaters R. J. Wawak, M. M. Wimmer, and H. A. Scheraga* Baker Laboratory of Chemistry, Cornell University, Ithaca, New York 14853-1 301 (Received: December 6, 1991; In Final Form: February 26, 1992) In this paper, we investigate the behavior of the recently developed diffusion equation method (DEM) for water clusters. The DEM is a global optimization method that involves the deformation of the potential energy surface to remove local minima of high energy. The deformed potential energy surface correspondsto the solution F(x,t) of the diffusion equation with the original energy function f(x) as the initial condition, i.e., F(x,O) =fix). Starting from the few local minima that survived for large values of t ( f can be interpreted here as time), consecutive local ” i z a t i o n s are carried out for gradually decreasing times. The trajectories of these time-reversing procedures lead to deep minima at time t = 0. Hopefully, one of them is the global minimum. The interactions between water molecules are more complex than those for argon clusters (Lennard-Jones potential) treated earlier by the DEM. The well-known MCY potential was used here to describe the water pair interactions. The diffusion equation was solved analytically in the space of the Cartesian coordinatesof all of the atoms, and then a subset of this space (for which the bond lengths and bond angles of the HzO molecules were fixed) was explored. The DEM WPS applied to systems containing up to n = 8 water molecules. For n I5, the method led to the global minima; for n = 8, we obtained a structure insignificantly different in energy and geometry from the global-minimum structure. For n > 5, the method failed to locate the global minima, but in these cases we found a strong indication that the minima given by the DEM are definitely wider than the known global minima. This result is similar to one obtained by Kostrowicki and Scheraga in an application of the DEM to oligopeptides. Keeping in mind the ultimate intended application of the DEM to the global optimization of polypeptides, two general conclusions about the method result from the present work (i) wider minima have a preference over narrower minima to survive in the deformation process; (ii) the trajectories in the timereversing procedures may split into two or more branches at some time t , and an improper decision as to which branch to follow may lead to a minimum of higher value. An efficient method to detect and examine those regions along the trajectories should be designed. Introduction various &entiflc problems such the low-temperature physics of spin g l w , l the m i a m p i c description of prebiotic evolution,2 To whom correspondence should be addressed.
0022-3654/92/2096-5 138$03.00/0
pattern recogniti~n,~ and protein folding4 frequently require the identification of the global minimum among numerous local minima of a relevant multivariable function. Different methods have been developed4 to attack this multiple-minima problem, including the recently introduced diffusion equation method 0 1992 American Chemical Society
The Journal of Physical Chemistry, Vol. 96, No. 12, 1992 5139
Diffusion Equation Method of Global Optimization (DEM).S This approach involves the deformation of the potential energy surface in order to remove the local minima of higher energy (see also ref 6). The original potential energy function is taken as the initial condition for solving the diffusion equation, and the deformation of the original potential energy surface corresponds to the solution of this equation. The DEM has been applied successfully to many standard test functions such as the Goldstein-Price, Hartman, Shekel, and Griewank function^.^^' It was also used to find the global-minimumstructures of argon clusters of up to 55 atoms with Lennard-Jones pairwise interactioms In this paper, we investigate the performance of the DEM in the more complex system of water clusters. Research on clusters has expanded recently because of the very efficient combination of molecular-beam techniques with absorption and mass and it is possible to study the dependence of the structural properties of clusters as a function of their size. Infrared predissociation spectra of water clusters of various sizes"J2 gave clear evidence of a crossover from the behavior of a system dominated by the lowest energy structure to liquidlike behavior (superposition of many possible structures) at cluster sizes somewhere between 6 and 19 molecules. Only comparatively small water clusters can be described by their global minimum-energy structures. Larger water clusters are best represented by a superposition of many thermally accessible local energy minima. Therefore, we restrict our study to clusters consisting of eight or fewer water molecules.
Diffusion Equation Method The diffusion equation (sometimes also called the heat equation) has the formI3 AxF(x,t) = aF(x,t)/at (1) with the initial condition: F(x,O) =Ax), where
t is greater than zero and is interpreted as time, and x = (xI, x2, ..., x,) E R". The functionflx) is the initial potential energy function, and xl, x2, ..., x, are independent variables that describe the positions of all atoms in the system under consideration. For example, for argon clusters x denotes the Cartesian coordinates of all argon atoms; for chain molecules with fixed bond lengths and bond angles, xlrx2, ...,x, denote all dihedral angles along the chain. In other systems, x may consist of Cartesian as well as angular coordinates. If we define #(x) as
~ ( x =) [ 2 ( ~ t ) I / ~exp[-llxl12/4t] ]-~ (2) then the solution to the diffusion equation withfix) as the initial condition can be represented as the convolution of the functions Ax) and #(XI: F(XJ) =flx)*g;"(x) =
s,
fly) F ( x - Y) dY
(3)
It should be noted that if t approaches zero, then # converges to the Dirac delta function in Rm,so that F(x,O) =Ax). Intuitively,fix) can be thought of as a distribution of a physical quantity such as concentration or temperature in R". Then, the diffusion process described by the formula for F(x,t) will usually lead to a smoother distribution of that quantity. The number of local minima of the function F(x,f)should decrease with increasing t. The philosophy of the DEM is based on the following assumption: only significant minima (especially the global one) will survive during the diffusion process. Therefore, the multipleminima problem is significantly simplified for the function F(x,t) at t = to (where to is large enough) and can be solved with traditional numerical methods. In the ideal case, only one local minimum remains (as for argon clusterss). Once all those minima are known, a kind of "time-reversing procedure" has to be applied because the positions of the surviving minima are usually not preserved during the diffusion process: the time is reduced stepwise starting from to. For each time tk, a local minimization is per-
formed for the function F(X,fk) by starting with the position of the local minimum reached in the minimization for time tk-1. The procedure is repeated until time t = 0 is reached. For further details, we refer to the original papers on the DEM.5-7*s The crucial point in the DEM is that the diffusion equation has to be solved analytically because a numerical solution would require a number of grid points which grows exponentially with the dimension of the problem. In many practical cases (e.g., chain molecules with fixed bond lengths and bond angles, or water clusters with rigid water molecules), it is extremely difficult or even impossible to solve the diffusion equation analytically in the space of independent variables as described above. This means that the application of the DEM in the aforementioned form is not possible. In those cases, the following modification of the DEM is introduced: (a) we treat the potential energy function in terms of the Cartesian coordinates of all the atoms in the model; (b) we solve the diffusion equation in these variables analytically (this means that all atoms are treated as if they were free, and the diffusion process evolves in the whole space of all Cartesian coordinates); (c) a subset Q of the space of all Cartesian coordinates is defined by those points that correspond to fixed bond lengths and bond angles. We introduce independent variables 6 = (t1, ..., &) where each t1,..., t k may take on any real value, and we define a transformationX Rk Q. More precisely, each of the values E Rkcorresponds to a point on Q and for any point w E Q there exists at least one € E Rksuch that X(€)= w. In practical cases, tl,..., t k might be Cartesian as well as angular coordinates. The local minimizations are carried out in these variables, i.e., we minimize the function F [ X ( € ) , f ]with , t: as the variable. Hence, we take into account the function F(x,f)restricted to points x belonging to the subset Q. The choice of tl,.,., Ik and a precise definition of the subset in the case of water clusters with rigid water molecules will be given in the section Choice of Coordinate System.
-
Technical Details We apply the DEM to clusters of water molecules of various sizes and use the MCY potentiall4*I5to describe the pairwise interactions between water molecules. MCY Potential. We use the MCY potential for the following reasons: (a) the potential is well-known and provides comparatively good results for the physical properties of bulk water;15(b) it is rather simple in form and therefore covenient for the investigation with the DEM; (c) extensive minimization computations on water clusters of up to eight molecules have been carried out with the MCY potential using a simulated annealing technique;16 hence, our results may be compared with the energy minima found in this earlier investigation. The MCY potential is a "true dimer" potential" in the sense that it is derived from a fit of a simple potential function to the results of an ab initio CI quantum mechanical calculation (including only the intermolecular electron correlati~nsl~) on pairs of water molecules. The water molecule is represented by a rigid four-point model consisting of HI, H2, M, and 0. In this model, an auxiliary point M carries the negative charge of the oxygen; the positive charges are located on points HIand H2. The geometry is indicated in Figure l. The MCY potential energy of a cluster of size n is
f = I C hj
(4)
SKjSn
where i and j refer to H 2 0 molecules, andfj (treated now as a pair potential) is a sum of all elements of the following matrix:l5
HI
H2
M
O
fm fm
f;lo+fHo
H,
fHH
fHH
H2
JHH
JHH
M
fm
fm 0 fAo+ JHo 0 foo
o -f&+fHo
fm
f;lO+fHO
9
(5)
5140 The Journal of Physical Chemistry, Vol. 96, No. 12, 1992
Wawak et al.
where
and q2 = 170.9389 A kcal/mol UHH
bHH = 2.760 844
= 666.3373 kcal/mol
b+HO= 2.961 895 A-I
aho = 1455.427 kcal/mol aio
Figure 1. Sketch of the MCY water model, including the virtual point M.IJ CbH = 0.9572 A; 0-M = 0.2677 A; H U H bond angle = 10 4 . 5 O .
biO = 2.233 264 A-'
= 273.5954 kcal/mol
am = 1088213.2 kcal/mol
boo = 5.152712 A-l
(7)
Only potential energy is considered in this problem. Therefore, in the absence of kinetic energy (i.e., at zero temperature), we simply consider various static configurationsof the system, with no question of "evaporation" of the cluster being involved. With the MCY potential, there are no local minima for clusters with larger intermolecular distances. Truncated MCY Potential. It is clear that the MCY potential is an unbounded function. To avoid singularities in the (+r-') repulsive part (which are regions of high values of the function f) and to avoid an overlap of H and A4 arising from the (-r-l) terms (which is an unphysical global minimum of the MCY potential and would cause numerical problems), we truncated all positive terms in eqs 6 at r = 1.5 A and all negative terms at r = 0.01 A. This means thatfHH(r) = fHH(S+), s+ = 1.5, for r < 1.5 A,fiO(r) = f i o ( s - ) , s- = 0.01, for r < 0.01 A, and so on. The value of the limt 1.5 A was taken reasonably large to decrease the influence of very high values of the function$ On the other hand, this number had to be significantly lower than the smallest expected atom-atom distance (the H...O distance of a hydrogen bond, which equals approximately 1.8 A for water clusterslet6). Then, the truncated MCY potential coincides with the original MCY potential in all regions of physical interest from the point of view of global optimization. The truncated MCY potential with s+ = 1.5 A and s- = 0.01 A [denoted by MCY(s+,s-)] can be represented by eqs 4 and 5, but now fHH
= EL(q*;s+) + EX(aw&wi;s+) fHM
= EL(-2q2;s-)
f Go = EX(a&,,b+H,;s+) fMM
f GO = EX(-aio,bio;s-)
(8)
= EL(4q2;s+)
foo = EX(ac,c,,boo;s+) where EL and EX denote electrostatic and exponential terms, respectively, and plr
EL@;s)(r) =
pls
r
5
(9)
r 5 s
r
!M= cos p cos y
-cos
p sin y
-sin j 3
a sin y - sin a sin j3 cos y
cos a cos y
+
p cos y
sin a cos y
- cos
a sin y
t
cos a sin
1
sin a sin p sin y
-sin a cos p
a sin p sin y
cos a cos p
(11)
used to calculate the Cartesian coordinates, in the main frame, of all atoms in the structure from the independent variables. This is how we define the transformation X mentioned in the Diffusion Equation Method section. The subset Q is obviously the set of Cartesian coordinates of all atoms in the structure that corresponds to all possible arrangements of rigid water molecules. Minimizrtioa procedure. Since the local minimizations in the time-reversing procedure usually start close to the local minima which we intend to reach at that particular time, we used the well-established Gill-Murray algorithm'*because of its quadratic convergence in those regions. It utilizes both the gradient and the Hessian of the function to be minimized. The solution to the diffusionequation (in Cartesian coordinates) with the MCY(s-,s+) potential as the initial condition, can be represented (in terms of the Fourier-Poisson integral)* with the help of eqs 4, 5, and 1-10, but now, instead of EL(p;s) and EX(a,b;s), we should write the convolutions EL(p;s)%&" and EX(a,b;r)*&". Further details are explained in Appendix A. We again recall that this solution, F(x,t), is calculated under the assumption that there are no intramolecular constraints for the water molecules, i.e., the atoms are free and unconstrained by bonds. On the contrary, the local minimizations are carried out in the independent variables
€ = (XI ,Y',Z, ,a1,19',Y,,...,X",Y",z",a",~n,?")
s
EX(a,b;s)(r) = [ae+
x axis and the hydrogen atoms are confined to the xy plane. The coordinates of all atoms in the main-frame system are defined by specifying the position of the local coordinate system. The position of the local coordinate system in the main-frame system is defined by the three Cartesian coordinates of its origin and three (Eulerian type) rotation angles (a,&?) for successive rotations around the local x, y, and z axes. Our independent variables, written as I = (X~,Y~~Z~,~I,BI,YI, ...,Xny"~",~~,~",?n), together with rotation matrices AI,..., An,each of the form given in eq 11, are
s
Choice of Coordinate System. Since the water molecules in our model are rigid, we need six independent variables to describe the position of each molecule uniquely. We proceed as follows: each water molecule is assigned to a local coordinate system so that the oxygen atom is at the origin, the point M lies on the positive
which means that now all the constraints of fixed bond lengths and bond angles have to be taken into account. Consequently, the fmt and second derivatives which are used by the Gill-Murray algorithm must be calculated in these variables. The analytical formulas for the gradient are somewhat messy but straightforward and are not given explicitly here. The Hessian is computed numerically from the gradient. The energy function for water clusters is a sum of atom-atom interactions where the atoms belong to different water molecules (see eqs 4 and 5 ) . Hence, every mixed derivative of second order (over the variables related to different water molecules io andj,) is given by the second
Diffusion Equation Method of Global Optimization derivative of the single term& over these variables, i.e., all second derivatives for pairs of molecules other than (iojo)are zero. Because of this fact, the numerical computation of the Hessian can be carried out efficiently and with high precision. Since translation and rotation of the whole system do not change its potential energy, the corresponding six degrees of freedom lead to numerical problems during the minimization. To avoid these degeneracies, the oxygen atoms of the first, second, and third molecules are placed at the origin of the coordinate system, on its x axis and in the x-y plane respectively (for clusters of four or more molecules; otherwise we fm only the first molecule in the coordinate system in the standard position, i.e., xl = y , = zI= al = 81 = y, = 0). Coustraint on Cluster Dimensions. After some test computations, we observed that the clusters ‘exploded”. This means that the diameter of the cluster, with the molecules interacting according to the transformed potential F(x,t), grows with time t. In other words, there exist local minima of the function F[X(€),t] at large t that correspond to large equilibrium distances. To avoid this effect, we introduced an auxiliary attractive function u(r) that allowed us to limit the expansion of the cluster. This function was added only to 0-0 interactions and was treated in the same way as the physical interaction terms. It equals zero for r < 100 A (it does not change the intermolecular interactions within physically reasonable distances), and introduces an attraction at r > 100 A; Le., u(r) = (u,/r)abgwith some u,,b,, > 0. This shape of u was chosen because of convenience in further calculations, and the coefficients a,, 6, were taken arbitrarily. Numerical experiments showed that with u, = 100 A kcal/mol and b, = 0.01 Awl,the diameter of the clusters never exceeded 55 A for any time t ; this maximum diameter was reached when t was close to 20.0, and the diameter decreased monotonically with increasing t for larger times because of the influence of the function u(r). This means that the influence of the unphysical interactions related to the function u(r) becomes too strong in those regions and, for this reason, we chose to = 20.0 for clusters of all sizes. Local Minim at tcrTo find all the local minima of the function FIX((),to] (Fincludes the transformed interactions coming from the function u(r), as well as the transformed truncated MCY potential) for time to = 20.0, we carried out two runs, each of 10 local minimizations, starting from random structures for each cluster size ranging from two to eight molecules. The random structures were chosen under the assumption that their diameter did not exceed 100 A. In the first run, we obtained one to three different local minima, depending on the size of the cluster. In the second run, only one additional local minimum was found for n equal 6,7,and 8. This indicates that we had found all, or almost all, local minima at to. Behavior during Time-Reversal. The next task was to carry out the time-reversing procedures for all those local minima at to. The time steps used for the time-reversal should not be too large to avoid a jump to a neighboring trajectory. On the other hand, small time steps require a larger number of local minimizations. To proceed safely and economically, we need a tool to detect all “dangerous” regions along the trajectories considered. In our calculations, we used a positive function u of time t related to a given trajectory which shows how fast the structure changes with time and enables us to decide on the sizes of the time steps. To define u(r), we recognize six different types of intermolecular atomatom distances: H-H, 0-0,M-M, MU-0,H--0, H-wM (the point M is treated here as an atom). For each of these six classes, we order all distances according to their magnitude. For structures related to times t and f - At, we calculate the absolute differences of the corresponding ordered pairs of values and take the maximum of these differences. The function u(r) is the sum of these six values (corresponding to the six classes mentioned above), divided by the current time step Ar. Hence, u ( t ) is in a sense a “derivative” of the geometry of the structure over time t and measures how fast this structure changes. Any sharp peak in u ( t ) at a given time t indicates that the structure changes significantly faster than for the neighboring times. Since it is rather arbitrary to define what “sharp” means, we designate any
The Journal of Physical Chemistry, Vol. 96, No. 12, 1992 5141 TABLE I: Comwted Prowrties of Water Clusters 2 3 4 5 6 7
1 2 2 2 4 2
1 1 1 1 1 2
8
4
2
-5.868 -15.838 -26.280 -34.363 -42.156 -5 1.442 -50.644 -62.156 -56.241
-2.93 -5.28 -6.51 -6.87 -1.03 -1.35 -1.23
2A 3A 4A SA 6A 1A
-1.17
8A 8B
-1.03
nr
nr -6.51 -6.81 -1.25 -1.68
4A 5A 6GL 7GL
-8.41
8GL
-1.68 -8.41
1GL 8GL
1B
Additiorial Results 7 8
-50.609 -61.039 -64.415
-1.23 -8.38 -8.05
7C 8C 8D
“Number of water molecules in cluster. bNumber of different configurations found in two runs (10 local minimizations in each run) starting from random structures at time to = 20.0. cNumber of structures at time t = 0 obtained after executing the time-reversing proccdure. dEnergy of the structure finally obtained (in kcal/mol). LEnergy of the structure finally obtained per water molecule (in kcal/mol). ’Symbol indicating the image of the corresponding structure shown in Figures 2-6. ZGlobal minimum energy per water molecule as reported in ref 16 (in kcal/mol). nr = not reported. Symbol indicating the image of the corresponding global minimum energy structure shown in Figures 2-6. peak that is at least twice as high as the neighboring values of as ‘sharpn. If the function u ( t ) does not exhibit any sharp peak in a certain region of t, we say that u ( t ) is smooth in this region. Since from our experience the largest changes in the function F(x,t) usually occur when the time t is close to zero, it is desirable to adopt smaller time steps in this region of t. We apply the following time schedule: u(t)
At = 1.0 for
1.0 < t I 20.0
At = 0.1 for 0.1 < t I 1.0 At = 0.01 for 0.01 < t I 0.1 At = 0.001 for 0.0 I t 5 0 . 0 1
and proceed according to the following rule: As long as the function u(t) is smooth, we assume that we move along the mathematical trajectory demanded by the diffusion equation, but we repeat the computation with 10 times smaller time steps in the time regions where peaks appear. The Occurrence of a sharp peak might indicate that we jumped to a different trajectory because of too large time steps. In the case when our trajectory has to be corrected and we move along a new path, we proceed again according to the rule outlined above. This method to check whether we stay within the proper trajectory turned out to be safe and efficient. To prove this, we repeated some computations with a time step which was even 100 times smaller than the initial time step and found no change in the trajectories. Nevertheless, this method does not detect ‘smooth bifurcations”, i.e., regions where the trajectory splits in two without being indicated by a sudden change in u ( t ) . To detect a bifurcation, one would probably have to use a kind of a second or third ‘derivative” of geometry over time t, with very small Ar along the whole trajectory. Such a procedure would be extremely time consuming. Further details and related examples are presented in Appendix B. Numerical Results and Discussion The results obtained by the procedure described above are summarized in the first part of Table I. We see that the DEM failed to find the global minima as reported in ref 16 for clusters of six, seven, and eight water molecules. Since our approach does not provide any information about the existence of smooth bifurcations, it is possible that their occurrence is partially responsible for those unsatisfying results. In the regions where the trajectory of a minimum splits smoothly into two (or more) trajectories with
Wawak et al.
5142 The Journal of Physical Chemistry, Vol. 96, No. 12, 1992 2A
3A
4A
SA
7A
Y 7GL
Figure 2. Global minimum-energy structures for n = 2-5 (2A-5A). Hydrogen bonds are drawn as continuous lines for H-0 distances 12.3
A.
6A
Fgrre 4. Comparison of the DEM structures (7A-7C) with that of ref 16 (7GL) for n = 7. Hydrogen bonds are drawn as continuous lines for He-0 distances 12.3 A.
6GL
8B Figure 3. Comparison of the DEM structure (6A) with that of ref 16 (6GL) for n = 6. Hydrogen bonds are drawn as continuous lines for H-0 distances 12.3 A.
decreasing time, the path that is chosen by the reversing procedure depends only on slight differences in the structure. These differences may be of the order of numerical inaccuracies. Therefore, we may lose the proper branch of the trajectory leading to the global minimum during the reversing procedure. Indeed, we observed that a very slight perturbation of the structure (of a magnitude of some A in interatomic distances), at a time t l , in a few cases resulted in a completely different minimum at a time t2 < tl;i.e., two almost identical configurations at t1led to different configurations at t2 in the time reversals. Such a phenomenon, which is not detected by a sharp peak in u(t) between times tl and t2,and therefore called a smooth bifurcation, leads to a new trajectory and usually results in a new minimum at f = 0. Of course, the function u(t) was also computed for this new trajectory in order to ensure that no sudden changes of the structure occur. For clusters of seven and eight molecules, smooth bifurcations were accidentally detected, and different minima were found at t = 0. These additional results are given in the second part of Table I, but we emphasize that they were not found by using a systematic approach. Our work demonstrates that the problem of smooth bifurcations is of crucial importance in applications of the DEM. For example, for clusters of eight molecules, we obtained surprisingly good results by a closer, but in some sense an accidental, investigation of such bifurcations. A systematic and not too time-consuming
7
Figure 5. Illustration of the DEM structures (8A and 8B) for n = 8, obtained by the systematictime-reversing procedure. Hydrogen bonds are drawn as continuous lines for H-0 distances 12.3 A.
approach to this problem is very desirable and is presently under investigation. The reason for the failure to find the global minima might possibly be due to an insufficient random search at to = 20.0. Therefore, we carried out test computations starting from the known global-minimum structures of ref 16 at t = 0 and increasing the time stepwise until t = to = 20.0. We found that, in all cases, the trajectory obtained by this "timeadvancingprocedure" jumped into one of the trajectories coming from the local minima found earlier by the random search. This may serve as an indication that the random search at to was sufficient. The computed structures of the clusters are shown in Figures 2-6. No structures were reported in ref 16 for n = 2 and 3. For the structure corresponding to the global energy minimum of the water dimer, we found the same result as that obtained by Clementi and co-workers15(and shown in Figure 2). For the water trimer, an energy minimum search19utilizing the MCY potential
Diffusion Equation Method of Global Optimization 8C
The Journal of Physical Chemistry, Vol. 96, No. 12, 1992 5143
TABLE II: Shape of Minima symbolb 6A 6GL 7A 7B 7c 7GL 8A 8B 8C 8D 8GL
rf 6 7
8
FE/nC -7.03 -7.25 -7.35 -7.23 -7.23 -7.68 -7.77 -7.03 -8.38 -8.05 -8.41
geom a$ 2.2 x loz8 1.9 x 1035 1.7 X 1.8 x 1037 3.1 x 1038 1.3 x 1 0 4 3 2.0 x 1048 9.6 x 1033 1.5 x 1053 2.3 x 1053 1.9 x 1054
'Number of water molecules in cluster. bSymbol indicating the structure (and its image in Figures 2-6). CEnergyof the structure per water molecule (in kcal/mol). dGeometricalaverage of the products of the eigenvalues of the Hessian at the minimum.
8GL
Figure 6. Comparison of the DEM structures (8C and 8D) for n = 8, obtained by perturbation of the structures along the time-reversing trajectory, with that of ref 16 (8GL). Hydrogen bonds are drawn as continuous lines for H.-O distances 52.3 A.
resulted in a structure with the same energy as the cyclic water cluster shown in Figure 2. The finding that the lowest energy form of the trimer is a cyclic structure is also confirmed by an earlier study by Clementi and co-workers20using the predecessor version of the MCY potential. Our computed structures for n = 4 and 5 (Figure 2) agree with those of ref 16. Our structure for n = 6 differs somewhat in energy but significantly in appearance from the global minimum energy structure of ref 16 (Figure 3); similarly for n = 7 (Figure 4). For n = 8, our structures and energies for 8A and 8B (Figure 5 ) differ significantly, but those for 8C and 8D (Figure 6) are very similar to the global minimum energy structure of ref 16. It is obvious that the chance for a minimum to survive during the diffusion process depends not only on the depth of this minimum but also on the shape of the surrounding potential surface. If the values of the function under consideration are low in this neighborhood, it is more likely that the minimum will survive. Consequently, wider minima have better chances to win than narrower ones with a similar depth. For a better understanding of our results, we estimated the width of all minima representing interesting structures. The product of all eigenvalues of the Hessian,calculated at the minimum, can serve as such an estimate. Unfortunately, since the Hessian is calculated in the independent variables Q = (xiyi ,zi ,T I ,...,Xn,Yn,znrLTnrBn,Yn), this product depends on the orientation of the cluster in the main coordinate system. This means that two minima representing the same geometry (thus, in this sense, equivalent) may produce different products of eigenvalues because of different orientations of the clusters. For example, if one of the angles 0 equals *90°, then the orientation of the corresponding water molecule depends only y (see matrix 1 l), and one of the eigenvalues of the Hessian on a i is zero; i.e., we gain one degree of freedom. Therefore, we selected 50 different randomly rotated positions of each minimum-energy cluster around the origin of the coordinate system. In the rotated structures, the oxygen atom of the first water molecule remained
0
0
9
10
II
12
t Figure 7. Parts of the plots of v ( t ) computed for one of the trajectories with n = 7: (-) time steps A f = 1.0; (---) time steps At = 0.1.
at the origin of the coordinate system. The oxygen atoms of the second and third molecules were placed on the rotated x axis and in the rotated xy plane, respectively. We then calculated the product of the eigenvalues of the 50 related Hessians and accepted the geometrical average of those products as a measure of the width of the minimum. The higher this average is, the sharper is the minimum. (Unexpected zero eigenvalues coming from the mentioned singularity of the matrix 1 1 may sometimes cause serious problems in local minimizations; a zero eigenvalue means that the function is flat in that degree of freedom, and a cannot handle this situation. In such cases, it was necessary to rotate the cluster in the main-frame system in order to find an orientation of the cluster for which all B values of the constituent water molecules were significantly different from f90°,) The results for the shape of the minima under consideration are presented in Table 11. The table shows that in all cases in which the DEM found minima significantly higher in energy than the higher energy "aare much wider than the global "a, the corresponding global minima. This means that in these cases the global minimum structures are much less flexible than the structures found with the help of the DEM (it is also in full agreement with the intuition derived from the images of those structures given in Figures 3-6). The fact that wide minima are preferred, compared with narrow ones, has also been observed in the application of the DEM to oligopeptides.21
Timing All the computations were carried out on an IBM 3090 supercomputer, using one processor. As expected, the numerical expense of our algorithm is only of the order of O ( d ) (in the case of clusters with painvise interactions,many terms in the exprespions for the gradient and Hessian equal zero, and their computation could be avoided by careful coding of the algorithm). To carry
5144 The Journal of Physical Chemistry, Vol. 96, No. 12, 1992
Wawak et al. Ax) =
400
0
t
12
1
14
16
20
I 0, where g ( x ) = [2(ar)1/2]-mexp[- 11~11~/4r]
(A4)
I
