Simulation of Graphene Nanoribbon Aggregation and Its Mediation by

Article pubs.acs.org/JPCB

Simulation of Graphene Nanoribbon Aggregation and Its Mediation by Edge Decoration Jonathan D. Saathoff and Paulette Clancy* School of Chemical and Biomolecular Engineering, Cornell University, Ithaca, New York 14853, United States S Supporting Information *

ABSTRACT: Large polyaromatic molecules, including synthetic graphene nanoribbons (GNRs), are the subject of considerable interest for a variety of electronic applications. For GNRs in particular, functional groups can be bonded along the ribbon edges to modify their dispersibility, self-assembly behavior, and electronic properties. However, these side chains are usually chosen in a “trial and error” fashion, without an underlying molecular-scale picture of the conformations they will adopt in solution and the resulting influence of such structures on macroscopically observable phenomena, particularly aggregation. In this study, we use molecular dynamics (MD) to predict the behavior of various side chains in different solvents as a means to understand how this influences aggregate morphologies and binding energies. Specifically, oligomeric PEG and n-alkoxy chains of varying lengths and grafting densities are examined in vacuum, water, and N-methylpyrrolidone. Examining the binding energies and side chain dispositions that occur with different sets of parameters allows us to suggest a combination of these variables that will minimize aggregational tendencies for the GNRs. The results underscore the value of molecular-scale computational techniques to understand the aggregational tendencies of 2D materials and guide the design of future polyaromatic edge modifications.

■

polyaromatic organic molecules.17,18 In addition, ligands can be used on GNRs for fine-tuning electronic properties.19 The flexibility in design of the GNR edges potentially improves the functionality of the material.20 However, it also adds a complicating factor, as the edge-modifying side chains need to be designed carefully to acquire the right balance of physical properties. Currently, side chains are chosen using a combination of physical intuition and trial-and-error. For example, the side chains used to disperse the GNRs in Narita et al.’s work were chosen because they were known to have successfully dispersed hexi-peri-hexabenzocoronenes, and several side chains were tested with the latter compound.12,15 Computationally screening possible edge configurations could make the selection process more efficient and provide design rules as well as the understanding required to make fine adjustments to the material’s properties. For example, methods such as molecular dynamics (MD) could be used to “synthesize” and test new chemistries far easier and faster than in experiments and in which a highly detailed atomistic representation can be built to provide a deeper understanding of the system. A Monte Carlo (MC) technique could have been used, but it would not have highlighted the impact of kinetics on the system. Hence, we decided to use MD for our studies. MD has been used previously to elaborate on the results of experimental studies of PAHs and GNRs.12,18 Of particular relevance to our studies, MD have been used to understand

INTRODUCTION Graphene nanoribbons (GNRs) are graphene sheets whose nanoscale widths make them suitable for electronic applications because of the opening of a band gap. GNR widths below 1.5 nm have been predicted to exhibit band gaps above 1.0 eV, comparable to Si.1,2 Very narrow GNRs would have band gaps that would theoretically provide suitable on/off ratios for use in logic devices.3,4 There are many methods to create GNRs in the laboratory, including the use of lithography to shape graphene sheets,5,6 the opening of carbon nanotubes through plasma or oxidative etching,7,8 and even the use of DNA as a template for CVD.9 In addition, GNRs can also be produced by a chemical synthesis route in which a highly conjugated polymer is oxidized to form a ribbon.10−14 This method can be used to produce GNRs with atomically precise edge configurations and small widths of only a couple of nanometers. The small, easily controlled geometry of these GNRs would provide consistent band gaps suitable for digital logic applications, which set them apart from the previously mentioned fabrication methods. These synthetic GNRs are also capable of having identical, regularly spaced ligands appended to the sides as a result of the fabrication process.11,12 These ligands are typically short, oligomeric chains and can be used to mediate the known (unwanted) tendency of GNRs and other polyaromatic hydrocarbons (PAHs) to aggregate, as well as aiding in selfassembly.12,15,16 Müllen and co-workers have made significant progress in developing GNRs and improving the solubility of hexi-peri-hexabenzocoronenes (HP-HBCs).12,15 There are also many other examples of using ligands to aid in dispersing large, © XXXX American Chemical Society

Received: October 9, 2014 Revised: December 22, 2014

A

DOI: 10.1021/jp510203j J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

the literature,28,29 we set the bending rigidity of the modeled graphene used in this paper to 34.3 kcal/mol to match existing ab initio calculations because the procedure used to calculate this value was most straightfoward.30 To do so, we adjusted the OPLS-AA force field by altering the dihedral coefficient of the C−C−C−C torsional term from 7.7 to 5.05 kcal/mol. This alteration was needed to obtain more correct binding energies for layered structures and those in aggregated structures with high grafting densities. Further details about this parametrization can be found in the Supporting Information. The SPC/E model was used to represent water molecules,31 and the SHAKE algorithm was applied to constrain the bonds and angles in each water molecule.32 In all the MD simulations, 2 fs time steps were used. Nosé− Hoover thermostats and barostats maintained the temperature and pressure at 300 K and 1 atm, respectively. In vacuum, an NVT ensemble (constant number of atoms, constant volume, and constant temperature) was employed, with a cutoff for intermolecular interactions set at 15 Å. For the solvent calculations, an NPT ensemble (constant number of atoms, constant pressure, and constant temperature) was used, with a 10 Å cutoff and a tail correction for van der Waals interactions. In addition, a particle−particle, particle−mesh approximation was used to accurately account for long-range Coulombic interactions.33 Periodic boundary conditions were used in all directions. The GNRs spanned the length of the simulation cell, thus connecting across the periodic boundary, in order to allow us to calculate per-length energy values without the need to consider the effects of the GNR tips (the ends of the ribbons). Figure 1A and Figure 1B show two views of the same simulation box: Figure 1A shows a cross-section of a GNR, while Figure 1B shows the length of the GNR within the periodic box. The GNR is connected to images of itself across the periodic boundaries at the left and right sides of Figure 1B. The box length was chosen to be 4.2 nm long in order to

how different side chains aid in the dispersion of PAHs. For example, Konatham and Striolo used MD to study the dispersibility of large PAHs, or “nanographenes”.21 In their study, they examined to what extent branched and linear alkane ligands disperse nanographenes and then later applied this knowledge to understand the dispersion of nanographenes in polymer matrices to modify thermal conductivities.22,23 In contrast to past MD studies, this work will examine a more diverse set of side configurations that might be used to modify GNRs. Here, potential energies of binding will be calculated rather than free energies, since they are less computationally expensive and permit a post-simulation deconstruction of the relative contributions from side chains, GNRs, and solvent. This sort of computational analysis is relatively fast, which has the advantage of affording us the opportunity to simulate and compare many more different GNR side chain configurations and side chain/solvent systems than could readily be studied with free energy approaches. As a result of the broader scope that our approach allows, we are more able to provide design rules to guide side chain choices in future studies. While entropies cannot be calculated through these simulations, the binding energies are equivalent to enthalpies, which provide valuable information on how GNRs within an aggregate interact with each other and the solvent. We envision that methods such as these can play an important, “triage” first step to identify more promising side chain choices. This smaller number of more promising choices can then be subjected to greater scrutiny by detailed and expensive free energy calculations that will account for entropic effects and, ultimately, tested experimentally. In this paper, we will examine the behavior of GNRs with PEG and n-alkoxy side chains of various lengths and grafting densities and observe their tendency to aggregate in vacuum, water (a poor graphene-dispersing solvent), and N-methylpyrrolidone (NMP, an effective graphene dispersing solvent).24,25 This study will be split into two sections. First, we will describe the disposition of various side chain configurations around the GNRs. Second, we will ascertain aggregate morphologies from steered molecular dynamics (SMD) calculations and calculate their binding energies. This will provide information on the relative ability of the different side chains to affect binding and link this to the side chain dispositions from the first task.

■

METHODS The LAMMPS molecular dynamics package developed by Sandia was used in each of the simulations described below.26 The OPLS-AA (all-atom) force field developed by Jorgensen et al. provided the intermolecular and intramolecular interactions.27 However, we made a few alterations to the force field to better describe interlayer distances, interlayer attractions, and bending rigidity of the graphene, as well as the partial charges on NMP to create the correct dipole moment. These changes were necessary, as the OPLS-AA force field was originally designed to represent small aromatic molecules, such as benzene, which we judged to be insufficiently accurate to describe graphene nanoribbon interactions. To ensure the most accurate graphene−(NMP) solvent and graphene−graphene interactions, we used the partial charges for the NMP solvent and the graphene van der Waals parameters used by Shih et al. to study the dispersion of graphene sheets in several different solvents.25 The results in this study had already been validated against experiment. While there are other estimates available in

Figure 1. Several images representing simulation set-up details. (A) and (B) are cross-section and side views, respectively, of a typical simulation box. The GNR is shown in cyan with blue side chains and surrounded by solvent. (C) shows the two side chain chemistries used here: PEG chains and n-alkoxy chains. (D) shows the three grafting densities used. Substitution sites for the side chains are shown in red. B


Article


Figure 2. Cross-sectional (left) and top (right) views of GNRs with side chain lengths varying from 6 to 18 atoms (bottom axis) for both n-alkoxy (top) and PEG (bottom) side chains for DG = 1. The GNRs in the top views are connected to images of themselves across periodic boundaries on the top and bottom.

aspect of aggregation behavior. In addition to the lack of or difference in solvent, there are three variables under consideration here: chain chemistry (n-alkoxy or PEG), grafting density (DG = 1, 2, or 4), and chain length. Figure 1 shows a typical solvent box used in this study as well as descriptions of the chain chemistries and grafting densities that were applied to the simulated GNRs. Unless otherwise mentioned, we used chain lengths containing 6, 12, and 18 “chain atoms”, with a “chain atom” being either a linked carbon or oxygen atom. The width of the GNRs simulated in this paper is 1.5 nm. This width combined with the armchair configuration would give the hydrogen-terminated GNR a band gap similar to silicon (1 eV),1,2 making it a reasonable GNR width for device applications. Single GNRs in Vacuum. In order to examine GNRs and side chains without the complicating effects of solvent, we studied the behavior of GNRs in vacuum with the lowest grafting density, DG = 1, different chain lengths, and both chemistries (n-alkoxy and PEG). In vacuum, a 1 ns simulation was sufficient to allow us to observe the conformational disposition of the side chains in vacuum. Figure 2 shows representative images of both cross-sectional and top−down views of the GNRs we studied. Both PEG chains and n-alkoxy chains performed similarly. The side chains prefer to adsorb on either face of the GNR, but this requires the side chains to be able to bend. Side chains that are six atoms long (or less) are unable to bend sufficiently to make significant contact with the GNR. Longer chains are able to wrap around the GNR. GNRs with side chains nine atoms long were included to highlight this change. Chains that were wrapped around the GNR desorbed very infrequently. Single GNRs in Solvent. Next, individual GNRs were simulated in both water and NMP solvents. The purposes of these calculations were to understand the disposition of side chains with different edge configurations and solvents, as well as to obtain the single ribbon energies needed to calculate ΔEbind. These solvents were chosen because they are both used frequently for graphene and CNTs, but they are expected to interact very differently with the graphene material. NMP is known to be effective at dispersing large polyaromatic compounds, like graphene,24,38 while water is very poor at dispersing graphene, requiring surfactants or charged function-

prevent individual chains from interacting with themselves while simultaneously minimizing artificial tension or compression enforced by the box dimensions. The length of the GNR was parallel to the y dimension of the box, which was held frozen, while the x and z dimensions were permitted to deform via the NPT barostat. In reality, GNRs dispersed in solution that are longer than their persistence length are expected to bend and coil.34,35 However, this would be very cumbersome to calculate, especially in solvent (as discussed in later sections), because these bent or coiled GNRs can be hundreds of nanometers long. Instead, placing GNRs across periodic boundaries, similar to the approach taken in modeling carbon nanotubes,36,37 allows us to calculate per-length values that can be applied to the overall ribbon. Most energies reported in this paper are given as “per length” quantities. For example binding energies were calculated, defined as the energy change that occurs when two ribbons come into contact with each other. These energies are calculated by eq 1: ΔE bind =

(E2 + ES) − 2(E1) L

(1)

where ΔEbind is the binding energy, E2 is the energy of a solvent box with two aggregated GNRs, ES is the energy of a solvent box with no GNRs, E1 is the energy of a solvent box with one GNR, and L is the length of the GNR. In each case, the number of solvent molecules in a solvent box is equal. We tested the assumption that ΔEbind was not lengthdependent by simulating hydrogen-terminated, periodic GNRs with lengths of 2.5, 4.2, and 5.9 nm for 12 ns each. The error in ΔEbind between the simulations for different GNR lengths was found to be ±0.49 kcal/(mol·nm), which was small in comparison to typical estimated errors of ΔEbind (±2.0 kcal/ (mol·nm)) and the magnitudes of ΔEbind relevant to this study (10 kcal/(mol·nm) or greater).

■

RESULTS

MD Simulations of Single Ribbons. Our first study examined single GNRs with various side chain configurations in vacuum, NMP, and water. This permitted us to understand the dispositions of the side chains in solvent, which is an important C


Article


Figure 3. Spatial distribution maps of the cross sections of OC11H23 (n-alkoxy) in water (top) and NMP (bottom) at all three grafting densities (increasing from DG = 1 on the left to DG = 4 on the right). The color represents the relative probability of finding a backbone atom at a point on the map, with white corresponding to the site most visited by a chain atom and dark red the least visited. The location of the GNR itself is shown in cyan.

Figure 4. Spatial distribution maps of the cross sections of PEG4 in water (top) and NMP (bottom) for all three grafting densities (key as in Figure 3).

D


Article

The Journal of Physical Chemistry B alities.39−41 With both solvents being polar, and using nonpolar alkane-tipped n-alkoxy chains and polar PEG chains, it was expected that the n-alkoxy chains should interact less favorably in both solvents than the PEG chains. This means we have different levels of solvent “quality” for four different combinations of side chains and solvents. For example, using water and n-alkoxy chains would represent a poor solvent for both the GNRs and side chains. Overall then, the specific choice of solvent and side chains provides a greater diversity of outcomes than in the vacuum case. All simulations were equilibrated for 2 ns, following which, data were collected over an additional 12 ns. This length of time for data collection was sufficient to reduce sampling error to approximately ±2.0 kcal/ (mol·nm). To describe the position of the side chains for the duration of the simulation, spatial distribution maps (commonly known as “heat maps”) were generated showing the relative probability of finding a chain atom at a particular point in space in the x and z dimensions. We only recorded the locations of carbon or oxygen chain atoms, excluding the oxygen atom directly bonded to the GNR, which shows little motion. Since the behaviors of the different chain lengths were comparable, results are only shown for one representative chain length, here 12-atom chains. Corresponding spatial distribution maps for other cases are included in the Supporting Information. As expected, GNRs with n-alkoxy chains in water show very similar behavior to those in vacuum calculations, as shown in Figure 3. Here, the chains are tightly adsorbed to the surface of the GNR in preference to extending into the poor solvent. However, their behavior is very different in NMP. In this case, most of the side chains point away from the GNR. The grafting density also affects the behavior of the side chains. In NMP solvent, at DG = 1, the chains sample a large section of space around the edge of the GNR, and some side chains periodically adsorb onto the GNR surface. At DG = 2, the side chains are roughly parallel with the GNR faces as the side chains begin interacting with each other. At DG = 4, because of steric congestion, the side chains are oriented in many directions, including several adsorbed on the GNR. While the PEG chains appended to GNRs in water appear to be similar to the vacuum case in Figure 4, details of the adsorption are somewhat different. Although PEG is known to be soluble in water, the chains were observed to adsorb onto the GNR. This behavior has been seen in other work.42 The distribution of the PEG chains is much more diffuse when NMP is the solvent. Here, the side chains prefer to remain along the edges, with a few chains adsorbing onto the GNR. MD Simulations of Two-Ribbon Aggregates. Having observed how side chains behave on single GNRs, we investigated how side chains affected aggregate binding energies and aggregate morphologies. Again, results for vacuum calculations are first shown to describe the system without the presence of solvent. Next, we describe SMD calculations that were used to generate different aggregate morphologies of each GNR variant. We calculated binding energies of these different variants and determined the contributions of side chains and GNRs to these binding energies. This provided a means to detect where changes in binding energy originated within the system: the GNRs or the side chains. MD Studies in Vacuum. In a previous section, we showed that side chains in vacuum and in water tended to wrap around the GNRs. These side chains can have a significant impact on the final aggregate morphology by preventing the GNRs from

coming into direct contact with each other, and hence, the side chains mediate aggregation. Examples of aggregates with 0, 1, and 2 layers of side chains separating the two GNRs are shown in Figure 5 in vacuum. The number of layers of side chains that separate two GNRs in a given aggregate will be called its “layer number” in the rest of this paper.

Figure 5. Cross-sectional views of aggregates formed by GNRs with OC17H35 side chains for DG = 2. The three images show example morphologies of aggregates with varying numbers of side chains sandwiched between the ribbons.

Wrapped side chains rarely desorbed from either face of the GNR. Side chains between the GNRs are highly kinetically trapped because they are effectively adsorbed onto two GNRs, which constrict the side chains between them. This meant that many simulations had to be run for each side chain configuration in order to calculate the likelihood that layered structures would form. Such simulations were initialized by first simulating isolated single GNRs in vacuum for 20 ps, then quickly drawing them together by a retracting harmonic spring. The final aggregates were then allowed to equilibrate for 6 ns to permit the side chains to relax. For each side chain configuration, 0−20 side chains were constrained to remain between the GNRs, while the rest were constrained to remain on the outside. These constraints were released after the GNRs aggregated. A more detailed description of the procedure is given in the Supporting Information. The results of these simulations are tallied in Table 1. The nalkoxy and PEG chains behaved almost identically. In vacuum, their particular chemistry does not play a big role in their dynamics, with the exception of the added flexibility of PEG chains to permit some layered structures to form with six-atomlong chains. Longer side chains can more easily fill the space between the GNRs, while side chains at higher grafting densities have a greater chance of remaining between the GNRs because of the greater number of side chains and the steric congestion along the edges. Solution Studies via a Steered MD Method. In order to observe layered aggregates in solution, we used a steered molecular dynamics (SMD) simulation approach.43 Single ribbons in solvent, discussed in the previous section, were replicated and placed 40 Å away from each other in a face-toface orientation in the solvent. The two GNRs were then drawn together at a pulling rate of 2 Å/ns constrained to approach in a face-to-face manner so that multiple layered structures could be observed in a single run. Different aggregated structures were found by locating local energy minima in the interactions of the E


Article

The Journal of Physical Chemistry B Table 1. Population of Layered Structures with Different Side Chain Configurations no. of layered structures (0-layered no./1-layered no./2-layered no.)a n-alkoxy chain length (chain-atom no.)

PEG chain length (chain-atom no.)

DG

6

12

18

6

12

18

1 2 4

(21/0/0) (21/0/0) (8/13/0)

(7/14/0) (7/14/0) (4/12/4)

(5/9/6) (3/9/7) (2/8/10)

(19/2/0) (16/5/0) (5/16/0)

(7/13/0) (7/14/0) (2/14/2)

(4/11/4) (3/10/4) (2/8/9)

a Twenty-one simulations were run in total for each side chain configuration. Aggregate structures that did not clearly fit into one of the three categories were not included in the tally.

demonstrating that both solvents are relatively poor for the nalkoxy chains and not the PEG chains. Of the variables we investigated, we found that changes in the grafting density have the most significant impact on the overall binding enthalpy. At DG = 1 in both Figures 7 and 8, there are very little differences in binding energy with increasing chain length. The side chains are spaced far enough apart from each other in this case that they cannot interact with each other. At the other end of the spectrum, at DG = 4, there is a significant positive contribution to the binding enthalpy from the graphene component of the system. Here, there is a great steric congestion along the edges of the GNRs due to the crowded side chains. The GNRs in this case must contort in order to come into contact with each other. Also, because of the contortion, a smaller area of the GNRs can come into contact with each other than in lower grafting density cases. The more complex trends found in Figure 8 reflect the insolubility of the n-alkoxy chains in combination with the tendency of side chains to wrap around the GNRs in water, which makes interpretation more difficult. In part, this is due to kinetic trapping of the side chains wrapped around the GNRs, which causes a significant amount of “noise” in the data. Some of this is evident in Figure 8E and Figure 8F, where within the same aggregate there are side chains that are wrapped and disordered and others that are well-ordered, much like the case of NMP in Figure 7. Binding Energies for Single-Layered Structures. Aggregates separated by one layer of side chains can be important where kinetic trapping makes it difficult to achieve the zero-layered aggregate or when the binding energy of the single-layered aggregate is lower than the zero-layered case, making it more stable. We observed aggregates to form with multiple separated side chains, as in the vacuum, but only single-layered structures could be consistently resolved with SMD. The energies were calculated as in the previous section: equilibrating for 1 ns and collecting data for 12 ns. But unlike the vacuum case, only one simulation was performed for each side chain configuration. GNRs with long side chains (12- to 18-atoms long) formed layered structures in water. However, in NMP, only OC17H35 side chains with DG = 4 produced a layered structure. Steric congestion along the edges, combined with the length of the side chains, caused a significant number of side chains to adsorb onto the GNR, making it possible for a layered structure to occur (Figure 3). Figure 9 gives the binding energy for each of the layered structures as well as binding energies of their nonlayered counterparts. The energy profile is complicated for the n-alkoxy case in water (Figure 9A) because of the combination of the nalkoxy’s insolubility and propensity for wrapping, which causes kinetically trapped states as discussed in the previous section. The n-alkoxy-decorated GNRs in NMP and the PEG-decorated

two ribbons (including interactions of the GNRs and the side chains) and correlated to visualizations of the aggregates that formed (Figure 6). The peaks between minima indicated

Figure 6. Ribbon−ribbon interactions in an SMD simulation. The two minima correspond to aggregates where the GNRs are in contact with each other (at a separation of 3.6 Å) or where they are separated by a single layer of side chains (at a separation of 7.5 Å).

distances where the side chains were being forced out from between the GNRs. In some cases, because of sluggish kinetics, the chains between the GNRs were not pushed out by the SMD procedure by the end of the simulation. Therefore, a force was temporarily added to the separating chains to “evict” them from the surrounding GNRs before the energies were collected. Finally, binding energies of the “zero-layer” and “single-layer” aggregates were calculated. Binding Energies for Zero-Layered Structures. Using eq 1 in the Methods section, we calculated total binding energies for GNRs with the side chains in direct contact with each other in NMP and water solvents; see Figures 7 and 8. These figures also show how the total energies can be broken down into components arising from the GNR and the side chains. The more negative the energy value, the stronger the binding (and the greater is the tendency for aggregation). We also present the conformations of the side chains corresponding to key points in these graphs. It is apparent that the side chains can adopt a wide variety of morphologies and that these morphologies have a direct bearing on the binding energy. There are a couple of obvious trends that reflect the reasoning behind choosing particular variables to be tested. First, the binding energies are generally lower for GNRs in water than NMP, reflecting that GNRs in NMP are more weakly bound. Also, the n-alkoxy side chains generally have a negative contribution to the binding enthalpy whereas this is not the case with PEG side chains, further F


Article


Figure 7. Potential energy per unit length for a variety of side chain configurations as a function of chain length in an NMP solvent. Colored lines show the effect of different side chain chemistries and grafting densities (solid lines for n-alkoxy and dashed lines for PEG, while the colors red, blue, and green represent DG values of 1, 2, and 4, respectively). The total binding energies are shown in the top left figure; the contributions to this total energy from the GNRs and the side chains are shown in the middle left and middle right plots, respectively. In addition, several points of interest are labeled in the total binding energy plot, and the respective morphologies of the aggregates are shown on the bottom.

side chains in water show interesting behavior. In both cases, the single-layered structures have lower binding enthalpies than their corresponding zero-layered structures at DG = 4. At this grafting density, the GNRs are highly contorted in the zerolayered state, and the side chains are also forced to interact with each other relatively unfavorably. Maintaining a single layer of side chains between the GNRs allows both the GNRs and side

chains to relax, stabilizing the single-layered structures. A more detailed breakdown of the components to the total enthalpy and renderings of selected structures can be found in the Supporting Information G


Article


Figure 8. Potential energy per unit length for a variety of side chain configurations as a function of chain length in a water solvent (key as in Figure 7). Analogous to Figure 7, this plot encapsulates the effect of all the variables studied in this paper for the behavior of GNRs in water.

■

DISCUSSION

makes generalizations difficult, we were able to determine which side chain characteristics are more effective than others at decreasing the binding energy of 2D ribbons like GNRs. With the discrimination that this has provided for side chain selection, we can apply it to screen quickly for side chain choices that can be synthesized in a lab. Overall, we summarize our observations as follows: What is the importance of each of the variables studied here? Our results show that the overall net effect of all the variables

We found that there were many nuanced ways in which side chain conformations are affected by their intrinsic characteristics of length, chemistry, and density and their interaction with their environment, especially the choice of solvent. We observed a rich variety of side chain conformations (Figures 3 and 4 make this apparent) and were able to relate them back to their effect on the binding energy of the nanoribbons. While this complexity of ribbon, side chain, and solvent interactions H


Article


large enough effect that it can even overcome the presence of a poor solvent. The steric crowding at such densities made it necessary for the GNRs to deform significantly if aggregated, making binding significantly less favorable. GNRs with highly branched side chains, which have been shown experimentally to be dispersible, could behave analogously to these GNRs with high grafting density.12 How does the ability of side chains to wrap around the GNRs impact aggregation? In vacuum and poor solvents (like water), side chains prefer to wrap around the GNRs. In some cases, this leads to aggregates where the side chains formed multiple layers between the GNRs. However, these aggregates invariably possessed a higher energy than systems without intervening side chain layers in either solvent, indicating they were only metastable states (see Figure 9). There are, however, cases in which side chain wrapping can be helpful. An interesting case in point occurred only for shorter PEG chains in a water solvent at low to intermediate densities, i.e., DG < 4. Here, the layered structures had a much higher (more positive) binding energy than the nonlayered case. This indicates that, in this particular case, the PEG chains behaved as a protective layer around the GNRs, preventing them from strongly binding together. However, using this general strategy to alleviate aggregation is problematic for several reasons. First, the binding energy of the layered structure is still significantly low, indicating that the binding is still strong. Second, the aggregate with the PEG side chain layer is metastable. If the barrier between the layered and nonlayered cases is small, the existence of this metastable structure should have little impact on aggregation. However, if the barrier is large, the zero-layer aggregate might not be easily obtainable and the single-layer case would make a large impact on the nature of aggregation. Metastable protective layers have been cited in the literature for dispersing graphene. For example, Shih et al. indicated that some of the most successful solvents at dispersing graphene form metastable solvent layers that slow the rate of aggregation.25 While further investigation of these layered structures is outside the scope of this study, it would be interesting to study them further to see what effects they have on aggregation, despite their failure in reducing the binding strength sufficiently to prevent aggregation in this study. Studies such as these could be used to discriminate among a wide range of possible solvent/GNR systems. In this study, it appears that GNRs with PEG chains in NMP with DG = 4 are the most viable for preventing aggregation. The PEG chains interact favorably with the solvent and do not form highly ordered low-energy domains along the sides. The NMP solvent interacts more favorably with the GNRs than water, making the strength of binding less strong. In addition, the high grafting density forces the GNRs to contort greatly when they come into contact with each other, again greatly decreasing the binding energy.

Figure 9. Comparison of binding energies between GNRs with and without intervening layers of side chains in good and poor solvents. (A) shows results for n-alkoxy chains; (B) shows results for PEG chains. Binding energies of GNRs with side chains that are 12 atoms long in H2O are shown in red. Those with side chains that are 18 atoms long in H2O are shown in blue, and those with side chains that are 18 atoms long in NMP are shown in green. Single-layered and zero-layered structures are represented by dashed and continuous lines, respectively.

studied on the binding energy is a delicate balance of multiple factors, making the outcome too difficult to have predicted a priori but uncovered by the comprehensive study here. These results will allow future studies to focus first on more important variables, especially grafting density. For example, side chains only appear to make significant impacts on the enthalpy when DG > 1. Also, layered aggregates are more likely to form with longer chains and higher grafting density. What is the effect of the grafting density? Despite the intricate relationship between the tested variables to affect the binding energy of the system, the grafting density had the single most significant impact on the binding energy. This was first apparent for single ribbons; Figure 4 showed the strong influence of grafting density on side chain interactions. The spatial distribution maps in Figures 3 and 4 showed that steric congestion along the edges forces side chains into more spreadout conformations. But more importantly, the summative plots in Figures 7 and 8 show that it is only the highest grafting density with a suitable chemistry (PEG chains) that leads to the positive binding energy that encourages dispersibility. This is a

■

CONCLUSIONS As mentioned in the conclusion of a recent review paper,44 solution-processed graphene-based materials are still relatively new and their capabilities are only now starting to be understood. We have provided much-needed fundamental understanding of GNR−solvent−ligand interactions that offers considerable assistance to researchers seeking to tailor the chemistry of the ligand, its length, and density, in combination with the choice of solvent. This is clearly relevant to the issue of I


Article


collaborators, Prof. W. R. Dichtel (Chemistry and Chemical Biology, Cornell University, NY) and Prof. Y. L. Loo (Chemical and Biological Engineering, Princeton University, NJ).

GNR dispersibility, which motivated our study, but the impact of this work is far more broadly based, as it speaks to the way that edge modifications can modify the properties of GNRs and, it might be reasonable to suggest, of other 2D materials. We have shown that simulations that probe the energetics of systems, such as those described in this study, can be used to quickly screen for side chain choices that are more likely to prevent aggregation in 2D systems like GNRs. These simulations can be run at a much faster rate than either physical synthesis in the laboratory or a computationally expensive free energy calculation that provides access to both energetic and entropic effects. We recommend a two-step procedure for side chain selection: fast energetic simulations (as here) to down-select among candidates to provide a small more tractable number of candidate side chain designs for subsequent free energy calculations. Such calculations will be the focus of a subsequent publication. We have provided a qualitative understanding of the link between side chain morphology and binding energy. Most importantly, we have provided an understanding of the relative importance of the major design variables that comprise the selection of the most appropriate side chain choice to minimize aggregation. While these individual recommendations (e.g., high density, long chains) might be intuitively reasonable a posteriori, it would have been impossible to predict in advance how the combination of these contributions (density, chemistry, solvent, and chain length) would lead to an overall recommendation of the best choice for the design of side chains for a given goal (here, minimizing aggregation). The results of the simulations provided in this paper provide a relatively fast computational route to predict side chain design/selection in other systems. This will be particularly helpful in cases where the nature of the interactions between different parts of the system (2D materials, side chains, and solvent) is more complicated. For example, these simulations could be used to design dispersible systems with mixed side chain chemistry where intuition is unlikely to be as clear. Computational tools such as these could be used to design systems where using intuition and experience alone would prove to be a formidable or uncertain task.

■

■

(1) Barone, V.; Hod, O.; Scuseria, G. E. Electronic Structure and Stability of Semiconducting Graphene Nanoribbons. Nano Lett. 2006, 6 (12), 2748−2754. (2) Son, Y. W.; Cohen, M. L.; Louie, S. G. Energy Gaps in Graphene Nanoribbons. Phys. Rev. Lett. 2006, 97, 216803. (3) Novoselov, K. S.; Fal’ko, V. I.; Colombo, L.; Gellert, P. R.; Schwab, M. G.; Kim, K. A Roadmap for Graphene. Nature 2012, 490, 192−200. (4) Schwierz, F. Graphene Transistors. Nat. Nanotechnol. 2010, 5, 487−497. (5) Chen, Z.; Lin, Y. M.; Rooks, M. J.; Avouris, P. Graphene NanoRibbon Electronics. Physica E 2007, 40, 228−232. (6) Han, M. Y.; Ö zyillmaz, B.; Zhang, Y.; Kim, P. Energy Band-Gap Engineering of Graphene Nanoribbons. Phys. Rev. Lett. 2007, 98, 206805. (7) Kosynkin, D. V.; Higginbotham, A. L.; Sinitskii, A.; Lomeda, J. R.; Dimiev, A.; Price, B. K.; Tour, J. M. Longitudinal Unzipping of Carbon Nanotubes to Form Graphene. Nature 2009, 458, 872−876. (8) Jiao, L.; Zhang, L.; Wang, X.; Diankov, G.; Dai, H. Narrow Graphene Nanoribbons from Carbon Nanotubes. Nature 2009, 458, 877−880. (9) Sokolov, A. N.; Yap, F. L.; Liu, N.; Kim, K.; Ci, L.; Johnson, O. B.; Wang, H.; Vosgueritchian, M.; Koh, A. L.; Chen, J.; et al. Direct Growth of Aligned Graphitic Nanoribbons from a DNA Template by Chemical Vapour Deposition. Nature Commun. 2013, 4, 1−8. (10) Cai, J.; Ruffieux, P.; Jaafar, R.; Bieri, M.; Braun, T.; Blankenburg, S.; Muoth, M.; Seitsonen, A. P.; Saleh, M.; Feng, X.; et al. Atomically Precise Bottom-Up Fabrication of Graphene Nanoribbons. Nature 2010, 466, 470−473. (11) Yang, X.; Dou, X.; Rouhanipour, A.; Zhi, L.; Räder, H. J.; Müllen, K. Two-Dimensional Graphene Nanoribbons. J. Am. Chem. Soc. 2008, 130, 4216−4217. (12) Narita, A.; Feng, X.; Hernandez, Y.; Jensen, S. A.; Bonn, M.; Yang, H.; Verzhbitskiy, I. A.; Casiraghi, C.; Hansen, M. R.; Koch, A. H. R.; et al. Synthesis of Structurally Well-Defined and Liquid-Processable Graphene Nanoribbons. Nat. Chem. 2014, 6, 126−132. (13) Vo, T. H.; Shekhirev, M.; Kunkel, D. A.; Francois, O.; Guinel, M. J. F.; Enders, A.; Sinitskii, A. Bottom-Up Solution Synthesis of Narrow Nitrogen-Doped Graphene Nanoribbons. Chem. Commun. 2014, 50, 4172−4174. (14) Vo, T. H.; Shekhirev, M.; Kunkel, D. A.; Morton, M. D.; Berglund, E.; Kong, L.; Wilson, P. M.; Dowben, P. A.; Enders, A.; Sinitskii, A. Large-Scale Solution Synthesis of Narrow Graphene Nanoribbons. Nat. Commun. 2014, 5, 3189. (15) Kastler, M.; Pisula, W.; Wasserfallen, D.; Pakula, T.; Müllen, K. Influence of Alkyl Substituents on the Solution- and SurfaceOrganization of Hexa-peri-hexabenzocoronenes. J. Am. Chem. Soc. 2005, 127, 4286−4296. (16) Ito, S.; Wehmeier, M.; Brand, J. D.; Kübel, C.; Epsch, R.; Rabe, J. P.; Müllen, K. Synthesis and Self-Assembly of Functionalized Hexaperi-hexabenzocoronenes. Chem.Eur. J. 2000, 6, 4327−4342. (17) Xiao, S.; Kang, S. J.; Wu, Y.; Ahn, S.; Kim, J. B.; Loo, Y. L.; Siegrist, T.; Steigerwald, M. L.; Li, H.; Nuckolls, C. Supersized Contorted Aromatics. Chem. Sci. 2013, 4, 2018−2023. (18) Yan, X.; Li, B.; Li, L. S. Colloidal Graphene Quantum Dots with Well-Defined Structures. Acc. Chem. Res. 2013, 46, 2254−2262. (19) Jippo, H.; Ohfuchi, M. First-Principles Study of Edge-Modified Armchair Graphene Nanoribbons. J. Appl. Phys. 2013, 113, 183715. (20) Hartley, C. S. Nanoribbons from the Bottom-Up. Nat. Chem. 2014, 6, 91−92. (21) Konatham, D.; Striolo, A. Molecular Design of Stable Graphene Nanosheets Dispersions. Nano Lett. 2008, 8, 4630−4641.

ASSOCIATED CONTENT

S Supporting Information *

GNR force field parameterization, simulation methods for aggregation in vacuum, all single GNR spatial distribution maps, full SMD results in solvent, and binding energy comparisons between select layered and nonlayered structures. This material is available free of charge via the Internet at http://pubs.acs.org.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*Phone: 607-255-8656. Fax: 607 255-9166. E-mail: Paulette. [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge the financial support of the National Science Foundation Nanoelectronics Beyond 2020 (NEB) program, Award CHE-1124754, and a gift from the NRI (Gift 2011-NE-2205GB). We also acknowledge the importance of the guidance and many helpful discussions of our NEB program J


Article

The Journal of Physical Chemistry B (22) Konatham, D.; Bui, K. N. D.; Papavassiliou, D. V.; Striolo, A. Simulation Insights into Thermally Conductive Graphene-Based Nanocomposites. Mol. Phys. 2011, 109, 97−111. (23) Konatham, D.; Striolo, A. Thermal Boundary Resistance at the Graphene−Oil Interface. Appl. Phys. Lett. 2009, 95, 163105. (24) Hernandez, Y.; Nicolosi, V.; Lotya, M.; Blighe, F. M.; Sun, Z.; De, S.; McGovern, I. T.; Holland, B.; Byrne, M.; Gun’Ko, Y. R.; et al. High-Yield Production of Graphene by Liquid-Phase Exfoliation of Graphite. Nat. Nanotechnol. 2008, 3, 563−568. (25) Shih, C. J.; Lin, S.; Strano, M. S.; Blankschtein, D. Understanding the Stabilization of Liquid-Phase-Exfoliated Graphene in Polar Solvents: Molecular Dynamics Simulations and Kinetic Theory of Colloid Aggregation. J. Am. Chem. Soc. 2010, 132, 14638− 14648. (26) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (27) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (28) Arroyo, M.; Belytschko, T. Finite Crystal Elasticity of Carbon Nanotubes Based on the Exponential Cauchy−Born Rule. Phys. Rev. B 2004, 69, 115415. (29) Lu, Q.; Arroyo, M.; Huang, R. Elastic Bending Modulus of Monolayer Graphene. J. Phys. D: Appl. Phys. 2009, 42, 102002. (30) Kudin, K. N.; Scuseria, G. E.; Yakobson, B. I.; C2F, B. S.; Nanoshell, C. Elasticity from ab Initio Calculations. Phys. Rev. B 2001, 64, 235406. (31) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (32) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (33) Hockney, R. W.; Eastwood, J. W. Computer Simulation Using Particles; Taylor & Francis Group: New York, 1988; pp 22−23. (34) Pang, A. L. J.; Sorkin, V.; Zhang, Y. W.; Srolovitz, D. J. SelfAssembly of Free-Standing Graphene Nano-Ribbons. Phys. Lett. A 2012, 376, 973−977. (35) Xu, Z.; Buehler, M. J. Geometry Controls Conformation of Graphene Sheets: Membranes, Ribbons, and Scrolls. ACS Nano 2010, 4, 3869−3876. (36) Uddin, N. M.; Capaldi, F. M.; Farouk, B. Molecular Dynamics Simulations of Carbon Nanotube Dispersions in Water: Effects of Nanotube Length, Diameter, Chirality, and Surfactant Structures. Comput. Mater. Sci. 2012, 53, 133−144. (37) Walther, J. H.; Jaffe, R. L.; Kotsalis, E. M.; Werder, T.; Halicioglu, T.; Koumoutsakos, P. Hydrophobic Hydration of C60 and Carbon Nanotubes in Water. Carbon 2004, 42, 1185−1194. (38) Hughes, J. M.; Hernandez, Y.; Aherne, D.; Doessel, L.; Müllen, K.; Moreton, B.; White, T. W.; Partridge, C.; Costantini, G.; Shmeliov, A.; et al. High Quality Dispersions of Hexabenzocoronene in Organic Solvents. J. Am. Chem. Soc. 2012, 134, 12168−12179. (39) Li, D.; Müller, M. B.; Gilje, S.; Kaner, R. B.; Wallace, G. G. Processable Aqueous Dispersions of Graphene Nanosheets. Nat. Nanotechnol. 2008, 3, 101−105. (40) Lotya, M.; Hernandez, Y.; King, P. J.; Smith, R. J.; Nicolosi, V.; Karlsson, L. S.; Blighe, F. M.; De, S.; Wang, Z.; McGovern, I. T.; et al. Liquid Phase Production of Graphene by Exfoliation of Graphite in Surfactant/Water Solutions. J. Am. Chem. Soc. 2009, 131, 3611−3620. (41) Green, A. A.; Hersam, M. C. Solution Phase Production of Graphene with Controlled Thickness via Density Differentiation. Nano Lett. 2009, 9, 4031−4036. (42) Xu, L.; Yang, X. Molecular Dynamics Simulation of Adsorption of Pyrene-Polyethylene Glycol onto Graphene. J. Colloid Interface Sci. 2014, 418, 66−73. (43) Park, S.; Schulten, K. Calculating Potentials of Mean Force from Steered Molecular Dynamics Simulations. J. Chem. Phys. 2004, 120, 5946−5961.

(44) Cheng, C.; Li, D. Solvated Graphenes: An Emerging Class of Functional Soft Materials. Adv. Mater. 2013, 25, 13−30.

K


Simulation of Graphene Nanoribbon Aggregation and Its Mediation by

Recommend Documents