Molecular Dynamics Simulations of the Melting of a Hexane

The potential is still isotropic and is only a function of the separation of the ... The effect of the isotropy of the nonbonded potential on the melt...
0 downloads 0 Views 316KB Size
Download PDF
+

+

Langmuir 1996, 12, 1557-1565

1557

Molecular Dynamics Simulations of the Melting of a Hexane Monolayer: Isotropic versus Anisotropic Force Fields G. H. Peters Chemistry Department III, H. C. Ørsted Institute, University of Copenhagen, Universitetsparken 5, DK-2100 Copenhagen Ø, Denmark, and Novo-Nordisk A/S, Novo Alle´ 1, DK-2880 Bagsvaerd, Denmark

D. J. Tildesley* Chemistry Department, Imperial College of Science and Technology, London SW7 2AY, United Kingdom Received April 17, 1995. In Final Form: November 3, 1995X Molecular dynamics simulations of a hexane monolayer on the basal-plane surface of graphite have been performed to investigate the effect of an anisotropic force field on the melting process. In two sets of simulations, which are carried out for a number of temperatures ranging from the solid to the fluid state of the monolayer, the molecules are described either by a skeletal model, where the interaction sites are represented by a “united atom” model, or by a model where the interaction sites are shifted away from the mass site (the “anisotropic united atom” model). Independent of the model, the low temperature configuration is an orientationally ordered herringbone structure, which on heating undergoes an orientational phase transition to a rectangular-centered structure where the molecules tend to align along one direction. The solid subsequently melts at approximately 175 K. However, the anisotropic force field promotes a perpendicular orientation of the backbone of the hexane molecules and, in contrast to the “united atom” model, molecules exhibit a smaller tilt and a lower percentage of gauche defects at melting. This is compensated by increased librational motion in the plane of the substrate.

1. Introduction In recent years there has been considerable interest in the study of the structural and thermodynamic properties of adsorbed monolayers by computer simulation.1-3 When used in combination with diffraction and calorimetric measurements, molecular dynamics (MD) has proved to be a particularly powerful technique. Monolayers of rare gases, especially krypton, on graphite have been investigated extensively4,5 and these simulations have been extended to complicated adsorbates, such as small linear6,7 polar,8 and tetrahedral molecules.9,10 Recently, long-chain molecules, and alkanes in particular, have been the subject of considerable theoretical and experimental research,11-29 due in part to the widespread application of scanning X Abstract published in Advance ACS Abstracts, February 1, 1996.

(1) Steele, W. A. Surf. Sci. 1973, 36, 317. (2) Nicholson, D.; Parsonage, N. G. Computer Simulation and Statistical Mechanics of Adsorption; Academic Press: New York, 1982. (3) Taub, H. In NATO Advanced Study Institute, Ser. C; Long, G. J., Grandjean, F., Eds.; Luxer: Dordrect, 1988; Vol. 228, p 467. (4) Abraham, F. F.; Rudge, W. E.; Auerbach, D. J.; Kock, S. W. Phys. Rev. Lett. 1984, 52, 445. (5) Rowley, L. A.; Nicholson, D.; Parsonage, N. G. Mol. Phys. 1976, 31, 365. (6) Talbot, J.; Tildesley, D. J.; Steele, W. A. Mol. Phys. 1984, 51, 1331. (7) Joshi, Y. P.; Tildesley, D. J.; Ayres, J. S.; Thomas, R. K. Mol. Phys. 1988, 65, 991. (8) Carlos, Ruiz-Suarez, J.; Klein, M. L.; Moller, M. A.; Rowntree, P. A.; Scoles, G.; Xu, J. Phys. Rev. Lett. 1988, 61, 710. (9) Severin, E. S.; Tildesley, D. J. Mol. Phys. 1980, 41, 1401. (10) Moller, M. A.; Klein, M. L. J. Chem. Phys. 1989, 90, 1960. (11) Xia, T. K.; Landman, U. Science 1993, 261, 1310. (12) Peters, G. H.; Velasco, E. Mol. Phys., in press. (13) Velasco, E.; Peters, G. H. J. Chem. Phys., in press. (14) Hansen, F. Y.; Newton, J. C.; Taub, H. J. Chem. Phys. 1993, 98, 4128. (15) Hansen, F. Y.; Taub, H. Phys. Rev. Lett. 1992, 69, 652. (16) Leggetter, S.; Tildesley, D. J. Mol. Phys. 1989, 68, 519. (17) Crowell, A. D.; Steele, R. B. J. Chem. Phys. 1961, 34, 1347. (18) Lysek, M. J.; LaMadrid, M. A.; Day, P. K.; Goodstein, D. J. Phys. Rev. B 1993, 47, 7389.

tunneling microscopy and other proximal probe techniques to study these adsorbates. The structural and orientational ordering of these alkane molecules physisorbed from solution onto various substrates is a topic of considerable interest.21,22 There is a delicate balance between the energy associated with the flexibility of the chains, the molecule-surface energy, and the intermolecular interactions within the adsorbates which gives rise to a rich phase diagram and the appearance of new structures. From the simulation perspective, these calculations are numerically intensive and a number of simplifications have been introduced to produce trajectories long enough to extract significant structural and thermodynamic information. The most important of these simplifications has been the use of the “united atom” (UA) approximation, where the hydrogen atoms of the methyl and methylene groups are subsumed by the increased van der Waals radius of the carbon atoms creating an isotropic potential around each pseudoatom. For a linear alkane CnH2n+1, this approximation reduces the number of force sites to n and this significantly reduces the computation time which is of the order of n2 for generating the neighboring list and n for performing the sums over the short-range (19) Lal, M.; Spencer, D. J. Chem. Soc., Faraday Trans. 2 1974, 70, 910. (20) Osen, J. W.; Fain, Jr., S. C. Phys. Rev. B 1987, 36, 4074. (21) McGonigal, G. C.; Bernhardt, R. H.; Thomson, D. J. Appl. Phys. Lett. 1990, 57, 28. (22) Thibaudau, F.; Watel, G.; Cousty, J. Surf. Sci. Lett. 1993, 281, L303. (23) Newton, J. C.; Dennison, J. R.; Wang, S. K.; Wang, R.; Taub, H.; Courad, E.; Shechter, H. Bull. Am. Phys. Soc. 1987, 32, 467. (24) Newton, J. C. PhD Thesis, University of Missouri-Columbia, 1989. (25) Wang, R.; Danner, H. R.; Taub, H. In Ordering in two dimensions; Sinha, S. K., Ed.; Elsevier: New York, 1980; p 219. (26) Krim, J.; Suzanne, J.; Schechter, H.; Wang, R.; Taub, H. Surf. Sci. 1985, 162, 446. (27) Morishige, K. J. Chem. Phys. 1992, 97, 2084. (28) Shirazi, A. R. B.; Knorr, K. Mol. Phys. 1993, 78, 73. (29) Yeo, Y. H.; McGonigal, G. C.; Thomson, D. J. Langmuir 1993, 9, 649.

+

1558

+

Langmuir, Vol. 12, No. 6, 1996

intermolecular interactions. A second model, which represents a more accurate physical picture of the alkane chain, and at the same time is still computationally tractable, is the “anisotropic united atom model” (AUA).30-33 In contrast to the UA model, which does not reflect the tetrahedral symmetry of the methyl and methylene groups resulting from the sp3 hybridization, it models the hybridization by shifting the force site away from the carbon atom (or mass site) along the bisector of the HCH bond angle toward the hydrogen atoms. The potential is still isotropic and is only a function of the separation of the interaction sites but is anisotropic with respect to the nuclei of the carbon atoms. A full description of the application of this model for alkanes appears in ref 30. This model has been successfully applied in the computation of the equation of state for n-alkanes30,31 and polyethylene34 as well as in the simulations of Langmuir monolayers.35,36 Though recently37-39 it has been shown that the “united atom” model also predicts the equation of state of n-alkanes accurately at low pressures/densities, it appears that the anisotropic contributions are important in determining structural changes at high densities as observed in simulations of condensed Langmuir films, where packing effects predominately determine the phase transitions. In ref 35, it was observed that the characteristic features of phase transitions, i.e. the onset of translational diffusion and the formation of gauche defects, are less pronounced when using the UA model for the segments in the alkane chains. It has also been observed that the AUA approach yields a better description of the crystal structure of polyethylene.34 For such dense systems, the packing is sensitive to the parameters and functional forms of the potential and similar effects might be expected for the melting of adsorbed monolayers at full coverage. A limited number of studies of phase transitions for adsorbed flexible molecules, have appeared in the literature so far.10,12-16,20,25-29 The effect of the isotropy of the nonbonded potential on the melting behavior of these layers is still unknown. An excellent system for such a study is n-hexane adsorbed on graphite, which has been extensively investigated. Experimental work3,23-26 and computer simulations using the UA approach12-14 have revealed the existence of a herringbone structure at low temperature, with a uniaxial commensurate-incommensurate transition observed on increasing the coverage. On heating, the layer undergoes a first-order melting transition at ≈175 K. The paper is arranged as follows. In the following section we describe the potential models. Section 3 contains the details of the simulations. In Section 4, we discuss the results obtained in this simulation study and compare the properties of the two models. Our conclusions are presented in Section 5. 2. The Model The potential energy surface of the adsorbed monolayer is described by the sum of surface contributions originating from the graphite substrate and the interactions between (30) Toxvaerd, S. J. Chem. Phys. 1990, 93, 4290. (31) Padilla, P.; Toxvaerd, S. J. Chem. Phys. 1991, 94, 5650. (32) Reimers, J. R.; Watts, R. O. Mol. Phys. 1984, 52, 357. (33) Lustig, R.; Steele, W. A. Mol. Phys. 1988, 65, 475. (34) Pant, P. V. K.; Han, J.; Smith, G. D.; Boyd, R. H. J. Chem. Phys. 1993, 99, 597. (35) Karaborni, S.; Toxvaerd, S. J. Chem. Phys. 1992, 96, 5505. (36) Adolf, D. B.; Tildesley, D. J.; Pinches, M. R. S.; Kingdon, J. B.; Madden, T.; Clark, A. Langmuir 1995, 11, 237. (37) Siepmann, J. I.; Karaborni, S.; Smit, B. Nature 1993, 365, 330. (38) Siepmann, J. I.; Karaborni, S.; Klein, M. L. J. Phys. Chem. 1994, 98, 6675. (39) Smit, B.; Karaborni, S.; Siepmann, J. I. J. Chem. Phys. 1995, 102, 2126.

Peters and Tildesley

the long-chain molecules. The surface contribution is described by a static external potential, originally proposed by Steele,1 and is an expansion in terms of a rapidly convergent Fourier series over the reciprocal lattice vectors of the graphite structure. The surface potential contributions to each pseudoatoms, located at b r ) (x,y,z), in the hexane molecules is given by ∞

Usurf(r b) )

∞

∑R

Uo(zR) +

∞

∑R m)1 ∑ Um(zR)fm(x,y)

(1)

In this description, the substrate is divided up into a series of layers, R, each of which interacts with the atoms of the hexane molecules. The surface potential is a sum of two terms. The first term accounts for dispersion and overlap forces between the adsorbate atom and the layers, assumed uniform, whereas the second term represents the corrugation of the layers. The function fm(x,y) characterizes the symmetry of the substrate layers. For a detailed description of the potential and its parameters the reader is referred to ref 1. The sums appearing in eq 1 were truncated appropriately. For the first term the maximum R was 20, and as in earlier simulations,16,40 only one term in the sums over R and m in the second term were included in the computation of the surface corrugation; this approximation has been shown to be sufficiently accurate41 in the study of the adsorption of various molecules on graphite. The Lennard-Jones (LJ) interaction parameters for the carbon atoms of the substrate are ss/kB ) 28 K and σss ) 3.4 Å,1 where kB is the Boltzmann constant. The Lorentz-Berthelot combining rules were applied to determine the LJ parameters for the interactions between the carbon atoms of the alkyl chains and the graphite substrate.17 The model adopted for the hexane molecules consists of intramolecular and intermolecular terms. The intraand intermolecular potential functions are n-1

b) ) Uintra(r

∑ i)1

n-2

Ub +

∑ i)1

n-3

Uθ +

∑ i)1

n-4

Uφ +

N

b) ) Uinter(r

ULJ,intra ∑ i)1

(2)

N

ULJ,inter + ∑Usurf ∑ i)1 i)1

(3)

where Ub, Uθ, Uφ, ULJ, and Usurf are bond-length extension, bond-angle distortion, torsional, Lennard-Jones, and surface energies, respectively. The three-body potential, Uθ, is modeled as a harmonic function (eq 4)

Uθ ) (1/2)kθ(cos θ - cos θ0)2

(4)

where kθ is the vibrational force constant and θ0 is the equilibrium value of the bond angle. The numerical values are kθ ) 460 kJ/mol/rad2 and θ0 ) 111.9°.42 The four-body, torsional potential, Uφ, is modeled as a polynomial expansion (eq 5) in cos φ i)5

Uφ )

ai cosi φ ∑ i)0

(5)

φ in eq 5 is the dihedral angle. Uφ has minima at φ ) 0 and φ ) (2π/3 corresponding to the trans and gauche states, respectively. The coefficients, ai* ) ai/kB/103 are: (40) Hentschke, R.; Schu¨rmann, B. L.; Rabe, J. P. J. Chem. Phys. 1992, 96, 6213. (41) Steele, A.; Vernov, A. V.; Tildesley, D. J. Carbon 1987, 25, 7. (42) Egbert, E.; Berendsen, H. J. C. J. Chem. Phys. 1988, 89, 3718.

+

+

MD and Adsorption of Flexible Molecules

Langmuir, Vol. 12, No. 6, 1996 1559

a0* ) 1.03776, a1* ) 2.42607, a2* ) 0.08164, a3* ) -3.12946, a4* ) -0.16328, and a5* ) -0.25273.31 The intramolecular and intermolecular van der Waals interactions are described by a sum of LJ interactions between pairs of pseudoatoms

[( ) ( ) ]

ULJ(rij) ) 4

σ rij

12

-

σ rij

6

(6)

where σ and  are the collision diameter and the well depth of the potential, respectively. The potential is truncated at 3.332σ, and for the UA model the numerical values of σ and /kB are 3.923 Å and 72 K, respectively.43 For the AUA model, the values for the LJ parameters are σ ) 3.527 Å and /kB ) 120 K for the methyl groups and /kB ) 80 K for methylene groups.30 The interaction site of the anisotropic potential is shifted by 0.32 Å from the center of mass of the methyl and methylene units.30 Bond lengths are fixed at the equilibrium distance of the covalent C-C bond by applying constraints, using Toxvaerd’s method44 which has been described recently in ref 45. 3. Simulation Details The MD simulations were performed in a rectangular cell with periodic boundary conditions applied in the x and y directions, along the graphite plane. To avoid any loss of molecules by desorption (a rare event at the relatively low temperatures studied) a reflecting wall was situated at a distance of z ) 50 Å from the graphite substrate. The velocity of the center of mass of the molecules was reversed along the z direction, normal to the substrate, whenever a molecule crossed that plane. The simulation cell dimensions in the x and y directions had the symmetry of the basal plane of graphite and was chosen to be as square as possible to minimize the effects of the periodic boundary conditions. The size of the simulation cell is 4x(3)agm Å by 2agn Å, where ag ) 2.46 Å is the lattice spacing of the underlying graphite lattice. m and n are integers having the numerical values of m ) 4 and n ) 14. The x and y dimensions of the cell are 68.1735 and 68.88 Å, respectively. The simulation was started at a low temperature, where the hexane molecules completely cover the substrate forming a herringbone structure23,26 with a density of 42 Å2 molecule-1. A completely commensurate structure consistent with periodic boundary conditions was obtained by placing 112 hexane molecules in the simulation cell. A detailed sketch of the commensurate unit cell is given in Figure 3 of ref 26. All reported simulations were performed at a surface coverage corresponding to a full monolayer. The equations of motion were integrated using a leapfrog algorithm.46,47 The temperature was controlled with a Nose-Hoover thermostat,48 which couples the system to a heat bath at the required temperature. Simulations were performed at a series of temperatures, from T ) 75 to 250 K. Initially the molecules were placed in the commensurate herringbone lattice in the orientation suggested by elastic neutron diffraction experiments.26 The molecules were in the all-trans configuration with the molecular plane parallel to the substrate. The atoms were assigned random velocities taken from a Gaussian distribution at (43) Ryckaert, J. P.; Bellemans, A. Faraday Discuss. Chem. Soc. 1978, 66, 95. (44) Toxvaerd, S. J. Chem. Phys. 1987, 87, 6140. (45) Peters, G. H.; Toxvaerd, S.; Svendsen, A.; Olsen, O. H. J. Chem. Phys. 1994, 100, 5996. (46) Toxvaerd, S. Mol. Phys. 1994, 72, 159. (47) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1989. (48) Hoover, W. G. Molecular dynamics; Springer: Berlin, 1987.

the required temperature. The thermostat equilibrated the temperature within a few thousand steps. A typical run for a given temperature consisted of 50 000 equilibration steps followed by a production period of 200 000 MD steps, during which the thermodynamic and structural properties were averaged. The equilibration was monitored by measuring the atomic and molecular temperatures. When both temperatures oscillate around the same value, the kinetic energy is uniformly distributed among the various degrees of freedom, and the system is at equilibrium. The long sampling times, in excess of 200 ps, were required to ensure the accurate calculation of the orientational and conformational distributions close to the phase transition. A series of runs were performed by progressively increasing the temperature and using the final configuration of one run to begin the next simulation. 4. Results In this section we present results for the thermodynamic and structural properties of both models. Figure 1 shows snapshots of the monolayer at temperatures of 125, 138, 156, and 188 K, which are selected to highlight structural changes observed in the system during heating. Snapshots on the left and the right panels show configurations taken from the simulations with the UA or AUA model, respectively. Figure 1a shows a low-temperature configuration at 125 K. Molecules are arranged in a herringbone (HB) packing, which is in registry with the substrate. On heating, the layer undergoes an orientational transition to an incommensurate phase which also has a rectangular-centered unit cell with the molecules oriented along a common direction; this is the RC structure. Figure 1c (left-band side) at 156 K shows domains of the RC structure in coexistence with isotropic (I) fluid-like domains. These oriented domains show a shorter coherence length as the temperature is increased. The final configuration at 188 K, shown in Figure 1d, shows the fluid at a temperature just above the melting transition. For both models at all the temperatures studied, the promotion of molecules to the second layer was not significant. This behavior is confirmed by the density profiles calculated along the surface normal (not shown). We begin by defining a number of useful order parameters for this system. The commensurate order parameter defines the extent of the registry of the adsorbate with the surface N

1

OPcomm )

cos(ri‚ga)〉 ∑ i)1

〈

N

(7)

where ri is the position of the centre of mass of the ith hexane molecule in the plane of the surface, ga is the reciprocal lattice vector of the translational lattice of the adsorbate i.e.

g1 )

(

(

)

π ,0 3ag

-π π , 6ag x3a

g2 )

(8)

)

(9)

g

where ag is the size of the graphite lattice vector. The herringbone order parameter is defined as

OPherr )

N

1 N

sin(2φi)(-1)k〉 ∑ i)1

〈

(10)

+

1560

+

Langmuir, Vol. 12, No. 6, 1996

Peters and Tildesley

Figure 1. Snapshots of distinct configurations taken at (a) 125 K, (b) 138 K, (c) 156 K and (d) 188 K (from the top). The left and right panels show configurations obtained with the UA and AUA model, respectively.

+

+

MD and Adsorption of Flexible Molecules

Langmuir, Vol. 12, No. 6, 1996 1561

Figure 2. Commensurate order parameters, OPcomm, for the AUA (4) and UA (b) models; g1 (- - -), g2 (s).

Figure 3. Herringbone order parameters, OPherr, for the AUA (4) and UA (b) models.

where φi is the angle between the long axis of the hexane molecule, defined to be the eigenvector of the moment inertia tensor corresponding to the smallest eigenvalue, and the glideline of the herringbone structure. The glideline is a line of translation-reflection symmetry and it is clear from Figure 1a that in this case, the glideline points along the box-fixed y-axis. The factor (-1)k is a phase factor and k is an integer describing the row number. In a perfect herringbone structure φi ) π/4, and hence OPherr ) 1. If the molecules are oriented randomly in the surface plane then OPherr ) 0, in the thermodynamic limit. In a system of this size only one of the possible three degenerate herringbone domains is observed. The nematic order parameter, OPnem, is defined as

OPnem )

N

1

∑〈cos 2(φi - φdir)〉

(11)

Ni)1

where φdir is the angle between the director and the box fixed axis

[

φdir ) tan-1

N

/

N

sin 2φi ∑cos 2φi ∑ i)1 i)1

]

(12)

Considering the surface registry, both the UA and AUA models exhibit the 2 × 2x3(R16.1°) commensurate structure at 75 K. The order parameters for the UA models are higher than the those for the AUA model at all temperatures and the fluctuation in the order parameter is smaller for the UA model. The order parameters OPcomm(g1) and OPcomm(g2) are shown in Figure 2. Interestingly for the AUA model OPcomm(g1) is significantly higher than OPcomm(g2). One explanation is that the AUA model is producing a uniaxial phase, i.e. commensurate in only one direction. However it is clear from the error bars in Figure 2 that the OPcomm(g1) has a much higher fluctuation than OPcomm(g2). This is indicative of a larger fluctuation around the commensurate sites in the box-fixed ydirection. This fits in with the idea that in the AUA model the hexane molecules are oriented with their backbones perpendicular to the surface and that there is more room for an oscillation in the space-fixed y-direction. The commensurate order parameter falls sharply to zero between 151 and 156 K. It is not possible to make definitive statements about the order of the phase transition using such a small system size in a simulation

Figure 4. Nematic order parameters, OPnem, for the AUA (4) and UA (b) models.

but the behavior of the order parameter in Figure 2 is indicative of a first-order HB to RC transition. The HB order parameter is shown for both models as a function of temperature in Figure 3. OPherr is close to 0.8 at low temperatures, indicating that the stable phase has a herringbone structure in which the angle between the molecular axis and the glideline is approximately 60°. The order parameter only attains the value 1 in the perfect herringbone lattice when φi ) 45°. The order parameter falls to zero at the orientational phase transition, which occurs between 150 and 156.25 K. At the same time, the change in OPnem, shown in Figure 4, indicates a transition to the RC phase, where molecules point along a director. The low value of OPnem of 0.4 can be explained by the presence of fluid-like domains in the layer, which are in coexistence with the RC phase. At temperatures less than 150 K, the layer is ordered, but with an OPnem of -0.6, which does not correspond to a nematic phase but to the herringbone structure. At temperatures higher than 175 K the RC phase disappears resulting in a reduction of the order parameter to a value of order 1/x(N); this corresponds to an orientationally disordered fluid. The translational correlations can also be monitored by calculating the two-dimensional radial distribution function, defined in terms of the density correlations of molecules contained in a first slab of width ∆z ) 6.5 Å, centered at z ) 3.6 Å.

+

1562

+

Langmuir, Vol. 12, No. 6, 1996

Figure 5. (a) Radial distribution functions, g(r), at selected temperatures for the UA model. (b) Radial distribution functions, g(r), at selected temperatures for the AUA model.

A series of g(r)’s at different temperatures for the two models are shown in Figure 5. For both models, the layers form a solid-like structure up to a temperature of ≈ 150 K. Above this temperature, the radial distribution functions show only short-range translational order indicating the coexistence between solid-like structure and isotropic fluid. Above ≈175 K, the g(r)’s display fluid-like structure. For the UA model (Figure 5a), the positions of the first peaks in g(r) for the HB and RC phases reveal that the incommensurate RC phase is slightly expanded. Measurements of the unit cell from the configurations indicate

Peters and Tildesley

an expansion of ≈10% between the two phases. For the AUA model, the g(r)’s (Figure 5b) show significantly less solid structure. The broader peaks for the AUA model reflect the increased in-plane mobility of the adsorbate. The first peak in g(r) is located at larger separations than the corresponding peaks for the UA model. This is due to the anisotropic nature of the AUA approximation and is an indication that molecules orient with their molecular plane both parallel and perpendicular to the graphite surface. In both series of simulations, the RC structure is maintained over a narrow temperature range, and the monolayer melts between 175 and 181.3 K, where, for the AUA model, the peaks located at r/σ ≈ 2.8 and r/σ ≈ 3.3 are less pronounced. Similarly, for the UA model, peaks located at r/σ ≈ 3.6 and r/σ ≈ 4.8 vanish at ≈ 175 K. It is difficult to determine the precise phase diagram from the simulations at one coverage. The equilibrium RC phase may have a different density from the HB phase and the regions of disorder seen in the snapshots of individual configurations could correspond to dislocations resulting from a constant area simulation at too high a density. On the other hand, for Θ ) 1.0 at temperatures above 156 K, it is possible that the models are in a twophase region and we observe coexistence. The latter explanation is supported by experiment.14,15 In this case the melting is from a two-phase region of RC solid in coexistence with an isotropic liquid to a fluid phase. The peak at r/σ ) 2.3 becomes significantly smaller at the transition, and peaks at larger r disappear. The second and third peaks in the fluid g(r) reflect the persistence of a certain amount of orientational short-range order in the fluid. The orientational order of the molecules in the plane of the substrate was monitored by computing the distributions of the Euler angles φ and ψ. The former is defined as the angle between the long axes of the molecules and the x axis in the x-y plane, whereas the latter measures the orientation of the plane formed by the carbon backbone of the hexane molecules with respect to the substrate surface (roll-angle distribution). Figure 6 show the sequence of distributions, f(φ), over a range of temperatures for the UA and AUA model, respectively. Two peaks at 30° and 150° are present at low temperatures. These peaks are asymmetrical due to the anisotropic crystal field experienced by the molecules, packing effects favoring angles less than 30° and greater than 150°. On heating the asymmetry and width of the peaks increase and for the UA model the herringbone structure transforms to the RC structure at just above 150 K, where the peaks are located at 0° and 180°, respectively. However, precursor effects of this new phase are already observed at 150 K for the UA model and the AUA approach. This can be seen by the small tails in the distributions developing toward 0° and 180° for the UA model and the broadening of the peaks for the AUA approximation reflecting again the higher mobility of the molecules. The melting transition takes place between 175 and 181.3 K and the peaks disappear, giving rise to a uniform distribution representing a fluid monolayer. Further structural information can be deduced from the distribution of the roll angle, f(ψ), displayed in Figure 7a (UA) and Figure 7b (AUA). The distributions generally show four peaks at low temperature, located at 0°, 90°, 180°, and 270°. The angles at 0° and 180° are associated with the plane of the molecules parallel to the surface, whereas those at 90° and 270° correspond to a normal orientation (molecular plane perpendicular to the surface). For the UA model, the motion of the molecular plane of the molecules is complicated. At 75 K the molecules

+

MD and Adsorption of Flexible Molecules

Figure 6. (a) f(φ) distributions computed from the angle between the long axes of the molecules and the x-axis at different temperatures for the UA model. (b) f(φ) distributions computed from the angle between the long axes of the molecules and the x-axis at different temperatures for the AUA model.

slightly favor the parallel orientation. Snapshots of different configurations taken at this state point do not show any pattern or regularity in position of these normal or parallel orientations. With an increase in the temperature, the normal orientation is favored, as a result of packing effects. For the UA model as the temperature increases the molecules move away from the surface and are packed more tightly in the surface plane due to thermal expansion. Under this pressure, molecules roll out of the surface plane and peaks are seen in the distribution function at 90° and 270°. As the solid transforms into the RC phase there is more room on the surface and the

+

Langmuir, Vol. 12, No. 6, 1996 1563

Figure 7. (a) Roll-angle distributions, f(ψ), at various temperatures for the UA model. (b) Roll-angle distributions, f(ψ), at various temperatures for the AUA model.

molecules roll back into the plane. This phase transition takes place at approximately 154 K but the roll angle distribution at 150 K shown in Figure 7a shows this effect just below the transition. On the other hand, the AUA model (Figure 7b) favors the normal orientation of the molecular plane. Although, parallel orientations are detected at low temperatures, the number of molecules whose molecular planes are oriented normal to the surface is much higher than for the corresponding UA model. The structural distributions indicate that the film undergoes a complicated melting behavior, consisting of an orientational transition prior to melting. Experimentally phase transitions can be detected by anomalies in the heat capacity (CV ) dU/dT, where U is the internal energy). In the simulations, we defined the internal energy, ULJ*(T), as the sum of inter- and intramolecular

+

1564

+

Langmuir, Vol. 12, No. 6, 1996

Peters and Tildesley

Figure 8. Internal energies as a function of temperature for the AUA (4) and UA (b) models.

Figure 10. Gauche defects expressed in terms of percent trans for the AUA (4) and UA (b) models.

molecules have a higher mobility in the surface plane as indicated by the f(φ) distribution, Figure 6b. Similarly, the number of gauche defects, shown in Figure 10, is higher for the UA model than it is observed for the AUA approximation. It appears that the anisotropic nonbonded potential has the freedom to adsorb the increased kinetic energy by orienting the molecular plane of the molecules perpendicular to the surface and creating sufficient space in the substrate plane to increase the mobility. This feature is not observed in the UA approximation where the kinetic energy is adsorbed by creating a higher percentage of gauche defects and increasing the out of plane tilting. 5. Conclusions

Figure 9. Tilt order parameters, OPtilt, for the AUA (4) and UA (b) models.

LJ energies, since this yields the most precise location of the melting transition. The energy as a function of temperature is displayed in Figure 8, where the inflection point, associated with the melting point, is at around Tm ) 175 K. Within the statistical error, UA and AUA show the same melting temperature, which is in good agreement with recent experimental results15 based on structural measurements. The breadth of the simulated transition region is caused by the finite size of the system.49 Information about the degree of tilting of the hexane molecules can be deduced from the tilt order parameter, OPtilt, defined in terms of the second Legendre polynomial (P2).

OPtilt ) 〈P2(cos θi)〉 )

1

N

(3 cos2 θi - 1)〉 ∑ 2N i)1 〈

(13)

θi is angle between the surface normal, zˆ , and the long axis of molecule i. OPtilt as a function of temperature is shown in Figure 9. For temperatures higher than 120 K, the OPtilt increases monotonically with temperature, reflecting the tendency of the molecules to librate out of the surface plane due to thermal motion. However, for the AUA model, this effect is less pronounced and, instead, (49) Hill, T. L. Thermodynamics of small systems, Part II; Benjamin: New York, 1964.

We have performed a comparative study of the melting process of a hexane monolayer on graphite to elucidate the effect of the isotropy of the force field on the melting behavior. We have chosen a hexane monolayer, since this system has been studied extensively by Taub and coworkers.3,15 Earlier simulations conducted by Hansen et al.14 were in good agreement with experimental observed phases. These studies showed that at low temperature, the monolayer exhibits a commensurate herringbone structure which, on heating, melts in two stages. Firstly, at monolayer coverage the solid undergoes an orientational phase transition to form a rectangular centered phase in coexistence with an isotropic fluid. This solid is incommensurate and is expanded with respect to the herringbone structure. The coherence length of the nematic patches decreases and the solid finally melts to an isotropic liquid phase. In the present study, the methyl and methylene segments of the skeletal chains are modeled by using either the “united atom” or “anisotropic united atom” potential model. Both models give good agreement with the experimental melting temperature and show the experimentally observed phases. The thermodynamic properties (melting temperature) and the basic features of the phases are the same for both models. For the UA model the molecular plane can be parallel and perpendicular to the substrate plane and the molecules move between these orientation with increasing temperature to minimize the stearic repulsion within the monolayer. The AUA model promotes a perpendicular orientation of the molecular plane allowing significantly more in-plane mobility at all temperatures.

+

+

MD and Adsorption of Flexible Molecules

Thus in the case of the UA model, the melting transition involves not only positional and orientational changes but also changes in the distribution of the rotation around the long axis of the molecule. At low temperatures, the molecular planes are mostly parallel to the substrate, due to the larger specific area available to the molecules in this model. However, as the layer expands at higher temperatures, the melting for the UA model proceeds by creation of a larger number of gauche defects and additional tilting of the molecues out of the surface. In contrast for the AUA model the perpendicular orientation of molecular plane at low temperatures reduces steric overlap and increases in-plane mobility. This results in a broader orientational transition and the melting proceeds without an increase in conformational disorder or a significant change in the out of plane order. This situation is significantly different from that observed for smaller diatomic molecules on graphite, where at coverages of one monolayer there is significant second layer promotion on melting. Hexane is able to create enough space to melt in the first layer by orienting the molecular axis or by tilting the molecule out of the plane and conformationally disordering. One key result is that the UA model appears to give a more commensurate structure at lower temperatures, and there is a difference in the

Langmuir, Vol. 12, No. 6, 1996 1565

orientation of the molecular backbone and the overall tilt of the molecule in the commensurate phase. These differences could be explored by additional diffraction measurements which are capable of elucidating the precise orientation of the molecular backbone with respect to the surface. Unlike the three-dimensional alkane crystals, the differences between the UA and AUA models for the adsorbed hexane are small; both models give a good account of the experimental observations. The major difference between the monolayer and the bulk crystal is that on the surface molecules can tilt out of the surface plane to relieve the packing constraints. There is no such escape in the bulk crystal and the bulk properties appear to be more sensitive to the precise details of the intermolecular potential. Acknowledgment. G.H.P. thanks Novo Nordisk A/S for providing access to their computational resources to conduct this study and acknowledge support from an European Biotechnology Grant (No. BIOZ-CT93-5507). D.J.T. thanks the EPSRC for a grant (GR/J/74459) for computing equipment. LA950306X