6850
J . Phys. Chem. 1989, 93, 6850-6855
PS I particle^.^ Since room-temperature curve analyses have not been published for absorption spectra of larger PS I particles and the temperature-dependent component bandwidths influence the EET dynamics (cf. eq 4), we used the latter simulation in our calculation. In particular, five Gaussian components t,(u) were used, centered at 650.7, 659.6, 668.6, 675.5, and 683.6 nm with the bandwidths and amplitudes reported in ref 7. The corresponding degeneracies g, were set equal to 1, 1, 2, 2, and 2, respectively, for calculations of EET for single-chromophore excitations, and equal to unity for EET between exciton components. The fluorescence spectra were similarly modeled as Gaussian, with bandwidths identical with those of the corresponding absorption components and peak heights proportional to those in t,(u). The fluorescence peaks were arbitrarily shifted 3 nm to the red from the absorption peaks. The resulting l / e decay times in P(u,t) were calculated by using eq 5 for pump-probe wavelengths A, = 660,665,670,675, and 681 nm. These are compared with the experimental depolarization lifetimes 7 in Figure 4, where the theoretical lifetimes evaluated for A, = 675 nm have been normalized to the experimental average at that wavelength. The predicted trends resemble the experimental pattern, particularly at the shorter wavelengths where slower depolarization is observed. Excessive significance should not be attached to these phenomenological calculations, which overlook spatial and orientational factors which are critical to EET in photosynthetic antennae with known g e ~ m e t r y . ~ ' It - ~also ~ ignores the differences between the absorption spectra of PSI-200 and highly enriched particles7 It demonstrates nevertheless that the essential wavelength dependence in 7 can be mimicked by a simple spectral overlap calculation. Such a calculation was suggested by the qualitative anticorrelation observed between the depolarization lifetime and the PSI-200 absorption spectrum. We similarly observed wavelength variation of pump-probe depolarization lifetimes in BChl a-protein aggregates from Prosthecochloris a e s t ~ a r i i .It~was ~ assumed in that work that rapid relaxation proceeded to the longest wavelength spectral component prior to EET. Since the depolarization lifetime should then be independent of wavelength in a structurally homogeneous Scheer, H.; Sauer, P. Photochem. Photobiol. 1987,46,427. (52) Holzwarth, A. P.; Wendler, J.; Suter, G. W. Eiophys. J. 1987, 51, 1. (53) Causgrove, T. P.; Yang, S.; Struve, W. S. J . Phys. Chem. 1988, 92, ( 5 1 ) Sauer, K.;
6790.
solution of BChl a-protein, we conjectured that the slower depolarization at the shortest wavelength may arise from EET among uncomplexed BChl a chromophores. In light of the recent evidence in ref 8 and in the present work that excitation in PS I equilibrates among spectral forms (rather than accumulating in the longestwavelength species), the wavelength variation in T in the BChl a-protein may be due instead to differential EET rates among the different BChl a spectral forms (cf. eq 4). The major conclusions of this work may be summarized as follows: (1) At least three exponential decay components are required to describe the Chl a isotropic photobleaching signal between 660 and 68 1 nm in PSI-200. Triexponential fits yield lifetime components of 1.5 ps) which contributed significantly to the power spectrum in the region of interest. In fact, the double peak observed in each spectrum and recorded in Table 111 is probably mainly due to this noise. Nevertheless, these calculations are sufficiently precise to provide a useful test for a cavity potential. Frequencies for the static host lattice have been obtained from MD simulations using very large masses for the atoms in the host lattice. Comparison of the frequencies given in Table I11 supports the findings of Basu and Mountain:23 the potentials fitted to the static hydrate give rise to essentially the same frequencies as the static hydrate simulations, but these are very different from those arising from the dynamic hydrate. However, the situation is much more encouraging when one considers the potentials fitted to the thermally averaged energies. In all three cases there is now a range of frequencies that is comparable to that defined by the two M D frequencies. Indeed, for the lower temperatures the range of frequencies is nearly coincident with that found in the simulations, and although the frequencies from the potential fitted at 270 K are too high, they still do overlap with the MD frequencies. Thus, it appears that the cell theory can give a reasonable description of the dynamics of the guest provided that the cavity potential is properly averaged. In applying the cell theory it is usually assumed that the one cavity potential is suitable for all conditions, and we tested this for the frequencies by taking the fitted potentials for one temperature and using them to predict the methane frequencies for the other two temperatures. In each case it was found that the anharmonicity of the fitted potential caused the frequency to change by -2 cm-' over the temperature range 145-270 K; this is only half the value of -4 cm-' actually observed in the M D simulations. It is possible that this difference is merely an artifact of the statistical noise mentioned above. However, the temperature variations in the broad envelope of the power spectra did seem to remain consistent when the length of each simulation was varied, and so it is possible that this reflects real changes in the shape of the cavity potential as the temperature changes. Langmuir Constants. As mentioned in the Introduction, Langmuir constants (C(T)) play a central role in predicting the phase behavior of gas hydrates, and so it is important that any theory of gas hydrates should adequately describe their values. Although the cell theory usually calculates C(7') with a spherically averaged potential, there are a variety of different ways that this average may be obtained; for example, one could use an unweighted spherical averageZoor, in view of eq 1, an exponentially weighted spherical average. Since these procedures will generally lead to different descriptions of the cavity potential, it is more appropriate to describe the Langmuir constant as C(T,U),and it therefore becomes important to know how sensitive C is to variations in U. The fitted potentials described above provide a convenient means of assessing how a spherical average of Umay be performed: since they were chosen to indicate the limits of the angular variation in (U(r,@,@)), they should provide upper and lower bounds for Langmuir constants obtained from any orientationally averaged cavity potential. Also, it is desirable to use a single cavity potential to describe a range of conditions; but we have already shown that Uvaries with temperature. Again the fitted potentials provide a convenient means of assessing the sensitivity of C to these variations in the cavity potential, and so we have calculated C(270 K, U(T))using the cavity potential determined at temperatures of 0 (the static hydrate), 145, 200, and 270 K. It is stressed that this temperature variation relates only to the conditions used to define the cavity potential and is additional to the well-known temperature dependence of C described in eq l . From the results of these calculations, which are listed in Table IV, it is clear that C is very sensitive to both the anisotropy and
Rodger
TABLE IV: The Langmuir Constant at T = 270 K model
C, lod ms2 kg-'
exact calculation (eq 1) from radial potentials" 270 K 200 K 145 K 0 Kb
3.00 i 0.09 1.62 3.64 6.22 14.71
1.64 3.61 6.61 15.30
1.78 3.83 5.48 15.71
2.27 5.06 7.51 16.20
"Temperatures refer to the temperature of the MD simulation to which the radial potential was fitted. The four values for each temperature are four different radial cross sections of the cavity potential (see text). bStatic hydrate model.
the temperature used to define the cavity potential. The extremes of the angular variation define a range of values that is typically about 30% of the midpoint of this range, and although this is an overestimate of the variation expected between differently weighted rotational averages, it does indicate the need to chose the weighting for the angular average appropriately. To further illustrate this point, we calculated C from both an unweighted and an exponentially weighted spherical average of the cavity potential and found that the results differed by about 20%. The Langmuir constants are even more sensitive to the temperature dependence of the cavity potential, with an increase in the temperature used to define the cavity potential from 145 to 270 K producing a 4-fold decrease in the predicted C(270 K). We note from eq 2 that the pressure at which the hydrate will begin to form (the hydrate pressure) is inversely related to the Langmuir constant (through the fugacity of the gas), and so these variations in the calculated value of C(270 K) will lead to large differences in the predicted hydrate pressures. In order to assess the validity of the various estimates for C(T), we have used eq 1 to perform an exact calculation of the Langmuir constant, and the result is given in Table IV; the three-dimensional integral was evaluated with a Monte Carlo technique and then averaged by sampling every fifth step of the trajectory. Due to computational limitations, it was only possible to do this at a single temperature-taken to be 270 K. In view of the results for guest frequencies it is likely that the comparison would be more favorable for the cell theory at lower temperatures; however, 270 K is more representative of the temperatures involved in most practical applications of the cell theory. Comparison of the exact and model results reveals that the exact value lies well outside the ranges predicted from the fitted potentials. In particular, it is well outside the range of the potentials fitted to the simulations of the hydrate at T = 270 K, indicating that a proper specification of the rotational average is not sufficient to give a good description of the Langmuir constant. This discrepancy is probably related to the treatment of the thermal average involved in eq 1. For the fitted potentials (or any of the spherically averaged potentials) it is assumed that the average over the lattice motion can be taken before the exponential of the potential is calculated, whereas this order is reversed in the exact calculation (see eq 1). This interchange will tend to weight the potential more heavily in favor of high energies than it should and will therefore predict a C( 7') that is too low, as is observed. This interpretation is also consistent with the spherical cell model improving as the temperature decreases, since a lower temperature will also make the thermal motion of the host become less significant. Thus, in order to obtain a good estimate of the Langmuir constant, the spherically averaged potential should be defined by exp{-VL(r)/kT) = ( 1 /4n)J
1
(expl-U(r,@,4)/kTl) sin 8 d 4 d@ (6)
and this will usually differ from the average that was shown to give a reasonable description of the guest frequencies:
It is worth emphasizing that this study does not invalidate the cell theory; rather it points out that the evaluation of the cavity
Cavity Potential in Type I Gas Hydrates potential should be done with care, and especially with due consideration to the property being studied. Since different properties will imply different weighted averages of the fundamental intermolecular potentials, they should be expected to give rise to different effective cavity potentials. From an empirical point of view, this means that it is possible to derive experimental cavity potentials but that the resulting potentials may differ according to whether guest frequencies, Langmuir constants, or some other properties are used for the experimental input. Similarly, it is to be expected that the cavity potential will change with the thermodynamic conditions considered, with the hydrate pressures being particularly sensitive to this effect.
4. Discussion One of the major observations to come out of the last section is that the cavity potential exhibits a pronounced temperature dependence, and this could have important ramifications for the way in which hydrate pressures are predicted; in particular, it was pointed out in the last section that the Langmuir constant (and hence the prediction of hydrate pressures) is particularly sensitive to this effect. Consequently, it is important to be able to describe the variations in the cavity potential induced by changes in temperature. It might appear that the inclusion of such a temperature dependence will merely increase the flexibility with which the model can be used empirically and will not provide a fundamentally more accurate description of hydrate stability; however, the results presented above suggest that the cell model will not satisfactorily predict hydrate pressures over a range of temperatures without it. As mentioned in the last section, the temperature dependence appears to be primarily related to changes in the size of the cavity, and so it should be possible to model this trend by allowing the cavity radius, R , to become a function of T . Under the constant-volume conditions used in the present study, any increase in temperature causes an increase in thermal fluctuations in the positions of the water molecules; this effect may be pictured as a thickening of the cavity walls and hence leads to a decrease in the cavity radius. However, hydrates are rarely either found or studied under isochoric conditions, and so an increase in temperature should also be accompanied by thermal expansion of the entire lattice, hence leading to an increase in the average radius of the cavity. Thus, the cavity potential will depend on a balance between two opposing effects: an increase in the cavity radius due to thermal expansion of the entire lattice, and an effective contraction of the cavity due to thermal motion increasing the width of g ( r ) for the cavity. A simple model for these two effects can be derived by assuming that the water molecules are harmonically bound to their average position R,,, with the force constant along the radius of the cavity being b. Since the guest-host interactions will be dominated by the closest interactions, it seems reasonable to assign the radius of the cavity to be a t the inside turning point of the water molecule’s oscillations. Thus, taking the oscillations to involve an average energy of kT, the cavity radius will be given by
R(T) = f ’ ~ 3 ( T , T 0 ) R , , ( T 0 )( 2 k T / b ) 1 / 2 where R,,(To) is the average cavity radius at a specified temperature Toandfl T,To)represents the relative volume expansion of a gas hydrate with respect to the specified reference temperature. Note that Ray(To)is just the cavity radius obtained from the observed crystal structure of the hydrate at the reference temperature, whilefl T,To)is an experimentally known function2* that is expected to be independent of the type of gas molecule c ~ n s i d e r e d .In ~ ~principle, the remaining parameter in eq 8, the force constant b, could be obtained from a harmonic analysis of the force field arising in the water lattice; however, in view of the simplistic nature of this model, it is probably better to treat b as an empirical parameter that is characteristic of the water lattice (and hence independent of the nature of the guest molecule). The other important conclusion to be made from the results of section 3 is that the anisotropy of the cavity potential cannot be ignored. Although it is possible to find a spherically averaged
The Journal of Physical Chemistry, Vol. 93, No. 18, 1989 6855 potential that gives rise to the correct Langmuir constant (eq 6), the exponential weighting in eq 1 will ensure that the resulting potential should vary with temperature (in a manner that differs from that discussed above), and it is not likely to be the same as the empirical potential derived from other physical properties (compare the way the fitted potentials in section 3 gave a reasonable description of the guest frequencies but a poor description of the Langmuir constants). Hence, it is likely that the use of a spherically symmetric cavity potential will necessitate the use of different molecular parameters to describe the behavior of the guest molecules in the bulk gas and in the hydrate. In view of these inadequacies in the spherical approximation, it is worthwhile considering how to describe the cavity potential’s dependence on the direction in which the guest molecule is displaced, Q (=e,$). Recently, there have been some very promising developments in models of the anisotropies in intermolecular potential^,'^ and since these involve methods of describing deviations from spherical symmetry, they should also be applicable to the cavity potential. According to these methods, the orientation dependence of the interaction between two molecules can be described by taking a conventional isotropic potential function and allowing the parameters (e.g., E, u, and a for a Kihara potential) to become functions of R; spherical harmonic functions then provide a convenient means of describing the angular functional forms for these parameters. For the cavity potential this means that the orientational dependence can be introduced by taking the parameters in eq 4 (R, e, u, and a ) and expressing them as the sum of a small number of spherical harmonic functions. One major advantage of this approach is that it is readily extended to account for anisotropy in the guest as well: the guest potential parameters (e, u, and a ) can also be expressed as sum of spherical harmonic functions and then coupled with the cavity potential through the usual combination rules. In order to limit the complexity of the resulting potential, it should be possible to focus attention on just one of two of these quantities, and in view of the studies on intermolecular potential^,^^ it is likely that the hard core size, a, and the cavity radius, R , will provide the best candidates. In this case eq 4 can be used by letting a = a ( Q ) or R = R ( Q ) . It is anticipated that explicit functional forms for a ( Q )or R ( Q )can be obtained by fitting the resulting cavity potential to that found in the MD simulations, and work on this task is in hand.
5. Conclusions In this article we have used the technique of molecular dynamics computer simulation to provide an exact description of a model system that is very similar to type I methane hydrate. We have been particularly interested in the potential energy field that arises in the small cavities of type I hydrates and have used the simulations to test a number of approximations that are frequently made in describing this cavity potential. Three major findings have emerged from this study: (1) the presence of kinetic energy in the host lattice leads to distortions in the cavity and so causes even the small cavities-which, being dodecahedral in shape, should be very nearly spherical-to become significantly anisotropic; ( 2 ) the thermal fluctuations in the positions of the host molecules also mean that the size of the cavities, and hence the nature of the cavity potential, will change with temperature even under isochoric conditions; and (3) although explicit hydrogen atom-guest interactions do not qualitatively change the cavity potential, they do make a substantial quantitative difference, and their presence appears to enhance effects 1 and 2. Note that since we have considered a perfectly spherical guest molecule, the first observation is not due to distortions in the host lattice caused by the inclusion of an aspherical guest but is actually a property of the host lattice. Ways of extending the existing cavity potential models to include these effects have been suggested. Acknowledgment. The author thanks Dr.J. T. K. Tan for his many helpful discussions. Registry No. H 2 0 , 7732-18-5; D20,7789-20-0; CH4, 74-82-8. (35) Stone, A. J.; Price, S.L.J . Phys. Chem.
1988, 92, 3325.
