Simulation and Characterization of Tetracosane on Graphite

Article pubs.acs.org/JPCC

Simulation and Characterization of Tetracosane on Graphite: Molecular Dynamics Beyond the Monolayer M.W. Roth,† L. Firlej,‡,§ B. Kuchta,*,§,∥ M. J. Connolly,§ E. Maldonado,§,⊥ and C. Wexler§ †

Department of Physics, Geology and Engineering Technology, Northern Kentucky University, Highland Heights, Kentucky 41099, United States ‡ Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Université de Montpellier, F-34095 Montpellier, France ∥ Laboratoire MADIREL, Aix-Marseille Université-CNRS UMR7246, 13396 Marseille, France § Department of Physics and Astronomy, University of Missouri, Columbia, Missouri 65211,United States ⊥ Department of Physics, University of Northern Iowa, Cedar Falls, Iowa 50614,United States ABSTRACT: We present the results of extensive fully atomistic molecular dynamics (MD) simulations of tetracosane (C24H50) bilayer and trilayer systems adsorbed onto the basal plane of graphite. At low temperature, both layers of the bilayer exist in well-defined solid phases. With increasing temperature, the system exhibits separated smectic phases that eventually lead to melting. During this process, we observed a strong interlayer translational correlation and mobility between layers; however, the upper layer presents more intra- (chain) and intermolecular disorder because of a lack of confinement and a greater distance to the graphite substrate. Simulations of the perpendicular trilayer patch show that gauche defects provide the main mechanism for spreading of the bottom and outer perimeter of the patch in the solid, leading to the ultimate collapse of the patch with increasing temperature and formation of a flat (parallel) trilayer that melts at a higher temperature than the bilayer structure. The wide variety of structural order parameters, thermodynamic functions, and probability distributions we employed provide a clear picture of the roles of gauche defects, confinement, and interlayer correlation in the phases and phase transitions exhibited by these confined organic layers.

1. INTRODUCTION Nanostructured systems are ubiquitous because the interactions that make the molecules aggregate into larger systems operate on natural subnanometer length scales. Consequently, an understanding of nanoscale phenomena is crucial for the development of novel materials tailored at the nanoscale and having desired functionalities and technological importance. Examples include low-dimensional semiconductor nanostructures (quantum wires,1 quantum dots2), carbon nanostructures (graphene,3 nanotubes4), and variety of biological nanostructures.5 Studies of such systems are of both fundamental and applied importance because they are central to a better understanding and the ultimate rational design of larger systems in which nanoscale subsystems are integrated, perhaps by self-assembly and/or selfreplication.6 Many nanostructured, self-assembled systems are supported on surfaces, and the surface contribution to the interactions in an adsorbed system is never negligible. The presence of the substrate modifies the interparticle interactions and determines the self-organization mechanism of adlayers at the nanoscale. Therefore, specific conditions in adsorption can be used for tuning the structure and self-assembly of molecules toward a desired configuration.7,8 In all such situations, the characteristics © 2015 American Chemical Society

of the substrate are crucial in determining the properties of nanomolecular structures. The primary mechanism of the system’s construction is defined by the mono- and multilayer formation and structural evolution of these structures as a function of temperature. In previous studies,9−13 we investigated the structures and phase transitions of self-assembled monolayers of alkanes on a graphite surface. We found that the systems’ evolutions with temperature toward disordered states have many characteristics of conventional phase transitions. However, the lamellar structure of the adsorbed monolayer at low temperature induces the creation of interesting mesophases that make the phase diagram even richer. Most noteworthy is the presence of smectic and/or nematic phases bracketed by the solid and liquid.9−13 The behavior of alkane multilayers is even more complex. The diffusive motion in a fluid bilayer of tetracosane C24H50 (C24) deposited on Grafoil was thoroughly investigated using highresolution quasielastic neutron scattering.14,15 The experimental spectra indicated that the diffusive motion in the fluid bilayer is Received: October 3, 2015 Revised: December 15, 2015 Published: December 18, 2015 984

DOI: 10.1021/acs.jpcc.5b09677 J. Phys. Chem. C 2016, 120, 984−994

Article

The Journal of Physical Chemistry C

Figure 1. Initial configuration for the C24/gr bilayer (left) and trilayer patch (right) systems. Computational cells: Xbilayer = 17 × 4.26 Å, Ybilayer = 53 × 2.46 Å; Xtrilayer = 17 × 4.26 Å, Ytrilayer = 27 × 2.46 Å. Lattice spacings are distorted due to differing perspectives.

slower than that in the previously measured monolayer.16,17 The authors speculated that the large length/cross section aspect ratio of the C24 molecule might inhibit interlayer exchanges of molecules.14,15 Initial simulations conducted in the united-atom approximation18 indicated that the bilayer melting temperature is significantly higher than that of the monolayer and that the upper layer is much more disordered than the bottom layer and melts at a lower temperature. Structural investigations of trilayer C24 and C32 systems indicated that the first two layers are parallel to the graphite substrate and that the third layer is perpendicular, which motivated our choice for the perpendicular trilayer patch discussed later.19 The broader goals of this article are (1) to promote the comprehension of the construction mechanism of self-assembled molecular adlayers, which is indispensable for understanding phenomena at the nanoscale, and (2) to forward an understanding of the effects of various adsorption situations. Specifically, we study and discuss structural and thermodynamic information gained from computer simulations of adsorbed tetracosane layer structures beyond the monolayer. We utilize fully atomistic MD simulations to study bilayers and trilayers of tetracosane adsorbed on graphite. The choice of these two systems was motivated by experimental observations of C24 multilayers,14−16 in particular, the spreading of the perpendicular patch in the trilayer system, revealing the importance of heightened gauche defect formation with increasing temperature. We employ a confluence of structural order parameters, thermodynamic functions, and probability distributions to elucidate the roles of gauche defects, confinement, and interlayer correlations in the phases and phase transitions exhibited by these confined organic layers.

arranged in a rectangular-centered structure, with the two layers staggered with respect to each other (Figure 1, left). (ii) For the trilayer, a third perpendicular layer was added atop the bilayer configuration (Figure 1, right). Such initial configurations were chosen to (i) represent the experimentally observed structural domains of the adsorbed films14−17,19 and (ii) model the dynamics of the vertical layer, governed by the formations of dihedral (gauche) defects. 2.2. Interaction Model. All interaction potentials used in the simulations presented here are in the standard CHARMM22 format.21 They account for bonded, nonbonded, and electrostatic interactions in the system. The two-body C−C bond stretching was described by the harmonic potential

ustretch = k(l − l0)2

(1)

where k is the bond stiffness and l0 is the equilibrium bond length. In our simulations, the C−H bonds were held rigid using the SHAKEH algorithm within the molecular dynamics package used (NAMD2),22 so that a larger time step (typically 1 fs) could be used during molecular dynamics (MD) simulations; moreover, given the relatively large values of k, this internal degree of freedom is typically not relevant to the dynamics of interest in this study. Two terms related to an angle deformation have been considered: (i) a three-body bond-angle (θ) bending ubend(θ ) = kθ(θ − θ0)2

(2)

where θ0 is the bond-angle equilibrium value and kθ is the angular harmonic stiffness constant, and (ii) a four-body dihedral torsion of the form udihed = kd[1 + cos(nφd − δ)]

2. SIMULATION DETAILS 2.1. Model of the System. The initial configurations of the tetracosane bilayer and trilayer systems are presented in Figure 1, with snapshots of typical equilibrium configurations for the bilayer system in the solid, smectic, and melt presented in Figure 2. The graphite substrate is modeled in an all-atom fashion, as six identical sheets stacked in an A−B−A−B graphite pattern. The tetracosane molecular structure and charge distribution over the atoms were taken from the Brookhaven Protein Data Bank (PDB).20 Hydrogen atoms are explicitly included in the model of the hydrocarbons.9,10 The number of alkane molecules in the bilayer simulation box is Nbilayer = 128 (64 molecules in each layer); for thetrilayer, Ntrilayer = 185 (32 molecules in each flat layer and 121 molecules in the perpendicular patch). The initial (low-temperature) configurations were prepared as follows: (i) In the bilayer, the molecules were placed in parallel lamellae

(3)

where kd is the torsional stiffness constant, n is the multiplicity factor, φd is the dihedral angle, and δ is an angular offset. Nonbonded interactions comprise two-body atom−atom interactions for (i) atom pairs either in the same or different molecules and (ii) adatom−substrate pair interactions. In addition to Coulomb interactions between (Mulliken) partial charges, the Lennard-Jones (LJ) interaction between atoms i and j is represented by ⎡⎛ ⎞12 ⎛ r ⎞6 ⎤ r0 ⎢ ⎜ ⎟ ⎜ uLJ(rij) = εij⎢⎜ ⎟ − 2⎜ 0 ⎟⎟ ⎥⎥ r ⎝ rij ⎠ ⎦ ⎣⎝ ij ⎠

(4)

where εij is the potential well depth and σij is the collision diameter. Lorentz−Berthelot combining rules were utilized to 985


Article


Figure 2. (Top and middle) Snapshots of the top and bottom layers, respectively, at (left to right) 275 K (solid), 350 K (upper-layer smectic), and 400 K (both layers melted). Atoms colored red signify that they are significantly higher off the surface relative to their respective layer. (Bottom) Composite showing upper (green) and lower (blue) layers together.

Correctly accounting for the intramolecular terms in the interaction model is extremely important. In fact, the interplay between the bonded and nonbonded intramolecular interactions directly defines the overall stiffness of the molecules, which, in turn, has a crucial effect on the mechanism of melting.14,16,23 We used the scaled 1−4 exclusion policy implemented in the CHARMM force field: It assumes that all 1−2 and 1−3 pairs of atoms (separated by one or two bonds in a 1−2−3−4−5− atom chain) are excluded from calculations of nonbonded interactions (bonded interactions are all-inclusive). For 1−5 and longer

calculate potential parameters in cases involving mixed interactions εij =

εiεj ,

σij =

1 (σi + σj) 2

(5)

All LJ potentials were taken out to a pair separation of rij = 7.5 Å and then smoothly diminished to zero at rij = 10 Å. Coulomb interactions were calculated using the particle mesh Ewald (PME) summation method. The values of potential parameters are listed in Tables 1 and 2. 986


Article

The Journal of Physical Chemistry C Table 1. Bonded Potential Parameters Used in the Simulationsa parameter K [kcal/(mol Ǻ 2)] l0 (Ǻ ) kθ (kcal/mol rad2)

θ0 (deg)

kd (kcal/mol)

n (units)

δ (deg)

value 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)

Figure 3. Order parameters OP2 (squares), OPnem (triangles), and Φ4 (circles) plotted for upper (green solid symbols) and lower (blue open symbols) layers of the bilayer at various temperatures.

a

All values are in CHARMM22 format.21 methyl (CH3) groups are represented by C3, and methylene (CH2) groups by C2; see elsewhere for further discussion.9

Table 2. Nonbonded Potential Parameters Used in the Simulationsa parameter

value

ε (kcal/mol)

0.07 (Cgr) 0.055 (C3) 0.08 (C2) 0.022 (H) 1.9924 (Cgr) 2.06 (C3) 2.175 (C2) 1.32 (H) 0 (Cgr) −0.27 (C3) −0.18 (C2) 0.09 (H)

σ (Ǻ )

q (esu)

Figure 4. End-to-end molecular lengths for upper (green solid circles) and lower (blue open circles) layers of the bilayer at various temperatures.

Simulation boxes commensurate with a graphite underlayer with dimensions of (Xbilayer = 17 × 4.26 Å, Ybilayer = 53 × 2.46 Å) and (Xtrilayer = 17 × 4.26 Å, Ytrilayer = 27 × 2.46 Å) were utilized. Periodic boundary conditions were implemented in the xy plane, and free boundary conditions were applied in the perpendicular (z) direction in both systems. In the trilayer case, our intention was to simulate a system that is infinite in the xy plane for the lower two layers and a finite patch for the upper layer of the trilayer system (Figure 1). The time step in the simulations was 1 fs, with a velocity rescaling thermostat called every 10 fs. The relaxation of the initial conditions for all systems was carefully examined, and the stability of bonded and nonbonded energetics in the system was checked to ensure that no noticeable trends occurred over the time scale of a typical run. Thermalization runs of few nanoseconds were conducted before the production runs [(3−5) × 106 steps, giving a 3−5-ns time interval] over which thermal averages were taken. 2.4. Characterization Methods. 2.4.1. Order Parameters. To quantify and characterize various aspects of the systems’ structures, we employed an extensive collection of order parameters10,11 that monitor molecular orientation. OPn monitors the collective azimuthal orientation of molecules and is defined as

a

All values are in CHARMM22 format.21 Like-species interactions are shown; the parameters for mixed interactions are given by eq 5. The data for the graphite carbons are indicated by Cgr. See elsewhere for further discussion.9.

distances, both electrostatic and van der Waals (vdW) energies are fully taken into account. At intermediate distances, for atoms separated by three bonds (1−4 interactions), the explicit nonbonded part, already partially included in the dihedral angle torsion term, is scaled down. The exact value of the scaling factor (SF) is subject to debate (see elsewhere for a full discussion23). In this work, we assumed SF = 0.45, identified as an appropriate value for C24.10,23 2.3. Simulation Details. The MD simulations were performed using the parallel NAMD2 program22running on a Linux Beowulf cluster. The NVT canonical ensemble MD method was used to simulate the tetracosane/graphite system.

OPn = 987

1 Nm

Nm

∑ cos nφi i=1

(6) DOI: 10.1021/acs.jpcc.5b09677 J. Phys. Chem. C 2016, 120, 984−994

Article


Figure 5. End-to-end molecular length distributions plotted for the (top) upper layer and (bottom) lower layer of the bilayer for various temperatures, color-coded by phase. Temperatures (back to front): 100, 150, 200, 250, 275, 300, 325, 350, 375, and 400 K.

where Nm is the number of molecules in the sample, the brackets represent time averages, and φi is the angle between the ith molecule’s principal axis (along its long dimension) and the x axis of the simulation box (azimuthal angle). OP2 = 1 when all molecules are oriented along the x direction, OP2 = −1 when all molecules are oriented along the ±y direction, and OP2 vanishes if either all molecules are oriented along the bisector of the xy plane or the molecules are randomly oriented. The orientation along an arbitrary direction (stacking) is characterized by the nematic order parameter, OPnem10,11,24 OPnem =

1 Nm

OPnem = 1 if all molecules in the simulation are aligned (in any direction φdir, even if this direction fluctuates), and OPnem vanishes when the orientation of the molecules is random. OPnem takes characteristic intermediate values when stacked domains are present, such as in the case of the “herringbone” configurations often observed for short linear molecules.9,11 The tilt order parameter, OPtilt,11,12 characterizes collective out-of-plane orientation N

OPtilt =

Nm

∑ cos 2(φi − φdir) i=1

⎡ ∑Nm sin(2φ ) ⎤ 1 i ⎥ tan−1⎢ Ni = 1 ⎢⎣ ∑ m cos(2φ ) ⎦⎥ 2 i i=1

(9)

where θi is the angle between the smallest moment of inertia (long axis) of molecule i and the direction normal to the adsorption plane. OPtilt = −0.5 if each molecule backbone is parallel to the xy plane, and OPtilt = +1 if all molecules axes are perpendicular to it. OPtilt = +0.25 for either free rotors in θ or a random static distribution of molecular orientations with respect to the adsorption plane. The order parameters Φn are utilized10 to track structural n-fold symmetry of the center-of-mass lattice positions by

(7)

where the director angle φdir is defined in each configuration, by maximizing OPnem11 φdir =

m 1 ⟨∑ (3 cos2 θi − 1)⟩ 2Nm i = 1

(8) 988


Article


Figure 6. Dihedral angle distributions plotted for the (top) upper layer and (bottom) lower layer of the bilayer for various temperatures colorcoded by phase. Temperatures (back to front): 100, 150, 200, 250, 275, 300, 325, 350, 375, and 400 K.

Figure 7. Azimuthal molecular orientation angle distributions plotted for the (top) upper layer and (bottom) lower layer of the bilayer for various temperatures and color-coded by phase. Temperatures (back to front): 100, 150, 200, 250, 275, 300, 325, 350, 375, and 400 K. Note lower and broader peaks for the upper layer toward the hightemperature side of the smectic, signaling the increased presence of gauche defects.

measuring the average bond order within a plane layer and are defined as Φn =

1 Nb

Nb

∑ exp(inψk) k=1

(10)

in-plane disordering, and vertical mobility of the adlayer. Because the tetracosane−graphite system does not show a clear commensurate phase, we used the molecule−graphite energy to characterize collective interlayer migration only. We also used the thermal average of the Lennard-Jones (van der Waals) energy to analyze the ordering and fluctuations within each system layer, in particular, the dynamics of the perpendicular patch in the case of the trilayer system

The centers of mass of each nearest-neighbor pair is connected by a vector (“bond”) that has a particular azimuthal orientation ψk. The index k runs over the total number of nearest-neighbor bonds Nb in the adsorbed layer. Here, “neighbors” are defined to be bonds of given distances that fall within minima of the pair correlation function. The parameter n accounts for layer symmetry and gives the number of nearest neighbors in the ideal structure (e.g., for a triangular lattice, Φ6 ≈ 1; for a square lattice, Φ4 ≈ 1; and all Φn = 0, for an isotropic or random system). Most importantly for this work, Φ4 drops sharply when the system undergoes a transition to the smectic mesophase, where lamellar shifting is present and the rectangular-centered lattice seen in the solid distorts. 2.4.2. Energy Parameters. In addition to order parameters, thermal averages of interaction energies can be used to characterize the system. The average atom−graphite energy ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σ σ 4 ε ∑ iC⎢⎢⎜ iC ⎟ − ⎜ iC ⎟ ⎥⎥ r riC ⎠ ⎦ ⎝ ⎠ ⎝ ⎣ iC i=1

ULJ =

(12)

2.4.3. Distributions. To complete our analysis of the transformations that the system undergoes with increasing temperature, we also accumulated distributions throughout the simulation of various structural parameters. In fact, the order parameters involve averages and moments of the distributions, and therefore, they statistically wash out fine details of sometimes subtle structural changes the system undergoes. In particular, the end-to-end molecular length distributions give layer-by-layer insight into the dynamics of the system because they help characterize the degree of deformation of the molecules. The dihedral-angle distribution tracks the presence of dihedral

Np

Ugr =

⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σij σij ⎢ ∑ ∑ 4εij⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ r ⎝ rij ⎠ ⎦ i=1 j>i ⎣⎝ ij ⎠ N

(11)

where Np is the total number of adlayer particles (atoms) in the sample, can be used to characterize changes in commensurability, 989


Article


Figure 9. Intramolecular pair correlation functions gintra(r) for the (top) upper layer and (bottom) lower layer at selected temperatures.

bilayer system gradually become disordered. OP2, OPnem, and Φ4 in conjunction provide a consistent overview of this process (Figure 3). However, we observed marked differences between the behaviors of the two layers: 1) In the lower layer, a rectangular-centered (RC) solid phase persists until T ≈ 275 K, where Φ4 sharply drops, and a smectic order appears. (2) On the other hand, in the upper layer, the Φ4 decay is gradual, indicating a more progressive destruction of the crystalline order and formation of a smectic phase. Both layers remain orientationally ordered (OPnem and OP2 remain large) until T ≈ 325−350 K, where the system melts and becomes isotropic. Analysis of the average end-to-end distance (Figure 4) and corresponding distribution (Figure 5) at various temperatures brings more insight into the loss of structural order of the system. We observed a relatively sharp decrease in the average length of the molecules in the lower layer at T = 325−350 K. The upper layer is considerably more disordered from lower temperatures; its transition from the solid to the smectic phase is more gradual (Φ4, Figure 3). Similarly, the double peak observed at intermediate temperatures in the molecular end-to-end length distributions (Figure 5) for the upper layer indicates the appearance of gauche defects and deformation of the molecular backbone. At high temperature, the appearance of a broad, structureless peak is consistent with a lack of confinement, favoring massive formation of gauche defects. This is corroborated by the dihedralangle distributions (Figure 6), where the gauche defect formation

Figure 8. Intermolecular center-of-mass pair correlation functions ginter(r) for the (top) upper layer and (bottom) lower layer of the bilayer at selected temperatures. The inset in the lower panel gives a pictorial representation of diffusive motion supporting lamellar shifting and resulting in skewing of the pair correlation functions.

(gauche) defects. The azimuthal-angle distribution provides information on the molecular orientation of the system and is instrumental in characterizing the range of temperatures for the smectic mesophase. Molecular (center-of-mass) height distributions aid in the understanding of vertical fluctuations through gauche defects, interlayer migration, and patch collapse. We also employed two distinct pair correlation functions, namely, ginter(r) and gintra(r), defined as the probability of finding one molecule center of mass (inter) or atom (intra) at a distance r from another. These correlations help identify the nature of the various phases, from very ordered solidlike phase to more disordered (liquid crystal or liquidlike) phases.

3. RESULTS AND DISCUSSION To provide a visual guide for the discussion of results, we present selected snapshots of the bilayer system in Figure 2. Specifically, there are significant differences between layers in the degree of chain ordering, slight differences in the temperature of the solid− smectic transition, and a high degree of dynamical interlayer correlation after a common melting of the two layers. 3.1. Structural and Orientational Disorder, the SolidSmectic Transition, and the Smectic Phase. With increasing temperature both the lower (contact) and upper layers of the 990


Article

The Journal of Physical Chemistry C is indicated by the presence of peaks at ±120° from the main peak at 180°. These results demonstrate that, for the molecules in the lower (contact) layer, their strong interaction with the substrate and their confinement between the substrate and the upper layer stabilize the ordered phase, suppress the formation of gauche defects, and make the solid−smectic transition sharp. The smectic phase is characterized by the lamellar motion of molecules, which generally conserve their straight shape and parallel orientation until melting. This is illustrated by the behavior of the azimuthal-angle distribution of the long molecular axis (Figure 7): The distinct peaks signal orientational epitaxy, preserved in both layers until the melting transition. At intermediate temperatures, the skewing of the peaks in Figure 8 is indicative of a smectic mesophase: The molecules slide past one another, causing the shift of the lamellae and an increase in the intermolecular center-of-mass distances. The intramolecular pair correlation functions in Figure 9 confirm that the lack of confinement on the upper-layer results in increased chain disorder. 3.2. The Melting Transition and Molecular Order. The behavior of the system with respect to the graphite substrate can be partially characterized by examining the center-of-mass height distributions (Figure 10), which allows us to track how readily

lacking confinement, which allows for more fluctuations, resulting in shallower, broader curves. The upper layer shows also considerably more chain disorder prior to melting, as evidenced directly in the intramolecular pair correlation functions, (Figure 9) which are also consistent with Figures 4−6. The analysis of the average end-to-end distance (Figure 4) and associated distribution (Figure 5) at various temperatures brings more insight into the loss of structural order and reveals a sharp change of molecular length when the temperature increases. At melting the molecules in both layers disorder simultaneously and become more globular in shape. There is also considerable molecules migration between layers (Figures 2 and 10) involving a confluence of promotion into the upper layers and demotion into lower layers. The molecule-graphite van der Waals energies clearly show that the preferential demotion seen at high temperatures is energetically favorable (Figure 11).

Figure 11. (Top) Intermolecular Lennard-Jones energies for the separate layers, between layers, and for the entire bilayer system. (Bottom) Molecule−graphite Lennard-Jones energies for the separate layers as well as for the entire bilayer system at various temperatures. The lines are provided only as guides for the eye.

We also calculate the static structure factor, a quantity which can be directly connected to experiment: S(q) = 1 + 2πρ

Figure 10. Molecular center-of-mass height distributions for the (top) upper layer and (bottom) lower layer at selected temperatures. Distributions are tracked from their initial configuration; that is, the peaks at small/large z in the upper/lower panels correspond to demotion/promotion of molecules between the layers.

∫0

∞

sin(qr ) [g (r ) − 1]r dr qr

(13)

where ρ is the average density of molecules per unit area. The temperature dependence of S(q) confirms also the previous analysis of the behavior of the system, from solid through the melt. As temperature increases, the peaks of S(q) shift, inward (toward lower q) or outward, indicating structural changes in the system (Figure 12). The outward shift is consistent with disorder in the solid and indicates that the distances between molecules

molecules migrate between layers. There is almost no migration of molecules between the layers prior to melting. The asymmetry in the solid height distributions is a result of the upper layer 991


Article


range. Three different regimes of trilayer behavior can be quantitatively and qualitatively delineated: (i) the perpendicular patch, (ii) the collapse of the perpendicular patch to a parallel multilayer, and (iii) the melting of the parallel layered system into the droplet. The initial, low-temperature structure with the third layer introduced as a perpendicular patch (Figure 1) was defined according to experimental observations by Taub and coworkers.14−16,19 With increasing temperature, the bottom of the perpendicular patch bounds and spreads in concert with the formation of a slight concavity on the external boundary of the patch. This is accomplished by dihedral ratcheting through gauche defect formation at the ends of the molecules. The patch spreading is accompanied by highly correlated lamellar motion of the two layers beneath (Figure 13). As a consequence, a significant decrease in the end-to-end length is observed (Figures 14 and 15). The confined middle layer exhibits far fewer gauche defects than the second layer in the bilayer, confirming that

Figure 12. Static structure factors S(q) for the bilayer at (bottom to top) T = 100 K (solid), 275 K (upper layer smectic), 350 K (both layers smectic), and 400 K (both layers liquid). Blue lines indicate peaks that shift outward with increasing temperature, red lines are for peaks that shift inward, and the green line shows an invariant trough that does not shift.

are decreasing (together with decreasing end-end length) at solid-nematic transition. The inward shift is associated with increasing distances between molecules at melting. When comparing the behavior of the bilayer structure to that of the monolayer,10 the following conclusions can be drawn: In the bilayer system, the (lower) layer confined between the substrate and the upper layer shows less disorder than the monolayer. The upper layer shows more disorder, due to (i) the lack of confinement, (ii) the larger distance from the graphite surface, and (iii) the interaction with a dynamic underlayer. The chain disorder that is the driving mechanism of the solid− smectic transition occurs at slightly different temperatures for the two layers. As the molecular motion between layers is highly correlated, the layers melt at the same temperature. It becomes difficult to compare the behavior of the two systems past the point where the upper layer preferentially demotes, but it is clear that the disorder present in the upper layer has a dramatic effect on that of the lower layer. 3.3. The Trilayer System. Figure 13 shows snapshots of all layers of a trilayer system, over the 100−350 K temperature

Figure 14. Atomic height distributions for each layer in the trilayer at various temperatures. Bottom layer, black; middle layer, blue; top patch, magenta.

Figure 13. Snapshots of the (top to bottom) top, middle, and lower layers at (left to right) T = 100, 200, 250, 300, and 350 K. Atoms colored red signify that they are significantly higher off the surface relative to their respective layer. 992


Article


both lower layers show fewer gauche defects than the monolayer at the same temperature. The melting of the system is initiated by local (in-plane) deformation of the intramolecular bonds. The strong interlayer correlations cause the melting temperatures for the two layers in the bilayer to match, even if, globally, the melting is smeared out over a range of temperatures. For the trilayer patch, although the dihedral defects are partially suppressed in the lower two layers, their lamellar motion is the main mechanism for spreading of the perpendicular patch bottom through ratcheting as well as its eventual collapse. Our simulations do capture the characteristics of system’s phases and important dynamics that characterize experimentally observed phase transitions. To reach this level of agreement when studying the systems showing very slow, diffusive behavior, relatively long (∼10-ns) simulations times were required.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Acknowledgment is made to the donors of The American Chemical Society Petroleum Research Fund (PRF43277-B5) and the University of Missouri Research Board and Research Council. The authors are grateful for useful discussions with Haskell Taub.

Figure 15. (Top) Average end-to-end molecular lengths for each layer in the trilayer at various temperatures and (bottom) tilt order parameter OPtilt for each layer in the trilayer at various temperatures.

■

confinement has more of an influence on the structure of the layer than the distance from the graphite substrate. As temperature increases, the patch remains vertical (Figures 13−15) until T ≈ 250 K and then collapses, forming a flat (parallel) multilayer structure (Figures 13−15). Melting into a droplet takes place at ∼300 K. Notice that, at 350 K, the middle layer expands, pushing the top layer upward (Figures 13 and 14), which occurs in concert with pronounced interlayer migration when melting occurs and is signaled in the height distributions as well as the tilt order parameter (Figures 14 and 15).

REFERENCES

(1) Tans, S. J.; Devoret, M. H.; Dai, H.; Thess, A.; Smalley, R. E.; Geerligs, L. J.; Dekker, C. Individual Single-Wall Carbon Nanotubes as Quantum Wires. Nature 1997, 386, 474−477. (2) Alivisatos, A. P. Semiconductor Clusters, Nanocrystals, and Quantum Dots. Science 1996, 271, 933−937. (3) Geim, A. K.; Novoselov, K. S. The Rise of Graphene. Nat. Mater. 2007, 6, 183−191. (4) Baughman, R. H.; Zakhidov, A. A.; de Heer, W. A. Carbon NanotubesThe Route Toward Applications. Science 2002, 297, 787− 792. (5) Niemeyer, C. M.; Mirkin, C. A. Nanobiotechnology: Concepts, Applications and Perspectives; John Wiley & Sons: New York, 2006. (6) Barth, J. V.; Costantini, G.; Kern, K. Engineering Atomic and Molecular Nanostructures at Surfaces. Nature 2005, 437, 671−679. (7) Schull, G.; Douillard, L.; Fiorini-Debuisschert, C.; Charra, F.; Mathevet, F.; Kreher, D.; Attias, A.-J. Selectivity of Single-Molecule Dynamics in 2D Molecular Sieves. Adv. Mater. 2006, 18, 2954−2957. (8) Arrigoni, C.; Schull, G.; Bléger, D.; Douillard, L.; FioriniDebuisschert, C.; Mathevet, F.; Kreher, D.; Attias, A.-J.; Charra, F. Structure and Epitaxial Registry on Graphite of a Series of Nanoporous Self-Assembled Molecular Monolayers. J. Phys. Chem. Lett. 2010, 1, 190−194. (9) Connolly, M. J.; Roth, M. W.; Gray, P. A.; Wexler, C. Explicit Hydrogen Molecular Dynamics Simulations of Hexane Deposited onto Graphite at Various Coverages. Langmuir 2008, 24, 3228−3234. (10) Firlej, L.; Kuchta, B.; Roth, M. W.; Connolly, M. J.; Wexler, C. Structural and Phase Properties of Tetracosane (C24H50) Monolayers Adsorbed on Graphite: An Explicit Hydrogen Molecular Dynamics Study. Langmuir 2008, 24, 12392−12397. (11) Roth, M. W.; Pint, C. L.; Wexler, C. Phase Transitions in Hexane Monolayers Physisorbed onto Graphite. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 155427.

4. CONCLUSIONS We have discussed here the structural evolution of alkane C24H50 multilayer systems adsorbed on graphite. Our calculations revealed, for the first time, the details of the nanoscale mechanism of layer formation in these systems. The applied large-scale molecular dynamics simulations allowed us to disentangle the roles of intramolecular and intermolecular degrees of freedom in the mechanism of the structural transformation. We show that the main driving force of the structural transformations of the multilayer system is the deformation of the alkane chain. However, each layer behaves in a different way, which is a typical property at the nanoscale. In fact, the local dynamic properties of each layer are correlated as a result of interactions with the neighboring layers and confinement effects. Both factors determine the amount of gauche defects. Gauche defects, when induced in the first parallel layer, are responsible for driving the solid−nematic transition. In the bilayer, the lower layer exhibits fewer defects than the monolayer at the same temperature, and the upper layer exhibits more. In the trilayer, 993


Article

The Journal of Physical Chemistry C (12) Pint, C. L.; Roth, M. W.; Wexler, C. Behavior of Hexane on Graphite at Near-Monolayer Densities: Molecular Dynamics study. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 085422. (13) Wexler, C.; Firlej, L.; Kuchta, B.; Roth, M. W. Melting of Hexane Monolayers Adsorbed on Graphite: The Role of Domains and Defect Formation. Langmuir 2009, 25, 6596−6598. (14) Hansen, G. Y.; Herwig, K. W.; Matthies, B.; Taub, H. Intramolecular and Lattice Melting in n-Alkane Monolayers: An Analog of Melting in Lipid Bilayers. Phys. Rev. Lett. 1999, 83, 2362. (15) Diama, A.; Simpson, M.; Taub, H.; Hansen, F. Y.; Dimeo, R.; Neumann, D.; Herwig, K. W.; Volkmann, U. G. Molecular Diffusive Motion in a Bilayer Fluid of Tetracosane Adsorbed on Graphite. Bull. Am. Phys. Soc. 2004, 49, 1345. (16) Hansen, F. Y.; Criswell, L.; Fuhrmann, D.; Herwig, K. W.; Diama, A.; Dimeo, R. M.; Neumann, D. A.; Volkmann, U. G.; Taub, H. Intramolecular Diffusive Motion in Alkane Monolayers Studied by High-Resolution Quasielastic Neutron Scattering and Molecular Dynamics Simulations. Phys. Rev. Lett. 2004, 92, 046103. (17) Matthies B., Diffraction studies of n-alkane fluids adsorbed on graphite. Ph.D. Dissertation, University of Missouri, Columbia, MO, 1999. (18) Pint C. L. Molecular simulation of melting in tetracosane (C24H50) monolayers and bilayers on graphite. 2006, arXiv.org e-Print archive. http://arxiv.org/ftp/cond-mat/papers/0602/0602478.pdf (accessed Nov 17, 2015). (19) Mo, H.; Taub, H.; Volkmann, U. G.; Pino, M.; Ehrlich, S. N.; Hansen, F. Y.; Lu, E.; Miceli, P. A novel growth mode for alkane films on a SiO2 surface. Chem. Phys. Lett. 2003, 377, 99−105. (20) Berman, H.; Henrick, K.; Nakamura, H. Announcing the Worldwide Protein Data Bank. Nat. Struct. Biol. 2003, 10, 980. See also http://www.wwpdb.org/ (accessed Nov 17, 2015). (21) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4, 187. See also http://www.charmm.org (accessed Nov 17, 2015). (22) Kale, L.; Skeel, R.; Brunner, R.; Bhandarkar, M.; Gursoy, A.; Krawetz, N.; Phillips, J.; Shinozaki, A.; Varadarajan, K.; Schulten, K. NAMD2: Greater Scalability for Parallel Molecular Dynamics. J. Comput. Phys. 1999, 151, 283. See also http://www.ks.uiuc.edu/ Research/namd/ (accessed Nov 17, 2015). (23) Firlej, L.; Kuchta, B.; Roth, M. W.; Wexler, C. Molecular Simulations of Intermediate and Long Alkanes Adsorbed on Graphite: Tuning of Non-Bond Interactions. J. Mol. Model. 2011, 17, 811−816. (24) Peters, G. H. The Molecular Dynamics Simulations of the Melting of a Hexane Bilayer. Surf. Sci. 1996, 347, 169.

994


Simulation and Characterization of Tetracosane on Graphite

Recommend Documents