Ind. Eng. Chem. Res. 1994,33, 957-964
967
Discontinuous Molecular Dynamics Simulation of Hydrogen-Bonding Systems Jin-Xing Liu,+Thomas L. Bowman 11, and J. Richard Elliott, Jr.*J Chemical Engineering Department, The University of Akron, Akron, Ohio 44325-3906
An efficient discontinuous molecular dynamics simulation algorithm is developed to simulate hydrogen-bonding systems. The efficiency of CPU time in the order of N log2 N is stressed by comparisons to several other methods and algorithms. To evaluate this algorithm, special attention is given to the simulation of thermodynamic and transport properties as well as oligomer distribution of the hydrogen-bonding water system by use of DMD-1B and -4Bsite models. The results are compared to the results of Monte Carlo simulations and Wertheim’s first-order theory. The overall agreement is satisfactory within a standard deviation.
Introduction Hydrogen-bonding and associating mixtures represent a broad class of extremely important yet challenging thermodynamic systems. Some examples of current research topics which basically encompass studies of molecular association would include polymer blending and protein dynamics. These topics may be more modern than ethanol + water distillation, but the problem is essentially the same: How do the strong attractive interactions of association compete with the strong core repulsions of nearby nonassociating atoms? Several promising methods of addressing this problem have been advanced recently. These methods range from molecular simulations to cluster diagrammatic expansions to classical chemical reaction theories. A number of recent papers have treated engineering aspects of cluster expansion theories for the thermodynamics of association (cf. Chapman et al., 1990; Panayiotou and Sanchez, 1991; Suresh and Elliott et. al., 1990, 1992), and these investigations have been very fruitful in characterizing complex solutions via fundamental physical relations. Despite the progress of these theories, however, there are a number of questions which only molecular simulations can answer. For instance, testing whether one theory for the effect of steric hindrance is more accurate than another requires detailed information about associations relative to position on the chain that experimental methods like spectroscopy cannot provide. To probe the details of evolution of association networks and their dynamic equilibria, there is no substitute for molecular simulation. General methods of molecular simulation have advanced rapidly in recent years, and it is now possible to purchase sophisticated software to probe the molecular details of complex phenomena. Molecular simulation capability has been credited substantially in the development of new pharmaceuticals, agrochemicals, materials, catalysts, and polymers. For simulating mixtures of large numbers of molecules, the algorithms applied are generally similar to the CHARMM technique (McCammon et al., 1977). The approach adapted in CHARMM is to assign point charges at various locations throughout the molecule. With enough point charges, it is possible to represent both the longrange polar forces and the short-range association forces around the molecule. These forces characterize the + E-mail:
[email protected]. t E-mail:
[email protected].
0888-588519412633-0957$04.50/0
equations of motion, and in principle, all of the molecular dynamics and thermodynamics can be solved via relatively simple integration algorithms. There are two problems with this approach, however, both relating to the simple inverse distance dependency of the potential energy of Coulombic point charges. First, hydrogen-bonding interactions are short-range and highly specific relative to orientation. Using long-range potentials to represent short-range interactions causes a need for a large number of point charges. Since the computation time of such an algorithm increases as W ,adding charge sites leads to a significant practical penalty. Generally, a compromise is reached among accuracy for the short-range properties, accuracy for the long-range properties, and computational tractability. Second, the range of a Coulombic potential is so long that it tends to be significant even at boundaries where it should be zero. For the simulation to proceed properly, this necessitates implementation of a corrective measure such as an Ewald sum. As a rule of thumb, applying an Ewald sum causes a simulation to run about 10 times slower. This additional slowness contributes another significant practical penalty when this common approach is applied. An alternative method of molecular simulation for associating systems has been pursued by Gubbins and coworkers (cf. Jackson et al., 1988) and Nezbeda and coworkers (cf. Kolafa and Nezbeda, 1987,1991). They have represented potentials of associatingmoleculesin a manner originally suggested by Andersen (1973) for focusing on the short-range nature of association. In this approach, the part of the potential which leads to association is represented by a small, short-range, attractive site located off-center from a relatively large repulsive site (Figure 1). These kinds of potentials form the basis for the cluster expansion theories developed by Wertheim (1984). Wertheim’s motivation for studying this kind of potential was that a number of simplifications could be applied in developing an analytic theory for the thermodynamics arising from these potentials. In fact, this theory has been the source of many of the recent advances in analytic theories. These potentials also have the advantage of concentrating directly on the short-range nature of association. A single site can represent the essential information in many cases. And no Ewald summation is required because the potential is short-range. But there is one significant problem with this approach. Because the potential function is discontinuous, the algorithm for simulating the dynamics of these molecules is entirely 0 1994 American Chemical Society
958 Ind. Eng. Chem. Res., Vol. 33, No. 4, 1994
A
B
C
Figure 1. Forms of potential models of simulating short-range association: (A) J o s h etal. (1987), r, = 0.550, L = 0.258;(B) Kolafa et al. (1987), r. = 0.150, L = 0.45~;(C) Jackson et al. (1988), r, = 0.575s,0 = i27.5O. The shaded areas indicate attractive sites. The open arean indicate repulsive sites. rc = radius of attractive site. L = distance between the repulsive site and attractive site in the same molecule.
different from the algorithm for continuous potentials. In fact, all simulations for these potentials to date have analyzed only the equilibrium properties via Monte Carlo simulation. Monte Carlo simulation is based on energy changes resulting from random displacements of molecules. Thus, the discontinuous potentials have little impact on the algorithm for Monte Carlo simulation, hut the dynamic information that might he helpful in understanding polymer diffusion or protein folding is lost. In the work presented here, we show how the benefits of short-range association potentials can he incorporated into a dynamic simulation through development of a specialized algorithm. The resulting algorithm offers a very large computational advantage over conventional dynamic simulation algorithms. The essence of this algorithm is that it takes advantage of the fact that discontinuouspotentials remain constant except when two specific sites collide, and these collision events may he scheduled in a manner which takes full advantage of computer data structures. The result is that simulation of the dynamics of large systems with long relaxation times should now he considerably more practical than previously envisioned. The implementation of a dynamic simulation for bonding systems runs into several difficulties that are not encountered in Monte Carlosimulation. These difficulties arise because Monte Carlo simulation is free to address only changes in configurational energy in random fashion whereas dynamic simulation must also address changes in kinetic energy as constrained by Newton's laws of motion. The objective of this paper is to address these difficulties systematicallyand show how they can he resolved to obtain an efficient and accurate simulation method. The presentation is organized as follows. The immediately following section reviews the background on dynamic simulations of discontinuous potentials leading to Rapaport's (1980) algorithm. The next section describes the difficultiesthat arise in extending the algorithm to bonding molecules and presents an empirical analysis of how these difficulties may he addressed. The last two sections provide illustrations of how the algorithm can he applied. All the simulations are conducted on microcanonical ensembles (i.e., NVE ensembles). These are intended primarily for illustrative purposes whereas future reports will focus on generating comprehensive results for large numbers of systems. Nevertheless, the sections on equilibrium results present new data for a model of water as well as a comparison to results from Monte Carlo simulations at the same conditions. The final section presents results for dynamic simulations that cannot he obtained in any other way, hut they are very limited in scope. Through the illustrative information about kinetics, transport, and thermodynamics that can he obtained from this method, it is hoped that a generalengineeringaudience can appreciate the benefit of this technique.
Y
-
O
100
200 300 400 Site Number
500
600
Figure 2. CPU time M site numbers on IBM RISC/6000-320H for hard sphere to 1ns of molecular time.
Previous Dynamic Simulations with Discontinuous Potentials Molecular simulations with discontinuous potentials are not new. As a matter of fact, the earliest dynamic simnlations were simulations of hard spheres, the simplest discontinuous potential (Alder and Wainright, 1959,1960). Thepurposeofthissectionandthenextistohrieflyoutline the technology to date so that the context of the extension presented here can he understood. Those who are familiar with this family of simulation methods may wish to move directly to the section on bonding molecules. A number of other investigators worked on simulations of hard spheres, with the algorithm eventually evolving to the status described by Wood (1975). This approach to simulation took on a new life, however, when Rapaport (1979) applied the method to simulating chain molecules. The key to Rapaport's success was to allow the bonds between sites on the same chain to vibrate within a square well. This made it possible to apply an algorithmvirtually identical to that of Alder and Wainright to the simulation of polymer-like chains. This approach has been extended by Chapela and co-workers (1984), and by Denlinger and Hall (1990). One important development that has rarely been cited, however, is Rapaport's development of a highly efficient algorithm for these kinds of simulations (Rapaport, 1980). Denlinger and Hall, for instance, applied an algorithm that was essentially the same as that of Alder and Wainright. A commonly availahleversionof the Alder and Wainright algorithm appears in the standard text on molecular Simulations by Allen and Tildesley (1987). To emphasize the significanceof Rapaport's algorithm (DMD for discontinuous MD), we have prepared Figure 2 as a direct comparison of Rapaport's algorithm (DMD) to that of Allen and Tildesley (A&T). Consideration of the nature of the algorithms indicates that the computation time required to simulate one nanosecond of the real time for the A&T algorithm increases as tCpU= 1.5
+ 0 . 0 0 3 p (min)
(1)
whereas Rapaport's algorithm (DMD) varies as
As simulations of novel systems like interfaces and large molecules require large numbers of sites, it is clear that this advantage is very significant. As another example,
1 P
0
500
1000
1500
I
4
a
!
*
I
I
I
0.00
0.01
0.02
0.03
2000
Site Number
Figure 3. CPU time vs site numbers on IBM RISC/6000-320H for butane to 0.02 ns of molecular time. For SHAKE algorithm, the time step of 2 fs is used.
we illustrate the computational efficiency of Rapaport’s algorithm (DMD) for chain molecules vs the SHAKE algorithm, which is a common technique for simulations of chain molecules with continuous potentials (Figure 3). The potential models for butane were constrained to bond lengths of L = 0 . 4 ~and bond angles of l l O o in both algorithms. On the basis of these figures it is clear that Rapaport’s (1980)algorithm (DMD)offers significant advantages that should be pursued. The problem is that no simulations have been performed by this method that could provide some insight into the complex equilibria and dynamics of bonding systems. In the presentation below, we discuss how this extension can be made and some of the difficulties that are encountered. We illustrate one methodology by which these difficulties may be resolved and develop a complete algorithm which can be generally applied to provide efficient insights into many engineering problems.
Fundamentals of the DMD Algorithm To understand how such improvements in efficiency are possible, we must review the methods of simulating discontinuous potentials vs the methods for continuous potentials. In simulating continuous potentials, a fixed time step is generally used for all the atoms in the simulation. At the beginning of the time step, the forces acting on each site must be summed over all the sites in the simulation, a summation of length N. Since this summation must be performed for all N sites, the total time spent in this summation is of order W . When sites are bound together as in chain molecules, one must apply a method of constraining the bond lengths and bond angles. This may be accomplished either by incorporating steep harmonic potentials between each bonded pair or by applying some constrained dynamics technique, such as the SHAKE method or the method of Edberg et al. (1986). In either case, the computation time required to simulate bonded sites increases by about an order of magnitude relative to that for unbonded sites. Now consider simulation via discontinuous potentials. This approach basically amounts to searching for future collisions and scheduling them. The forces between the atoms not colliding do not change because the potential is flat there. Those atoms simply get closer together or farther apart. The colliding atoms exchange momentum, and their future prospective collisions must be computed
I 0.04
Figure 4. Time steps vs reciprocal of site numbers for hard sphere molecules. Curve fit: At = 151(1/N) ( N = site number).
and scheduled. The obvious thing to do at the time of the collision is to update the positions of all the interaction sites in accordance with the time since the last collision, but that would require a loop of order N when only the coordinates of two sites have changed in any significant way. One trick then is to compute false positions where the two sites just colliding would have been at some previous big update of all positions such that their velocities and positions are correct after the collision. Then all future events will be correct, and that is the primary concern. Another trick is to break the simulation box into many small cells, each about the size of the largest repulsive site. Then cell crossings are tracked in much the same way as collisions and the only possible collisions, are with sites in adjacent cells. Therefore, the search for future collisions need not proceed over all sites in the fluid. A similar cell structure may be applied to simulations of continuous potentials, and may reduce their computation times to order N. But, the range of a potential like the Lennard-Jones potential is much larger than that of a hard sphere potential so that the cells must also be very large. Generally, the sites included in summing the forces in a continuous potential model must be 1-2 orders of magnitude more than in a discontinuous potential model, with proportional penalties in computation time. As an example, in simulating liquid butane at 0.52 g/cm3,roughly 500 sites must be summed in the force calculation via a Lennard-Jones potential whereas only 30 sites must be checked for future collisions via a hard chain potential. The advantages of false positions and cell structures were known before Rapaport presented his algorithm (cf. Erpenpeck and Wood, 1977). What Rapaport contributed was a sophisticated’algorithm for scheduling collisions. In previous algorithms, identifying the next event required a search of roughly order I W 2 . For example, the A&T algorithm performs a linear search over all N sites to identify the next event. Rapaport developed an algorithm for keeping track of events that implements the data structure of a binary tree (cf. Aho et a1.,1974). Tracking events with a binary tree requires a computational effort of order logz(N). This means that all computations per collision require a time of order logp(N). The only strong dependency on N for the Rapaport algorithm is then the fact that the time between collisions is proportional to N-1, as shown in Figure 4. Therefore, simulating a fixed amount of time, like a nanosecond, requires more collisions for a large system than for a small system.
960 Ind. Eng. Chem. Res., Vol. 33, No. 4,1994
Table 1. DMD-B Simulation Results vs J o s h ' s Result at the Same Conditions. Josh's results DMD results with r. = 0.550, with r. = 0.550,
L = 0.25 0.005s 0.5 0.2454
L =0.25~ 0.5 DMD-4s rite model
DMD-1B site model =
Donor or Acceptor
0
= Donor;
(11
-
Acceptor
I21
Figure 5. (1) DMD-1B site model. CT is the diameter of the hard sphere home site,L is the distance between the center ofhard sphere home site and the center of bondable site (donor or acceptor); r. is the diameter of bondable site. (2) DMD-4B site model. c is the diameter of the hard sphere home site; L is the distance between the center of hard sphere home site and the center of acceptor site, the distance between the center of the hard sphere home site and the center of donor is 012; r, is the diameter of donor or acceptor.
Applying DMD to Bonding Molecules In principle, the extension of the DMD algorithm to associatingmolecules is straightforward. The square-well association sites must simply be added to the model and the equations of motion integrated. The square-well associations are different from the hard sphere repulsions in that they may "bond", when two sites are attracted to each other, or "dissociate", when two square-wells pull away from each other, or "bounce", when the association sites try to pull away from each other hut their momenta are not high enough to overcome the association energy. But incorporating these considerations is straightforward. What islessstraightforwardis exactlyhow theassociation sites should be modeled. That is, what values should be assigned to their distance from the center of the repulsive site? What size should be assigned to the sites? How do these considerations impact the dynamics of hydrogen bonding which may be known from other sources? And, a very tricky question, what mass should be assigned to the proton acceptor association sites? One might expect that the mass of the acceptor site should he like that of an electron pair, which is practically negligible, but then the acceptor site could not influence the motions of the atoms to which it is attached. To answer these questions, we have undertaken several comparative studies of the various possibilities for each question. The answers to these questions lead to the development of a specific algorithm and model for the molecular interactions which we refer to as DMD for bonding (DMD-B). The 1B site and 4B site models are shown in Figure 5. Effects of Site Size a n d Distance from the Center of the Repulsive Site. As the bonding sites become larger, they tend to be more disperse in the way they represent the association energy. This trend suggests several questions regarding the best representation of the association sites. As a result of Wertheim's firsborder theory, discussed by Jackson et al. (1988) in this context, the precise shape of the bonding sites is relatively unimportant to theequilibrium resultsas long as theenergy of bonding and the volume available for bonding are maintained. These quantities are characterized by the depth of the association well 6 and the quantity Kab defined 2
Ksb = %[In(
rc + 2L 7 ) ( 6 r :
+ 18r;L
- 24Ls)
+ (r, +
v
0.2454
(&b)
2.69 19 0.4429
elk (K) Z = PIpkT hydrogen bond number X
2100 2.62 0.031 20 1.95
2700
10'
0 - W bond angle (deg) temperature (K)
* 0.31 * 0.0199 * *
180 18.08 308 16.2 0.108
298
time step (fs)
X A = mole fraction of acceptor sites; v = ( d 6 ) ~ 9 9Z; = compressibility; L = hydrogen bond energy. Table 2. Effects of Potential Well-Width on Simulation Results in DMD-B Simulation of 1B Site Model. potential well-width 0.10s 0.05~ 0.01.7
*
av lifetime of hydrogen
(K.b) X 10'
0.52 0.071 0.31 0.026 0.23 0.0052 188.9 135.9 43.97
hydrogen bond "umber temperature (K)
19 5.21 299*43.2
*
284:
X A= 0.5;r, = 0.15~;the center of well = 0.4750; = 0.2454; elk = 2700 K molecular time I 0.6 na. where o is the diameter of hard sphere site, r, is the diameter of bondable site, and L is the distance between the center of a hard sphere site and the center of a bondable site. It should be noted that the vibrating bonds of the
DMD-Balgorithmcausethevalue0fK.btochangethrough the changes in bond length. This means that the value of K.bneeds to he computed for several cases: (1)theaverage (&,) over all sites, (2) the average (&) over bonded sites, and (3) the average (&b) for unhonded sites. By considering all three cases, we can gain some insight into the similaritiesand differencesbetween the various models. Table 1showsthat theaverage (K&) in DMD-B simulation may he slightly different from that in Monte Carlo simulation even though L and re are the same in both simulations. That is because that the formation of a hydrogen bond stretches the bond length between the hard sphere site and the donor or acceptor in DMD-B. The simulation shows that the Kabvariation ranges from 0.24 X lo-'$ to 0.73 X 10-W for the large site model (rc = 0.550) at the well-width of 0.010 and from 0.171 X l(r8 to 0.29 X 10-W for the small site model (rc= 0.150) a t the same well-width (Table 2). Two inferences may be drawn from these observations. First, the large bonding site model is much more sensitive to the variability introduced by bond vibrations. It seems that for the DMD-B model in which the bond vibrations are involved the equilibrium results are affected to some extent by the precise shape of the bonding sites. The small site model is more appropriate for the DMD-B algorithm because the vibrations have a smaller impact on the bonding volume. Second, by a comparison of the extent of bonding and compressibility in all these models, it seems that the equilibrium properties are not all that sensitive to vibrations in the value of (Kab). So some variability in this parameter is permissible. The formula of &b for the 4B site model is slightly different since the bond length between proton donor, or acceptor, and the big home site is different from that for the 1B site:
Ind. Eng. Chem. Res., Vol. 33, No. 4, 1994 961
where L is the distance between the center of the hard sphere site and the center of acceptor site. The meanings of rc and u are as same as those in eq 3. The angle between hydrogen-bonding sites is known to be near 180'. As shown in Table 1,the angle varies over a range of 160-200' with a mean of 180'. One possible reason why the value varies so widely might be that the high density provides a driving force away from linear bonding in favor of better packing. But when we used different densities to execute the 4B site model simulation discussed below, the results did not show a significant difference in bond angles. This implies that this type of potential model is such that the 180' angle is not strictly enforced as long as the bonding sites can overlap to some degree. Effect of Well-Widthbetween Covalently Bonded Sites. The well-width between covalently bonded sites is analogous to the steepness of the well between bonded sites in a continuous potential simulation. By this well being incorporated, the bonded sites are constrained to remain within a certain distance from each other, this distance approximating the bond length. Allowing this bond distance to vibrate within its well instead of being constrained to an exact value greatly simplifies the algorithm. But the impact of this looseness should be carefully considered when the potential model is designed and the computation speed is optimized. When the potential is designed, it is important that the loose constraint on the bond length not have a significantly negative impact on the nature of the bonding volume. This impact can be straightforwardly minimized by maintaining a very narrow well-width. When the computation speed is optimized, however, it is important to provide as wide a well as feasible, because a narrow well leads to a large number of tedious intramolecular collisions and a short time step. The examples considered below provide some insight into this mini-max problem. In the first example we consider a single association site model (Figure 5). Table 2 shows the effect of varying the well-width from 0 . 1 ~ to 0 . 0 1 ~on the thermodynamic properties, the time step, and the hydrogen bond lifetime while the smaller site model is adopted. Note that the wider well-width leads to a significantly longer time step and longer lifetime of the hydrogen bond in this 1B site model. The longer bonding lifetime for the looser bond is consistent with the availability of more degrees of freedom to spread the bonding energy, as developed in the theory of unimolecular reactions. The sensitivity to the well-width means that the well width must be carefully studied in building models of reaction kinetics. We have applied the well-width of 0.01~ in most calculations below, and we find reasonable agreement with the available spectroscopic data for water. Progressing to the narrower well-widths shows that the standard deviations of the hydrogen-bond number in the simulation box as well as the thermodynamic properties are significantly reduced when the well-width decreases to 0.01~.The variation in (Kab) suggeststhat the bond lengths of the small site model tend to be stretched when allowed to vibrate. This trend is confirmed in Table 3, which shows that unbonded donors and acceptors tend to stretch their bonds as well as the
hydrogen-bonded species. From Table 3 it is clear that the variability in the bond length is very small when the well-width is 0.01~. Similar considerations can be applied to the kind of 4B site model (Figure 5) that one may apply to water. Table 4 shows the variations in properties for the water model relative to variations of the well-width. Once again the standard deviations are broader when the well is wider. One important note is that the extent of association (as reflected by fractions of monomer and dimer) increases with the increase of the well-width for the 4B site model. It should be noted that when the well-width is 0.01u, the result for the oligomer distribution matches the Monte Carlo results of Kolafa and Nezbeda (1987) for similar geometric and energy parameters. The oligomer distribution may provide a more sensitive indication of the optimal well-width. Another important phenomenon is that the well-width weakly affects the average lifetime of the hydrogen bond for the 4B site model, as indicated in Table 4. This is different from that of the 1B site model. In the 1B site model, there is only the vibration between bondable site and hard sphere home site so that the wider well-width results in lower vibrating frequency and longer hydrogenbond life, but in the 4B site model, the situation is different. There are not only vibrations between bondable sites and the hard sphere site but also intramolecular collisions between bondable sites. So relaxing the well-width does not affect the hydrogen-bond lifetime as much as in 1B site model. The final result for the hydrogen-bonding lifetime is something that should be measurable spectroscopically. The closest measurements we have encountered are measurements of rotational time correlations by Klier (1973) and Kerstel et al. (1992). Klier's measurements indicate a lifetime of about 1ps, but it is not clear that they are measuring exactly what we call a bonding lifetime. I t is reassuring that these values are on the same order of magnitude, however. Estimations based on treating the bond energy as an excitation energy also lead to similar values. Effect of Acceptor Mass. One problem in a dynamic simulation that does not arise in Monte Carlo simulations is the assignment of a value for the mass of each site in the simulation model. This does not present a problem in conventional MD simulations because the mass of each atom provides an obvious choice. When the association sites are broken down into proton donors and proton acceptors, however, it is not entirely obvious what mass should be assigned to the proton acceptors. These acceptor sites represent regions of enhanced electron density which can draw donors in to form a hydrogen bond. A Monte Carlo simulation assigns only energy to these associations and randomly varies the positions of all the sites in search of the overall minimum in free energy. Therefore, the question of mass never arises. In a dynamic simulation, however, the equation for change in site velocity for an attractive collision is scaled by the reciprocal mass of the site, so a mass of lo4 atomic mass units (amu) for an acceptor would lead to extreme elevations of the velocities of the acceptor sites. These higher velocities would then lead to a large number of short time step collisions and extremely tedious calculations of collisions that have very little impact on the progress of the simulation. Table 5 shows some results with different masses of acceptor. There are some changes in the simulation results, especially for collision time step and the lifetime of the hydrogen bond, when the mass of acceptor changes from 1.0 to 0.1 amu. The trends here suggest that reducing the mass to 0.1 amu
962 Ind. Eng. Chem. Res., Vol. 33, No. 4, 1994 Table 3. Bond Length Comparison (After 1 ns Molecular Time) to that in Initial Configuration in the 4B Site Model Which Can Be Used To Simulate the Behavior of Water.
potential well-width initial config 0.1OU
0.05~ 0.01u a9
av length of intramolecular O-D bond (not hydrogen bound) 0.500~ 0.519 f 0.0032 0.505 0.0006 0.5002 f 0.00003
= 0.26; c/k = 2014 K;
av length of intramolecular 0-D bond (hydrogen bound)
*
o
0.526 f 0.0013 0.507 f 0.0015 0.5004 f 0.00004
= diameter of oxygen.
Table 4. DMD-B Simulation Results of 4B Site Model of Water at Different Potential Well-Widths (Molecular Time = 1 ns, Average Time = 0.5 ns). 0.1OO 0.050 0.01O well-width 0.0038 time step (fs) 0.017 0.012 4.47 5.26 av lifetime of hydrogen 6.03 bound (ps) 2.22 Z = P/pkT 2.13 2.22 332.6 f 45.2 315.4 f 43.5 299.0 f 8.0 temperature (K) 180 f 25.1 O-D-0 bond angle (deg) 180 f 29.9 180 f 28.1 distribution of I-mer 0.245 0.280 0.280 10.124 0.160 0.167 20.119 0.082 0.102 30.060 0.070 40.089 0.059 0.045 0.073 50.042 0.045 0.046 60.297 270.230 0.385 0 Initial temperature = 158.5 K; 7 = 0.26; elk = 2014 K; r, = 0.150; acceptor mass = 1.0 amu.
Table 5. DMD-B Simulation Results of 4B Site Model with Different Acceptor Masses (Well-Width = 0 . 0 1 ~ ) ~ mass of acceptor (amu) 1.0 0.5 0.1 0.0038 0.0036 0.0018 time step (fs) av hydrogen bond 6.03 5.60 2.46 lifetime (pa) 2.23 f 0.11 2.27 f 0.09 Z = PIpkT 2.22 f 0.07 299.0 f 8.0 299.5 f 10.5 302.7 f 9.4 temperature (K) 180 f 25.0 180 24.8 O-D-0 bond angle (dee) 180 f 25.1 distribution of I-her -. 10.280 f 0.030 0.273 f 0.030 0.287 f 0.044 20.163 0.167 0.161 0.119 0.126 0.091 30.089 0.080 0.077 40.073 0.082 0.056 50.042 0.041 0.039 6270.230 0.237 0.287
*
av length of intramolecular O-A bond (not hydrogen bound) 0.4500 0.470 f 0.025 0.457 f 0.0013 0.4503 f 0.00003
av length of intramolecular O-A bond (hydrogen bound)
-. 120
1
100
1
0.477 f 0.0007 0.461 f 0.0010 0.4504 O.ooOo2
I
-
2
I'
-1
80
I
---
I
1g/crni
0t -20 ' 20
1
40
,
,
I
80 100 120 Number of Molecules 50
i 1
140
Figure 6. Comparisons of CHARMM (version 3.0) vs DMD-B on Silicon Graphics 4D/240 for 1 ps simulated time. CHARMM calculations are based on ST2 model of water. DMD-B calculations are based on the 4B site model of water.
parameterization for comparison to the DMD-B method. It should be noted that Ewald sums were included in the CHARMM model because of the long-range nature of the Coulombic forces. Figure 6 shows that the DMD-B algorithm retains its benefits in computational efficiency even for hydrogen-bonding systems. In closing,we should note that Ding and Goddard (1992) have presented an algorithm which also claims W dependency. We would remark, however, that the W dependency is only apparent for systems larger than 15 000 atoms.
a Initial temperature = 158 K; t-k = 2014; 9 = 0.26; r, = 0.150; molecular time = 1.0 ns; average time = 0.1 ns.
Analysis of Association Trends over a Range of Bond Energies and Densities
leads to rapid vibrations that decrease the bonding lifetime in a way that may be misleading. Many of the bonds broken by the rapid vibrations are quickly reformed such that the discontinuous potential is not accurately portraying bond formation. On the basis of these results we conclude that an acceptor mass of 1.0 amu represents a reasonable value and, being equal to the mass of proton, does simplify the data structure of initial configuration. Comparison of the DMD-B Algorithm to Conventional Methods. Software for hydrogen-bonding polyatomic molecules is rapidly becoming a standard tool in the pharmaceutical and specialty chemical industries, and the various methods are very similar in their algorithms since all are based on continuous potentials with point charges or multipoles to mimic the hydrogen-bonding effect. The principal source of debate centers on the parameterization of the locations and strengths of the Coulombic contributions. It is unlikely that any parameterization will accurately mimic all of the results of bonding and polarity at any time in the near future. Therefore, we have somewhat arbitrarily selected one representative form of software based on the CHARMM
Having developed an efficient tool for investigating association effects, we would now like to apply it to illustrate what insight may be gained relevant to predictions that can be made about dynamic and thermodynamic properties over a range of conditions. Throughout these computations, we neglect the disperse attraction energy that could be modeled by the usual centered square-well potential with a range of 1 . 5 ~ These . disperse attractions can be treated by the methods presented here, but the longer range of the potential does lead to some reduction in computational efficiency. Since the roles of disperse attractions have been studied previously and are already understood to some degree, it is more interesting to focus specifically on the roles of association in this analysis. By varying bond energies and densities for modeled water, we can gain some insight into the role of association in altering the diffusion coefficients as well as the structure and thermodynamics. Table 6 shows the results for simulations of modeled water at two different densities and two different association energies. Note that the extent of association at low energy is very small at both values of density. There are barely any trimers forming
Ind. Eng. Chem. Res., Vol. 33, No. 4, 1994 963 Table 6. DMD-B Simulation Results of 4B Site Model at Different Densities and Different Hydrogen Bond Energies. time step (fa) av lifetime of hydrogen bond (pa) 2 = P/pkT temperature (K) 0-D-0 bond angle (deg) fraction of monomer computed by Wertheim theory distribution of I-mer 123456617(I
7 0.130; tlk 0.003 88 6.3 1.271 f 0.101 300.5 f 9.2 180 24.38 0.5876
2014 K
0.544 f 0.035 0.226 0.110 0.071 0.035 0.035 0.010 0.004
7 = 0.260; c/k
7 = 0.130; elk = 1007 K 0.003 91 1.5 1.525 f 0.036 302.1 1.4 180 f 15.78 0.9788
= 2014 K
7 0.260; c/k = 1007 K 0.003 82 2.3 2.360 f 0.084 300.7 f 2.6 180 f 23.18 0.9391
0.003 77 6.0 2.226 f 0.074 299.0 f 8.0 180 25.08 0.2908
*
*
0.980 i 0.020 0.019 0.001
*
0.280 f 0.030 0.167 0.119 0.089 0.073 0.073 0.042 0.230
0.922 0.033 0.068 0.009 0.001
The well-width is 0 . 0 1 ~mass ; of acceptor is 1 amu; rc = 0 . 1 5 ~molecular ; time = 1 ns; average time = 0.15 ns.
at low density and low energy. The lifetime of the hydrogen bond is also much shorter at low energy than at higher energy, as one would expect. Table 6 also indicates that the lifetime of a hydrogen bond is a strong function of bonding energy but is a weak function of density. I t seems that the hydrogen-bond clusters appear to be bigger if the density is higher, which is reflected by the monomer fractions and the oligomer distributions in Table 6. The monomer fractions in all these cases can match the data calculated by Wertheim’s theory within one standard deviation, as shown in Table 6. It is noted that the narrower the deviation of the temperature, the less difference between the theoretical results and the simulation results on fraction of monomer.
1.oq
I
I
I
I
0 - 7)=0.065 ’1=0065
- 7)=0.130 7)=0 130 - ~7 = 0 .195 195 - 0=0260 0=0.260
0.8 A
s >
N-
0.6
=: h
4 v
? 0.4
h
0
v
>
0.2
0.0
Analysis of Effect of Hydrogen Bonds on Self-Diffusivity One of the big advantages of MD vs MC simulation is that transport properties may also be derived. Here we investigate the velocity autocorrelation and self-diffusivity of hard sphere and 4B site models by use of the DMD-B algorithm. We have used the autocorrelation approach to the diffusivity instead of the displacement for two reasons. First, it is much easier to calculate since the velocities do not need to be corrected for false positioning. Second, our simulation times are much longer than typically practiced in 1970, so the difficulty of treating the hydrodynamic tail is minimized. The results are shown in Figures 7 and 8. For hard spheres, the diffusion coefficients match those of Alder et al. (1970),thus substantiating the validity of the algorithm. Similar to that in the hard sphere fluid, the lower the density the longer the time of velocity correlation before the uncorrelated line is reached. As shown in Figure 8, hydrogen bonding reduces the selfdiffusivity of the fluid as expected. This is because the hydrogen-bondingchains and clusters restrict the diffusion of the molecules. Figure 8 shows, however, that the hydrogen bonding affects the diffusivity more at low density than at high density.
Conclusions We have developed and demonstrated a molecular simulation algorithm for discontinuous molecular dynamics of bonding systems (DMD-B). The advantage of a molecular dynamics algorithm over the previously available Monte Carlo algorithm is that dynamic properties like bonding lifetimes, reaction rates, and diffusivity may be addressed. The disadvantage is that the fundamental issues involving the impact of bond vibrations and the
00
0.3
15
0.6 09 1.2 Correlation Time ( p s )
Figure 7. Velocity correlation of 4B site model of water with hydrogen bonds at 307 K and different densities.
1
005
I
I
I
1
010
015
020
025
‘I
Figure 8. Comparison of diffusion coefficientswith (1) hard sphere model and (2) 4B site model of water with HB. 7 = (46)pZ; p is molecular number density.
mass of the proton acceptor site must be resolved. We have undertaken a demonstrative analysis which shows that accurate representations can be obtained with small sites of re = 0 . 1 5 and ~ bond lengths of 0.475-0.5u,while proton acceptor masses of 1.0 amu are maintained if vibration well-widths of 0.01~are maintained. The simulation results for equilibrium properties agree, within stated tolerances, with previous MC results regarding extents of association and distributions of oligomer
964 Ind. Eng. Chem. Res., Vol. 33,No. 4,1994
size.With regard to dynamic properties, the simulations indicate a hydrogen-bonding lifetime of 6 ps (Table 4) for simulated water at 9 = 0.26 and 300 K, which is similar in magnitude to the value reported by Klier (1973). The role of hydrogen bonding in altering diffusivity is to decrease the diffusivity more strongly at low densities. Through these demonstrations we believe we have substantiated the validity of this method of molecular simulation and the kinds of novel insights that it can provide. It may not be as efficient as some biased MC methods for obtaining the equilibrium properties, but it provides an extremely efficient method of simulating the dynamic properties.
Acknowledgment This material is based upon work supported by the National Science Foundation under Grant No. CTS 9110285. The government has certain rights in this material.
Literature Cited Aho, A. V.; Hopcroft, J. E.; Ullman, J. D. The Design and Analysis of Computer Algorithms;Addison-WesleyPublishing Co: Reading, MA, 1974;p 115. Alder, B. J.; Wainwright, T. E. Studies in Molecular Dynamics. I. General Method. J. Chem. Phys. 1969,31,459. Alder, B. J.; Wainwright, T. E. Studies in Molecular Dynamics. 11. Behavior of a Small Number of Elastic Spheres. J. Chem. Phys. 1960,33,1439. Alder, B. J.; Gass, D. M.; Wainwright, T. E. Studies in Molecular Dynamics. VIII. The Transport Coefficients for a Hard-Sphere Fluid. J. Chem. Phys. 1970,53,3813. Allen, P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford Science Publication: Oxford, U.K., 1987;Chapters 3 and 5. Andersen, H. C. Cluster Expansions for Hydrogen-bonded fluids. I. Molecular Association in Dilute Gases. J. Chem. Phys. 1973,59, 4714. Chapela, G. A,; Martinez-Casas, S. E.; Alejandre, J. Molecular Dynamics for Discontinuous Potentials. I. General Method and Simulation of Hard Polyatomic Molecules. Mol. Phys. 1984,53, 139. Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radoaz, M. New Reference Equation of State for Associating Liquids. Znd. Eng. Chem. Res. 1990,29,1709. Denlinger, M. A.; Hall, C. K. Molecular-dynamics SimulationResults for the Pressure of Hard-chain Fluids. Mol. Phys. 1990,71,No. 3, 541.
Ding, H. Q.; Karasawa, N.; Goddard, W. A,, 111. Atomic Level Simulations on a Million Particles: The Cell Multipole Method for Coulomb and London Nonbond Interactions. J. Chem. Phys. 1992,97,4309. Edberg, R.;Evans, D. J.; Morriss, G. P. Constrained Molecular Dynamics: Simulations of Liquid Alkanes with a New Algorithm. J. Chem. Phys. 1986,84,6933. Elliott, J. R.,Jr.; Suresh, S. J.; Donohue, M. D. A Simple Equation of State for Nonspherical and Associating Molecules. Znd. Eng. Chem. Res. 1990,29,1476. Erpenpeck, J.J.; Wood, W. W. Modern Theoretical Chemistry;Berne, B. J., Ed.; Plenum: New York, 1977;Vol. 6B,p 1. Jackson, G.; Chapman, W. G.; Gubbins, K. E. Phase Equilibria of Associating Fluids Spherical Molecules with Multiple Bonding Sites. Mol. Phys. 1988,65,1. Joslin, C. G.; Gray, C. G.; Chapman, W. G.; Gubbins, K. E. Theory and Simulation of Associating Liquid Mixtures. 11. Mol. Phys. 1987,62,843. Kerstel, E. R. T.; Meyer, H.; Lehmann, K. K.; Scolee, G. The Rotationally Resolved 1.5Hm Spectrum of the HCN-HF Hydrogenbonded Complex. J. Chem. Phys. 1992,97,8896. Klier, K. Potential Correlation Motion in Liquid and Adsorbed Water. J. Chem. Phys. 1973,58,737. Kolafa, J.; Nezbeda, 1.Monte Carlo simulations on Primitive Models of Water and Methanol. Mol. Phys. 1987,61,161. Kolafa, J.; Nezbeda, I. Primitive Models of Associated Liquids. Equation of State, Liquid-gas Phase Transition, and Percolation Threshold. Mol. Phys. 1991,72,777. McCammon, J. A.; Gelin, B. R.;Karplus, M. Dynamics of Folded Proteins. Nature 1977,267,585. Panayiotou, C.; Sanchez, I. C. Hydrogen Bonding in Fluids: An Equation-of-state Approach. J. Phys. Chem. 1991,95,lOO90. Rapaport, D. C. Molecular Dynamics Study of a Polymer Chain in Solution. J. Chem. Phys. 1979,71,3299. Rapaport, D. C. The Event Scheduling Problem in Molecular Dynamics Simulation. J. Comput. Phys. 1980,34,184. Suresh, S. J.; Elliott, J. R.,Jr. Multiphaee Equilibrium Analysis via a Generalized Equation of State for Associating Mixtures. Znd. Eng. Chem. Res. 1992,31,2783. Wertheim, M. S. Fluids with Highly Directional Attractive Forces. I. Statistical Thermodynamics. J. Stat. Phys. 198$,36,19. Wood, W. W. Computer Studies on Fluid Systems of Hard-core Particles. Fundamental Problems in Statistical Mechanics ZZfi North-Holland/American Elsevier: Amsterdam, New York, 1975; p 331.
Received for review July 8, 1993 Revised manuscript received January 10,1994 Accepted January 24, 1994. Abstract published in Advance ACS Abstracts, March 1, 1994.
