Principles Governing Control of Aggregation and Dispersion of Graphene and Graphene Oxide in Polymer Melts

Controlling the structure of graphene and graphene oxide (GO) phases is vitally important for any of its widespread intended applications: highly ordered arrangements of nanoparticles are needed for thin‐film or membrane applications of GO, dispersed nanoparticles for composite materials, and 3D porous arrangements for hydrogels. By combining coarse‐grained molecular dynamics and newly developed accurate models of GO, the driving forces that lead to the various morphologies are resolved. Two hydrophilic polymers, poly(ethylene glycol) (PEG) and poly(vinyl alcohol) (PVA), are used to illustrate the thermodynamically stable morphologies of GO and relevant dispersion mechanisms. GO self‐assembly can be controlled by changing the degree of oxidation, varying from fully aggregated over graphitic domains to intercalated assemblies with polymer bilayers between sheets. The long‐term stability of a dispersion is extremely important for many commercial applications of GO composites. For any degree of oxidation, GO does not disperse in PVA as a thermodynamic equilibrium product, whereas in PEG dispersions are only thermodynamically stable for highly oxidized GO. These findings—validated against the extensive literature on GO systems in organic solvents—furnish quantitative explanations for the empirically unpredictable aggregation characteristics of GO and provide computational methods to design directed synthesis routes for diverse self‐assemblies and applications.


DOI: 10.1002/adma.202003213
the characterization of GO so uncertain, relating the structural or compositional properties of GO flakes to their final morphology in composites is a challenging task. Molecular simulation offers a powerful method by which to study these systems, whereby we can control the composition and structure of GO, examine thermodynamically stable morphologies of GO, and view dispersion mechanisms with a high level of fidelity. Our simulations have the ultimate aim to provide computational methods to design synthesis routes for different self-assemblies [3] and applications. [4] In this study, multiscale modelling is used, combining coarse-grained (CG) and atomistic molecular dynamics, to examine the effects of chemical composition (flake oxygen content and its distribution of oxygen) on the properties of graphene and GO composites (morphologies, interlayer distance). As we have shown in our study of polymer intercalation in clay-platelets, CG simulations can approach length and time scales of micrometers and microseconds, respectively, required to observe the intercalation and dispersion processes. [5,6] CG simulations include molecular specificity and, using our recently introduced accurate models of GO, we capture the realistic domains of oxidation and graphitic domains seen experimentally on GO flakes [7] (for more details, see Computational Methods and Supporting Information). More fundamentally, the work described in this article is based on the GraFF forcefield, which we introduced for graphene 2 years ago. [8] We consider four different oxidation states of GO flakes, varying from zero oxidation (graphene) to C:O atomic ratios of 10.0, 5.0, and 2.5. Typical values of GO oxidation found experimentally lie between 4 and 2, with reduced grapheneoxide or modified oxidation methods with reducing reagents greater than 4. [9,10] GO flakes of diameter 10 nm were chosen. Although this is small compared to experimental sizes, it is large enough to contain typical sizes of graphitic domains at low oxidation states (experimentally reported sp 2 cluster sizes range between 2-18 nm [11][12][13] ). Herein, we consider the aggregation tendency of GO flakes, which are related to the size of the graphitic domains, unlike elastic properties, which have a critical length before reinforcement occurs due to "shear lag effects," estimated to be of the order 3 μm for monolayers of exfoliated graphene. [14] In addition, the small percentage of Controlling the structure of graphene and graphene oxide (GO) phases is vitally important for any of its widespread intended applications: highly ordered arrangements of nanoparticles are needed for thin-film or membrane applications of GO, dispersed nanoparticles for composite materials, and 3D porous arrangements for hydrogels. By combining coarse-grained molecular dynamics and newly developed accurate models of GO, the driving forces that lead to the various morphologies are resolved. Two hydrophilic polymers, poly(ethylene glycol) (PEG) and poly(vinyl alcohol) (PVA), are used to illustrate the thermodynamically stable morphologies of GO and relevant dispersion mechanisms. GO self-assembly can be controlled by changing the degree of oxidation, varying from fully aggregated over graphitic domains to intercalated assemblies with polymer bilayers between sheets. The longterm stability of a dispersion is extremely important for many commercial applications of GO composites. For any degree of oxidation, GO does not disperse in PVA as a thermodynamic equilibrium product, whereas in PEG dispersions are only thermodynamically stable for highly oxidized GO. These findings-validated against the extensive literature on GO systems in organic solvents-furnish quantitative explanations for the empirically unpredictable aggregation characteristics of GO and provide computational methods to design directed synthesis routes for diverse self-assemblies and applications.
GO can vary significantly depending on the method used to create it and the extent of oxidation, leading to a wide distribution of flake sizes and composition of functional groups, depending on the production method. [1] The composition may also change over time. [2] It is of no surprise, then, that with edge sites compared to sites on the basal surface (5%) ensures that edge effects are not important in our conclusions and that our results can be extrapolated to larger sizes of GO flakes. Each coarse-grained simulation contains five hexagonal GO flakes with a volume fraction of approximately 3.5 %. For the purposes of sampling and uncertainty quantification, we use ensembles of simulations for each oxidation state consisting of eight replicas starting from a randomly dispersed configuration, and eight from a stacked/aggregated configuration (see Supporting Information). The flakes and initial structures are chosen in Figure 1. Figure 1b also shows the distinct domains of sp 2 and sp 3 phases as created by our GO model builder. [7] Each replica contains a different initial velocity distribution drawn randomly from a Maxwell-Boltzmann distribution, and is simulated using a simulated annealing routine by heating to a high temperature (500 K) followed by cooling down to room temperature (see Computational Methods). This procedure ensures as far as possible that we probe thermodynamic equilibrium properties. [15] A variety of morphologies are found. To categorize the morphologies and to subsequently examine their relationship to the flake composition, we created a morphology definition based on the distribution of nearest neighbor distances for CG flake atoms. The atomic nearest-neighbour distances (d) for each flake CG atom to atoms in other flakes is subdivided into three categories: i) d ≤ 8 Å, the atom is in an "aggregated" environment; ii) 8 Å ≤ d < 18 Å, the atom is an "intercalated" environment; and iii) d ≥ 18 Å, the atom is in an "unbound" environment. The percentage of atoms in each environment is used to define the overall flake morphology. The four flake morphologies are designated as "intercalated," "aggregated," "partially aggregated," and "dispersed"; the relationship to their constituent atomistic environments is shown in Figure 2e, with specification of the morphologies in Figure 1. Starting from a dispersed initial condition (i.e., at t = 0, % intercalated atoms = % aggregated atoms = 0), the most common flake morphology for poly(vinyl alcohol) (PVA) systems is intercalated, irrespective of oxidation state (Figure 2a). The exception is graphene, which is predominantly aggregated, and the system with a C:O ratio of 10.0, which is a mixture of intercalated and partial aggregation. The diagonal line in Figure 2 indicates a perfectly "assembled" structure, that is, all atoms are in either an aggregated or intercalated environment. PVA-graphene (oxide) systems will preferentially self-assemble, with a tendency toward forming intercalated morphologies as the degree of oxidation increases.
Intercalated morphologies are less favored for poly(ethylene glycol) (PEG), and the self-assembled structures contain only aggregated moieties (Figure 2c). The higher the oxidation state of the flake, the greater the tendency for dispersion. Highly oxidized flakes remain dispersed, with only a small amount of intercalation and aggregation, but for flakes with a C:O ratio greater than or equal to 5.0, the flakes assemble over the graphitic regions. Note that this is at a higher oxidation state than the aggregation seen in the PVA systems (at C:O ratio of 10.0).
In the case of graphene systems, for both PVA and PEG polymer, the flakes are predominantly aggregated, the majority of flakes directly interacting with each other via attractive van der Waals interactions, with no polymer resident between the flakes. For GO with a C:O ratio of 2.5 in PVA, the GO flakes lie directly above each other with polymer between the flakes (illustrated in Figure 3b). The density profile of CG atoms perpendicular to the surface in Figure 3c Figure S19, Supporting Information. Note that points may lie on top of each other; for example, graphene in PVA starting from a dispersed initial configuration (blue dots in (a)) comprises 80% aggregated flakes (see Figure S19, Supporting Information). e) Illustration of morphology classification as a function of percentage of atoms in aggregated and intercalated environments. The two simulation snapshots shown are for the PVA-dispersed C:O ratio of 10.0 system. One simulation has formed an intercalated assembled structure (i), while another has aggregated via the graphitic domains on the flake (ii). Polymer molecules are not displayed to aid visualization.
It is well known that PVA forms a hydrogen-bonded network with GO hydroxyl groups. [16] This highly attractive interaction between the polymer and GO oxygen functionality leads to a very dense layer of immobile polymer on each flake, compared to bulk PVA molecules, as shown by the density profile in Figure 3c. When these "effective particles" of GO flake and polymer drift close to each other, the matching of densities produces an effective polymer-mediated attractive force, and a stack of assembled intercalated GO flakes is the result. For PEG, the amount of hydrogen bonding is reduced as it can only act as a hydrogen bonding acceptor, as opposed to PVA which can be both an acceptor and donator, as well as hydrogen bond with itself. PEG polymers on the surface of the GO flake are therefore less tightly bound (Figure 3f), and the polymermediated force between flakes is consequently less. A snapshot of the dispersed GO with a C:O ratio of 2.5 in PEG is shown in Figure 3e.
Our simulations also provide insight into aggregation via graphitic domains. The PVA polymer C:O ratio of 10.0 systems reside in both the intercalated and the partially aggregated regions in Figure 2a. GO sheets have aggregated via the graphitic domains on the surface of the flake in some replicas, while in other replicas, an intercalated structure has self-assembled (see Figure 2-i,ii). Therefore, when the graphitic domains are large enough (in this case, for diameters >8 nm), aggregation is possible in PVA. However, as intercalated structures are also observed in some replicas, the energy of these structures must be comparable. For low oxidation state GO in PEG polymer, aggregation via the graphitic domains occurs for the system of C:O ratio 5.0 leading to partially aggregated morphologies, while in system of C:O ratio 10.0, fully aggregated morphologies arise, the oxidized regions of the flake aggregating along with the graphitic domains.
These findings demonstrate the greater interaction of PVA polymer molecules with the hydroxyl groups on the flakes' surfaces compared to PEG molecules, thereby making the intercalated morphology much more favorable than the aggregated one for PVA GO systems, even via graphitic domains. As a consequence, aggregation occurs for smaller graphitic domains in PEG than for PVA.
In Figure 2b, we show the categories of assembly manifested by the PVA system when the flakes are initially in a stacked configuration. These are all predominantly still assembled, comprising a mixture of aggregated and intercalated arrangements. Visualizing the morphologies formed, we observe these structures are aggregated over areas containing a large number of hydroxyl groups on the GO surface for high oxidation states and graphitic regions for low oxidation states. However, for high oxidation states, we also observe fully intercalated interflake galleries, as well as partially intercalated ones. These simulations illustrate the large kinetic barriers to full intercalation, as the process requires the breaking of hydrogen bonds between flakes before the polymer can move into the interflake gallery. Where the number of interflake hydrogen-bonds is large, this is expected to be a slow process.
Experimentally, the self-assembled intercalated GO structure is also observed in the dispersion behavior of GO in ethanol, which is the monomer unit of PVA. When GO powder is mixed with ethanol (or methanol), the solvent molecules intercalate into the inter-flake gallery, leading to an expansion in flake separation. A wide range of separations has been reported, from approximately 9 to 16 Å, [2,17] the latter corresponding to the bilayer of polymer observed in our PVA simulations. The zeta potential of GO in ethanol is low, indicating the GO flakes are relatively uncharged, [18] analogous to the neutral GO flakes in our simulations. Similarly, most GO flakes produced using the Hummers method [19] have the C:O ratios between 2 and 3, corresponding to our high oxidation state GO flakes. [20] Due to PVA and PEG having high melting points close to their decomposition temperatures, PVA and PEG with GO have been processed mainly from aqueous solutions, which will have different GO structures from those observed in our simulations due to oxidation of the graphene edge groups. The intercalation of the equivalent monomer solvent is therefore the ideal experimental analogue with which to compare our computed structures, and due to the similar structures found, we can be confident that the assembled morphology is the thermodynamically stable state for highly oxidized GO in PVA. The corollary to our conclusion, therefore, is that where GO sheets are found to be dispersed in PVA during experiment, this is a metastable state arising from the processing conditions (for example, an initial aqueous dispersion of GO, followed by sonication, etc). For PEG systems, we can compare with the dispersion characteristics of GO in diethyl ether; we find that the solubility is three times larger than for ethanol, [17,21] although still lower than other aprotic polar solvents, such as dimethylformamide. [22] This is consistent with our findings: There is less thermodynamic driving force for the self-assembled intercalated structure compared to PVA, although a small amount of aggregation can be discerned.
Our discovery of the low-energy assembled structure of GO in PVA also explains the observation that relatively high GO percentage weight composites (5 % by weight, compared to our 3.5%) in PVA exhibit intercalated structures of comparable separation. [23] In addition, lower concentrations of GO/ PVA composites formed from sonication of aqueous GO with dissolved PVA manifest some restacking in field-emission scanning electron microscopy (FESEM) images into groups containing a few sheets. [24] The hydrogen bonded structure observed in our simulations further accounts for the processes that occur during the creation of highly ordered PVA and GO thin films using layer-by-layer self-assembly [24] and vacuum-assisted self-assembly/flocculation. [25] Despite the very different experimental setup-in the latter case, an aqueous dispersion of GO and PVA, with an external force removing the water molecules through filtration via a membrane-the final thin film has a similar interlayer spacing to our self-assembled PVA GO flakes. There has been some debate in the literature as to whether the formation of this highly-ordered, bricks-and-mortar structure is a result of polymer-mediated interaction forces, or due to the external forces. [26] Our work shows that, if the GO flakes are relatively uncharged, the hydrogen bonds between the GO flakes and the polymer drive the self-assembly of the final highly ordered structure. For PEG-based GO thin films, however, this driving force is less and the final structure is likely to be more driven by the external pressure. [25] Our study demonstrates that aggregation over graphitic domains at lower oxidation states is possible, but the critical size of graphitic versus oxidized domains is dependent on the polymer and its interaction with the hydroxyl groups on the flake surface. The degree of aggregation and the assembly of aggregated domains into a supramolecular structure as a function of oxidation state have not been experimentally examined in organic media, but Krishnamoorthy et al. synthesized GO in aqueous solution at a variety of different oxidation states using a modified Hummers method [11] and recorded the variation in aggregation using X-ray diffraction. They found two peaks corresponding to both graphitic aggregation and intercalated water at intermediate oxidation values. Our work shows that simulation can predict whether the thermodynamically stable state involves aggregation, intercalation, or a combination of the two. This insight could in future be employed as a first step in controlling the cooperative behavior of single GO flakes for specific applications.
Understanding the relationship between flake composition and structure is critical for optimizing the performance of graphene (oxide)-based nanocomposites. From our findings, we can specify a set of general rules relating the composition of the flakes to the morphologies formed, when under thermodynamic control, as follows: Rule 1: Large graphitic domains tend to aggregate regardless of the surroundings. The critical size of graphitic domains that produces aggregation depends on the polymer chemistry, however. Rule 2: Graphene-oxide domains will form hydrogen bonds wherever possible. Rule 3: Polymers with the capacity to form strong hydrogen bonds will produce intercalated morphologies. Rule 4: Polymers that are only hydrogen bond acceptors have a lower propensity for forming intercalated morphologies and a greater tendency to aggregate over graphitic domains.
These rules show that graphene-oxide can be regarded as "amphiphilic," with self-assembly through aggregation over the graphitic domains, or through hydrogen-bonded interactions with polymer over the oxidized regions. [27] We conclude that the formation of liquid crystal phases of GO sheets in the dilute or semi-dilute regime should be possible via these interactions, and is not purely due to excluded volume effects. [28] Tuning the balance between the different types of interaction through altering the amount of oxidation of the GO flake or the polymeric medium creates GO materials with desired ordered microstructures, [29] without the need for highly disruptive processing methods such as sonication, which can fracture the flakes. [30] These multiscale simulations furnish a major improvement in our understanding of graphene and GO aggregation and dispersion. Future work will seek to further understand the role of surface functionalization and how the topology of oxidized and unoxidized domains leads to "patchy particles," which could be designed to self-assemble into precise and predictable extended structures. This will require extension of the multiscale simulation scheme to include carboxylic acid sites on the edges of graphene flakes, which are important in determining the dispersion in other media such as water, when such groups are charged at normal pH and wherein electrostatic repulsion is expected to be the dominant factor in creating dispersed morphologies. [31]

Computational Methods
Graphene and GO possess a planar structure. GO flakes contains functional groups on the basal plane and free edges. To simplify the coarse-grained GO models, only oxidation in the form of hydroxyl groups on the basal plane was considered in order to maintain focus on the comparison of oxidation level without variation in oxidation type.
In addition to drawing on the new GraFF potential parameterization, [8] the present work was built upon a program that generates realistic atomistic representations of GO flakes, which used quantum mechanical simulations of reactive intermediates [32] to determine the progress of oxidation on the flake surface. [7] This algorithm recreated the two-phase nature of oxidized and unoxidized graphene domains on GO flakes observed in microscopy experiments, and is a significant improvement on the previous standard of an uncorrelated random distribution of oxidized functional groups, based on the Lerf-Klinowski model. [33] Coarse-Grained Potentials: The multiscale scheme used coarse-grained molecular dynamics simulations parameterized from short timescale classical atomistic molecular dynamics simulations to create interaction CG-MD parameters, where the coarse-grained particles represent groupings of several atoms. [6] To construct the coarse-grained model of GO, a hybrid atomistic-coarse grained approach was used. First, sp 2 carbon atoms of the GO framework had a 1:1 mapping to CG atoms, and interact with other sp 2 carbon atoms using the OPLS potential, [34] in the same manner as atomistic graphene (called CG type C). All remaining atoms were grouped together as follows: a united atom approach was used to generate the edges of the graphene flakes by mapping the edge carbon and hydrogen atoms into a single CG bead (type CH); each hydroxyl group on the graphene surface was represented by a single CG bead that maps to the sp 3 carbon, the oxygen, and the hydrogen atoms (type COH). Interaction potentials not defined by the OPLS forcefield (both bonded and non-bonded) were computed using the iterative Boltzmann inversion (IBI) scheme (see Supporting Information for more details).
To create a CG representation of each polymer, a mapping of each monomer unit of the polymer to a single polymer CG bead was used. The procedure to create the CG interaction parameters for PEG and PVA using the IBI scheme has been described in previous papers. [5,6] In this study, PEG and PVA polymers containing 100 monomer units were used. To represent the interactions between the CG polymer atoms and the graphene surface, interaction potentials were constructed that recreated the atomistic density profiles perpendicular for the graphene surface, calculated from fully atomistic simulations of the graphene-polymer interface using the OPLS forcefield. To construct the CG interaction potentials for the COH CG bead with the C CG bead and with itself, it was matched to atomistic potential of mean force calculations for a single COH group interacting with a pristine graphene flake and a graphene flake containing a COH group, respectively. Further details are available in the Supporting Information.

Coarse-Grained Model Systems:
To construct the initial graphene/GO polymer systems, the GO flakes were created through the algorithm mentioned above. [7] These flakes were then either 1) randomly placed in the simulation box, or 2) started in a stacked configuration with a 6 Å separation for GO models and 3.5 Å for graphene models (see Figure 1). The initial dimensions of the simulation cells were 20 nm × 20 nm × 20 nm, and polymer molecules were subsequently added around the flakes with atomistic resolution at a density of 0.7 g mL -1 . The atomistic systems were subsequently relaxed for 2 ns using NpT conditions and mapped to their CG equivalent. The CG simulations were performed in the NVT ensemble, with the relaxed atomistic lattice parameters lying between 17.5 and 18.5 nm. The different systems are listed in Table S3 of the Supporting Information.
Ensembles and Averages: To increase the sampling, 16 replicas of each system were created, with 8 starting from the randomly dispersed configuration, and 8 from the stacked configuration. These replicas were given a different random seed for the generation of atomic velocities drawn from a Maxwell-Boltzmann distribution. There are, therefore, 16 various replicas for each of the 4 oxidation states (graphene, GO = 2.5, 5.0, and 10.0), leading to 128 simulations in total. Replica simulation building and subsequent analysis used the VECMA toolkit (www.vecma-toolkit.eu/toolkit). All CG simulations were performed using a simulated annealing routine to overcome energy barriers. The simulations were simulated at 500 K for at least 80 ns, and cooled down to 300 K over 30 ns. Due to the softness of the coarse-grained potentials, the dynamics were often much faster than in the all-atom description. In the Supporting Information, this speed-up through a comparison of mean-squared displacement at the atomistic and coarsegrained levels was estimated, with a factor of approximately 3 for PEG and 13 for PVA-based systems, resulting in total rescaled timescales of between 330 and 1430 ns. All subsequent analyses were performed on the configurations generated at 300 K. All simulations were performed using the LAMMPS molecular dynamics code, [35,36] using a RESPA multitimestepping integrator, with 5.0 fs for the non-bonded interactions, and 1.25 or 0.75 fs for the bonded interactions. More details are in the Supporting Information.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.