3228
Langmuir 2008, 24, 3228-3234
Explicit Hydrogen Molecular Dynamics Simulations of Hexane Deposited onto Graphite at Various Coverages M. J. Connolly,† M. W. Roth,*,† Paul A. Gray,‡ and Carlos Wexler§ UniVersity of Northern Iowa, Department of Physics, Cedar Falls, Iowa 50614, UniVersity of Northern Iowa, Department of Computer Science, Cedar Falls, Iowa 50614, and UniVersity of Missouri, Department of Physics and Astronomy, Columbia, Missouri 65211 ReceiVed October 2, 2007. In Final Form: December 21, 2007 We present results of molecular dynamics (MD) computer simulations of hexane (C6H14) adlayers physisorbed onto a graphite substrate for coverages in the range 0.5 e F e 1 monolayers. The hexane molecules are simulated with explicit hydrogens, and the graphite substrate is modeled as an all-atom structure having six graphene layers. At coverages above about F = 0.9 the low-temperature herringbone solid loses its orientational order at T1 ) 140 ( 3 K. At F ) 0.878, the system presents vacancy patches and T1 decreases to ca. 100 K. As coverage decreases further, the vacancy patches become larger and by F ) 0.614 the solid is a connected network of randomly oriented islands and there is no global herringbone order-disorder transition. In all cases we observe a weak nematic mespohase. The melting temperature for our explicit-hydrogen model is T2 ) 160 ( 3 K and falls to ca. 145 K by F ) 0.614 (somewhat lower than seen in experiment). The dynamics seen in the fully atomistic model agree well with experiment, as the molecules remain overall flat on the substrate in the solid phase and do not show anomalous tilting behavior at any phase transition observed in earlier simulations in the unified atom (UA) approximation. Energetics and structural parameters also are more reasonable and, collectively, the results from the simulations in this work demonstrate that the explicit-hydrogen model of hexane is substantially more realistic than the UA approximation.
I. Introduction Because of its utility, stability, and geometry, much experimental and theoretical work has been completed on systems involving graphite.1,2 Hexane on graphite has been studied experimentally3-5 and computationally.6-13 Experimentally, uniaxially incommensurate (UI) or commensurate herringbone (HB) phases are seen at low temperatures (depending on coverage), which transition into a rectangular solid/liquid coexistence region, melting finally at temperatures ca. 175 K.3-5 At near-monolayer coverages, the melting temperature remains fairly constant, and as the coverage decreases to about F ) 0.5, the melting temperature drops to about 150 K.4 Until recently, most computer simulations of hexane on graphite utilized molecular dynamics (MD) methods and employ the united atom (UA) approximation. In the UA approximation, methyl * To whom correspondence should be addressed. E-mail:
[email protected]. † Department of Physics, University of Northern Iowa. ‡ Department of Computer Science, University of Northern Iowa. § University of Missouri, Columbia. (1) Bruch, L. W.; Cole, M. W. Zarembam, E. Physical Adsorption: Forces and Phenomena; Oxford University Press: Oxford, 1997. (2) Shrimpton, N. D.; Cole, M. W.; Steele, W. A.; Chan, M. H. W. Rare Gases on Graphite Surface Properties of Layered Structures; Benedek, G., Ed.; Kluwer: Amsterdam, 1992. (3) Krim, J.; Suzanne, J.; Shechter, H.; Wang, R.; Taub, H. Surf. Sci. 1985, 162, 446-451. (4) Newton, J. C. Ph.D. dissertation, University of Missouri-Columbia, 1989. (5) Taub, H. NATO AdVanced Study Institutes, Series C: Mathematical and Physical Sciences; Long, G. J., Grandjean, F., Eds.; Kluwer: Dordrecht, The Netherlands, 1988; Vol. 228, pp 467-497. (6) Hansen, F. Y.; Taub, H. Phys. ReV. Lett. 1992, 69, 652-655. (7) Hansen, F. Y.; Newton, J. C.; Taub H. J. Chem. Phys. 1993, 98, 41284141. (8) Velasco, E.; Peters, G. H. J. Chem. Phys. 1995, 102, 1098-1099. (9) Peters, G. H.; Tildesley, D. J. Langmuir 1996, 12, 1557-1565. (10) Peters, G. H. Surf. Sci. 1996, 347, 169-181. (11) Herwig, K. W.; Wu, Z.; Dai, P.; Taub, H.; Hansen, F. Y. J. Chem. Phys. 1997, 107, 5186-5196. (12) Roth, M. W.; Pint, C. L.; Wexler, C. Phys. ReV. B 2005, 71, 155427155439. (13) Pint, C. L.; Roth, M. W.; Wexler, C. Phys. ReV. B 2006, 73, 8542285431.
(CH3) and methylene (CH2) pseudoatoms replace the respective real functional groups in a molecule. The UA approximation saves significantly on the computational effort but has significant shortfalls. The most significant are (i) the lack of in-plane space occupation due to the missing hydrogens in the UA model, (ii) the anisotropy introduced by the terminal hydrogens which is averaged out in the UA model, (iii) the effect of the interaction of the hydrogens with the graphite substrate which can be substantially different than that of the UA model, and (iv) the UA model significantly underestimates the moment of inertia of the molecule. Such UA simulations have provided a framework for advancing our understanding of physisorbed alkanes and are capable of reproducing the melting temperature at completion (F ) 1) fairly accurately. In the most recent studies,12,13 at F ) 1 a solid herringbone phase persists until a transition temperature T1 = 130 K. Then, there is a transition to an orientationally ordered nematic mesophase up until T2 = 172 K above which there is an isotropic liquid. In the near-monolayer range13 (0.875 e F e 1.05) a uniaxially incommensurate herringbone (UI-HB) solid is present. The simulations in the UA approximation, however, have some serious pitfalls: for example, at moderate-to-high coverages the adsorbate molecules are more prominently rolled on their side perpendicular to the surface of the substrate, which is in significant disagreement with experiment.4-7 Since the solid-nematic transition temperature is very sensitive to coverage, an inaccurate description of molecular rolling may result in erroneous characterization of both low- and intermediate-temperature phases. Because of the three main limitations of the UA model mentioned earlier, we expect that including explicitely the hydrogens in the MD simulations will have a significant impact on the simulation results. First off, since the hydrogens occupy space through their collision diameters, including them stifles the molecules’ ability to order and stack as seen in the UA nematic
10.1021/la703040a CCC: $40.75 © 2008 American Chemical Society Published on Web 03/07/2008
MD Simulations of Hexane on Graphite
Langmuir, Vol. 24, No. 7, 2008 3229 Table 1. Dimensions (X,Y) of the Graphite Substrate for the Various Coverages Used in This Work and Corresponding Integer Valuesa F (monolayer)
X (Å)
Y (Å)
nx ) X/4.26
ny ) Y/2.46
1.000 0.966 0.933 0.878 0.614 0.509
68.16 68.16 68.16 68.16 85.20 93.76
68.88 71.34 73.80 78.72 86.10 98.40
16 16 16 16 20 22
28 29 30 32 35 40
a
Figure 1. Perspective snapshot (allowing distortion of atomic sizes) for the initial C6/gr herringbone configuration at area density F ) 1. Adjacent graphene sheets alternate between purple and gray; hexane carbons are blue and hydrogens are red.
mesophase. Specifically, the UA model incorporates a spherically averaged collision diameter for methyl and methylene pseudoatoms, but in the explicit-hydrogen simulations, the presence of the hydrogens gives rise to an effectively anisotropic methyl/ methylene collision diameter. As a result there are directions along which the methyl/methylene collision diameter is on the order of the sum of those for carbon and hydrogen, and this is how the hydrogens occupy space and prevent close stacking. Second, since the methyl and methylene constituents of the hexane molecule are anisotropic, including the hydrogens can help better represent molecule-substrate interactions, stifling the amount of rolling and tilting seen in the UA model especially at higher coverages. Since stacking, rolling, and tilting are significant transition mechanisms in the UA simulations, the results of the explicit-hydrogen model may differ significantly from the UA picture, even though the UA model was meant to represent the same physics. The main purposes of this work are (i) to simulate the hexane/graphite (C6/gr) system at various coverages with explicit hydrogens on the adsorbate molecules and (ii) to compare the results to experiment and the results of UA simulations. This allows for a more realistic representation of various physical aspects of the C6/gr system which in turn will give more meaningful and relevant insight into its dynamics and behavior at phase transitions. II. Computational Method A. About the Simulation Program: NAMD.14 Since a major objective of this work is to simulate hexane with explicit hydrogens, each molecule will now contain more than three times the number of atoms it did in the UA approximation. Moreover, since the hydrogen atoms are so light, their simulation usually requires a smaller time step. As a result, computational time cost increases considerably. For our simulations we utilized the NAMD code14sa parallelized MD simulation package which has been carefully and thoughtfully developed and validated for different systems such as nucleic acids15 and lipid bilayers.16 The optimizations and parallelization of NAMD permitted us to offset the extra cost of the hydrogen inclusion by running our simulations in parallel computer clusters. In addition, we have written pre- and postprocessors in order to generate system specific input files and to reduce the resulting output files, respectively. (14) Kale, L. et al. J. Comp. Phys. 1999, 151, 283-312. See also http:// www.ks.uiuc.edu/Research/namd/. (15) Jha, S.; Coveney, P. V.; Laughton, C. A. J. Comp. Chem. 2005, 26, 1617-1627. (16) Benz, R. W.; Castro-Roma´n, F.; Tobias, D. J.; White, S. H. Biophys. J. 2005, 88, 805-817.
The number of hexane molecules is kept constant (N ) 112).
B. Simulation Setup. The hexane molecule definition in our study was obtained from the Brookhaven Protein Data Bank (PDB).17 A constant N ) 112 hexane molecule system is used for each coverage in an initial herringbone configuration. Because NAMD currently does not include an analytical expression for adatom-substrate interactions, the substrate must be modeled in an all-atom fashion. The graphite is modeled as six identical stacked graphene sheets, stacked in the known (-A-B-A-B-) pattern in the z direction. At the temperatures considered in this work, the dynamics of the graphite does not contribute appreciably to those of the adlayer and so the graphite is static, effectively serving as a lattice of interaction sites for the adsorbed layer. A snapshot of the system’s initial configuration, including the graphite substrate and a F ) 1 hexane layer, is shown in Figure 1. The (x,y) dimensions of the graphite sheets are varied for each coverage studied, respecting the constraint that the computational cell must remain commensurate with the graphite substrate. At nearmonolayer densities, the system is in a uniaxially incommensurate phase and so, in order to mirror such symmetry in the simulations, only the Y parameter of the computational cell is varied for the four highest densities studied here. For all densities studied, the initial configuration was a herringbone lattice expended uniformly in the y direction that was allowed to equilibrate. At densities where vacancies begin to form, other initial configurations, such as a patch in the center of the computational cell were utilized and the results obtained were not appreciably different from those with the expended configuration. Appropriate substrate dimensions are shown in Table 1. All simulations in this study are constant molecule number, coverage, and temperature (N ) 112, F, and T, respectively). Periodic boundary conditions are used in the (x,y) plane and free-boundary conditions are applied in the vertical z direction. To maintain a constant temperature, velocity rescaling is utilized. The time step for all simulations is 1 fs, and each simulation ran for 250 000 equilibration steps and then averages are taken over the next 500 000750 000 time steps. C. Interaction Potentials. All particle-particle interactions in the simulations presented here are in the standard CHARMM22 format.18 The internal bonded degrees of freedom consider interactions between atoms within the same molecule created by a chemical bond. The first internal molecular degree of freedom considered is two-body C-C bond stretching, whose force arises from the harmonic potential ustretch ) k(l - l0)2
(1)
Here k is the bond stiffness and l0 is the equilibrium bond length. The C-H bonds were held rigid so that as large a time step as possible could be used and also because such an internal degree of freedom is probably unimportant to the dynamics of interest in this study. The three-body bond angle bending is considered (17) http://www.wwpdb.org/. (18) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; States, S.; Swaminathan, S.; Karplus, M. J. Comp. Chem. 1983, 4, 187-217. See also http:// www.charmm.org.
3230 Langmuir, Vol. 24, No. 7, 2008
Connolly et al.
ubend(θ) ) kθ(θ - θ0)2
(2)
Table 2. Bonded Potential Parameters Used in the Simulationsa parameter
Here θ is the bond angle, θ0 is its equilibrium value, and kθ is the angular harmonic stiffness constant. Four-body dihedral torsion of the form udihed ) kd{1 + cos(nφd - δ)}
K (kcal/mol ) l0 (Å) kθ (kcal/mol rad2)
[( ) ( ) ] r0 rij
12
-2
r0 rij
kqiqj rij
(6)
A. Unit Coverage (F ) 1). The F ) 1 low-temperature solid is a commensurate herringbone structure where almost all molecules lie flat to the surface. This can be seen visually in Figure 2. As temperature is increased for the system, various order parameters and physical quantities are calculated in order to quantify and understand its behavior and changes.13 The herringbone order parameter OPherr gives a measure of the orientational order of the molecular axis with respect to the graphite substrate. Here the “molecular axis” is along a molecule’s smallest principal moment of inertia. OPherr is defined by
〈∑ Nm
i)1
(-1)j sin (2φi)
All values are in CHARMM2218 format; methyl (CH3) groups are represented by “C3” and methylene groups by “C2”. Table 3. Nonbonded Potential Parameters Used in the Simulationsa parameter
value -0.055 (C3) -0.08 (C2) -0.022 (H) -0.07 (C) 2.06 (C3) 2.175 (C2) 1.32 (H) 1.9924 (C) -0.27 (C3) -0.18 (C2) 0.09 (H) 0 (C)
(kcal/mol)
σ (Å)
III. Results and Discussion
Nm
δ (deg)
(5)
which arises due to the internal charge redistribution between carbons and hydrogens (even though hexane is not polar, there is a small charge distribution). The electrostatic interactions are calculated by the well-known particle mesh Ewald (PME) summation technique and involves pairs on the same, as well as other, molecules. The nonbonded potential parameters are given in Table 3.
OPherr t
n (units)
a
are utilized to calculate potential parameters in the cases involving mixed interactions. All Lennard-Jones lattice sums are taken out to a pair separation rij ) 7.5 Å and then smoothly diminished to zero at a separation of rij ) 10 Å. Coulomb interactions are also included
1
kd (kcal/mol)
(4)
1 σij ) (σi + σj) 2
uC )
θ0 (deg)
6
The rij-12 term results in a repulsive force that becomes considerable when electron clouds of interacting atoms overlap. Ultimately, such repulsion has its origin in the Pauli Exclusion Principle. The rij-6 dispersion term represents an attractive force and has its origin in the van der Waals forces arising from fluctuating dipolar interactions. In eq 4, ij is the potential well depth and σij are the well depth and collision diameter, respectively, for interactions involving particles i and j. Lorentz-Bertholot combining rules ij ) xij,
222.5 1.53 58.00 (C3-C2-C2) 58.35 (C2-C2-C2) 36.00 (H-C3-H) 36.00 (H-C2-H) 115.0 (C3-C2-C2) 113.6 (C2-C2-C2) 115.0 (H-C3-H) 115.00 (H-C2-H) 0.15 (C3-C2-C2-C2) 0.15 (C2-C2-C2-C2) 0.195 (H-C2-C2-H) 0.16 (H-C3-C2-H) 1 (C3-C2-C2-C2) 1 (C2-C2-C2-C2) 3 (H-C2-C2-H) 3 (H-C3-C2-H) 0 (C3-C2-C2-C2) 0 (C2-C2-C2-C2) 0 (H-C2-C2-H) 0 (H-C2-C2-H)
(3)
is also included, where kd is the torsional stiffness constant, n is the multiplicity, φd is the dihedral angle, and δ is an angular offset. This expression for torsional energy, and hence the description of gauche defect formation, differs from that used in previous work6-10,12,13 because in this case the hydrogens also contribute through methyl torsion. Still, in our simulations gauche defects are calculated as rotations around bonds in the carbon backbone of the molecule. The bonded potential parameters used in the simulations are given in Table 2. Nonbonded interactions include two-body C-C and C-H interactions for atom pairs either in different molecules and also adatom-substrate pair interactions. The first type of nonbonded interaction is represented by a modified the Lennard-Jones potential uLJ(rij) ) ij
value
Å2
〉
(7)
where φi is the angle between the axis of greatest moment of inertia of molecule i and the x axis. The parity of the integer j
q (esu)
a All values are in CHARMM22 format.18 Like-species interactions are shown here, and parameters for mixed interactions are given by eq 5. Here data for the graphite carbons are represented with “C”.
accounts for the different orientational herringbone sublattices. OPherr is equal to unity if all molecular axes are oriented at 45° and 135° for each sublattice, respectively, and drops as a result of thermal fluctuations to zero for random in-plane orientations. Since the F ) 1 static herringbone lattice for the systems studied here has angles at about 38° and 153°, OPherr has a zerotemperature limiting value of about 0.866 and drops from there. At low temperatures, the large value of OPherr (Figure 3) along with the sharp peaks in azimuthal angle probability distributions P(φ) (Figure 4) characterize a strong herringbone solid. The microscopic roll angle ψ is defined as
Ψ ) cos-1
{
}
[(b r j+1 - b r j) × (b r j-1 - b r j)]‚zˆ |(b r j+1 - b r j) × ( b r j-1 - b r j)|
(8)
Here b rj is the position vector for atom j and the expression for ψ involves the two other atoms that comprise the triplet making up the bond angle. This roll angle takes on a value of 0° when the plane consisting of the three molecules in the bond is parallel to the graphite substrate and takes on a value of 90° when the plane is perpendicular to the substrate, and it gives an idea of the degree of rolling of segments of the molecules in the overlayer.
MD Simulations of Hexane on Graphite
Langmuir, Vol. 24, No. 7, 2008 3231
Figure 4. Azimuthal angle probability distributions P(φ) for various coverages in the solid (black), at the onset of orientational order loss (blue), at the onset of nematic order loss (green), and in the liquid (purple). The plot for F ) 0.614 has no color because neither OPHerr or OPNem show signatures of transitionns. Perspective angles for various coverages are different in order to better visualize features at each temperature.
Figure 2. Snapshots of the system at F ) 1 for the low-temperature solid (T ) 100 K), at the onset of herringbone orientational order loss (T ) 140 K), at the onset of nematic order loss (T ) 162 K), and in the isotropic liquid (T ) 180 K). Color scheme is the same as in Figure 1.
Figure 3. Herringbone order parameter OPHerr (left) and nematic order parameter OPnem (right) as functions of temperature for coverages F ) 1 (black), 0.966 (red), 0.933 (green), 0.878 (blue), and 0.614 (purple).
If the molecules have their carbon backbone flat against the substrate, ψ ) 0° or 180° and if they are “on their side,” ψ ) 90°. Figure 5 shows the microscopic roll angle probability distributions; it is clear that most of them are flat at low temperature, consistent with experimental diffraction results.4,6,7 As temperature is increased, it is evident from the data presented that the herringbone solid loses its orientational order at around T1 ) 145 ( 2 K. Note however, that all the data presented suggest that considerable orientational subsists even after the herringbone disappears. The nematic order parameter shown in Figure 3 illustrates that the system indeed is in a mesophase, which persists up until melting into the isotropic fluid at about T2 ) 160 ( 3 K. The nematic order parameter gives a measure of the long-range orientational order (and thereby the nematic character) of the system. It is given by12,13
OPNem t
1 Nm
〈∑ Nm
i)1
〉
cos 2(φi - φdir)
(9)
Figure 5. Microscopic roll angle probability distributions P(ψ) for various coverages in the same format as Figure 4.
where φi is the angle that the axis of molecule i makes with the x axis and φdir is a director which is essentially the average orientation of the system.13 In UA calculations13 the herringbone phase persists up until about T1 ) 145 K and then changes to a strong nematic, which in turn persists until about T2 ) 170 K when the system melts into the isotropic liquid. The nematic is still present in the fully atomistic simulations, but the degree of order as indicated by the value of the nematic order parameter is substantially less than that seen in the UA case, as seen by the absence of a sharp increase in OPNem after herringbone order is lost but and in the absence of high peaks in the azimuthal angle probability distributions in the region T1 < T < T2. The substantially weakened nematic phase in the all-hydrogen model is more consistent with experiment3-7 in that no true nematic phase is observed. Rather a lattice of rectangularly centered islands mobile in the fluid is inferred from the scattering data. Pair correlation functions and static structure factors in Figure 6 confirm the two transition temperatures. Still, there is no signature seen in experimental results for the herringbone-nematic transition at T1, but rather the melting transition at T2 ≈ 172 K, which is a few Kelvin higher than calculated in this work.
3232 Langmuir, Vol. 24, No. 7, 2008
Connolly et al.
Figure 6. Pair correlation function (left) and static structure factor (right) at selected temperatures for F ) 1.
Figure 7. Average Lennard-Jones interaction as a function of temperature for F ) 1.
In order to understand why the nematic mesophase is suppressed in the explicit-hydrogen calculations, we also examine the average adatom-adatom Lennard-Jones interaction as a function of temperature
〈∑ ∑ [( ) ( ) ]〉 σij
N
U)
4ij
i)1 j>i
rij
12
-
σij rij
6
(10)
which are shown for unit coverage in Figure 7. The herringbone to nematic transition was accompanied by a sharp drop in the molecule-molecule interaction energy in the UA model12,13 (an increase in |U|). In contrast, here we see a gradual reduction in |U|. The relative suppression of the nematic phase brought about by including the hydrogens has to do in one sense with the hydrogen atoms occupying more in-plane space. In the UA model the anisotropy of the hydrogens are averaged out and the interaction centers (pseudoatoms) have the opportunity to sit very near their equilibrium positions when oriented and stacked as in the nematic phase. In the explicit-hydrogen picture, however, the hydrogens are weakly interacting centers but since they are concentrated at specific spots they are able to prevent the entire molecules from getting so close together and sampling the steep repulsive part of the Lennard-Jones interaction potential. Hence, the carbons are farther apart than in the UA model and their nonbonded interaction energy does not drop at T ) T1. Hence, the explicit-hydrogen molecules are simply not allowed to stack and pack closely together; hence, the Lennard-Jones interactions are weaker. In addition, the anisotropic interaction of the hydrogens with the substrate is important, which is evidenced by the more reasonable rolling behavior afforded in the allhydrogen model.
Figure 8. Snapshots at three submonolayer coverages. Color scheme is the same as in Figure 1; the graphite substrate for the lower two coverages is removed for clarity but the approximate size of the computational cells are shown. To emphasize topological features, the cells are distorted so the areas do not visually scale. Table 4. Positions of the Peaks in the Azimuthal Angle Probability Function P(O) for Various Coverages in the Low-Temperature Solid F (monolayer)
peak positions (deg)
1.00 0.966 0.933 0.878 0.614
28, 153 29, 151 31,149 30, 150 0, 35, 61, 104, 153, 172
B. Partial Coverages (F < 1). As the coverage is decreased slightly below unity, the herringbone transition at T1 shows very little change. It is interesting to note that, based on the data in Figure 3, F ) 0.966 actually retains order to slightly higher temperatures than F ) 1. This is because at F ) 1 there are vertical atomic fluctuations and ultimately about 1% layer promotion (Figure 2)which creates in-plane room which in turn lowers the transition temperatures. At F ) 0.966 the spreading pressure is somewhat lowered and vertical fluctuations are not as prevalent. Hence, there is less in-plane mobility and the transition temperatures are not lowered. Moreover, nematic mesophase exists over a very short temperature range at F ) 0.966, consistent with the conclusion that preventing close stacking of the molecules weakens the mesophase. However, at F ) 0.878, T1 decreases dramatically. This decrease is directly related to the onset of vacancy and vacancy patch formation in the systemsnot a UI phase any more. Inspection of Figure 8 (left) confirms the presence of vacancy patches and oriented domains. In fact, by F ) 0.614 the topology is a connected network and is very similar to that seen at, F ) 0.509, shown in Figure 8 (right). Below coverages of about 0.5, the system breaks up into individual patches, as seen in Figure 8. It is interesting to note that, although more prevalent at F ) 0.878, all submonolayer systems shown have static gauche defects usually present at grain boundaries. While the herringbone phase subsists, the peaks of the azimuthal angle distribution are not exactly at 30° and 150° but do show a collection of individual peaks representing many differently oriented domainssa coverage dependence in reasonable agreement with experiment.4,6,7 Table 4 shows some related numerical data. The interpretation that the connected network existing at F ) 0.509 and 0.614 is comprised of oriented domains is supported not only visually in Figure 8 but also by the appearance of the azimuthal angle probability distributions shown in Figure 4. Inspection of the dihedral angle distributions in Figure 9 suggests that in this model the herringbone order-disorder transition depends on the presence of gauche defects, because in the highest
MD Simulations of Hexane on Graphite
Langmuir, Vol. 24, No. 7, 2008 3233
Figure 11. Same as Figure 6 but for F ) 0.614. Figure 9. Dihedral angle probability distributions P(φd) in the same format as Figure 4.
Figure 10. Same as Figure 6 but for F ) 0.878. Figure 12. End-to-end H-H length as functions of temperature at various coverages. Format is the same as Figure 3.
densities defects are not present until OPHerr starts to drop but at lower densities they are present in the solid where defects are evident. In fact the presence of gauche defects at lower temperatures at lower densities is surprising and suggests that hexane is not so rigid as previously believed. Moreover, the presence of gauche defects dynamically reduce the time-averaged Lennard-Jones nonbonding energy so as to promote the weakness of the nematic phase with respect to the UA model. For the highest densities examined, the behavior of OPNem in Figure 3 shows considerable insensitivity of the melting temperature T2 with decreasing coverage until at about F ) 0.9. Since the presence of vacancies and oriented domains wash out angular order parmeters, we depend on the pair correlation function and static structure factors (shown in Figures 10 and 11) to help determine the location of T2. It is evident that T2 is somewhere between 130 and 140 K, again in reasonable agreement with experimental values of between 150 and 160 K.4 It is interesting to note that dihedral angle probability distributions in Figure 9 show the presence of considerable gauche defects at temperatures that depend on coverage. In fact, inspection of the left panel of Figure 8 confirms the presence of gauche defects in the low-temperature submonolayer solid. The presence of such defects can be seen in the snapshots as molecules whose backbones are not straight and in all the submonolayer cases are likely found in the region between two or more grain boundaries. Furthermore, it can be seen that at low temperature the gauche defects allow the molecule to conform to the graphite substrate. Such results strongly suggest that the hydrogens interact strongly with the graphite substrate vis-a-vis lateral corrugation and the
gauche defects in fact promote the domain formation by relieving grain boundary stress. So the gauche peaks for F ) 0.878 are static in nature. Then, they disappear as grain boundary stress is relieved and reappear dynamically at higher temperatures. Surprisingly, the end-to-end H-H lengths in Figure 12 seem to show something contrarysthat, the molecules begin to bend at nearly the same temperature for all coverages. Inspection of the end-to-end length distribution (not shown here) shows a single peak as well. Such results may be understood in that since, at low temperatures the static submonolayer gauche defects cause the molecule to conform the substrate, they do so in such a way that the H-H distances are the same. This is contrary to the dynamic higher-temperature gauche defects, where methyl torsion is present. All the dihedral and end-to-end distance results, coupled with the data in Table 4, suggest that as coverage increases up to the monolayer there is some strain placed on the molecules in the herringbone solid which is relieved as coverage is decreased but at the cost of grain boundary formation. Moreover, this also explains the insensitivity of the transition temperatures with decreasing coverage near F ) 1, which is anomalous in simple physisorbed systems. In turn, this suggests that true monolayer completion for this interesting system may not occur at the value previously and usually deemed F ) 1 in this system.
IV. Conclusions The work presented here allows us to make several important conclusions about simulating hexane on graphite with explicit hydrogens:
3234 Langmuir, Vol. 24, No. 7, 2008
(i) Since the hydrogens represent concentrated interaction centers but are averaged out in the UA model, explicit-hydrogen molecules are prevented from coming so close together as they do in UA simulations. Hence, the carbon backbones are prevented from aligning and stacking closely, as seen in the UA model, which in turn weakens the explicit-hydrogen nematic. Gauche defects also promote a dynamic reduction in strength of the adatom-adatom energy, which further weakens the nematic. (ii) The herringbone solid-to-nematic phase transition is fairly insensitive to coverage near completion and begins to become affected when defects (vacancy islands; percolated networks) are present in the system. This is in disagreement with previous work with the UA model,13 where the same transition temperature (T1) increases with decreasing coverage and vice versa. Since the molecules in the UA model roll in the solid and have to tilt in order to form the nematic, decreasing stress through decreasing density lessens the likelihood for such dynamics to be present. The explicit hydrogen molecules do not significantly roll or tilt prior to the transition. (iii) The loss of herringbone order at higher densities is directly related to the presence of gauche defects in the explicit-hydrogen molecules. In previous work with the UA model,13 gauche defects became prominent throughout the nematic but not so much before.
Connolly et al.
(iv) As coverage is decreased down below about F ) 0.9 some gauche defects are present in the solid but the end-to-end H-H distance drops at the same temperature regardless of coverage, suggesting that the hydrogens are able to notch into the graphite substrate. (v) There appears to be strain present in the higher-coverage systems which is relaxed at lower densities. The strain present helps explain the insensitivity of the transition temperatures at the highest densities examined. (vi) The spreading pressure observed in both the UA and explicit-hydrogen simulations is comparable. The fact that it is not zero at F ) 1 indicates that the point of monolayer completion based on simulations must be further investigated, and work with larger systems is highly desirable. Acknowledgment. Acknowledgment is made to the Donors of The American Chemical Society Petroleum Research Fund (PRF43277-B5) and the University of Missouri Research Board for the support of this research. This material is based upon work supported in part by the Department of Energy under Award No. DE-FG02-07ER46411. The authors acknowledge fruitful discussions with Haskell Taub and Flemming Hansen. LA703040A
