Computer Simulation of the Gas/Liquid Surface BY GUSTAVO A. CHAPELA AND GRAHAM SAVILLE Department of Chemical Engineering and Chemical Technology, Imperial College, London, S.W.7 AND JOHN S. ROWLINSON Physical Chemistry Laboratory, South Parks Road, Oxford Received 13th January, 1975 The gas/liquid surface of a system of 255 Lennard-Jones (12,6) molecules has been simulated by Monte Carlo sequences at three reduced temperatures which span most of the liquid range. The results at the two higher temperatures show a monotonic profile of density as a function of height, but those at a temperature close to the triple point suggest that the liquid just below the surface can form layers; that is, that local density is an oscillatory function of height. This phenomenon is compared with similar previous results and it is concluded that it may well be the consequence of the constraints applied to the simulated system rather than an inherent property of a free liquid surface.Preliminary densities for a mixture have also been calculated. Molecular theories of surface tension and adsorption relate these macroscopic properties rigorously to the singlet and pair distribution functions of molecules and the variation of these functions across the gas/liquid interface. They provide, however, only approximate methods for the a priori calculation of the distribution functions, and one way of checking these approximations is to simulate the behaviour of a system by molecular dynamic or Monte Carlo techniques. Such methods have been used successfully to study the molecular distribution functions in bulk liquids for nearly twenty years but the study of surfaces is more recent and more difficult.The most obvious difficulty is the choice of a realistic constraint which fixes a plane surface in space without disturbing its shape or density profile. In the real world this constraint is provided by the gravitational potential, but this is ineffectual in a system containing 100-1000 molecules, which is as large as can be handled on a computer. Over a distance of ten molecular diameters the change in the earth’s gravitational potential is negligible compared with the intermolecular potential energy. The constraints that have been chosen for computer simulation have been either an artificially strong gravitational potential confined, perhaps, to one part of the system, or the application to the “ bottom ” of the cell of an external intermolecular potential which is chosen to represent, as realistically as is convenient, the field provided by the bulk liquid. An additional constraint, whose distorting effect is hard to estimate, is that imposed by the use of cyclic boundary conditions in the two directions, x and y , parallel to the surface.If the simulated cell has a width L in the x and y directions then this constraint requires that the surface be horizontal on a square grid of spacing L. It is known that a free liquid surface oscillates vertically to give surface waves of a length comparable with the wavelength of visible light (“ ripplons ”), and there is no reason to suppose that there are not shorter waves of 22G.A . CHAPELA, G. SAVILLE A N D J. S. ROWLENSON 23 all lengths down to that of the molecular spacing. The systems simulated on the computer suppress all oscillations whose wavelength is longer than L. In view of these difficulties it is, perhaps, not surprising that the simulations of similar systems, but with different constraints, have not always led to the same conclusions. In particular it is uncertain whether, or to what degree, there is a " layering " in the liquid just below the surface. This term is used to describe an oscillatory variation of density in the z, or vertical, direction with a period of about one molecular diameter. The results in this paper throw some light on this point, without settling it conclusively, and will be supplemented by further results, given at the meeting.MONTE CARL0 SIMULATION The present study is a Monte Carlo simulation of a system of 255 molecules between each pair of which there is a Lennard-Jones (1 2,6) potential u(r), characterised by the usual parameters of a collision diameter o [that is, u(a) = 01 and a depth E (that is, umi, = -8). The system is confined to a rectangular prism 50 square on its base and 250 in the vertical, or z, direction. Cyclic coordinates are used in the x and y directions and all potentials terminated at r = 2.50 [where u(r) = -0.016~1 so that a molecule cannot interact directly with its image. Potentials are set positive infinite for r < 0.80 to avoid computational difficulties. Below the base of the cell, z = 0, there is a simple representation of the bulk liquid obtained by assuming that each molecule occupies all positions at random.From this continuum representation there is emitted into the cell (z > 0) a Lennard-Jones (9,3) potential given by,l u(z) = p ~ ~ [ ( ~ ~ - ~ ~ ] (0.80 < z < 2.50) where the subscript w denotes " wall " and where ew/8 = (10/3)* O ~ / O = (2/15)*. The reduced density of the fluid, p, is (Na3/V). The top wall of the cell (z = 250) is similarly represented but, because of the lower density of the gas, this representation differs little from a plane infinitely repulsive wall at z = (25.0-0.8)~. Three Monte Carlo sequences have been generated, at reduced temperatures of z = kT/e of 0.701, 0.918 and 1.127. These temperatures are, respectively, close to that of the triple point, at a moderately high vapour pressure, and fairly close to the critical point.If argon is conventionally represented by a (12,6) potential with Elk = 119.8 K, then these temperatures correspond to 84, 110 and 135 K. The sequence at z = 0.701 was started at a randomly chosen configuration whose overall volume was chosen to match the expected liquid density. The first 1.5 x lo6 con- figurations were used to establish equilibrium, the approach to which was monitored by recording the total configuration energy of the fluid, the liquid and gas wall energies, and the density profile, p(z). The same properties and the surface tension were then calculated at each temperature during sequences of equilibrium states of from 3 to 6 x lo6 steps.The statistical average which gives the surface tension is where A is the area of the cell (250~) and the prime denotes differentiation. By symmetry, this sum vanishes in the isotropic fluids and so the whole contribution to y comes from the surface layer. To avoid the effects of the wall the sum was taken only over those molecules for whose centres z 2 60.24 COMPUTER SIMULATION OF THE GAS/LIQUID SURFACE A sequence for an equimolar mixture was generated at a temperature equivalent to 2, = 0.918 and Zb = 0.701 ; that is, for two molecules for which &,,/ebb = 0.763. The diameters were assumed to be equal, cr, = bab = bbb, and the cross-energy to be given by Berthelot’s rule &,b = (&,,Ebb)*. The walls at top and bottom were taken to be composed of a single substance whose properties were those calculated from the van der Waals 1-fluid appr~ximation.~ RESULTS The primary results are the density profiles, p(z), for each sequence. Those for the single component system are shown in fig.1-3. Near the bottom of the cell, z = 0, the density has a strong maximum caused by the packing of the first layer of 210 FIG. 1.-The density profile, pas a function of the height, z at a temperature of T = 0.701 for a r m of 5.6 x 106 configurations. z i o FIG. 2.-The density profile, p as a function of the height, z at a temperature of 7 = 0.918 for a run of 5.6 x 106 configurations.G . A. CHAPELA, G . SAVILLE AND J . S . ROWLINSON 25 molecules in the liquid against the flat repulsive wall potential of eqn (1).This peak is followed by one or two others which resemble the usual oscillations of g(r), the pair distribution functions in bulk liquid. At the two higher temperatures these oscillations f ~ ~ " " " ~ 1 1 ' 1 ' ' ' 5 10 15 4 0 FIG. 3.-The density profile, p as a function of the height, z at a temperature of T = 1.127 for a run of 3.5 x lo6 configurations. disappear beyond about z = 30 and the bulk liquid is then of constant, if slightly "noisy", density until the surface is reached. This is marked by a smooth and apparently symmetrical fall of density from that of the liquid to that of the gas. At z = 1.127 this fall is spread over about 5a, at z = 0.918 over about 40, and at z = 0.701 over about 3 0 . The density profile at the lowest temperature differs qualitatively from the other two in that it shows pronounced layers out to a distance of 100, a distance which would be presumed to be beyond the range of influence of the wall.The density of the bulk phases is readily found by averaging over the horizontal The discussion of these is deferred to the next section. 1 I I I 1 I I t 0.3 0.5 0.7 0.9 P FIG. 4.-The density of the bulk phases as found by simulation of systems of one&hase by Hansen and Verlet l 1 (triangles), and by simulation of systems with an interface by Liu lo (O), by Lee et aL9 (O), and in this work (0).26 COMPUTER SIMULATION OF THE GAS/LIQUID SURFACE (or oscillatory) portions of p(z). These densities are compared in fig. 4 with those obtained by others for a (12,6) fluid from simulations on systems both with and without a gas/liquid surface. The agreement is good.The surface tension is subject to greater statistical error and, as is shown in fig. 5, the agreement between different sets of computer results is satisfactory but not as good as for density. The first results for the mixture (1.2 x lo6 configurations) show that the molecule of lower E moves preferentially to the surface. The value of z at which p has fallen half-way from its value in the liquid to that in the gas is greater for the component of lower energy by about 0.50. These first results show no layering of the sum of the two densities, [p,(z) + pb(z)], but a degree of out-of-phase oscillation between them. Since, however, the run started from a randomly uniform local composition this partial layering may be a consequence of the initial movement of one component to the surface.Lattice models do not show this partial la~ering.~. A longer run is needed to resolve this point and it is hoped to present further evidence at the meeting. 2;o- 1:5- b4 * x $ 1.0- 0s - \ ”\ 0 \ 0 I I I I I I I I 0.6 0.8 I .o 1.2 1.4 7 FIG. 5.-The surface tension y as a function of temperature from the simulations of Liu lo (O), Lee et aL9 (0) and in this work (0). The line is the result of the perturbation calculations of Lee et al. THE EVIDENCE FOR LAYERS The first simulation to show strong layering was a molecular dynamic sequence with 200 molecules, made by Croxton and Ferrier.6 Since this was for a two- dimensional array it is uncertain evidence of the behaviour of real liquids, but Croxton ’ has since argued that liquid metals, in particular, might be expected to have a layered structure.A molecular dynamic run by Opitz * on 300 Lennard-Jones molecules in three dimensions at (apparently) z = 1.03 shows a bulk liquid which might have feeble layers of a spacing of about 2.50. However, the volume of the bulk liquid is not sufficient to show conclusively that the weak oscillation is not random. The evidence in fig. 3 of this work (T = 1.127) is that there are no oscillations of such long wave- length at these high temperatures.G, A. CHAPELA, G. SAVILLE AND J . S. ROWLINSON 27 There have been two recent Monte Carlo simulations of Lennard-Jones systems which should be closely comparable with the work reported here. Lee, Barker and Pound placed 1000 molecules in a rectangular prism whose base was 6.540 square and whose height was 400.The liquid was confined to the centre of the prism by a strong " gravitational '' potential which dissuaded the molecules from entering the parts of the prism near the ends, that is z < 20 and z > 380. At a temperature of z = 0.7035 they obtained two surfaces whose density profile at the surface resembled that of fig. 1. However below both surfaces the liquid was strongly layered; in each case there were about 8 maxima and minima in p(z), the separation of the maxima was close to CT and the peak-to-trough height was 0.2~-0.3~. Liu lo chose a system of 129 molecules without a gravitational field, and, as in this work, confined the liquid to one end of the prism by an intermolecular potential.His was less realistic, however, since it consisted of a square-well of arbitrary depth and width. He made runs at z = 0.75,0.90, 1.10 and 1.30, and claimed to find a " striking layered structure " at the two lower temperatures. It is, however, hard to see anything but random noise in his broad histograms, and the peak-to-trough heights of these are no more than 0.05~. The results of Optiz, of Liu and those shown in fig. 2 and 3 suggest that layering is not present at high reduced temperatures. If it is a true property of a gaslliquid surface, that is, one not induced by the constraints necessary for computer simulation, then it occurs only at low reduced temperatures. The three results at the lowest temperatures (Lee et al.z = 0.7035 ; Liu, z = 0.750 ; and fig. 1, z = 0.701) are all close to the triple point, which has been put at z = 0.68 i- 0.02 by Hansen and Verlet l 1 after extensive simulation of the bulk liquid and solid. In no case however were the x and y dimensions of the cell an integral multiple of the crystal lattice spacing, and, at least in our case, even the lowest layer of molecules appeared to be arranged without long-range order in the x-y plane. Moreover our liquid density, p = 0.81, is appreciably below Hansen and Verlet's estimate of 0.85+_0.01 for the freezing liquid. We conclude that the layering is not likely to be the onset of crystallization, but is a property of the constrained fluid systems under study. It is hard on this evidence to settle the question of whether or not the layers would persist in the absence of the constraints, but there are some indications that they would not.Thus our results and those of Lee et al. show oscillations which differ in amplitude by a factor of 2. We apply the constraint to the bulk liquid and the amplitudes decreases slowly as z increases. The last recognizable peak in fig. 1 is at z = 9.80 whilst the centre of the interface is at z = 12.80. Lee et ul. apply the constraint to the gas side of the interface and so ensure that their surface is even more closely confined to a plane than is imposed by the cyclic coordinates in the x and y directions. Their oscillations appear to start in the interface, to reach a maximum at a depth of about 40 below it and then to die away in the bulk liquid. The position and size of the oscillations appear therefore to depend on how the constraints are applied and, in particular, to depend on how strongly the surface is constrained to local planarity.It is possible that if all restriction on planarity could be lifted then the layers would disappear. This hypothesis could be tested by re- peating the runs in cells with substantially larger x and y dimensions. We are indebted to Dr. N. G. Parsonage for advice on the Monte Carlo simulation of surfaces, to Dr. J. A. Barker for an early copy of his paper, and G. A. C . thanks CONACYT (Mexico) for their support.28 COMPUTER SIMULATION OF THE CAS/LIQUID SURFACE ' D. M. Young and A. D. Crowell, Physical Adsorption of Gases (Butterworth, London, 1962), * F. P. Buff, 2. Elektruchem., 1952, 56, 3 1 1 . chap. 1. T. W. Leland, J. S. Rowlinson and G. A. Sather, Trans Faruduy Suc., 1968,64, 1447. S . Ono and S. Kondo, Handbuch der Physik (Springer, Berlin, 1960), vol. 10, p. 262 et seq. J. E. Lane and C. H. T. Johnson, Austral. J. Chem., 1967,20, 61 1. C. A. Croxton and R. P. Ferrier, J. Php. C, 1971, 4, 2447. ' C. A. Croxton, Liquid State Physics (Cambridge U.P., 1974), chap. 4. * A. C. L. Opitz, Plzys. Letters, 1974, 474 439. J. K. Lee, J. A. Barker and G. M. Pound, J. Chem. Phys., 1974, 60, 1976. J. P. Hansen and L. Verlet, Phys. Rev., 1469, 684, 151. lo K. S. Liu, J. Chem. Phys., 1974, 60,4226.