743
J. Phys. Chem. 1985,89,143-747
Simulation Study of Melting for Diatomic Molecules wlth Quadrupolar Interactions Vinayak Kabadi and W. A. Steele* Department of Chemistry, The Pennsylvania State University, University Park, Pennsylvania 16802 (Received: July 5, 1984; In Final Form: September 10, 1984)
Computer simulations for model homonuclear diatomic molecules are reported for a fixed density at temperatures ranging from the solid through the melting transition into the liquid. The effect of varying quadrupole moment upon the melting parameters was studied and it was found that large changes were produced when the quadrupole moment was increased from zero to a finite value. In addition to the thermodynamic properties, it is shown that structure-sensitive functions such as the pair correlation show striking changes upon melting. It is concluded that the temperature, the entropy, and the energy of melting are quite sensitive to the quadrupole moment.
Introduction In recent years, molecular dynamic computer simulations have been used extensively to study the configurational and dynamic properties of fluids and solids of diatomic molecules with and without quadrupole-quadrupole interactions. The system most studied has been nitrogen. Simulations of the equilibrium and dynamic properties of liquid nitrogen show good agreement with experimental data for thermodynamic properties, structure factors, and self-diffusion Numerous studies of the a,0, and y phases of solid nitrogenelo have also been reported. The general conclusion has been that the long-ranged quadrupolar interactions affect the properties of the solid phase much more than they do the liquid phase. Neither the equilibrium nor the dynamic properties of solid nitrogen are reproduced correctly by potentials including only atom-atom interactions between the molecules. In fact, the a-0 phase transitions are reproduced more accurately by a spherical Lennard-Jones potential plus ideal quadrupolar interactions9J0 than by the simple diatomic models. A few studies of the liquid phase of diatomic systems of longer molecules (reduced bond length L* = 0.425 to 0.793) interacting by atom-atom potentials with and without quadrupolar interactions”*l2have been published. To our knowledge, no computer simulations have been performed on the solid phase or the solid-liquid transitions of these systems. English and Venables13have evaluated the minimum potential energy solid phases of a number of diatomic systems, using an atom-atom potential with variable molecular length and quadrupolar interaction. The crystal cohesive energies for different possible lattice structures were very close in many cases so that comparisons with the lattice structures of real systems were not entirely successful. In this work, we have studied the solid phase and the solid-liquid transition of a system of diatomic molecules of fixed length (L* = 0.547) at constant density p* = 0.625. This length parameter is equal to that used by Streett and TildesleyI2 and corresponds approximately to that for a chlorine molecule. In addition to evaluating equilibrium properties including pair correlation functions, energies, and pressures as temperature was varied through the melting region, a study of the effect of quadrupolar interaction upon melting was undertaken. Thus, we compare
,*e
(1) Barojas, J.; Levesque, D.; Quentrec, B. Phys. Rev. A 1973, 7 , 1092. (2) Cheung, P. S. Y.; Powles, J. G. Mol. Phys., 1975, 30, 921. (3) Cheung, P. S. Y.; Powles, J. G. Mol. Phys. 1976, 32, 1383. (4) Weis, J. J.; Klein, M. L. J . Chem. Phys. 1975, 63, 2869. (5) Klein, M. L.; Weis, J. J. J. Chem. Phys. 1977, 67, 217. (6) Klein, M. L.; Levesque, D.; Weis, J. J. Phys. Reu. B 1980, 21, 5785. (7) Klein, M. L.; Levesque, D.; Weis, J. J., J. Chem. Phys. 1981, 74, 2566. (8) Quentrec, B. Phys. Rev. A 1975, 12, 282. (9) Mandell, M. J. J . Chem. Phys. 1974, 60, 1432. (10) Mandell, M. J. J . Chem. Phys. 1974, 60,4880. (11) Singer, K.; Taylor, A,; Singer, J. V. L. Mol. Phys. 1977, 33, 1757. (12) Streett, W. B.; Tildesley, D. J. Proc. R. SOC.London, Ser. A 1977, 43-
JJJ, LJ7.
(13) English, C. A.; Venables, J. A. Proc. R . SOC.London, Ser. A . 1974, 340, 57.
0022-3654/85/2089-0743$01.50/0
results obtained for nonquadrupolar molecules with those for = 1.00. molecules having quadrupole moment Q* (= Q/(CU~)~/’) Pronounced changes were observed, which is somewhat surprising in view of the fact that molecular shape and density were chosen to be the same in both sets of simulations.
The Potential and the Simulation Method In this, as in many other simulations, the potential energy of two homonuclear diatomic molecules is taken to be the sum of four interactions between pairs of atoms not on the same molecule:
where pij is the atom-atom separation. For homonuclear linear molecules, the quadrupole term is the first nonvanishing term in the expansion of the electrostatic potential. The energy of interaction of two point quadrupoles can be written14 3 Q’ - 5(C0S2Oi Lpa = -[1
4ri;
4
+ cos2 ej) + 2(sin ei sin Oj cos 4ij -
COS
8i COS
Oj)2
- 15 cos2 8i
COS’ Oj]
(2)
where rij is the separation distance and Oi, Oj and +ij are the polar and azimuthal angles which describe the relative orientations of two axially symmetric linear molecules. The simulation method used is nearly the same as that of Cheung and Powles2 and Streett and Tildesley.12 The only significant difference is that q ~ a t e r n i o n s ’ ~instead J ~ of Euler angles are used to represent the orientations of molecules. Simulations have been carried out with either 108 or 256 molecules in a cube with periodic boundary conditions. Runs were normally started from a face-centered-cubic (fcc) a-N2structure. In some cases around the phase transition region, runs were started from the final configurations at the next lower or higher temperature and, during the equilibration stage, the temperature was either increased or decreased gradually until the required value was reached. After equilibration, the simulations were carried out for 6000-1 2000 steps during which averages were taken to obtain equilibrium properties. Time steps in the range of 0.0002 to 0.0004 reduced units were used, smaller time steps being necessary in the phase transition region. A cutoff distance of 2 . 5 ~to 3u was used for the calculations of forces and torques. Corrections were made for the contributions to pressure and energy due to long-ranged interactions, assuming that there is no correlation between pairs of molecules separated by distances greater than the cutoff distance. (14) (15) (16) (17)
Buckingham, A. D. Q. Rev. Chem. SOC.Lond. 1959, 13, 183. Evans, D. J. Mol. Phys. 1977, 34, 317. Evans, D. J.; Murad, S.Mol. Phys. 1977, 34, 327. Steele, W. A. J . Chem. Phys. 1963, 39, 3197.
0 1985 American Chemical Society
744
The Journal of Physical Chemistry, Vol. 89, No. 5, 1985
Kabadi and Steele
TABLE I: Thermodynamic Properties Simulated for Diatomics with L * = 0.547, Q * = 1.00, p* = 0.625" 104 x At*
mean
mean
square
P*
E*
forced
square torqued
2
-6.80 -0.75 0.306 2.33
-20.42 -18.83 -18.53 -17.86
969.1 2999 3238 3792
62.93 166.3 180.5 203.1
8000 11000 7000
4 4 4
7.52 8.31 9.98
-15.79 -15.57 -15.11
4878 5305 6361
251.9 272.2 323.3
12000 10000 10000
3 3 3
11.33 12.83 14.38
-14.73 -14.39 -13.96
7382 8459 9674
368.4 421.4 480.6
T*
N,
0.595' 1.30 1.386 1.502'
3300 6000 10000 10000
6 4 2
1.598' 1.692' 1.918' 2.09 2.308 2.530
"N, is the number of time steps after equilibrations, and At* is the length of each time step in reduced units. 'Temperature very close to the transition temperature. 'These runs used 256 molecules. All the other runs carried out with 108 molecules. dReduced units are used: force X a/ 2.0) had all the properties of a liquid. For these systems, the rotational time scale is much smaller than the translational one so that orientational melting could occur much more rapidly than translational. Since the results for a couple of state points in the temperature range of the second phase showed a slow drift of energy and pressure, the simulations in the range 1S O < P < 2.0 were rerun with 256 molecules. The new runs were also performed over much longer times (6 to 10 reduced units) and a simulation a t a given temperature was started from the final equilibrium conf&yrations at the next lower or higher temperature. The results were different from those of the previous runs and indicated that the orientationally disordered "solid" phase was unstable relative to a liquid phase. We further checked the results at P = 1.598 by first starting from solid configurations and then reproducing the same results obtained starting from liquid configurations. Summaries of the results for the equilibrium thermodynamic properties are given in Tables I and I1 for Q* = 1.0 and 0.0, respectively. Mean square forces and torques were also computed and are listed.
-4
0.6
0.8
I
I
I
I
1.0
1.2
1.4
1.6
15
1"
Figure 2. Same as Figure 1, but for Q* = 0. The start of melting is less certain in this case because of a scarcity of results in this region, but the vertical line again indicates the estimated value.
In contrast to the usual constant pressure transition, melting at constant density occurs over a range of temperatures. Since the simulations are made at a relatively small number of state points, our estimates of melting P,p*, and thus of Ap* and M* are somewhat uncertain. Values for these parameters were taken from the plots of Figures 1 and 2. Of course, the melting P, p* for p* = 0.625 are those where the low-temperature curves first exhibit a break in slope;melting energy and pressure changes were estimated from the differences between the extrapolated highand low-temperature portions of the curves. Striking changes in melting temperature and energy are observed when a quadrupolar term is included in the interaction potential for this model diatomic molecule. Consequently, several other calculations were undertaken in an attempt to better understand this phenomenon. In addition to structural information such as that given by pair correlation functions, we calculated the entropy of melting hs* = M * / P = 0.4 and 1.4 for Q* = 0.0 and 1.0, respectively. Since the estimated uncertainty in these numbers is -lo%, it is evident
Simulation of Diatomic Molecule Melting I
I
I
I
I
o* = 0.0
I
I
o*= 1.0
o'=o.o
I
0" = 1.0
r"
Figure 4. Center-center (angle-average) correlation functions are shown; numbering system and Tz values are the same as for Figure 3.
-
*.
m
in
I
D
I
a*:
r*
Figure 3. Site-site correlation functions for a number of simulated systems are plotted here as a function of reduced separation distance r*. Three curves are given for Q* = 0, two for the liquid (dashed lines) at T* = 0.90 and 1.40 (numbered 2 and 3) and one for the solid (solid line) at Tz = 0.60 (number 1). For Q* = 1, the two liquid-statecurves are at Tz = 1.60 and 2.31 (6 and 7). Two solid-statecurves are also shown, for Tz = 0.60 and 1.30 (4 and 5). that the inclusion of quadrupole moment has a significant effect upon both AE* and M*of melting. Because both the dynamical and the structural results (see below) indicate that order in rotational as well as translational degrees of freedom is lost, we have estimated separately the quadrupolar and the site-site contributions to the total potential energy of these phases. One question to be answered is: can one ascribe the large increase in AE* when Q* is increased entirely to quadrupolar energy, or is the site-site portion of the total also affected (via changes in the structure)? The answer can be found in the entries of Table I11 which show that the melting change in the site-site part of the energy is almost three times as large for the quadrupolar molecules as for the nonquadrupolar. Also note that the quadrupolar energy for the liquid just above its melting point is far from the zero value expected for an orientationally disordered phase. In fact, only half of the average electrostatic energy is lost on melting. Not only does the quadrupolar interaction have a significant effect upon the structure of the solid phase near its melting point, the quadrupolar differences in the melting entropy as well as the component parts of the melting energy appear to be related to the structural alterations upon melting for the quadrupolar vs. nonquadrupolar systems. To substantiate this, we now show simulated structural functions. Structure As is well-known, there are several ways of presenting information concerning molecular pair correlations in dense phases. One can evaluate various point-point correlation functions such as g&), the site-site function (with sites taken to be coincident with the atoms of the molecules); g=(r), the center-of-symmetry
-0.5
1.0
-
-1.0-
-
1
. I
5
0
u r*
Figure 5. A low-order coefficient in the spherical harmonic expansion of the center-center pair correlation function is shown here for the same states as in Figure 3.
function; or g&), the site-to-center-of-symmetrycorrelation. Each of these may be viewed as an angle-averaged term in the expansion of the full intermolecular correlation function
where the Y/,,,(wi)are spherical harmonics with arguments equal to the orientation angles of the axis of molecule i. It is often convenient to take the centers of symmetry as origins,'' which would make gooo(rij) gE(r),but other choices are possible and, in some respects, preferable. Some of the results obtained for the pair correlation functions are shown in Figures 3-8. In order to better understand the significance of these curves, we should recollect that it is widely believed that the structural properties of simple liquids are primarily determined by the harshly varying repulsive portions of the intermolecular potential functions. Consequently, one anticipates that moderate changes in temperature at fixed density will not significantly affect the structure, nor will the introduction of multipolar electrostatic energies. Few systematic computer simulation studies of liquid-state structure for variable temperature at fixed density have been reported, nor have there been extensive structure simulations across the melting transition. In comparing the curves of Figures 3-8, especially the sets with different Q*,
746 The Journal of Physical Chemistry, Vol. 89, No. 5, 1985 I
Kabadi and Steele
I
-
-
*.
*-
I
0
0
0 t
N P
-80
2
I
Figure 6. A spherical harmonic coefficient of the pair correlation function having quadrupolar symmetry is shown here for the same states as in Figure 3. I
I
0':o.o
I
I
O*.l.O
-
*,--
N
P
I
I
I
2 r
Figure 7. A spherical harmonic coefficient exhibiting correlations in azimuthal angle is shown here for the states of Figure 3.
there is one complicating feature to keep in mind which arises from the lower melting point for Q* = 0; namely, most of the reduced temperatures for this case are nearly a factor of two smaller than those for Q* = 1. (The simulations of Streett and Tildesley12allow one to make comparisons at constant p*, T* but In order to facilitate comparison, solid-state simvariable ulations for Q* = 0 and for Q* = 1 were run at a single T* = 0.6. Angle-averaged pair correlation functions are shown in Figures 3 and 4 for the site-site and the center-center coordinate systems, respectively. The curves for the solid phase show the expected strong, long-ranged correlations. The peak positions and areas occur at the correct distances and yield the correct coordination numbers for the fcc a-N2 structure. After melting, much of the detailed structure disappears, as expected. Most of the differences between the solid-state gs,(r) for Q* = 0 and Q* = 1 may be ascribed to the difference in T*. The four liquid-state g&) shown for the two systems are interesting more for their insensitivity to T* and Q*than for their variation. This feature is also present in the work of Streett and Tildesley'2 who found that changing Q* at fixed T* = 2.20, p* = 0.524 gave almost no shift in the average site-site potential energy (AE,,* = 0.03 f 0.03, for Q*
e*.)
=()--re*=1).
I
I
2
I
3
I
I
I
I
2
3
r*
r*
I
I
Figure 8. A higher-order nonquadrupolar symmetry coefficient is shown for the states of Figure 3.
It is concluded that the repulsive part of the force is indeed playing a dominant role in determining these curves. Also, the similarity between the solid-state g,, curves at T* = 0.6 for the two Q* values leads one to the conclusion that the electrostatic energies have no strong effect on the structure of the solid. One of the more unsatisfactory features of this (and many other solid state simulation studies) is the great difficulty encountered in determining which of the many possible crystal packings corresponds to the most stable solid. It is extremely difficult (i.e., impractical) to simulate freezing of a sample of 256 molecules starting from a disordered liquid-state configuration, and indeed only a few such simulations have ever been reported in the lite r a t ~ r e . ' ~ ,Furthermore, '~ transition from the initial choice of crystal structure in a constant p, E simulation to a form with lower free energy is equally difficult. Even in the absence of such a proof, there are two reasons for believing that the a-N2 structure is the most stable solid for the molecules studied here. First, English and Venables13 have shown that this phase has the lowest potential energy of the various crystal structures that they treated, for the values of L* and Q*used in here. Secondly, we did show that an initial configuration of randomly oriented molecules on the fcc lattice relaxed into the a - N 2 structure at a temperature less than the melting point. Neither of these arguments is definitive and constant p,Tsimulations are now under way in an attempt to obtain additional confirmation for this conclusion. Returning to the pair correlation functions, we now consider for the behavior of the spherical harmonic coefficients gIltm(r*) the center-center origin. In contrast to the site-site coefficient, the glrm(r)show striking changes in the r* region corresponding to the first peak in gooo(r).Figures 4-8 indicate that the molecules penetrate much more in the liquid than in the solid, from a 1.2 in the solid down to r* = 1.O in minimum approach of r* the liquid. At the same time, the shapes of the "first peak" in the gll,,,,(r*)are significantly altered upon melting and, furthermore, the detailed nature of the alteration appears to depend somewhat upon the quadrupole moment. It is straightforward to calculate the positions and areas of the first peaks in the glF,(r) for the rigid a-N2 structure at p* = 0.625. One finds a nearest-neighbor distance r* = 1.3 13, which corresponds nicely to the maxima found in the first peaks in gooo(r*) for all systems studied. The spherical harmonics are also readily (18) Mandell, M.J.; McTague, J. P.; Rahman, A. J . Chem. Phys. 1976, 64, 3699.
(19) Mandell, M. J.; McTague, J. P.; Rahman, A. J. Chem. Phys. 1977, 66, 3070.
(20) See, for instance: Gray, C. G.; Gubbins, K. E., 'Theory of Molecular Liquids" Oxford University Press: London, 1984.
Simulation of Diatomic Molecule Melting evaluated. For example, consider Y2,0(Oi)(= (3/2 cos' Oi 1 / 2 ) ( 5 / 4 ~ ) 1 / 2 )where , Oi is the angle between the molecular axis (oriented along a diagonal of the cubic unit cell) and the intermolecular vector to one of the 12 nearest neighbors, for r* N 1.3. One quickly finds that g200(r*)vanishes when a sum over nearest neighbors is performed. Evidently, the small peaks found for the solids at r* near 1.3 are due to translation-libration displacements from the 0 K rigid structure. This argument is supported by the two Q* = 1 curves for the solid, which show that the first pair of small peaks decrease as the temperature is lowered. Furthermore, these peaks increase greatly on melting because of the great increase in librational disorder associated with this change. A comparison of the liquid-state curves indicates that peak height is more sensitive to Q* than to T* for these systems. For gzzl(r*),the sum of the spherical harmonics for the nearest neighbors in the solid vanishes because c#J~~,the difference in azimuthal angles, is 90° for all nearest-neighbor pairs. However, a large nonzero value is obtained for g220(r*)at the nearestneighbor distance. As the solids are melted, a dramatic change occurs in these spherical harmonic coefficients (and others not shown). In the first place, the penetration to shorter distances observed in gooo(r*)has the effect of producing a significant peak for 11" # OOO. Intuitively (and numerically), this in the gIpm(r*) is associated with configurations where both molecules are roughly perpendicular to the line of centers. For molecular pairs in this range of distance, one expects a negative peak in g200and positive peaks in gZz0and in g m ; gzzl would be expected to be zero for such configurations and we see only a small peak. At larger 1.3, the calculated gIrm indicate that a distances where r* significant remnant of the a-N2 solid-state structure persists in the liquid (another indication that this is the most stable solid for these molecules). Finally, note that the variations in peak height with T* and Q* at constant p* are small for the liquid glrm. In the case of the g22m.this observation is somewhat surprising since the quadrupolar energy has the correct symmetry to couple strongly with these components of the pair correlation function. Although it is hard to find simulation studies of molecular liquids at fmed p* and variable P,many of the phenomena noted here can also be detected in less systematic simulations reported previously.20 Furthermore, theoretical treatments such as RAM21,22have been proposed which are at least in qualitative agreement with our curves. The feature of this work which is most novel is the observed difference between solid- and liquid-state gIrm(r*) at the nearest-neighbor separations (in center-center coordinates), computed for systems where the site-site pair cor(21) For instance, see: Smith, W. R.; Nezbeda, I. In "Molecular-Based study of Fluids"; Haile, J. M., Mansoori, G. A., Eds.; American Chemical Society, Washington, DC, 1983; Adv. Chem. Ser. No. 204. (22) Nezbeda, I.; Smith, W. R.; Labik, S.J . Chem. Phys. 1983, 79,6242.
The Journal of Physical Chemistry, Vol. 89, No. 5, 1985 147 relations for the same range of separations hardly change on melting.
Conclusions In general, the use of periodic boundary conditions tends to stabilize the solid phase somewhat, so our estimates of melting temperature, etc. should be regarded as upper limits of the actual values for these model potentials. Note, however, that the 256particle system was definitely better behaved than the 108-particle case with respect to loss of translational order above the melting point. Furthermore, both systems exhibited reversible orientational disordering, beginning at the melting temperature for the 256particle simulation. That is, systems that were initially ordered translationally but disordered orientationally readily became ordered at temperatures below this melting point. We speculate that the orientational disordering process may couple with the translational disordering to give phase change properties that are less sensitive to system size than those for spherically symmetric potentials. Simulations of larger systems, plus studies of more nearly spherical potentials (Le., smaller L*), would be of interest in this respect. It is clear that melting in both the quadrupolar and the nonquadrupolar diatomic systems is accompanied by significant alterations in structure. The thermodynamic changes on melting, and particularly AE*,can be understood in terms of these structural alterations. The decrease in the height of the first peak in gss(r)shown in Figure 3 gives rise to a significant site-site contribution to AE*. For the quadrupolar systems, the decrease in the first (negative) peak in the gZh and the appearance of large 1.0 (shown in Figures 6 and 7 for m = 0 and m peaks at r* = f 1, respectively) give a significant quadrupolar contribution to AE*,and to Ap*, since the quadrupolar pressure is just equal to (5~*/3)E*,,~. Nevertheless, the quadrupolar parts of the energy and pressure are quite significant in the liquid just above the melting point and actually decay quite slowly with increasing temperature at fixed p * . Although a number of interesting theoretical papers dealing with melting of monatomic solids have appeared r e ~ e n t l y , 2there ~~~ appears to be no quantitative treatment of the melting of molecular solids. It is hoped that the results presented here will stimulate work on this problem. Acknowledgment. The authors gratefully acknowledge the support provided by the National Science Foundation for this work. (23) March, N . H. Mol. Phys. 1983,46,657. Weeks, J. D.; Broughton, J. Q.J. Chem. Phys. 1983, 78, 4197. (24) Bagchi, B.;Cerjan, C.; Rice, S.A. J. Chem. Phys. 1983, 79, 5595, 6222. (25) Ramakrishnan, T. V.; Yussouf, M. Phys. Reu. B 1979, 29, 2775. Ramakrishnan, T. V. Phys. Rev. Lei?. 1982, 48, 541. Haymet, A. D. J.; Oxtoby, D. J . Chem. Phys. 1981, 74, 2559.
