7280
J. Phys. Chem. 1990, 94, 7280-7288
Monte Carlo Calculation of Orientationally Anisotropic Pair Distributions and Energy Transfer in a Model Monolayer Mengtao Het and Paul Siders* Department of Chemistry, University of Minnesota-Duluth, Duluth, Minnesota 5581 2 (Received: June 28, 1989; In Final Form: January 31, 1990)
The NPT (number of particles, pressure, temperature) Monte Carlo method has been used to calculate the equation of state and orientation- and distance-dependent pair distributions for a model monolayer of hard “UF0”s. A UFO is defined to be the intersection volume of two identical spheres whose centers are more than one and less than two radii apart. As pressure is increased from zero, the UFOs rapidly tip perpendicular to the plane and then more gradually align. By use of the calculated equilibrium configurations, and Forster’s orientation-dependent dipole-dipole energy-transfer rate constant, the bimolecular energy-transfer rate constant was calculated as a function of pressure and density. The macroscopic rate constant is sensitive to the extent of orientational ordering in the model monolayer.
Introduction We have used the Monte Carlo method to investigate dipoledipole energy transfer in a model one-component monolayer of oblate molecules. Our model particle is a hard spheroidlike “UFO”. A UFO is the intersection space of two identical spheres. (See Figures 1 and 17.) The UFOs are allowed to rotate and translate freely, except that their centers are confined to a plane. The elongation of the UFO is a / b = 6 in the system we have studied. The effect of different elongation has been studied elsewhere, for hard spheroids on a plane.’,* The simple geometry of UFOs leads to an analytic overlap criterion, presented in Appendix A. This model was chosen to represent, albeit crudely, a monolayer of large aromatic molecules such as chlorophylls, porphyrins, or polycyclic aromatics. The UFOs represent the repulsive interactions between such large, flat molecules. We have studied the thermodynamic behavior of the model monolayer, pair correlations in the monolayer, and the effect of these equilibrium properties on the kinetics of a simple bimolecular reaction within the monolayer. As described below in the Theory section, we used the constant-NPT Monte Carlo method to produce equilibrium configurations for model monolayers of 32 and 128 UFOs, with periodic boundary conditions. From these configurations, we have calculated the full angle-dependent single-particle and pair correlation functions. The pair correlation functions depend on four orientation variables, as well as on the distance between UFO centers. The angle dependencies of these correlations are discussed in the Results and Discussion section and are compared to correlations in related model^.^-'^ Phase behavior of the model is discussed in terms of the pressure-density isotherm and nematic order parameters. A high-density discontinuous freezing transition is observed, roughly in agreement with the results of Frenkel et al., for spheroids in three dimensions.’*” We observe a lower density transition to a nematic phase, much as observed by Vieillard-Baron’* and Frenkel and EppengaI3 for hard ellipses in two dimensions. As a direct application of the calculated pair correlation functions, we have calculated the rate of dipole-dipole energy transfer within the monolayer, using Forster’s orientation-dependent energy-transfer rate constant.I4 Energy transfer between substituted porphyrins in a monolayer has been studied experimentally by Tweet et al.,15*16 Gonen et al.,” and Frackowiak and Szurkowski.18 The work of Gonen suggests that energy-transfer fluorescence quenching is sensitive to an orientational phase transition in the monolayer. We have found that orientational ordering substantially affects energy transfer in our model system.
* Author to whom correspondence should be addressed. ‘Present address: Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis M N 55455. 0022-3654 f 90 f 2094-7280$02.50f 0
At low density the rate constant is near zero and the reaction appears second order. At intermediate densities, where average UFO separations are such that the Forster mechanism is most applicable, the reaction shows third-order density dependence. The order returns to two when the monolayer freezes. Reaction orders for surface reactions have been found to differ from those for bulk reactions in other model systems, also. The origin and significance of such differences are discussed, for example, by Silverburg and B e n - S h a ~ l and ’ ~ Kopelman.20 Theory Monte Carlo Calculation. The UFO hard-particle geometry is defined in the Introduction and discussed further in Appendix A. The length of the major (C,) symmetry axis is 26, and the length of the perpendicular C, axis is 2a. The centers of the UFOs lie in the ( x , y ) plane of the laboratory coordinate system. Orientation of a UFO (relative to a fixed coordinate system in the plane) is described by Euler angles: a,/3, and y.21 The intermolecular potential is independent of the third Euler angle y, so equilibrium orientation distributions are constant with respect to ( 1 ) Wesemann, J.; Qin, Lihong; Siders, P. Langmuir 1989, 5, 1358. (2) S!ders, P. Mol. Phys. 1989. 68, 1001-1013. (3) Lippel, P. H.; Wilson, R. J.; Miller, M. D.; Woll. Ch.; Chiang, S. Phys. Rev. Lett. 1989, 62, 171-174. (4) Baum. R. Chem. Eng. News 1989, Feb 6, 27-28. ( 5 ) Ghazi, F.; Rigby, M . Mol. Phys. 1987, 62, 1103-1 1 IO. (6) Nezbeda, I . Czech. J . Phys. 1980, 30, 601-612. (7) Frenkel, D.; Mulder, B. M . Mol. Phys. 1985, 55, 1171-1192. (8) Perram, J. W.; Wertheim, M. S . Chem. Phys. Lett. 1984, 105, 277-280. (9) Ward, D. A.; Lado, F. Mol. Phys. 1988, 63, 623-638. (IO) Cummings, P.; Nezbeda, I.; Morriss, G . Mol. Phys. 1981, 43, 1471-1475. ( 1 1) Frenkel, D.; Mulder, B. M.; McTague, J. P. Phys. Reu. Lett. 1984, 52, 287-290. (12) Vieillard-Baron, J. Mol. Phys. 1974, 28, 809-818; Chem. Phys. 1972, 56, 4729. (13) Frenkel, D.; Eppenga, R. Phys. Rev. 1985, 31, 1776-1787. (14) Forster, T. Discuss. Faraday SOC.1959. 27, 7-17. (15) Tweet, A. G.; Bellamy, W. D.; Gaines, G. L., Jr. J . Chem. Phys. 1964, 4 1 , 2068-2077. (16) Gaines, G. L., Jr.; Tweet, A. G.;Bellamy, W. D. J . Chem. Phys. 1%5, 42, 2193-2199. (17) Gonen, 0.; Levanon, H.; Patterson, L. K. Isr. J . Chem. 1981, 21, 271-276. (18) Frackowiak, D.; Szurkowski, J . Mol. Cryst. Liq. Cryst. 1984, 1 1 1 , 359-392. (19) Silverberg, M.; Ben-Shaul, A . J. Stat. Phys. 1988, 52, 1179-1 195. (20) Kopelman, R. Science 1988, 241, 1620-1626. (21) Arfken, G . Mathematical Methods for Physicists; Academic Press: New York, 1970; pp 178-180. A UFO is oriented as follows: Define a local (x,y,z) coordinate system on the UFO, with the z axis being the U F O s principal axis. Initially, the local coordinate system is aligned with the laboratory system. Rotate the local system through an angle a about the laboratory z axis. Then rotate through an angle B about the (newly rotated) local y axis. Finally, rotate through y about the local z axis. All rotations are in the right-hand sense.
0 1990 American Chemical Society
The Journal of Physical Chemistry, Vol. 94, No. 18, 1990 7281
Energy Transfer in a Model Monolayer
The configurational averages (density, order parameters, compressibility, etc.) reported have been examined for serial correlation. No values are reported for which significant linear correlation was found. Pair Distribution. It is at the level of mutual orientation of molecular pairs that orientational order affects bimolecular reaction rate. Pair order is described quantitatively by the pair distribution function and pair cbrrelation function. The location and orientation of the ith particle is given by Xi = (xi,y,,a,,fl,,Y,). We use single and pair probability distributions as defined by Zannoni (ref 27, Chapter 3, eqs 6 and 7), except that we use the independent variable /3 rather than cos 0. In eqs 2 and 3, Z N is
X
p’lvfl)
x
Figure 1. Two UFOs. Plane polar relative coordinates (r2,,B)and Euler angles are shown.
7. The Monte Carlo method of statistical mechanics is reviewed in several excellent books.22-24 This work used the Monte Carlo method at constant N P T (number of particles, pressure, and temperature). The NPT method yields the pressure-temperature equation of state directly and allows the size of the simulation region (planar area, in this simulation) to equilibrate. The intermolecular potential used is a “hard” potential, so the overlap criterion for UFOs presented in Appendix A is central to the simulation. To reduce the number of pair tested for overlap, a neighbor table was maintained. Random numbers were generated with an algorithm published by Whichmann and Hillz5 The basic cell into which N particles are placed is a rectangle. Both the height and width of the rectangle fluctuate freely. Particle displacements are attempted in the x and y directions, and then rotations are attempted. Rotational moves are made in the variables ( a , cos P). The acceptance ratios for cell size, displacement, and rotation were each maintained at about 30% following the suggestions of Kolafa.26 We find that, at sufficiently high density, the thermodynamically stable phase is ordered in both a and 0. We refer to that ordered state as “nematic”. We use the single-particle orientation distributionfla, cos P) to locate a nematic director and to calculate orientational order p a r a m e t e r ~ fandfi2. ~~ We use the Q-matrix m e t h ~ d ’ ~to . ~locate ’ the director. Because of the symmetry of our model system, when a director exists it lies in the (x,y) plane. The orientational order parameters f22 and fi0 defined by Zannoni2’ are given in eq 1. In eq 1, anem is the value of a along
= ( N / Z N )Jexp[-U(IXNJ)/ksT](fisin
p’2’(X,J,)
d22 = (3/8)’/2(1
- cos2 8)
which the nematic director lies. (The value of anem is calculated from a configuration of finite size, so even in the absence of a ordering, a value of a,,,,, is defined. In an a disordered state, a,, is merely an average a.) The p a r a m e t e r ~ fandf22 ~~ are complex. In this paper, we report only their norms. The parameter Ilf2011 reflects only ordering in 6 and is zero in the absence of 0 ordering. The parameter llfizll reflects both P and a ordering and is zero unless both a and fl are ordered.
fi
J
(N(N - 1 ) / ZN) exp [-WIXN), / kr3 TI ( sin i= 1
Bi) l u N 3 l
(3)
the configurational partition function for N particles. In the absence of pair correlations, the two-particle distribution is the product of the two one-particle distributions. The pair correlation function g‘2)is the ratio of the two-particle distribution to the product of one-particle distributions. The deviation of d2)from unity quantifies pair correlations.
d2’(XlJ2)= p’2’(Xlx2)”’’)(X1)
p’I)(X2)1
(4)
In the case of particles that are not orientable, the distribution functions depend only on positions, and not on orientations. In such a case d2)defined above reduces to the radial distribution function g‘2)(r21), with r21being the distance between the centers of particles 1 and 2. In the model considered in this paper, of course, orientations must be considered explicitly so g‘2)cannot be so completely simplified. We need not consider the full set of variables (XIJ2J.Because our model is translationally invariant, we can use relative coordinates for the position of the center of UFO 2 relative to UFO 1. Denote the positional and orientational parts of Xi by ii = (xi+) and Qi = (ai,Pi,yi);-Xi = (ii,Qi). The position of particle 2 relative to particle 1 is r21 = i2 - i l , and fi2)(X1J2)= pC1)(il)-
P ’ 2 V 2 1 ,Ql $2).
For the pair probability distribution, we may measure a2relative to a I . Let azI = a I - a2. Moreover, f i 2 ) does not depend on (yI,y2).Denote the probability densities for a I ,yl, and y2 as follows: f(aI).f(yI),f(y2). Note that W ( i I=) p andfly,) = fly2) = 1/(27r). The location of a1is arbitrary, because the entire monolayer is invariant with respect to rotation, soflal) = 1/(27r). = ~p’~’(~21,~21~Pz,P1)/(27r)~
(5)
Equation 6 defines a reduced pair distribution function G. G(72IP21432,Pl)
= p’2’(~21,~21,P2,PI)/
bfo(a21)f o ( P 2 )
fo(81)l (6)
The distribution f u n c t i ~ n s f ~ ( a randfo(&) ~~) are the low-density ideal-gas limiting distributions for aZ1and Pi, normalized to unity. sin Pi(so that dPi = 1) andfo(a21) That is,f0(Pi) = = 1/(27r). The reduced pair distribution G is closely related to the pair correlation function d2), as in eq 7. For an isotropic fluid, G is
sfo(&)
f ( a 2 , )f ( P 2 ) (22) Valleau, J. P.; Whittington, S. G.Modern Theoretical Chemistry 5A: Statistical Mechanics; Berne, B. J., Ed.; Plenum Press: New York, 1977; Chapter 4. Valleau, J. P.; Torrie, G.M. Ibid. Chapter 5 . (23) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (24) Binder, K. Monte Carlo Methods in Statistical Physics, 2nd ed. Binder, K., Ed.;Topics in Current Physics 7; Springer: Berlin, 1977; Chapter 1. (25) Whichmann, B.; Hill, D. Byte 1987, March, 27-28. ( 2 6 ) Kolafa, J. Mol. Phys. 1988, 63, 559-579. (27) Zannoni. C. The Molecular Physics of Liquid Crystals; Luckhurst, G. R., Gray, G.W., Eds.; Academic Press: New York, 1979; Chapter 3, pp 51-83.
(2)
=
fl2’(X1J2)
do2 = (3/2) cos2 /3 - 1 / 2
Bi)lc”,I
1-1
flPJ
g‘2’(~21,a21,P2,PI) = fo(a21) f O ( P 2 )
fo(P1)
G(~zI,a21rfl21fll)
(7)
identical with g‘2). We chose to calculate G rather than d2) because G represents total pair orientational order, while g(*) reflects only orientational order in excess of orientational order in the single-particle distributions. Moreover, g‘*)is difficult to calculate numerically in highly ordered fluids, because the single-variable orientation distributions f ( a Z Iand ) f(&) are then sharply peaked. The normalization of G,for a distribution of N particles, is N - 1. The long-range (i.e., large 11iZl11 behavior ) of G, in the absence of long-range order, is given in eq 8. The
7282 The Journal of Physical Chemistry, Vol. 94, No. 18, 19‘90
lim
ll~2lll--
G ( ~ ~ I . C Y ~ I ,=P ~8~f(azl)f(P1)f(P2)/(sin ,PI) PI sin Pz)
(8) reduced pair distribution G is calculated on a grid of values of its independent variables {r21,8-al,a21,~2,~l). The relative position vector 721is represented in plane polar coordinates rZI= 1172111and 6. The angle 6 is measured relative to a I : 6 - a l . We have chosen a coarse uniform grid, as follow< A(6-al) = Am21 = A& = A81 = n/10 and Ar21 = b. We have found 104-105 configurations sufficient to reduce relative uncertainties in G to 3-7%. Pair correlations have been calculated by other authors in terms of projections on orthogonal functions: spherical harmonics and related functions. We have not used such an expansion technique because the expansions have been found to be only slowly convergent, especially near contact of two particles. The region near contact is especially important to us, because of our intent to calculate bimolecular reaction rates. The relationship of the pair probability distribution fi2)to the probability of occupancy of each “box” in r21,8-al,a21,P2,Pll space, P$$k, is given in eq 9. The function fibox is simply the fraction
He and Siders usual e x p r e ~ s i o n . ~ *We . ~ ~define the average unimolecular rate constant to be Pi,in eq 15. The constant kuniis the ensemble
rate/area = kUniN(N - I)/area
(16)
average of k,, taking a UFO donor and acceptor pair as a unimolecular complex. The theory of microscopic rates of bimolecular reactions in liquids is in general difficult. We have chosen to study the rate of dipole-dipole energy transfer. Energy transfer is important in many real monolayer and bilayer systems that display orientational order. The porphyrin derivatives involved in some real energy-transfer systems are indeed flat, disklike molecules. The geometry of their intermolecular interaction potentials may approximately be represented by the UFO geometry for which we have discussed equilibrium distributions above. In our model monolayer, the energy donor and acceptor have identical interaction potentials. The donor may be imagined to differ from the fi2)(721421?P2,PI)= acceptor by metal substitution, an attached substituent, or simply We use the unimolecular energy-transfer rate constant ( N - ~ ) ~ ~ ( ~ 2 l ~ ~ - ~ l , ~ 2 l ~ P * ~ P I ) / [ ~ ( ~ 2 1 ) A(9) ~ 2 I A P 2excitation. APII of Forster,14 because it is simple and has been widely used and of all pairs that fall into the box in which the argument of P!!; tested by others. The Forster rate constant is given in eq 17. In lies. The normalization is CP’& = 1. The factor AA(rzl)is the area of a single (6 - a , ) box in an annulus about particle 1 in which particle 2 may lie. The calculation of AA is described in Appendix C. We note that ( N - l ) / ( p A A ) = ( N - l)/[N(AA/A)] = l / ( U / A ) , where A is the average area of the Monte Carlo cell. eq 17, 7 is the mean lifetime of the energy donor, Rois a critical The reduced pair distribution is calculated in the form of eq 10. transfer distance, and rzl is the distance between the centers of the reactants. The orientation dependence of k, is in the factor G(r2i76-ai . 0 2 1 8 2 t P 1 ) = K, where 421is the angle between the transition dipole moment 8r~~(r21,6--(yl,~2l,P2,PI) vectors of molecules 1 and 2, while 42 and dl are the angles (10) [AA(r21)/AlAa21sin (P2) sin (PI) AP2APl between these respective vectors and the line joining the centers of molecules 1 and 2. The UFO model interaction potential does Reaction Rate Constant. The microscopic rates of many binot include any dipole moment. We add a dipole to every UFO molecular reactions are dependent upon the separation and mutual after calculating pair distributions, in order to model dipole-dipole orientation of the reacting molecules. The observable macroscopic energy transfer between UFOs. Any effects these dipoles would rate constant is the average of the microscopic rate over a range have had on density and on pair distributions are neglected. In of reactant orientations and separations. Reaction systems that each particle, a dipole is placed perpendicular to the UFO principal allow control of orientational order may yield rate constants symmetry axis and in the plane of the semimajor axis a . It is in sensitive to the extent of ordering. Chemical reaction, reactant this plane that the carbon skeleton of a molecule represented by orientation, and reactant translation all are dynamic processes. the UFO would lie. If all occur on similar time scales, then separating their effects The angles 421,42, and 41can be expressed in terms-of the Euler is difficult. However, the initial rates of fast reactions, and the angles (a21,P2,y2,/31,yIJ and the polar angle 8 - a l . That transrates of reactions slow compared to rotational or translational formation is discussed in Appendix D. The Euler angles y2 and diffusion, are approximately the average of the microscopic rey I express rotation about the principal symmetry axes of UFOs action rate over the static equilibrium (or near-equilibrium) 2 and 1, so the pair distribution fi2)is independent of y2 and y l . distribution of reactant separation and orientation. We consider Therefore, the k , to be averaged with f12)is itself first averaged in this section the role of equilibrium fluid structure (especially with respect to y2 and y l , as in eqs 18 and 19. This average can orientational order) on macroscopic rates of orientation-sensitive reactions. Let ku(721,R2,Rl)be the “unimolecular” microscopic rate constant for a chemical reaction. The rate constant k, is unimolecular in the sense that it applies to a particular reactant pair at specified separation and orientation. The macroscopic, observable reaction rate per area is given in eqs I 1 and 12. be done analytically. The resulting expression for L,,is given in rate/area = ( 1 / A ) l k , ( X l , X 2 )P(2)(Xl,X2)dXIdXz ( 1 1) Appendix D. The integral over r21must be evaluated from rZI = 26 to rZI= a,where 6 is the semiminor axis of the UFO. The = P S ~ , ( % I . Q ~ , Qfi2)(721,Q2,Ql) ~) dQ2 dQl (12) range of the pair distribution, however, is limited by the number of particles used in the Monte Carlo simulation. We evaluate the The rate per area equals p2kbr.so kbl is related to the pair probintegral numerically from 26 to 206. At large r Z l ,fi2)ceases to ability distribution Pcz) as in eq 13. depend on r Z 1so , all_the large-r2, radial dependence of the integrand is explicit in k,. We integrate analytically over r21from kbi = (1 / P I Sku(721,Q2*fil)fi2)(7213Q29Rl) d72l dQ2 dQ1 (1 3) 20b to m, taking fi2)to be independent of rZI(Le., a function of angles only, as in eq 8) over that range. The integrand is so near In the special case of an isotropic fluid, and a reaction dezero for rZI> 206 (for our choice of Ro) that this approximation , reduces pendent only upon reactant separation ( k , = l ~ & ~ ~eq) )13 is entirely adequate. to eq 14. Equation 14 is the two-dimensional analogue of the
4
kbi
= ~ku(r21)g(2)(r21)r21 dr21 (isotropic fluid and reaction)
(14)
(28) Keizer, J. Chem. Rev. 1987, 87, 167-180. (29) Keizer, J . J . Phys. Chem. 1982, 86, 5052-5067
The Journal of Physical Chemistry, Vol. 94, No. 18, 1990 7283
Energy Transfer in a Model Monolayer
"'i 10,
2 1
128 hard UFOs
0
,
,
o n , o - ,
,
0 00
0.2
04
0.6
,
1 08
PI4 Figure 2. Pressure-density isotherm for 128 UFOs;a / b = 6
Figure 4. Single-particle distribution, at P/p,kT = 0.20.
"em
Figure 5. Single-particle distribution, at P/p,kT = 4.70.
P/po
I
Figure 3. Orientational order parameters for 128 UFOs;a / b = 6 .
Results and Discussion Phase Behavior of UFOs on a Plane. The Monte Carlo calculated equilibrium densities and configurations are used directly to analyze the phase behavior of the system. The freezing transition from nematic to solid phases is clear in the graph of pressure versus density, Figure 2, and appears to be first order. The densities graphed are relative to close-packing density po. Close-packing density for the UFO system is derived in Appendix B. The UFO freezing is similar to the 5% density change that occurs when a monolayer of hard spheres freezes as P/pokT = 8.1 .j0 The nature of the transition from liquid to liquid crystal on a two-dimensional surface is an unsettled problem. Various researchers have found quite different results?J2 We have examined several features of the transition to a nematic phase in the UFO monolayer. The nematic order parameters llfi0ll and llf2211are plotted versus density in Figure 3. The parameter llfioll indicates p ordering. The UFOs rotate up out of the plane (Le., their Ps move to 3/2) at low densities, within the liquid phase. Because of this ordering, no completely isotropic liquid phase exists in this system, except at zero density. The parameter llfi211 is sensitive to a ordering and increases rapidly between p / p o = 0.4 and p / p o = 0.5. (The increase begins at lower density in the 32-particle system than in the 128-particle system.) This a ordering is due to alignment of the UFO long axes into a nematic state. No density discontinuity is apparent. If the a ordering is a thermodynamic phase transition, as opposed to simply an increase of order within a single phase, our data suggest that the transition is either continuous or weakly first order. The Single-Particle Distribution Function. The ordering of /3 and a is shown by the single-particle distribution functionflcos p,a-am). In the low-density reference state, cos p and (a- am,,,) are uniformly distributed. That is,f(cos b,a-anem) 1 / ( 4 ~ )as p 0. Figure 4 shows that p is ordered even at low density. At tends to 0 or T , as shown in Figure 5 . higher density a - anem The Reduced Pair Distribution Function C. For comparison with previous studies of related model systems, we calculated the
-
-
(30)Hoover, W . G.;Ree, F. H. J . Chem. Phys. 1968, 49, 3609-3617.
-
",/porT=Z74 F/p&T=078
Figure 6. Radial correlation function at three pressures.
100
80
0 0
5
10
15
20
Figure 7. Radial dependence of the reduced pair distribution function G , for two configurations, at P/p,kT = 4.70. T : 8, = p2 = */2; 0 - a, = 0, T ;a2,= a/2. Parallel: & = p2 = ~ / 2 0; - a, = 0, a;a2,= 0, ?r.
radial distribution function d2)(riI), which can be obtained by averaging eq 7 over angles. As shown in Figure 6 , the radial distribution function is sensitive to the surface pressure. The UFO radial distribution function lacks the sharp peaks seen for spheres at high density. It is difficult to compress UFOs into fixed lattice structures, even in the solid phase at p 7 0.7~0. Figure 7 shows G ( r i 1 )for "T" and "parallel" orientations. ( T has the parameters PI = Pi = ail = ~ / 2 6, - a I= 0; the parallel = 6 - a l = 0.) The radial parameters are PI = Pi = a/2,
He and Siders
7284 The Journal of Physical Chemistry, Vol. 94, No. 18, 1990
0
Figure 8. Average reduced pair distribution function ( G ) ,at P/pokT =
r/2;a21 = 6 -
Figure 11. Reduced pair distribution function G. 02 0; P/PokT = 2.74.
0.20.
Figure 9. Average reduced pair distribution function (G),at P/p,kT = 4.70.
Figure 12. Reduced pair distribution function G . CY, N
LV
f12
LV
p2 N r/2;6 -
N
r/2;6 -
0; P/p,kT = 0.20.
c-Zf
Figure 10. Reduced pair distribution function G . b2 a, u 0; P/p,kT = 0.20.
r/2;a2, = 6 -
distribution of prolate spherocylinders in T and parallel orientations ha5 been studied by others.5.6 The G functions for UFOs and spherocylinders in these two cases are similar. Both peak sharply near contact in the parallel orientation and have broader peaks at contact in the T orientation. Such similarity might be expected, in that a UFO behaves somewhat like a cylinder when the UFO "stands up" with p */2. Figures 8 and 9 show the dependence of { G ) on r21and 8 - a i . For these figures the dependence of G on PI, p2, and q1has been integrated out. Regardless of orientation, UFOs cannot approach closer than 2b, so G(r21 a ) . There are four spherical surfaces, and hence four spheres, to consider. Overlap is checked for by checking for intersection of each of the defining surfaces of one UFO (UFO A) with each of the defining surfaces of the other UFO (UFO B). Let the four surfaces be numbered 1,2,3, and 4, with surfaces 1 and 2 bounding UFO A and surfaces 3 and 4 bounding UFO B. Let Z,, 1 Ii I4, be the vector from the (arbitrary) origin to the center of the sphere on which surface i lies. The vector . if I I i A to the center of UFO A (B) is denoted FA ( i B )Clearly, I2b, then the UFOs must overlap, regardless of their orientations, and if I I i A > 20, then the UFOs do not overlap. Those simple situations are not considered further in this appendix. We begin by deriving an overlap criterion for surfaces 1 and 3. If no intersection is found between those two surfaces, then surfaces 1 and 4, then 2 and 3, and then 2 and 4, must be checked
-Energy -I ransrer in. a Moaei Monolayer . A
\
The Journal of Physical Chemistry, Vol. 94, No. 18, I990 7287
I . . .
I
for intersection. Define a local unit vector BZA along the symmetry axis of UFO A . Note that ?I = i A - ( R - b)OzA and 22 = i A + ( R - b)tZA. These vectors are sketched in Figure 18. A general vector 7 lies on surface 1 if and only if 117 - ?1112 = R2
(AI)
R - b I(7 - ZI)*DZA IR
(A2)
R - b I (i($) - Z3)ePz~I R
(AI 1) By inserting eq A10 into eq AI 1 and using the shift angle 6, it can be shown that the range of 4 that lies on surface 3 is eq A12.
and If the two surfaces intersect, their intersection must be part of the intersection of their defining spheres. If the two spheres do not intersect, then the surfaces do not intersect. Suppose the spheres do intersect: then their intersection is a circle. We parametrize the circle with the angle 4, measured from the plane containing Z3, and i A , (The case (Z3 - ?I)ll(iA - zl) is handled later in this appendix as special case 1, because the present procedure then fails to define pyA and CXA.)Construct a local coordinate system:
p, =
?3 - ?I - q1'
%A
(?A - ZI) x b, = - Il(?A - ZI) Ozll
-
The vectors and the angle 4 are sketched in Figure 19. In terms of 4, the intersection circle of the two spheres' surfaces is i(4) = (ZI Z3)/2 a(;, cos 4 + pY sin 4) (A4)
+
+
where uI
( R 2- IlZ3 - Z111'/4)'/'
('45)
Expression A4 may be substituted into condition A2 to determine the range of 4 over which the spheres' intersection circle and surface 1 of the UFO intersect. Substitution yields expression A6. R - b - llZ3 - ?lll(L?z'?ZA)/2 Icos 4 I PZAll
R - llZ3 - Zill(ez'8~~)/2 (A61 ullOz
;ZAIl
It can be shown that the right-hand limit is not less than unity. Moreover, if the surfaces overlap at any 4, then they overlap at 4 = 0. The upper limit can therefore be replaced with unity. Also, the lower limit cannot with physical significance be 2(R - b) and is null otherwise. Special case 3: PB X tc = 0 and CB.Cc = + 1. There is no overlap between surfaces 1 and 3 . Special case 4: ZB X = 0 and ZB.CC = -1. The J. range on surface 3 is ( 0 , 2 r ) if l1Z3 - Zlll > 2(R - b ) and is null otherwise. If the J. range is ( 0 , 2 r ) ,then we may of course set 6 = 0. Special case 5: 2A X 2, = 0 and CA.Cc = + I and P B X 0c = 0 and pB.bC = -1. Surfaces 1 and 3 overlap if and only if l1FB i A l l < 2b. Appendix B. UFO Close Packing
In order to calculate close-packing area, we chose the structure shown in Figure 20. In this structure /3 of the UFOs is r / 2 and a1- a2is 0 or r. Let the length of the semimajor axis of the UFO be u and that of the semiminor axis be b. Line BC is parallel to the major axes, and angle ACB equals r/2. Points A and B are the centers of the two spheres in contact. The spheres' radius is R = (a2 b 2 ) / ( 2 b ) .Length AB = 2R, AC = 2R - b, and ZZ = AB2 Close-packing density po = 1 / [ 2 b ( = ) ] = I / [2b(2a2 b2)1/2].If one sets u = 11 A and a / b = 6 , then po = 2.9 X 1O-Io mol/cm2.
+ z2. +
Appendix C. Calculation of AA at Small rZ1 To calculate the AA in eq 10, let a be the length of the semimajor axis and b the length of the semiminor axis and let Spro,
be the cross section of a UFO on the two-dimensional surface. At small pair separation, rZI 2a, the area available to particle 2 is actually less than [rzlArzlA(O-a1)], because particle 1 covers a portion [2r/A(O-al)]A,x,Iof the annulus. The effective area available to particle 2 is AA = [r21Ar21A(&al)] -A,,,,. The area AA is, of course, a function of r Z 1 .To calculate f i 2 ) and G, we use AA in place of [ r 2 1 A r 2 1 A ( O - ~Aexcl, l ) ] . for a particular radial box at rZl,is Spro,(rZl,P) divided by the number [2r/A(O-al)]of (0 - a l ) boxes in that annulus. Then 0
0 5 r215 26
r21A(B-al)Ar21- AeXcI 26 5 rZ15 a r21A(O-al)Ar21 a 5 rZ15 m
The statements apply only to the reference state that is orientationally uncorrelated. In the ideal-gas reference state we chose, two molecules do not interact even in the close-packed situation, except that they are governed by the hard-core potential.
J . Phys. Chem. 1990, 94, 7288-7293
7288
So what must be calculated is AA(rzl)in the region r2, c [2b,a]. We use Monte Carlo integration to calculate the irregular cross-section area Amd. The probability density of j3 in the isotropic ideal state is sin a. The averaged AA9 ( A A ), is easily computed as follows: ( A - 4 ) = r21A(6-al)Ar21- (A,,,,)
The area S,,. is easily calculated because of the simple UFO geometry. As shown in Figure 21, SProj always lies in the intersection of two identical circles. As shown in Figure 22, the distance llhlh211between the centers of these two circles is
llhlh211 =
a2 - b2 7 sin B
and the radius llhlGll of the circles on the two-dimensional surface is
((h,GII = [ ( a 2+ b2)’ - ( a 2- b2)2~os2 @I1/’/(2b) Appendix D. Average of k , with Respect to yI and y2
We use Forster’s rate constant averaged over rotations of UFOs about their symmetry axes. The rotation angles (also the third Euler angles) are y I and y2 for UFOs 1 and 2 . This appendix discusses explicitly the average presented in eq 19. As stated in eqs DI and D2,
(71.72)dependence of k, is related to the (yI.y2)average of K’ in eq D1. The orientation factor ( 2 )is calculated by transforming K’ from the coordinates of eq 17 (421,42,41) to Euler angles and integrating. The transformation and integration are done analytically. Define unit vectors PI and p2 in the symmetry planes of the UFOs. These represent the directions of transition dhole moments. In terms oiEuler angles cos a, cos p, cos y I - sin aI sin y, i = 1,2 (D3) p, = sin a, cos pi cos y,+ cos a, sin y I -sin Dl cos y ,
1
i
Let be the vector from the center of UFO 1 to the center of is denoted r 2 , , and T2, = (r21cos 8, r2, UFO 2. The norm of sin 6.0). Identities D4 allow ( K ~ to ) be expressed easily in terms of and the unit dipole moment vectors F2 and fi2. The result for ( K ’ ) is eq D5. cos
= pI-p2
cos
= pI+21/r21
cos
+2
= p2+21/r21
(D4) 3 2 ( ~= ~ 9) cos [4(8 - a l ) - 2(a2 - a , ) ] sin’ PI sin2 p2 6 cos [2(8 - a , ) - 2(a2 - .,)](I + cos’ PI) sin’ p2 6 cos [2(8 - a , ) - ( a 2- a , ) ] sin (2&) sin (2p2) 6 cos ( 2 8 - 2al) sin’ pl(l cos2 &) 2 cos (a2- aI)sin (2P,) sin (2P2) cos (2a2 - 2 a l ) sin’ PI sin’ p2 + 2[9 cos2 PI cos’ P2 + COS’ PI + COS’ /32 + 91 (D5)
+
+
We can compare the above result for ( K ~ )to Figure 4 of Gonen et al.” We set 8 - a , = 0, a2 = a I ,and P2 = 0,and plot ( K ’ ) versus PI in Figure 23. The least favorable orientation for energy transfer is seen to be PI = 58’. The (K’) plotted is averaged over y2 and yl, so it should lie between the two curves in Gonen’s Figure 4. Indeed it does.
A Simple Method To Estimate Free Energy from a Molecular Simulation: Renormalization on the Unit Interval B. Jayaramt and D. L. Beveridge* Department of Chemistry, Hall- Atwater Laboratories, Wesleyan University, Middletown, Connecticut 06457 (Received: January 3, 1990; In Final Form: April 4, 1990)
A theory to estimate the free energies of aqueous solutions expeditiously from a single conventional molecular simulation is proposed. Free energies are expressed as a function of the average internal energy and minimum energy by transforming the problem to the unit interval and evaluating much of the expression analytically. The theory is illustrated by applications to liquid water models and ionic solutes in water. The free energies of liquid water models (TIP4P, MCY-CI, and SPC/E) are in excellent accord with the results from other methods (average error is less than 3%). The free energies of ions (Li+, Na’, and CI-) in water are reasonably good (average error is less than 20%), considering the approximationinvolved. Applications to aqueous solutions of NaDNA in the presence of a simple salt indicate that the magnitude of the calculated free energies are consistent both qualitatively and quantitatively with polyelectrolyte theory for counterion condensation.
There is considerable current research interest in the numerical calculation of the free energy of chemical and biomolecular systems via molecular simu1ation.l Various methods such as thermodynamic integration, perturbation method, and potential of mean force calculations have been applied to this problem, all of which to date involve a series of simulations, each of which are
individually costly and time consuming. There is a critical need for a theoretically robust and computationally efficient algorithm for free energies in which the result can be obtained much more rapidly. We explore herein the possibility of obtaining accurate estimates of free energy from a single conventional simulation. The essence of the proposed procedure is a renormalization to the unit interval, whereby the configurational energy is trans-
Present address: Department of Chemistry, Indian Institute of Technology, Hauz Kaus, New Delhi, I10016 India.
1989, 18, 431.
Introduction
( I ) Beveridge, D. L.; DiCapua, F. M. Ann. Reu. Biophys. Biophys. Chem.
0022-3654/90/2094-7288%02.50/0 0 1990 American Chemical Society
