Computer Simulation of Gas-Liquid Surfaces : Molecular Fluids BY STEPHEN M. THOMPSON School of Chemical Engineering, Cornell University, Ithaca, New York 14853, U.S.A. Received 5th May, 1978 The molecular dynamics method is applied to the coriiputer simulation of the gas-liquid surface of molecular fluids. Two fluids (nitrogen and chlorine) are modelled with a site-site potential at three temperatures for each system. Surface tensions are in good agreement with experiment and show a significant size effect, while density against orientation profiles show orientational ordering in the surface zone at low temperatures. 1. INTRODUCTION Here we describe the results of molecular dynamics (MD) computer simulations of the gas-liquid surface of two molecular fluids, modelled on nitrogen and chlorine, at three temperatures from the triple point region upwards. Results are presented here for the density against orientation profiles and surface tensions, while results for selected non-equilibrium properties will appear elsewhere at a later date.Computer simulation was chosen as the tool used to obtain these results for two main reasons. First, a number of properties of surfaces are difficult to obtain experimentally, for example, the density against orientation profiles and pair distribution functions. Secondly the theoretical calculation of other surface properties for which defined experi- mental techniques exist, such as the surface tension, is very difficult unless one makes drastic simplying assumptions, as the theory usually involves experimentally inaccessible quantities.All these properties are obtainable in principle by performing computer simulation experiments. So far as we are aware, this is the first work of this type on molecular fluids,l although a body of literature exists on monatomic fluids.2-8 Calculation of the density against orientation profiles enables one to study the preferred orientations (if any) in the surface zone; the effect of molecular shape on surface tension can also be investigated. 2. MOLECULAR MODELS Each fluid was modelled by means of a site-site or atom-atom p~tential.~ Each molecule is assumed to consist of two interaction centres (commonly assumed posi- tionally coincident with the atomic nuclei). The intermolecular potential is then given by 2 2 ~ ( r l 2 ~ 0 1 ~ 0 2 ) = 2 2 uLJS (IY~") - rj(2)l).(2.1) i = l j = 1 Where ri(") is the position vector of site i on molecule k, rI2 is the centre to centre108 COMPUTER SIMULATION OF GAS-LIQUID SURFACES separation, cok = {Ok& are the polar angles specifying the orientation of molecule k, and ~Ljs(r) is a shifted Lennard-Jones (12, 6) interaction lo where uL,(r) is a Lennard-Jones (12,6) interaction uLJ(r) = 4&[(O/r)” - (O/r)6] (2.3) and Y, is the potential cut-off distance. Here E has the dimensions of energy and 0 TABLE 1 .-PAIR POTENTIAL PARAMETERS reference for fluid (&lk)lK alnm Lla parameters NZ 37.3 0 . 3 3 1 0 0.3292 11 Cl* 1 7 3 . 5 0 . 3 3 5 3 0.6080 1 2 has the dimensions of length. The values used are given in table 1. The cut-off distance was given in each case by rc = 2.50 + L.The second term in eqn (2.2) is the usual shifting that results from using a trun- cated p ~ t e n t i a l , ~ ~ ~ while the final term is introduced to ensure that both the potential and its derivatives are continuous at the cut-off. This leads to better dynamics (see below). The results are reduced by means of the parameters E~ and oR, which are defined as the well depth and collision diameter respectively of the spherically averaged pair potential : (2.4) <u(aR, co1c02))wlwz = 0 (2.5) <u(YR, co1c02))wlw2 == -&R (2.6) where rR is the location of the minimum, This enables the results to be directly compared. Values of the parameters These values were obtained by numerical integra- cR, O, and rR are given in table 2.tion. TABLE 2.-sPHERICALLY AVERAGED POTENTIAL PARAMETERS N 2 l3 97.04 0.3754 0 . 4 1 57 0.3292 ClZ 2 3 3 . 3 9 0.4568 0.4984 0.6080 Currently under investigation is the addition of point charges of magnitudes q to the sites and -2q to the centres of mass, where q is chosen to reproduce the known quad- rupole moments of N2 and Cl,.S . M. THOMPSON 109 3. SIMULATION Each system consists of 216 molecules confined to square prisms of dimensions x,, yL (=xL) and z,( > x,) with the usual periodic boundary conditions in the (xy) plane of the surface. The liquid was confined to the centre of the cell with vapour phases at either end. Thus an infinite film is simulated, with two surfaces that are stable without the inclusion of external force^.^-^*^ The cell widths were xL = 2r, (5.65840 and 6.2160 for N2 and C1, respectively) and the cell heights were 190 for both systems at the lower two temperatures and 220 at the highest temperature.The translational and angular velocities of the molecules, ui and mi, respectively, are chosen randomly subject to the constraints that the desired kinetic energy is obtained, this being properly divided between translations and rotations, and that the centre of mass of the system has no resultant translational or angular velocity. At the lowest tem- peratures the molecular centres of mass are assigned coordinates on an f.c.c. lattice structure of the appropriate liquid density, with empty vapour phases. An initial equilibration period of typically lo4 integration steps is rejected, during which the lattice structure melts and the density profile becomes stable.The runs at the higher temperatures are started from the last configuration of the run at the next lowest temperature. A vapour molecule which attempts to leave the cell through one of the end walls is returned to the cell by simply bouncing it off the wall. This procedure is thought to have negligible effect on the surface structure. The position of the centre of mass was monitored during the run, and moved on average a distance of 0.005~ per 10 000 steps. Initially the centre of mass is fixed in the centre (z,,, = zJ2) of the cell. The centre of mass motion was found by solving the equations Fi = M-'Fi (3.1) using a third-order predictor-corrector method. Here ri is the centre of mass vector of molecule i, M is the molecular mass, H i is the force on i due to all the other mole- cules, and the dots indicate differentiation with respect to time.The rotational motion was found by solving the equations hi, = I-'TiP (3.2) using a second order predictor-corrector method, where miP and Ti, are the principal angular velocity of and torque on molecule i respectively, and Zis its moment of inertia. Molecular orientations were expressed by means of quaternion parameter^,'^ giving singularity-free equations of rotational motion with better resultant energy conserva- tion. The equations of motion were integrated with a reduced time step of 0.0015 [ t * = t ( ~ M / a ' > ~ ] . s for N2 and C12, respectively. Total energy was found to oscillate around a fixed value with an amplitude of -0.25% an1 a period of 3000-5000 time steps.The use of a truncated but not shifted pair potential changes this behaviour to a slow monotonic drift. The use of the same algorithm to simulate a bulk phase gave total energy values exact to five significant figures over 2000 steps, thus showing the inhomogeneity of the system to be a cause of numerical instability. The calculations were performed on a PDP 11/70 computer. The computer programs were overlaid and occupied a maximum of 28.5 K 16-bit words of core. Various time-saving schemes were implemented. A run of 50 000 integration steps took 100-170 h, depending on the temperature. The method used to solve eqn (3.2) has been described previou~ly.'~ In real units this is 4.39 x and 3.53 x1 1 0 COMPUTER SIMULATION OF GAS-LIQUID SURFACES 4.PROPERTIES CALCULATED The density against orientation profile p(z, 0) gives the probability density for finding a molecule at height z with the molecular axis at an angle 8 to the z axis.15 It may be expanded as a sum of Legendre polynomials in cos 8: For symmetrical linear molecules pl(z) vanishes when I is odd. The total density profile (centres of mass) is pCOM(z) = [ p(z, 8) sin 8 d8. (4.2) Thus where dl0 is a Kronecker delta. The individual coefficients can be calculated during the simulation from the formulae 3.4) s = l i where the sum over s is over a sequence of Ns time steps, the sum over i is over those molecules with coordinates in the range z - Az < zf < z + Az, and A is the surface area.A value of Az = 0.050 was used in this work. The coefficients I = 0, 2 and 4 were calculated. The surface tension is calculated during the simulation by means of the equation where i and j are molecular subscripts and the r, depend implicitly on i and j . Note that the calculation applies for a truncated Lennard-Jones diatomic fluid, with the exception that the resulting shift in the coexisting densitiesS is not allowed for (see below). If the sum over i a n d j in eqn (4.5) is taken over all molecular pairs, then the result must be divided by two because there are two free surfaces. A correction for the missing tail of the potential was added in a similar manner to that explained previously.8 5. RESULTS The averages obtained during the simulations are given in table 3.The coexisting densities are plotted along with experimental values in fig. 1. Values obtained for a monatomic (1 2, 6) fluid modelling argon are plotted for comparison. Simulation values are low as has been observed previously:s this is due to the shift in the pair potential. While the values for N2 are similar to those for Ar, being only slightly higher, values for C12 are much higher, showing a large size effect. Total density profiles are plotted in fig. 2 and 3 for N2 at 66.5 K and C12 at 171.1 K, respectively. These are averages for the two surfaces in the box in each case; the profiles for the two surfaces were in almost exact agreement with each other. These two profiles are similar in their major details : the " structure " observable in the liquid region changed in form during the simulation and is not considered significant.The estimated error in liquid density in table 3 is based on the size of these oscillations.S. M . THOMPSON 111 Apart from the size of these oscillations being much smaller at higher temperatures, the other profiles are very similar to the ones shown here and are therefore not plotted. Here it is important to note that the real-time duration of the present simulations is less TABLE 3 .-AVERAGES OBTAINED DURING SIMULATIONS Na 42.75 x lo3 66.5 0.887 0.971 1.26 1.20 1.58 1.876 NI 50.0 x 103 77.2 0.803 0.918 0.89 0.94 1.68 2.194 N1 40.0 x lo3 96.6 0.649 0.811 0.46 0.48 3.36 1.755 CI, 50.0 X lo3 171.1 1.335 1.398 2.43 2.54 0.81 1.763 CI, 50.0 x lo3 195.8 1.272 1.350 2.30 2.24 0.99 1.763 C11 50.0 x lo3 218.8 1.211 1.306 2.09 1.95 1.10 1.763 52.0 zt0.035 50.07 h2.2 f0.025 50.04 52.6 zt0.017 50.04 k5.2 h0.066 *0.10 *5.9 f0.036 *0.10 56.9 50.029 f0.05 a Number of time steps used in run proper.Temperature. Reduced liquid density from simulation. * Reduced I Reduced surface tension from experiment. liquid density from experiment. g Reduced surface thickness. * Real time duration of simulation. Reduced surface tension from simulation. 0 a El 0 I . 0.5 0 . 7 '.* kT/ER O.' 1 . o Fig. 1 .-Coexisting liquid densities as a function of reduced temperature (a) Ar - 0 , sirn"; expt; (b) N2 - x , sim"; 0, expt; (c) C12 - x, sim"; 0 expt. than half that of those for Ar in ref. (8) and hence larger statistical error must be ex- pected. A surface thickness may be defined in terms of the gradient of the total density profile at the Gibbs dividing surface.8 Values obtained in this manner are112 COMPUTER SIMULATION OF GAS-LIQUID SURFACES 1 .0 - 0.75 "b" - h, - 5, 0.50- < 0.25 0.0 shown in fig. 4 together with argon values for comparison.* Values for N2 and Ar are very similar, but with an apparently more rapid divergence for the former at - - + 1 higher temperature. Thicknesses for C12 are smaller, which is reflected in the surface tension values (see below). In fig. 5 is plotted the second density harmonic coefficient p2(z) for C12 at 171 K. I /a FIG. 3.-Total density profile for C12 at 171.1 K. The density orientation profile constructed from this coefficient and po(z) is plotted in fig. 6. For N2 there was no evidence of any structure in these curves at any of the temperatures studied although the statistical noise was rather large.This was trueS . M. THOMPSON 4.0 3 . 0 - F -m 2.0 1 . 0 - 0.0 113 - - 0 O x x o I I L 0.7 0.8 0.9 1 .o 8 I * x 0 B FIG. 5.-Second density harmonic coefficient p2(z) for C12 at 171.1 K. The arrow locates the Gibbs dividing surface.114 COMPUTER SIMULATION OF GAS-LIQUID SURFACES 0.0 FIG. I I L 6.-Density against orientation profile po(z) + p2(z)P2(c0s 0) for C12 at 171.1 K. correspond to the asterisks in fig. 5 . to locates the Gibbs surface. The heights P I IS. M . THOMPSON 115 also for the p2(z) coefficients for C12 at 218.8 K and for the p4(z) coefficients for both systems at all temperatures. The p2(z) coefficient for C1, at 195.8 K was similar to that presented here, although of reduced amplitude.It also appeared to be shifted in towards the liquid phase relative to the Gibbs surface. The non-zero values in the bulk liquid region in fig. 5 are statistical noise : here the coefficient differed substanti- ally in both box halves and changed in form as the simulation progressed, while this was not so for the structure in the surface region. A pronounced (height dependent) tendency to adopt preferred orientations is indicated. In the liquid phase at the posi- tive maximum in p2(z), the molecules have a tendency to orient with their axes per- pendicular to the plane of the interface (@ = 0), while at the negative minimum (0.050 on the vapour side of the Gibbs surface) the molecules tend to orient with their axes parallel to the interfacial plane.These results are in qualitative agreement with the predictions of first order per- turbation theory using a spherical reference potential with anisotropic overlap forces.15 A quantiative comparison of simulation and perturbation theory for the potential employed in the simulation is in progress. We find little tendency for the near-spherical molecule N2 to adopt preferred orientations in the surface zone, while this tendency is pronounced for C1, at low temperatures. The surface tensions obtained are plotted in fig. 7 together with those of argon* for comparison. Agreement with experiment is good, although this is probably fortuitous in view of the coexisting densities obtained.Little effect of deviation from sphericity is shown by N,, while again the tendency is pronounced for Cl,. This is intuitively obvious. Probably the inclusion of multipolar forces will significantly alter the surface tension. The author thanks the S.R.C., the National Science Foundation and the Petroleum Research Fund (administered by the American Chemical Society) for support of this work. He also thanks Professors K. E. Gubbins and J. S. Rowlinson and Mr. S. Murad for many discussions. S. M. Thompson and K. E. Gubbins, presented at the Symposium on Computer Simulation of Bulk Matter from a Molecular Perspective, American Chemical Society Meeting, Anaheim, California, March 14, 1978. To appear in the ACS Symposium Series. G. A. Chapela, G. Saville and J. S. Rowlinson, Faraday Disc. Chem. Soc., 1975,59, 22. J. K. Lee, J. A. Barker and G. M. Pound, J . Chem. Phys., 1974,60, 1976. F. F. Abraham, D. E. Schreiber and J. A. Barker, J . Chem. Phys., 1975, 62, 1958. A. C. L. Opitz, Phys. Letters A , 1974, 47, 439. K. S. Liu, J . Chem. Phys., 1974, 60, 4226. M. Rao and D. Levesque, J . Chem. Phys., 1976, 65, 3233. 1133. J. R. Sweet and W. A. Steele, J . Chem. Phys., 1967, 47, 3029. lo W. B. Streett, D. J. Tildesley and G. Saville, as ref. (1). l1 P. S. Y. Cheung and J. G. Powles, Mol. Phys., 1975, 30, 921. l2 K. Singer, A. Taylor and J. V. L. Singer, Mol. Phys., 1977, 33, 1757. l3 S. M. Thompson, D. J. Tildesley and W. B. Streett, Mul. Phys., 1976, 32, 711. l4 D. J. Evans and S. Murad, Mol. Phys., 1977, 34, 327. l5 J. M. Haile, K. E. Gubbins and C. G. Gray, J . Chem. Phys., 1976, 64, 1852. l6 R. C. Reid, J. M. Prausnitz and T. K. Sherwood, The Properties of Gases and Liquids (McGraw- l7 J. J. Jasper, J . Phys. Chem. Ref. Data, 1972, 1, 841. * G. A. Chapela, G. Saville, S. M. Thompson and J. S. Rowlinson, J.C.S. Faraduy ZZ, 1977, 73, Hill, New York, 3rd edn, 1977).