Screening Effects on Hydrogen Bonding in Chain Molecular Fluids

Jul 3, 1996 - AIChE Journal 2010 , NA-NA. Article Options. PDF (320 KB) · PDF w/ Links (297 KB) · Full Text HTML · Abstract · Figures · References · C...
0 downloads 0 Views 320KB Size
Download PDF
Ind. Eng. Chem. Res. 1996, 35, 2369-2377

2369

GENERAL RESEARCH Screening Effects on Hydrogen Bonding in Chain Molecular Fluids: Thermodynamics and Kinetics Jin-Xing Liu† and J. Richard Elliott, Jr.*,† Chemical Engineering Department, The University of Akron, Akron, Ohio 44325-3906

Development of a special molecular dynamics algorithm (discontinuous molecular dynamics with bonding) was previously reported for simulating molecules via potentials with small squarewell potential bonding sites that mimic hydrogen bonding. The present study extends this approach to fluids composed of two new molecular models and focuses especially on the kinetics of bond associations and dissociations in addition to the equilibrium effects. The first molecular model comprises two repulsive sites and three bonding sites (simulated methanol). The second molecular model comprises five repulsive sites where bonding sites are placed at the end and the middle of the chain (a simulated 5-mer diamine). Equilibrium properties focus on the fractions of bonding sites that remain unbonded after extended equilibration. At moderate to high packing fractions, the equilibrium fraction of unbonded acceptors for the middle site of the 5-mer is roughly equal to that of the end site. The fluid structures around end sites and middle sites are analyzed in order to suggest how the equilibrium effects may be analytically predicted. The rates of bond association and dissociation are studied for each type of bonding site as a function of temperature, molecular packing fraction, position along the chain, and bonding energy. The kinetic data more strongly suggest the presence of a mild screening effect by showing a shift in the rate constant. At all densities, the bonding in the 5-mer is higher than expected, indicating a stabilization effect that was not previously anticipated. These data point the way toward a fundamental theory of the kinetics of hydrogen bonding. 1. Introduction The ultimate goal of synthetic chemistry is to produce a chemical structure that possesses precisely the properties desired by the engineer who designed it. The desired properties may range from improved compatibility for a polymer blend via copolymers to enhanced resistance to viral infection via modified proteins. The common need underlying all of this breadth is to understand the detailed interactions between segments of the various molecules well enough to manipulate their behavior to the engineer’s advantage. It will take time to generate the necessary information for a complete understanding; a series of systematic investigations will be required. The following paper contributes to this systematic procedure by investigating the competition between the attractive interactions of hydrogen bonding relative to the repulsive interactions of a chain molecule’s backbone. We begin by studying one relatively small molecule that might be modeled by a simpler and more traditional spherical approximation and by contrasting it to a slightly larger chain molecule which accentuates the need for a nonspherical development. One of the differences between chain molecules and small molecules is that intramolecular interactions play significant roles in dictating the behavior of the overall fluid. For hydrogen-bonding systems, an example of such an intramolecular interaction is the repulsive screening of prospective bonding sites by other sites on the same chain. By the term “screening” we refer to something more general than steric hindrance. When †

E-mail: [email protected]; dickelliott@ uakron.edu.

S0888-5885(95)00107-2 CCC: $12.00

screening occurs by design of a specific molecular arrangement, it is referred to as steric hindrance. Screening, on the other hand, refers generically to any manner of inhibiting intermolecular contacts. The purpose of the investigation described in this paper is to apply a newly developed method of molecular simulation (Liu et. al, 1994) to explore how these repulsive screening effects compete with bonding attractions in relation to the kinetics of bonding as well as the thermodynamics of bonding. Since significant components of rate phenomena may be attributed to screening effects, this study should provide insight that is fundamental to a broad range of kinetic issues. Our methodology derives from equilibrium studies of hydrogen bonding in chain molecular systems. Recent progress in the study of the equilibrium thermodynamics of hydrogen bonding has demonstrated the efficacy of modeling hydrogen bond formation by small squarewell attractive sites located near the surface of a large repulsive site (Figure 1a). Potential models like these form the basis of Wertheim’s (1984a,b; 1986a,b) theory, which has been successfully adapted to describe a broad range of complex thermodynamic phenomena in hydrogen-bonding solutions. Many molecules can be described by combining repulsive sites with bonding sites to mimic the geometry of the molecule of interest. The accuracy of Wertheim’s theory has been demonstrated through the developments of SAFT theory (Chapman et al., 1989, 1990) and ESD theory (Elliott et al., 1990; Puhala and Elliott, 1993). In contrast to the equilibrium studies, kinetic studies of hydrogen bonding remain relatively underdeveloped. During further developments related to the ESD work, © 1996 American Chemical Society

2370

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996

Figure 1. Illustrations of molecular models used to study bonding. The small sites represent bonding sites. The correspondence between shading and the functional group is given in (a) a simple repulsive sphere with three square-well bonding sites, one donor and two acceptors, which we refer to as “spheranol”; (b) simplified representation of methanol; (c) 5-mer with a bonding “head” and a repulsive “tail”; (d) current “periodic table” of DMD-B. Table 1. Parameters of the Simulated Models of Methanol and 5-Mer Chain diameter of home sites (σ) diameter of proton donor diameter of proton acceptor bond length between adjoining repulsive sites square-well width between repulsive sites bond length between home site and donor/acceptor square-well width between repulsive site and donor/acceptor intramolecular bond angle mass of home site (amu) mass of donor/acceptor site (amu)

0.3 nm 0.15σ 0.15σ 0.51σ 0.1σ 0.475σ 0.01σ 109°28′ 14.5 1.0

we recently extended methods of molecular dynamics simulation to be applicable to the type of hydrogenbonding model described above. The advantage of molecular dynamics simulation relative to the Monte Carlo methods that previously had been used (Jackson et al., 1988; Kolafa and Nezbeda, 1987) is the information gained about the rates of bonding and dissociation. The combination of kinetics and thermodynamics permits new insights to be gained for many hydrogenbonding molecules. In the following section we list and explain the molecular models that we use in the simulations. In section 3, we give a brief description of the new molecular dynamics technique that we employed for the simulations. Details of this method are available in a previous publication (Liu et al., 1994). The essential features of Wertheim’s theory are outlined in section 4 with a special emphasis on the role of fluid structure in the theory. The results of equilibrium effects and kinetic effects are described and discussed in section 5. In the conclusions, we summarize key findings with an emphasis on what is new and what we expect from future work. 2. Molecular Models We consider two specific molecular systems. Both are illustrated in Figure 1 and the parameters are listed in Table 1. Figure 1b shows a simple molecule that is representative of “methanol”. This model is slightly oversimplified in that the size of the oxygen atom is assumed to be identical to that of the methyl group and that the hydrogens on the methyl group have been lumped into the single repulsive sphere of the methyl group. These simplifications help to clarify the issues to be addressed below. In particular, they establish a “head” group that is capable of hydrogen bonding and

a “tail” group that is purely repulsive. Note that in methanol both of these sites occur at the ends of the molecule since the molecule is so small. Note also that the screening which can occur on end atoms is very limited. To study the role of screening vs bonding, we must have a longer molecule, and we must study how bonding with sites in the middle of the molecule contrasts with bonding at the ends. This purpose is served by the somewhat peculiar molecule illustrated in Figure 1c. This molecule has a head and a tail but also includes bonding sites in the middle of the 5-mer chain. We refer to this molecular model simply as “5-mer” although it could be considered as a representation of (1-ethyl, 1-methylamine) amine. Note that the acceptor site in the middle of the molecule is identical to the acceptor on the end except for the extent to which it is screened by neighboring repulsive sites. Thus, comparing bonding between these two sites provides a measure of the screening effect. The tail group is identical to the head group except for the absence of bonding sites. Thus the molecular distributions around the tail provide a measure of the distributions to be expected when the bonding energy approaches zero. A typical hydrogenbonding energy of AD/k ) 2013.5 K (equivalent to 4 kcal/ mol) was used for the simulations in both tables. This value for the energy of hydrogen bonding is similar to that adopted by Sindzingre and Klein in the molecular dynamics simulation of methanol (Sindzingre and Klein, 1991). The bond angles on the 5-mer were constrained to be ∼110° while the dihedral angles were freely rotating. The diameters of the bonding sites were set to be 0.15 times the diameter of the repulsive sites in all cases. For further details of the molecular models, we refer the reader to Table 1. 3. Discontinuous Molecular Dynamics Technique Discontinuous molecular dynamics (DMD) simulation basically amounts to searching for future collisions and scheduling them. The forces between the atoms not colliding do not change because their potential is flat. Those atoms simply get closer together or farther apart. The colliding atoms exchange momentum, and their future prospective collisions must be computed and scheduled. At the time of the collision the positions of all the interaction sites should be updated 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 of the modern DMD technique 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. 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. The binary tree structure is used in tracking events. In previous algorithms, identifying the next event required a search of roughly order N1/2, but tracking events with a binary tree only requires a computational effort of order log2(N). This means that all computations per collision require a time of order log2(N). This algorithm can be extended to chain molecules by applying square wells between covalently bonded

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996 2371 Table 2. DMD-B Simulation Results of Methanola packing fraction (η ) VmolF) av lifetime of HB (ps) av collision time step (fs) compressibility (Z ) p/FkT) fraction of monomer av temp (K) radial distrib of dumbbell as ref system [g(σ+)EE] radial distrib of eq hard sphere at same η [gHS(σ+)] 1-mer 2-mer 3-mer 4-mer 5-mer 6-mer 7-mer 8-mer or more

η ) 0.1 8.1 6.43 × 10-3 1.199 ( 0.019 0.7450 ( 0.0381 299.2 ( 1.5 0.6463 1.2963 0.745 0.227 0.028 0.000

η ) 0.2 8.5 6.32 × 10-3 1.766 ( 0.064 0.5100 ( 0.0331 303.1 ( 6.0 1.1414 1.7188 0.510 0.260 0.133 0.047 0.040 0.009 0.001 0.000

η ) 0.3 9.4 6.36 × 10-3 2.749 ( 0.048 0.3250 ( 0.0331 299.2 ( 5.2 2.1811 2.3470 0.325 0.213 0.158 0.090 0.051 0.062 0.041 0.060

η ) 0.4 9.5 6.23 × 10-3 4.646 ( 0.015 0.2150 ( 0.0341 300.4 ( 9.1 4.0889 3.3333 0.215 0.150 0.111 0.100 0.095 0.064 0.055 0.210

a The system is an NVE ensemble of 108 molecules with periodic boundaries. /k ) 2013.5 K. Repulsive site diameter σ ) 0.30 nm. Center of well between repulsive sites is at 0.51σ; well width is 0.1σ. Center of well between repulsive site and donor/acceptor is at 0.475σ; well width is 0.01σ. Molecular time of 1 ns.

Table 3. DMD-B Simulation Results of 5-Mer packing fraction (η ) VmolF) av lifetime of HB (ps) av collision time step (fs) compressibility (Z ) P/FkT) segment fractn of unbonded acceptor head segment middle segment av temp (K)

η ) 0.1 9.0 3.72 × 10-3 1.941 ( 0.070

η ) 0.25 6.0 3.60 × 10-3 3.636 ( 0.086

η ) 0.3 6.6 3.46 × 10-3 3.843 ( 0.068

η ) 0.4 7.6 3.45 × 10-3 9.059 ( 0.050

0.8959 ( 0.0131 0.8966 ( 0.0143 308.9 ( 2.7

0.5773 ( 0.0197 0.5925 ( 0.0227 308.5 ( 4.1

0.4168 ( 0.0212 0.4141 ( 0.0147 298.3 ( 4.4

0.1294 ( 0.0119 0.1250 ( 0.0168 302.6 ( 4.9

a The system is an NVE ensemble of 108 molecules with periodic boundaries. /k ) 2013.5 K. Repulsive site diameter σ ) 0.30 nm. Center of well between repulsive sites is at 0.51σ; well width is 0.1σ. Center of well between repulsive site and donor/acceptor is at 0.475σ; well width is 0.01σ. Molecular time of 4 ns. Average time of 2 ns.

sites. These square wells are centered on the covalent bond length and assigned a width within which the bond may vibrate. In our simulations, the covalent bond lengths and vibration wells between repulsive sites were fixed at 0.51 ( 0.05 times the repulsive site diameter; the bond lengths and vibration wells between repulsive site and bonding site (donor/acceptor) were fixed at 0.475 ( 0.005 times the repulsive site diameter. Because all sites are simply represented as colliding spheres, the calculation of pressure (or equivalently, the compressibility factor) reduces to a simple summation of momentum changes over all collisions through the virial theorem (cf. Dennlinger and Hall, 1990). It should be noted that pressure calculations are somewhat less straightforward for chain molecules via Monte Carlo simulation. The implementation of the binary tree in the DMD algorithm was originally developed by Rapaport (1980) and applied to chain molecules by Rapaport (1979). The extension of DMD to bonding molecules (DMDB) is straightforward. The square-well association sites must simply be added to the model and the equations of motion integrated. However, the square-well associations are different from the hard-sphere repulsions in that they may “bond”, when two sites are attracted to each other, “dissociate”, when two square wells pull away from each other, or “bounce”, when the association sites try to pull away from each other but their momenta are not high enough to overcome the association energy. The adaptation of the DMD-B methodology for thermodynamics and transport properties was previously demonstrated for spherical molecules (Liu et al., 1994). The extension of the algorithm to study chain molecules and the kinetics of hydrogen bonding as well as the thermodynamics are entirely consistent with the detailed development presented by Liu et al. (1994). One more note about simulation technique for these systems is in order. Note that the average lifetime of a

hydrogen bond is about 3-9 ps (Tables 2 and 3; also see Ladanyi et. al., 1993, and Matsumoto et al., 1990), and that the estimated temperatures exhibit significant fluctuations from the targeted values of 300 K even after 2 ns of simulated time. These observations are significant because many previous dynamic simulations have focused on time frames of ∼150 ps. Our results suggest that the strong energetics of hydrogen bonding lead to a much longer time constant for equilibration, because it takes longer for the bonding sites to sample prospective partners and find their optimal chemical potentials. When coupled with the inherently slower diffusion of chain molecules at high packing fraction, the dynamics of these systems should be expected to require much longer simulations than were common in previous studies. These observations underscore the importance of developing a highly efficient molecular simulation methodology as we move toward detailed analysis of more complex chemical systems. 4. Wertheim’s Theory for the Extent of Bonding at Equilibrium Wertheim’s (1984a,b; 1986a,b) first-order thermodynamic perturbation theory (TPT1) permits estimation of the extents of bonding for potential models like those described in section 2. Given the extent of bonding, the equilibrium thermodynamic properties can be efficiently and accurately estimated (cf. Puhala and Elliott, 1993). The TPT1 theory is distinct from chemical theory in that it provides a direct link between the microscopic potential model and the macroscopic thermodynamic properties. Chemical theory describes the extents of bonding through the parameters that characterize the equilibrium constant, but there is no connection with the potential function, except through Wertheim’s theory. The advantage of having a direct link between the potential model and the thermodynamics is that the

2372

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996

theory can be tested through molecular simulation without the introduction of adjustable parameters. In every test to date, the theory has proved to be remarkably accurate (cf. Ghonasgi and Chapman, 1993). The essence of Wertheim’s theory is to relate the extent of bonding to the geometric characteristics and energy of the bonding sites. To achieve this, the attractive interactions of the bonding sites are viewed as perturbations to a purely repulsive reference fluid. The probability of bonding is particularly influenced by the structure of the reference fluid (as reflected by the radial distribution function) in the vicinity of the bonding sites. The equations below outline the results of this approach, and we refer to the original literature for further details of the derivation. In general, the extent of bonding is related to the fluid structure and bonding by the following equation (Wertheim, 1984ab, 1986ab; Jackson, et al., 1988):

1

XA ) 1+

∑ D

(1a)

FXD∆AD

1

XD ) 1+

∑ A

(1b)

FXA∆AD

were XA is the fraction of unbonded acceptors, XD is fraction of unbonded donors, and F is number density of chain segments at the packing fraction of molecules. The analytical expression for a purely repulsive reference system can be written as (Jackson et al., 1988)

∆AD ) 4πgref(σ+)KADFAD

∫σr

KAD ) σ2

〈fAD(12)〉ω1,ω2 dr12 FAD

c+2d

FAD ) exp(AD/kT) - 1

(2) (3) (4)

where gref(σ+) is the contact value of the radial distribution function of the purely repulsive reference fluid, σ is the diameter of the repulsive site, rc is the diameter of the bonding site, d is the covalent bond length of the bonding site to the repulsive site, AD is the hydrogenbonding energy; k is the Boltzmann constant, T is temperature, and 〈fAD(12)〉ω1,ω2 is the angle average of the site-site Mayer f-function over all orientations of molecules 1 and 2. For the methanol model in Figure 1b, eqs 1a and 1b form a closed system that can be solved for XA and XD. The fractions of acceptor/donor sites unbonded become

-(1 + F∆AD) + x(1 + F∆AD)2 + 4F∆AD XD ) 2F∆AD

(5)

noting that

FXA∆AD ) F∆AD2XA ∑ A

and

FXD∆AD ) F∆ADXD ∑ D

for methanol because the two acceptors are identical by symmetry arguments and there is only one donor. The equations above illustrate the roles of the geometric factors and the fluid structure in a reasonably simple manner, but there is a subtle problem with

generalized implementation of Wertheim’s theory. As originally developed, the theory only applies to spherical molecules, for which the choice of a hard-sphere reference system is obvious. When one considers ways of adapting the theory to nonspherical molecules, one must choose between a hard-sphere repulsive reference system or a hard-chain reference system. The most obvious choice would seem to be the hard-chain reference system (i.e., gref ) gHC, where gHC is the site-site distribution function of the hard-chain fluid). It should be noted, however, that a key consideration is the relatively narrow region where the bonding is expected to occur. Chapman and co-workers (Ghonasgi and Chapman, 1993) have suggested that the bonding can be modeled by estimating the structure from an equivalent isolated sphere (EIS) at the same packing fraction. The EIS approximation is to set gref(σ+) ) gHS(σ+). To elucidate the best approximation, we consider both approximations in section 5 and compare them. Our calculations approximate the hard-sphere radial distribution function as described by Verlet and Weis (1972). For the 5-mer model in Figure 1c, it is necessary at the outset to make some assumptions about the equivalence between sites at the head of the molecule vs those in the middle. Otherwise, one obtains more unknown variables than equations. As our initial estimate for the theoretical treatment, we neglect screening effects and assume that the head sites are equivalent to the middle sites. With this assumption the fraction of unbonded acceptor sites is

XA )

-(1 + F∆AD) + x(1 + F∆AD)2 + 8F∆AD 4F∆AD

(6)

noting also for the 5-mer,

FXD∆AD ) 3F∆ADXD ∑ D

and

FXA∆AD ) F∆ADXA ∑ A

because each acceptor can bond with any of the three donors, but there is only one acceptor of the specified type on each molecule since we compare estimates for the head acceptor and the middle acceptor in the discussion below. 5. Results and Discussion Tables 2 and 3 summarize our results for methanol and 5-mer. Our results are discussed in three parts. In section 5.1 the equilibrium screening effects are illustrated through the monomer fraction for methanol and fraction of unbonded acceptor for 5-mer and by referring to the radial distributions of the corresponding reference fluids. The monomer fraction for methanol is defined as the fraction of all methanol molecules that are not bonded on any of the three bonding sites. The fraction of unbonded acceptors in the 5-mer is similar but compares only the acceptors at the head site to those at the middle site. It is necessary to independently recognize two segments in the 5-mer that might participate in bonds. Turning to section 5.2, the kinetics involved in the hydrogen bonding and dissociation at low density are explored in detail. The kinetic trends at low density reflect the essence of the screening effect. In section 5.3, the rate constants of hydrogen bonding at high density are discussed, explaining the temperature dependency of the forward rate constant at zero activation energy.

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996 2373

Figure 2. Fraction of monomer vs packing fraction of methanol at 300 K: Wertheim’s TPT1 model compared to DMD-B simulation results. The hydrogen bond energy is 2013.5 K.

Figure 3. Site-site radial distribution of methanol and hard dumbell fluids at 300 K and packing fractions of 0.1. The hydrogen bond energy is 2013.5 K for methanol.

5.1. Segment Fractions of Unbonded Acceptors at Different Temperatures and Different Packing Fractions. The equilibrium results for the 5-mer are somewhat more complicated than those for the methanol model. Therefore, we discuss the 5-mer after providing preliminary understanding from the explanation of the results for methanol. 5.1.1. Extent of Bonding in Methanol. The results for methanol are presented in Table 2. Results for a similar molecular model have been presented by Kolafa and Nezbeda (1987) based on Monte Carlo simulation. The results in Table 2 are comparable to those of Kolafa and Nezbeda. Figure 2 illustrates the trend of the monomer fraction with respect to density and includes the predictions based on Wertheim’s theory using the EIS approximation and also applying the site-site distribution function from the appropriate hard-chain reference system, a hard diatomic. Clearly, the EIS approximation provides an accurate estimate of the extent of bonding in this system. The next problem is to understand why it is so accurate for methanol. Figure 3 compares the distribution functions of the two prospective reference systems mentioned in section 4 and the site-site distribution function around the tail site in the methanol molecule. The spherically averaged

Figure 4. Segment fraction of unbonded acceptors vs molecular packing fraction for 5-mer at 308 K and hydrogen bond energy of 2013.5 K.

distribution function for the tail site (gTT) shows clear evidence of a screening effect in that the function is less than unity for r < 1.51σ. This reduction is a result of the covalently bonded atom blocking access to the central site over a significant range of the angles included in computing the spherical average. The distribution function for the tail site is practically identical to that of the diatomic molecule (gEE), indicating that the tail site is unaffected by the extent of hydrogen bonding at the head site. The distribution function of the hard sphere (gHS), on the other hand, rises to a value much greater than unity for r < 1.51σ. Considering that the spherically averaged distribution functions around the methanol sites most closely resemble those around the hard diatomic molecule, the higher accuracy of the EIS approximation may still seem to be a mystery. The answer to this mystery lies in the observation that bonding does not proceed via spherically averaged local structure. As discussed by Ghonasgi and Chapman (1993), bonding proceeds through the alignment of bonding sites over a narrow range of prospective bonding angles. Over this narrow range of bonding angles, the site-site distribution functions more closely resemble those of a hard sphere than those of a hard diatomic molecule. This is a subtle but fundamental insight that impacts much of the analysis of geometic effects in bonding situations. Since local compositions are generally based on spherically averaged distribution functions, it should be of special interest in studies of the effects of local compositions on reactivity (cf. Roberts et al., 1995). Based on this analysis of simulated methanol, we conclude that the extent of bonding is unaffected by screening in such a small molecule. To see any effect of intramolecular interactions, it is evidently necessary to consider larger molecules. The roles of intramolecular interactions become more apparent in the analysis of the 5-mer. 5.1.2. Extent of Bonding in a 5-Mer. General results for the 5-mer simulation are given in Table 3, and the overall trend with respect to packing fraction is illustrated in Figure 4. Figure 4 also shows results of the EIS approximation as well as a modified theory. Two significant observations can be drawn from these results: (1) the extent of association at the middle site is indistinguishable from that at the head site, and (2) the extent of association is much greater than expected from the EIS approximation. The first observation

2374

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996

Figure 5. Radial distributions of 5-mer chain molecules without hydrogen-bonding sites at 300 K and a packing fraction of 0.3.

Figure 6. Segment fraction of unbonded acceptors vs temperature for 5-mer at packing fraction of 0.25 and hydrogen bond energy of 2013.5 K.

indicates that the effect of enhanced screening at the middle site is either negligible or it is canceled by some other effect not previously anticipated. From the radial distributions of a reference 5-mer chain without bonding sites, we can see a big difference between end site (gHH) and middle site (gMM) (Figure 5). The differences between middle site and head site in Figure 4 do not reflect this screening effect to the extent indicated by Figure 5. It should be noted that Figure 5 is based on spherically averaged site-site distribution functions whereas it would be more appropriate to average over a narrower range of angles where bonding might occur. Nevertheless, the results in Figure 4 appear to discount the need of accounting for screening effects in this system. The second observation indicates that there is at least one unanticipated effect operating on this molecular system. It turns out that a single modification to the theoretical formalism is capable of explaining both of these observations. That is, there is apparently some stabilization of the bonded species that is not recognized in Wertheim’s formalism. Wertheim’s formalism simply describes the probability of species coming together; the probability of their staying together is assumed to be the same for all species. This assumption is reasonable if all bound species are identical, as in Wertheim’s original formulation for pure spherical species. But extensions of Wertheim’s theory beyond his original development must apparently account for any subtle differences between the systems originally considered and those of interest in extensions. In retrospect, the source of the stabilization in the 5-mer simulated here may seem obvious. The detailed results in Table 4 show that there is a strong tendency for species that are bound at, say, the head segment to

also be bound at the middle segment. Due consideration of the geometry of such a doubly bound species shows that it is a hexagonal ring, equivalent geometrically to dimer formation in a carboxylic acid. Such hexagonal rings have been recognized as stable for many years. What is interesting is that this stabilization is encountered through entirely classical mechanics; DMD-B does not recognize quantum effects. We can formulate a simple correction to the EIS approximation by applying a “stabilization factor” to the R values of species that are singly bound. The mass balances of species that are singly bound vs those that are doubly bound are fairly straightforward and yield the dashed curve in Figure 4 as the result for the average extent of association for all acceptor sites. A single value for the stabilization factor of 18 has been applied at all densities in Figure 4. By symmetry, doubly bound species are as equally bound at the head site as at the middle site and this largely explains the first observation listed in the preceding paragraph. Further clarification of the roles of screening and stabilization for singly bound species should await simulations of species with single bonding segments, because the fraction of singly bound species in the present study is too small to allow a statistically strong analysis. This clarification study is in progress and will be reported in the near future. The fraction of unbonded acceptors is also weakly correlated with the temperature, as illustrated in Figure 6. It may be noticed that the fraction of unbonded acceptors on the middle segment increases more than that of the head segment with the increase of temperature. This may be due to the rotational swing of the sites surrounding the middle segment increasing screening since the higher temperatures give rise to more

Table 4. Fractions of Unbonded Acceptor Sites Based on the Condition of the Comolecular Acceptor Site in Di-5-mer Simulationsa acceptor

comolecular acceptor

a

packing fractn of molecule (η)

Xm A

XhA

Xm A

0.1 0.25 0.3 0.4

0.902 0.825 0.665 0.303

0.901 0.803 0.669 0.256

0.856 0.277 0.235 0.131

not bonded

XhA

not bonded on middle segment Xm A

not bonded on head segment XhA

not bonded on either segment XA

0.854 0.248 0.238 0.106

0.897 0.593 0.414 0.153

0.896 0.577 0.417 0.129

0.808 0.476 0.277 0.0393

bonded

The hydrogen-bonding energy (/k) is 2013.5 K.

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996 2375

rapid rotation. The temperature effect also plays a role in the reaction rates as discussed further in the next sections. The advantage of the quantitative analysis presented here is that we can envision development of a predictive theory. Wertheim’s theory with the EIS approximation provides a firm basis for initial approximations. The presence of intramolecular repulsive sites leads to a screening effect that must depend on the number of repulsive sites in the chain. Quantifying the magnitude of this effect will involve straightforward extension of studies like the one presented here. A number of chain lengths must be studied, and the position of the bonding sites must be varied systematically with single bonding segments per molecule. Compilation of the data will yield a series of correction factors for the hard-sphere reference system. Finally, the roles of multiple intramolecular bonding sites must be characterized in a similar way. It would be especially interesting to include multisegment backbones which mimic protein backbones. The incorporation of these correction factors into Wertheim’s theory would provide a comprehensive theory capable of predicting the roles of hydrogen bonding that would substantially advance the current state of the art. 5.2. Screening Reflected in Rate Constants of Hydrogen Bonding at Low Density. The rates of hydrogen bond association and dissociation are easily evaluated in our simulations, but the rate constants require a detailed definition before proceeding to analysis of their trends. In order to define a rate constant, it is necessary to divide the rate by a concentration. This concentration must be specifically defined, and it must be recognized that the concentrations are fluctuating with time. We have chosen to define a time interval of 50 ps as the basis for our calculations of rate constants. We found that the overall concentrations changed smoothly over this interval, thus providing a reasonable basis for analyzing trends in the bonding rate phenomena. The rate constants from all the 50 ps intervals were averaged, and the standard deviation was computed at each set of conditions. The definition of the association rate constant, ka, on the basis of a 50 ps interval is then

ka ) no. of bonds formed/{NA50(10-12)[A][D]} where no. of bonds formed is a simple count of the number of bonds formed during a 50 ps interval. [A] is the concentration (mol/cm3) of acceptor sites not bonded averaged over that same 50 ps interval; [D] is the concentration (mol/cm3) of donor sites not bonded averaged over that same 50 ps interval; NA is Avogadro’s number. Similarly, for the dissociation rate constant, kd

kd ) no. of bonds broken/{NA50(10-12)[AD]} where [AD] is the concentration (mol/cm3) of bonded sites averaged over that same 50 ps interval. Analogous relationships for the 5-mer model simply recognize the distinction between a head segment and a middle segment. It should be noted that these definitions are in terms of concentrations, not activities. Traditionally, one might expect to relate the rate constants to each other via the equilibrium constant. However, the high densities and solution nonidealities to be expected for this

Figure 7. Association rate constant (ka) of hydrogen bond formation vs temperature for 5-mer at a packing fraction of 0.1. Table 5. Apparent Activation Energies (K) in Figures 7-10 (ln k vs Reciprocal Temperature) η

head segment

middle segment

0.1, dissociation 0.3, dissociation

-775 -1609

-473 -1779

system make it difficult to apply any simple relations. At this point, we prefer not to overinterpret the few results that we have generated to date. Therefore, we limit the discussion below to some of the more obvious trends, such as the difference between the head segment and the middle segment. These kinetic results do show similarities to traditional kinetic data and we hope to develop a more complete description of them as we perform further investigations. Figures 7-10 present rate constants vs temperature for the association and dissociation rates at two densities. Table 5 lists the slopes in Figures 7-10 as apparent activation energies. All lines in Figures 7-10 were generated from the simulation data using leastsquares regression methods. Figure 7 indicates that the rate constant on the head segment is higher than that on the middle segment. Here the effect of screening on the rate is much more significant than the effect on the extent of bonding discussed in the previous section. Figure 7 also shows that the forward rate constant changes with temperature. Since the activation energy of bonding is zero in principle, a hydrogen bond should form whenever the proton donor meets the proton acceptor. If the overall hydrogen-bonding process were controlled by the intrinsic reaction step, the lines in Figure 7 should be flat according to the Arrhenius rate law. However, the figure shows that the bonding rate constant increases with temperature, just as the diffusion coefficient is increasing. This observation would seem to be more consistent with the collision theory of rate constants than the Arrhenius perspective (cf. Moore, 1972). Figure 8 shows the correlation between the rate constant of dissociation and temperature at a packing fraction of 0.1. As shown in Table 5, the slopes of both lines are significantly less than the intrinsic activation energy for dissociation, given by the hydrogen bond energy in our models. In unimolecular dissociations, it is not unusual for the apparent activation energy to be slightly less than the known dissociation energy, but the reduction is not usually as significant as indicated in Figure 5. This is a result that requires further investigation. Similar to the association rate constant,

2376

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996

Figure 8. Dissociation rate constant (kd) for hydrogen bond breaking vs temperature for 5-mer at a packing fraction of 0.1 and bonding energy of 2013.5 K. The units for kd are s-1.

Figure 9. Association rate constant (ka) of hydrogen bond formation vs temperature for 5-mer at a packing fraction of 0.3.

there is a significant difference between the head segment and middle segment for dissociation at this density. 5.3. Screening Reflected in Rate Constants of Hydrogen Bonding at High Density. Figure 9 shows that the association rate constant at a moderately high density exhibits a mild temperature dependence similar to that in Figure 7. Once again, diffusion is likely to be playing a role in this behavior. The magnitude of the rate constant at high density is larger than that at low density. For the association rate constant, the difference between the middle segment and head segment is still evident but smaller in magnitude. Figure 10 shows that the apparent activation energies for dissociation are near to the hydrogen-bonding energy at this high density. Thus, the process of dissociation seems to be easier to understand at high density for these hydrogen-bonded species. The problem at high density, however, is that the dissociation rate constant shows much broader fluctuations. For the dissociation rate constant, the difference between the head segment and the middle segment at the higher density is much less significant. 6. Conclusions and Recommendations Based on several simulations for relatively simple molecules, we have been able to show several instances of intramolecular interactions affecting the intermolecular hydrogen bonding. When considering the equilibrium bonding, we have analyzed the site distribution

Figure 10. Dissociation rate constant (kd) for hydrogen bond breaking vs temperature for 5-mer at a packing fraction of 0.3 and bonding energy of 2013.5 K. The units for kd are s-1.

functions to show that the screening from the nearest neighboring atoms can be practically neglected, consistent with findings by Chapman and co-workers for molecules like methanol. For larger molecules, however, the intramolecular interactions cause somewhat complicated phenomena. At low density, the effect of screening is evident in the rate constants for association and dissociation. The similarity in the extent of bonding between the head site and middle site in the 5-mer suggested some influence which tended to cancel the screening in the association rate constants. Furthermore, the deviations from the EIS approximation suggested that the equilibrium extents of bonding were higher overall than would be expected from a screening effect. By invoking a correction factor accounting for the stabilization of a hexagonal ring when both sites in the 5-mer were bound, we could explain the higher extent of bonding and much of the cancellation of the screening effect. Consistent with these observations, we expect that the screening effect should be stronger for longer chains. Furthermore, the prospect for intramolecular bonding in longer chains brings up many possibilities for extending this work to analyze loop formation in cross-linking systems, as well as improving our understanding of screening vs bonding. To explore these possibilities, we are currently working on simulation of a 9-mer which we expect to show more dramatic screening effects at the middle site. This should prove to be a valuable complement to the 5-mer study in that it will indicate a trend in the screening vs chain length. Acknowledgment This material is based upon work supported by the National Science Foundation under Grant CTS 9110285. Furthermore, Walter Chapman was very helpful in clarifying the nature of the equivalent isolated sphere approximation. Nomenclature [A]:

concentration (mol/cm3) of acceptor sites not bonded averaged over 50 ps interval

[D]:

concentration (mol/cm3) of donor sites not bonded averaged over 50 ps interval

[AD]:

concentration (mol/cm3) of bonded sites averaged over 50 ps interval

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996 2377 NA:

Avogadro’s number

T:

temperature

XA:

fraction of unbonded bonding sites

d:

covalent bond length of chain repulsive sites

k:

Boltzmann constant

ka:

hydrogen bonding rate constant

k d:

dissociation rate constant

rc:

diameter of bonding site

gEE:

site-site distribution function for end sites on a hard diatomic molecule

gHC(σ+):

contact value of the hard-chain site-site distribution function

gHS(σ+):

contact value of the hard-sphere radial distribution function

gTT:

site-site distribution function for tail sites on a methanol molecule

〈fAD(12)〉ω1,ω2: angle average of the site-site Mayer ffunction over all orientations σ:

diameter of chain repulsive site

AD:

hydrogen-bonding energy

F:

number density of chain repulsive sites at the packing fraction of molecules

Literature Cited Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. SAFT: Equation of State Solution Model for Associating Fluids. Fluid Phase Equilib. 1989, 52, 31. Chapman, W. G.; Gubbins, K. E.; Jackson, G.; Radosz, M. New Reference Equation of State for Associating Liquids. Ind. Eng. Chem. Res. 1990, 29, 1709. Dennlinger, M. A.; Hall, C. K. Molecular Dynamics Simulation Results for the Pressure of Hard-Chain Fluids. Mol. Phys. 1990, 71, 541. Elliott, J. R., Jr.; Suresh, S. J.; Donohue, M. D. A Simple Equation of State for Nonspherical and Associating Molecules. Ind. Eng. Chem. Res. 1990, 29, 1476. Ghonasgi, D.; Chapman, W. G. Theory and Simulation for Associating Chain Fluids. Mol. Phys. 1993, 80, 161.

Jackson, G.; Chapman, W. G.; Gubbins, K. E. Phase Equilibria of Associating Fluids Spherical Molecules with Multiple Bonding Sites. Mol. Phys. 1988, 65, 1. Kolafa, J.; Nezbeda, I. Monte Carlo Simulations on Primitive Models of Water and Methanol. Mol. Phys. 1987, 61, 161. Ladanyi, B. M.; Skaf, M. S. Computer Simulation of HydrogenBonding Liquids. Annu. Rev. Phys. Chem. 1993, 44, 335. Liu, J. X.; Bowman, T. L., II; Elliott, J. R., Jr. Discontinuous Molecular Dynamics Simulation of Hydrogen-Bonding Systems. Ind. Eng. Chem. Res. 1994, 33, 957. Matsumoto, M.; Gubbins, K. E. Hydrogen Bonding in Liquid Methanol. J. Chem. Phys. 1990, 93, 1981. Moore, W. J. Physical Chemistry, 4th ed.; Prentice-Hall: Englewood Cliffs, NJ, 1972. Puhala, A. S.; Elliott, J. R., Jr. Correlation and Prediction of Binary Vapor-Liquid Equilibrium in Systems Containing Gases, Hydrocarbons, Alcohols, and Water. Ind. Eng. Chem. Res. 1993, 32, 3174. 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. Roberts, C. B.; Brennecke, J. F.; Chateauneuf, J. E. Solvation Effects on Reactions of Triplet Benzophenone in Supercritical Fluids. AIChE J. 1995, 41, 1306. Sindzingre, P.; Klein, M. L. A Molecular Dynamics Study of Methanol near the Liquid-glass Transition. J. Chem. Phys. 1991, 96, 4681. Vasudevan, V. J.; Elliott, Jr., J. R. Structure and Thermodynamics of Chain Molecular Fluid Mixtures. Mol. Phys. 1992, 75, 443. Verlet, L.; Weis, J. J. Molecular Motions in Simple Liquids. Faraday Symp. Chem. Soc. 1972, 3, 116. Wertheim, M. S. Fluids with Highly Directional Attractive Forces. I. Statistical Thermodynamics. J. Stat. Phys. 1984a, 35,19. Wertheim, M. S. Fluids with Highly Directional Attractive Forces. II. Thermodynamic Perturbation Theory and Integral Equations. J. Stat. Phys. 1984b, 35, 35. Wertheim, M. S. Fluids with Highly Directional Attractive Forces. III. Multiple Attraction Sites. J. Stat. Phys. 1986a, 42, 459. Wertheim, M. S. Fluids with Highly Directional Attractive Forces. IV. Equilibrium Polymerization. J. Stat. Phys. 1986b, 42, 477.

Received for review February 10, 1995 Accepted April 15, 1996X IE9501079

X Abstract published in Advance ACS Abstracts, June 1, 1996.