J. Phys. Chem. 1995, 99, 11805-11812
11805
ARTICLES Molecular Dynamics on Deformed Potential Energy Hypersurfaces Jaroskw Pillardy and Lucjan Piela* Quantum Chemistry hboratory, Department of Chemistry, University of Warsaw, Pasteura I , 02-093 Warsaw, Poland, and Interdisciplinary Center of Mathematical and Computational Modelling, University of Warsaw, Banacha 2, 02-093 Warsaw, Poland Received: December 4, 1994; In Final Form: April 16, 1995@ In the classical molecular dynamics for a complicated system the trajectory remains relatively close to any arbitrarily chosen starting configuration of the nuclei. This is because of a prohibitively large volume of the configurational space as well as of energy barriers between numerous potential energy basins. In the present paper the molecular dynamics is performed on a deformed potential energy hypersurface. The deformation is tunable, and when increased, the barriers become lower or even disappear; therefore, the system is able to explore much larger regions of the configurational space. The deformation is slowly released, and finally the molecular dynamics is continued on the undeformed potential energy hypersurface. As a result of such a procedure, some low-energy structures may be found in a more effective way. The method has been applied to the clusters of the N Lennard-Jones atoms, N = 5 , ..., 66, which is the most extensive systematic study.of the clusters published till now. The energies of the clusters are lower than those obtained by the classical molecular dynamics (at the same numerical effort) for the starting points being far away from the lowest energy region. Contrary to the molecular dynamics, they are independent of the starting point to a significant degree. The size of the most stable cluster tumed out to be a stepwise increasing function of N, which reflects building of new atomic shells. However, for N = 38 some qualitative change of the structure to the one characteristic for the infinite crystal occurs and the cluster size becomes abnormally low when compared to the clusters with N < 38 and N > 38.
1. Introduction The electronic ground-state energy of a molecule in the adiabatic approximation plays a role of the potential energy for the motion of the nuclei and is a complicated function of their configuration. A proper quantum-mechanical calculation of this function for medium- and large-size molecules is beyond the capabilities of current computers. The problem is that one needs to carry out such calculations for a number of points in the configurational space, and the number goes up exponentially with size of the molecule. Moreover, the energy is a nonconvex function of the nuclear coordinates, Le., there is usually more than one stable configuration. In fact, for large molecules there is myriad of stable conformations, and the question of which of them corresponds to the lowest energy is a very difficult and unsolved problem of theoretical chemistry. The larger the molecule the more critical becomes the question. In such a situation one is forced to abandon quantum mechanical approaches and to use semiempirical or empirical force fields, i.e., some approximate expressions for the electronic energy of a molecule as a function of its nuclear positions. There are many force fields in use in the literature, e.g., ECEPP,' AMBER,2and CHARMM.3 All of them take into account only qualitative and simple features of the proper potential energy function. To compare theory and experiment one has to choose among a vast number of stable configurations one that corresponds to the lowest free energy. The free energy is a global thermodynamical quantity-at a given temperature it is a number
* To whom correspondence should be addressed. @
Abstract published in Advance ACS Abstracts, July 1, 1995.
corresponding to the whole hypersurface (or even some manifold of them). In practice, in thermodynamics the time scale is important, and therefore, the free energy is usually connected to some limited region of the hypersurface, which the system is able to explore during a given time period at a given temperature. Therefore, the free energy of an isomer becomes a function of the corresponding region of the hypersurface. What counts in this region is how many low-energy states are localized in it. First of all, this depends not only on its lowest-energy well but also on how wide the well is. The wider the well, the more numerous are the accessible vibrational states and, therefore, the more important is the entropy factor in the free energy. The last question is not addressed in the present paper. In order to find the lowest-energy well for a given system (a flexible molecule or a cluster of atoms or molecules), a number of theoretical methods have been designed! In almost all of them the potential energy expression contains both the shortand long-range interactions, and for this reason the hypersurface becomes extremely complicated. Much simpler hypersurfaces result when only the electrostatic interactions are taken into a c ~ o u n t . ~Nevertheless, even in this case the hypersurface remains prohibitively complicated to be effectively explored. A way to single out the low-energy structures is to use the molecular dynamics (MD). At a given temperature the point representing the system under study is running on the hypersurface exploring regions of the potential energy values accessible for the temperature used. The lowest potential energy found on this trajectory may be used as an upper bound for the global minimum value. The problem is that the MD trajectory (although able to surmount some small potential energy barriers between basins) practically remains close to the starting point
0022-365419512099-11805$09.00/0 0 1995 American Chemical Society
11806 J. Phys. Chem., Vol. 99, No. 31, 1995 used. In other words, the system is trapped in some potential energy basin, and this basin may be far away from the basin of the most stable configuration. For example, in the MD simulations for biomolecules the starting point is usually chosen to be the X-ray determined geometry, and the system always stays relatively close to this conformation during a run of the molecular dynamics. In the present paper a deformation of the potential energy hypersurface followed by the MD is proposed in order to find low-energy structures of the system under study.
Pillardy and Piela 45
- ..-
15
,.I
25
P s
15
05
0 5
.l 5
2. Deformation of the Potential Energy Function The idea comes from the global optimization methods, where the global minimum of a nonconvex function V(x) is searched.6 The function is deformed into V(x,t) using a deformation parameter t > 0, such that V(x,O) = V(x). The deformation procedure is designed in such a way that, when t increases, any basin (potential energy valley) changes its size and shape and becomes more and more shallow. Some basins, especially the shallow ones, disappear and others grow at their expense. Continuing the process, for sufficiently large t > tmm, one hopefully gets a single basin of the global minimum of V(x), which has grown at the expense of all other basins. The position of this single minimum of V(x,tmax)is attainable by a gradient trajectory algorithm from any starting point. Then an algorithm called the reversing procedure is applied. In the reversing procedure, one minimizes V(x,tmax-6), 0 < 6 *: 1, with respect to x using as the starting point the minimum position of V(x,t,-). Thus, a new position x is obtained. Repeating this step many times and therefore gradually releasing the deformation by decreasing t, one finally gets x corresponding to the minimum of V(x,O) = V(x). Hopefully, after the reversing procedure one obtains the global minimum of V(x). The value of can be tmax estimated or calculated.’ A few deformations have been tested in the literature. One is the diffusion equation method (DEM):,’ where the potential energy is smoothed according to the heat equation as if the function V represented a temperature or concentration distribution and t denoted time. Another deformation tested was the so-called shift method (SM),8 where in a pairwise additive potential any two-body (distance dependent) interaction is shifted toward its origin (the shift being ?-dependent), thus also producing a smoothing. Both methods have been able to find the lowest-energy structures for some selected argon In both methods, however, the single minimum emerged at a collapse position of the whole system to a single point. Then, during the reversing procedure a “big bang” took place, and the system evolved usually to the global minimum of the energy. A drawback of this procedure is that the collapse is quite far from the global minimum configuration, whereas the main idea of the method was to have the single minimum already quite close to the global minimum position. For this reason another deformation procedure called the distance scaling method (DSM) is proposed in the present paper, which (contrary to the DEM and SM) does not change the position of the minimum of the Lennard-Jones atom-atom potential. The Lennard-Jones potential is now expressed with the scaled distance r‘ (the scaling parameter being t > 0):
+
r r*t r’ = l+t In the present paper
r* = ro (2) where ro is the distance corresponding to the energy minimum
.2 5
Figure 1. Deformed atom-atom Lennard-Jones potential, u[r‘(r,t)]of eq 3 as a function of r for t = 0, 1, 9;t = 0 means the undeformed potential, r* = ro.
of the usual Lennard-Jones potential,
The position of the minimum of u(r) remains at r = r* for both unscaled and scaled potentials for any value o f t . As it is seen from eq 1 the value of u[r‘(r*,t)]is independent o f t . Thus, by changing the deformation parameter t, one obtains a family of potentials u[r‘(r,t)],each member of the family going through the point (r*,u(r*)). The plot of u [ f (r,t)]calculated for a few values of t is given in Figure 1.
3. Molecular Dynamics on Deformed Potential Energy Hypersurfaces The main purpose of this section is to obtain shallow energy basins in order to allow the system to pass the barriers. Thus, the region of the space visited by the system will get larger, and possibly one may reach such regions of space which would be not visited at the MD with the original potential. This may facilitate finding low-energy regions that may be very far from the starting point. The deformation is applied by changing the deformation parameter t in the formula for the pairwise additive potential energy N
V(%O = C u [ / ( r , . , t ) l
(4)
i 0. (3) The procedure is the same as in point 1, but three starting points equivalent to the lowest-, a medium-, and the highestenergy minima have been chosen. (4) For a long MD run for a chosen t > 0, and then for each minimum found, a reversing procedure is applied (without MD). ( 5 ) For a long MD run starting for a chosen t > 0, the deformation is slightly released after each MD step (the reversing procedure combined with the MD). No cutoffs (of energy or distance) have been used in the present paper. The minimization procedure used at the beginning was the steepest-descent algorithm, followed by application of the conjugated gradient algorithm (until the minimum was found). There are a lot of other adjustable parameters possible for each of the described procedures. The temperature of the system may depend on the deformation parameter t or may be kept constant during the reversing procedure, and the system may be kept at a constant temperature during the MD run (“bath”) or not (pure classical dynamics). The initial starting point for the MD at t = tmaxmay be chosen as a cube or a line. Finally, a step and a starting point of the reversing procedure have to be chosen.
Pillardy and Piela We have tested most of the above mentioned procedures. Usually we have applied tmax =9 and a step 6t = 1, but we have also tried tmax= 16 and a ti series such as (9, 7, 5 , 3, 2, 1.5, 1, 0.5, 0) and many others. Especially successful was a linear dependence of the system temperature T as a function of t, e.g., T starting from 150 K for tmax= 9, and increasing with t, up to 240 K. During the MD run the system temperature was in most cases rescaled after a chosen number of steps (1/ 50 of the all MD steps), but also classical MD runs have been performed. Despite the relatively high temperature, no single case of “exploding” the clusters was observed. Results of the methodologies from points 1 to 4 after a suitable adjustment of the parameters described above were similar, but different energy minima have been often obtained. The most reliable was the procedure of point 2 with the linear dependence of T(t) and T being rescaled. The initial starting point was only marginally important for the MD carried out on the deformed hypersurfaces.
References and Notes (1) Nemethy, G.; Pottle. M. S.; Scheraga, H. A. J . Phys. Chem. 1983, 87, 1883. (2) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. J. Comput. Chem. 1986, 7, 230. (3) Brooks, B. R.; Brucoleri, R. E.; Olafson, B. D.; Slater, D. J; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (4) See a review article Piela, L.; Olszewski, K. A,; Pillardy, J. J. Mol. Struct.: THEOCHEM 1994, 308, 229. (5) Piela, L; Scheraga, H. A. Biopolymers 1987, 26, 533. (6) Piela, L; Kostrowicki, J.; Scheraga, H. A. J. Phys. Chem. 1989, 93, 3339. (7) Kostrowicki, J; Piela, L.; Cherayil, B. J.; Scheraga, H. A. J. Phys. Chem 1991, 95, 4113. (8) Pillardy, J; Olszewski, K. A.; Piela, L. J. Phys. Chem. 1992, 96, 4337. (9) Shalloway, D. In Recent Advances in Global Optimization; Princeton University Press: Princeton, NJ, 1992, p 433. (10) Pearlman, D. A.; Case, D. A.; Caldwell, J. C.; Seibel, G. L; Singh, U. C.; Weiner, P; Kollman, P. A. AMBER 4.0; University of California: San Francisco, 1991. (11) Mackay, A. L. Acta Crystallogr. 1962, 15, 916. (12) Hoare, M. R. Adv. Chem. Phys., 1979, 40,49. (13) Van de Waal, B. W. J. Chem. Phys. 1989, 90, 3407. (14) Wales, D. J. J. Chem. Phys. 1989, 91, 7002. JP9432197
