On the Consistency of the Exfoliation Free Energy of Graphenes by Molecular Simulations

Monolayer graphene is now produced at significant yields, by liquid phase exfoliation of graphites in solvents. This has increased the interest in molecular simulation studies to give new insights in the field. We use decoupling simulations to compute the exfoliation free energy of graphenes in a liquid environment. Starting from a bilayer graphene configuration, we decouple the Van der Waals interactions of a graphene monolayer in the presence of saline water. Then, we introduce the monolayer back into water by coupling its interactions with water molecules and ions. A different approach to compute the graphene exfoliation free energy is to use umbrella sampling. We apply umbrella sampling after pulling the graphene monolayer on the shear direction up to a distance from a bilayer. We show that the decoupling and umbrella methods give highly consistent free energy results for three bilayer graphene samples with different size. This strongly suggests that the systems in both methods remain closely in equilibrium as we move between the states before and after the exfoliation. Therefore, the amount of nonequilibrium work needed to peel the two layers apart is minimized efficiently.


Introduction
Graphenes are two-dimensional, single-layer carbon nanosheets with unprecedented physical, mechanical, optical and electronic properties [1][2][3][4]. They are classified with different layered materials like boron nitrides, metal oxides, dichalcogenides and the recently introduced class of metal organic framework nanosheets (MONs) [5][6][7]. Most of their fascinating properties are attributed to the atomic-scale thickness, the continuous 2D connectivity and the ultra-large specific surface area. Graphene nanosheets can be produced by liquid-phase exfoliation and dispersion of graphites [8][9][10][11][12]. This is achieved by chemical functionalization and sonication of graphitized matrices in the presence of certain solvents. Highly polar solvents can destroy the graphitic structure or leave defects and functional moieties on the dispersed layers. By definition, the defects and functional groups change the intrinsic atomic structure and the electrical and mechanical properties of graphene. Graphene oxide, for example, is an insulator rather than a semi metal, and therefore, it is conceptually different than graphene. The special properties of graphenes are also expected to change as the number of layers increases. This is important, because partial exfoliation may release multilayers and clusters instead of large scale of graphene monolayers. In addition, the exfoliated monolayers tend to aggregate into multilayer configurations within the solvent. Nowadays, several high-yield exfoliation methods have been reported that produce unoxidized and defect-free graphene. The material is used to produce graphene-based composites or films, a key requirement for applications such as thin-film transistors, conductive transparent electrodes, photovoltaics and biomedical implants [13][14][15]. The availability of high-quality graphene samples has increased the interest for explicit molecular simulations of the exfoliation processes and relevant applications [16][17][18][19].
The energy required to exfoliate graphene is balanced by the solvent-graphene interactions [9]. Solvent molecules enter the inner core of the graphite and break the Van der Waals forces between the layers. In order to simulate exfoliation, we configure paths either by pulling one layer at a distance from the remaining structure or by breaking the interactions of the layer in small changing steps [20,21]. If the steps are effectively small, the changes on the system are considered reversible, and the reversible work needed to transform one state into the other is equal to the free energy difference of the states before and after exfoliation. The free energy difference between two states, labeled A and B, is given by ∆F = F A − F B . ∆F is related to the ratio of the partition functions of the states [22]. If the states lie far apart in phase space, the estimation of this ratio is intractable [23,24]. The reversible path serves to overcome this hurdle. It provides a connection between the two states so that we can evaluate the changes in F using thermodynamic integration [25][26][27]. That is, we split the overall perturbation into small perturbation steps for which the phase spaces continuously overlap. In this respect, the ratio of the partition functions for the consecutive perturbation steps is calculated more easily [28,29].
We may configure various paths to transform state A into state B [30,31]. A common approach is to use a reaction coordinate to move reversibly from the A-like state to a B-like state. Another approach is to modify the system's Hamiltonian. In this case, we mix the energies or the parameters of states A and B, according to a decoupling parameter, λ, which varies from 0 to 1 [32][33][34]. In many circumstances, it is convenient to introduce an intermediate state (or states) labeled C and evaluate the free energy difference by subtracting the free energies with respect to the state C: [35][36][37][38][39].
In a previous work, we employed one method to compute the exfoliation free energy of graphenes [40]. In brief, we considered a normal and a shear reaction coordinate to exfoliate (dissociate) a graphene layer from a bilayer configuration in an aqueous solution. We computed the free energy differences using umbrella sampling in small steps along the exfoliation paths. We reported that the free energy difference was greater when the exfoliation was coordinated on the shear than the normal direction. This was attributed to the awareness that the shear exfoliation of graphenes is reversible whereas the normal exfoliation is not. Notably, this outcome is in accordance with experiment. For instance, ball milling utilizes high shear force to delaminate layered materials and generate 2D nanosheets [41,42]. In this regard, molecular simulations succeed to realize a proof of principle. That is, the slip between the layers takes place in the in-plane direction, under the effect of shear force to yield free-standing graphene monolayers.
In this work, we compute the exfoliation free energy of graphenes using decoupling molecular simulations. The simulations are performed in two stages. In the first stage, we decouple the Van der Waals interactions of a single layer starting from a bilayer graphene configuration, in the presence of saline water. In the second stage we decouple the interactions of the same layer starting from a single layer configuration. The exfoliation free energy is computed by subtracting the free energy differences between the two stages. We consider three bilayer graphene samples with small, medium and large size. Then, we compare the free energy estimates of the decoupling simulations with those of the umbrella sampling in which, the same graphene nanosheets dissociate on the shear direction in respect to the bilayer plane.

Materials and Methods
Solvation free energies were computed by simulation by decoupling a nanoparticle from the solvent by thermodynamic integration, using the idendity, where H is parameterized Hamiltonian, and λ is the decoupling parameter, which parameterizes the atomistic interactions between the solvent and the nanoparticle [43,44]. The coupled state (λ = 0) corresponds to a simulation where the nanoparticle is interacting fully with the solvent, whereas the uncoupled state (λ = 1) corresponds to a simulation where the nanoparticle is not interacting with the solvent. Considering the van der Waals interaction between two atoms i and j, the λ-dependent Lennard Jones expression takes the form, where r ij = |r i − r j | is the interatomic distance, ij and σ ij are the Lorentz-Berthelot combination parameters. We set α = 0.5 and n = 1. The term, aλ n , at the denominators makes the potential convergent as r → 0, for λ < 1. This is useful especially for advanced steps of the decoupling (0.7 < λ < 1) as it improves the sampling at distances close to σ ij from the atoms of the solute [32,45,46]. To configure the bilayer graphenes we considered two copies of a single layer graphene. The single layers were all-carbon sheets, having a honeycomb pattern of carbons and peripheral hydrogens. The layers were planar and square. We generated the coordinate files of the graphene layers using the BuildCstruct script (http://chembytes.wikidot.com/ buildcstruct, accessed on 28 July 2021) [47]. We built three graphene layers with edges 1.0 nm, 1.5 nm and 2.5 nm. Using these structures, we configured three bilayer graphenes with small, medium and large size, respectively. We placed the copied layers in a parallel orientation, at a distance 0.3 nm from each other. We modeled the sp 2 carbon atoms on the basis of the OPLSAA references of naphthalene and of aliphatic carbons [48][49][50][51]. We modeled hydrogen interactions using the interaction parameters of benzene hydrogens. The structures were uncharged. We used the GROMACS package, version 2018 to build the simulations [52]. All simulations were submitted to the high performance computing services of the Greek National Infrastructure for research and technology, GRNET-ARIS.
We placed the bilayer graphene at the center of a cubic box. We set the nearest edge of the cube at 1.5 nm from the outermost atom of the bilayer. We set the cut off radius to 1.4 nm to be consistent with the minimum image convention. We used periodic boundary conditions in all directions. We solvated the simulation box with simple point charge (SPC) water and we added 100 mM NaCl. Ions were interacting fully with the SPC molecules. They were not interacting electrostatically with the bilayer, because the graphenes were uncharged. Complementary, we performed equivalent simulations for single layer graphene samples with small, medium and large size. We followed the same simulation protocol as for the bilayer graphenes, only that we used a single layer instead of two. The dimensions of the simulation cubes for all the studied systems and the number of molecules described therein, are listed in Table 1.
Before the decoupling simulations, the system energy was relaxed using a steepest descent minimization, followed by an equilibration in the NPT ensemble over 1 ns. We used the Berendsen coupling to regulate the temperature at 310 K and the pressure at 1 bar. The final configuration of the last NPT run, was used as starting configuration for the decoupling simulations. The bilayer graphene and the solvent (water and ions) were coupled to separate coupling paths. We used the Nose-Hoover method for the temperature coupling and the Parinello-Rahman for the pressure. We labeled the two graphene layers using a different index. One layer served as an immobile reference on which, we applied position restraints on the atoms. The other layer was interacting with the reference layer and the ambient molecules on the basis of Equation (2). We used 31 λ steps evenly distributed in [0,1], with dλ = 0.032. Simulations for each of the 31 total λ steps ran over 10 nanoseconds. The thermodynamic integration was computed using the Bennett acceptance ratio (BAR) method [25].
We considered a two-stage process to compute the exfoliation free energy of graphenes. First we deleted a single graphene sheet from the bilayer configuration, then we solvated the graphene sheet back into the solvent. The free energy difference on the first stage was computed by the decoupling simulations of bilayer graphene samples. Likewise, the free energy difference on the second stage was computed by the decoupling simulations of single layer graphenes. The exfoliation free energy was estimated by subtracting the free energy difference on the second stage from that on the first stage.
A different technique to compute the exfoliation free energy of bilayer graphenes is to use umbrella sampling simulations. The umbrella sampling method was discussed in detail in a previous work [40]. There, we employed umbrella sampling to compute the binding free energies of the same, as in this work, three bilayer graphene samples. Although we used the term "binding" instead of "exfoliation" in the definition of the free energy, we refer to the same thermodynamic quantity. In brief, we configured a path in which a single graphene layer dissociated from a bilayer configuration (state A) to an arbitrary far distance (state B) in the solvent. This was achieved, by setting a force to pull the graphene along the coordinate of the path. The other layer (reference) remained at a fixed position. We performed umbrella sampling on a sequence of configuration points along the dissociation path. Each umbrella simulation output a probability distribution function. It was critical, that the neighboring umbrella distributions along the sequence of configurations were overlapping. If two consecutive umbrella distributions did not overlap, we sampled more configuration points until the new distributions bridged the non-overlapping gaps. In this respect, the spacings between the sampled configurations on the path from state A to state B, could be uneven or arbitrarily small. After completing a sufficient amount of umbrella sampling, we gathered the corresponding probability distributions and computed the free energies using the weighted histogram analysis method (WHAM) [53]. Table 1. Model properties of graphene layers with small, medium and large size, edge lengths of the simulation cubes and the number of solvent molecules (water and ions) contained in the cubes for the single and bilayer graphene configurations.

Graphene layer Small Medium Large
Side Length (nm)

Results
In Figure 1, we present configurations from the decoupling simulations performed on three bilayer graphene samples with small, medium and large size. We place the samples at the centre of a cubic box and solvate with saline water. In Figure 1, we decouple only the red graphene layer. The black layer interacts regularly with the environment and serves as a reference. We scale down the Van der Waals interactions of the decoupled graphene by increasing the parameter λ on the basis of Equation (2). At small λ values, the bilayer configurations are stable due to the adequately strong interlayer interactions. The layers preserve a parallel orientation and they can only twist at small angles against each other. With increasing λ, the Van der Waals interactions of the decoupled graphene decrease. This makes the decoupled graphene disconnect from the reference layer. As λ → 1, the decoupled graphene moves freely inside the box. The decoupling interactions become small enough, so that we may to observe configurations in which the two layers intersect. This is due to the term αλ n on the denominators of Equation (2), which makes the Van der Waals potential go to zero in a well-behaved manner with λ. The reference layer remains at a fixed position through the simulations, by imposing position restraints on the atoms. However, we can see in the snapshots of Figure 1 that when the two layers are apart, the reference layer is displaced from the box centre to the opposite direction of the decoupled graphene. This is because we change the periodicity on the representation of the simulated trajectories, setting the center of mass (COM) of both layers at the center of the simulation box.   Figure 2, shows the derivatives of the free energy as a function of the decoupling parameter, λ. Increasing λ from zero to one represents a gradual change from a system with full Van der Waals interactions, to a system where the decoupled graphene is not interacting with the environment. We use 31 λ steps evenly distributed in [0,1]. We consider both single layer and bilayer configurations. In the bilayer configurations only one layer is being decoupled. The free energy values are normalized by the area of the graphenes in order to highlight the finite-size effects of the calculations. Such effects are expected to be appreciable for small interfaces. As a result, small graphenes obtain greater free energy derivatives than the large. We also observe greater energy derivatives for the bilayer than the single layer configurations. This is due to the additional interactions of the second layer. At the low λ range, the energy difference is positive indicating the stability of the starting configurations. At advanced steps of the decoupling, λ > 0.7, the free energy difference is negative. At high λ range, the decoupling interactions are near zero so that the water molecules are allowed to sample the volumes occupied by the atomic sites of the decoupled graphene. This is the same reason why we may observe the layer intersections in Figure 1. The small interaction energy implies a rearrangement of the solvent phase that creates new, lower in energy configurations, making the energy decrease. The bilayer graphenes obtain a first negative energy difference at higher λ, than the single layer graphenes. This is because of the interactions of the reference layer that do not change upon decoupling.
Likewise, the negative energy differences of the bilayers are smaller than those of the single layer configurations. In Figure 3 we show the relative free energies of the decoupling simulations as a function of λ. The energies result from the summation of the energy derivatives shown in Figure 2, up to the value of λ. Using decoupling simulations we simulate the transformation of a single layer (state B) or of a bilayer (state A) into a system where one layer is removed (state C). At full decoupling (i.e., λ = 1) the relative free energy of the single and the bilayer graphenes is expressed by F B − F C and F A − F C , respectively. In the panels of Figure 3, the free energies converge nearly on the same value, regardless of the size of the layers. This is important as it approves that we used appropriate spacings dλ in the thermodynamic integration. The convergence of the free energy differences also confirms the statistical consistency of our calculations, since the contributing integrals have resulted from individual simulations.
The exfoliation free energy of the bilayer graphenes is given by the difference between the relative energies in the two panels in Figure 3, i.e., The free energies are plotted in the left-hand panel in Figure 4, as a function of the decoupling parameter, λ. The free energy is proportional to λ up to a range, where the energy presents a step-like increase. This step is attributed to the peak of the bilayer energy curves at 0.7 < λ < 0.9, in Figure 3. At higher values of λ, the energy values reach a plateau. At λ = 1 we compute comparable exfoliation free energies for the three systems. In general, it is difficult to obtain converged free energy estimates by subtracting the energies of different Hamiltonians. The descrepancies of the free energy are attributed to the different contributions of the solvent-solvent interactions to the overall energy. Solvent-solvent interactions can be negligible compared to the magnitude of the solvent-graphene or the graphene-graphene interactions. However, when the graphene layers are small, the contributions of the solvent-solvent interactions may affect the result. The exfoliation free energies along with the values and statistical variances of F B − F C and F A − F C for the single and bilayer configurations are listed in Table 2.
is the exfoliation free energy for the corresponding bilayer graphene samples computed using umbrella sampling in a previous work [40]. The energies are normalized by the area of the graphene layers, expressed in kJ mol −1 nm −2 . For the sake of comparison, we plot the exfoliation free energies of the bilayer graphenes computed with umbrella sampling simulations, in the right-hand panel of Figure 4. The umbrella simulations are detailed in a previous work [40]. The exfoliated graphene is pulled on the shear direction relevant to the initial plane of the bilayer. It is pulled over 0.5 ns using a pull rate 10 nm ns −1 , so that a final COM distance of 5 nm between the layers is achieved. At this distance, the layers do not interact with each other, so that we may assume the layer exfoliated. The free energies are expressed in Figure 4 as a function of the center of mass (COM) distance between the layers. The free energies obtain a plateau when the interlayer distance becomes greater than the lateral size of the graphenes. Pulling further the graphene, does not change the free energy of the system. The umbrella simulations give comparable free energies with the decoupling simulations for the same graphene samples. This means that the two methods, i.e., shear pulling umbrella sampling and decoupling of van der Waals interactions of the graphene layers, are similar in terms of reversibility.

Discussion
Free energy calculations are independent of the simulation protocol that is used. However, the level of precision does depend very much on the choice of the transformation path [54,55]. Different configurations having an absolute free energy F A should on application of an adiabatic transformation all end up with the same free energy F B . If this condition is not satisfied the transformation has been carried out too rapidly and the transformation strictly speaking is not adiabatic. The states, A and B, before and after exfoliation are separated by a free energy barrier. With decoupling simulations, we cross the barrier by introducing an intermediate state C. With umbrella simulations, we cross the barrier with more umbrella sampling. Nevertheless, there are two distinct types of barriers. The first type is is due to interactions of the perturbed degrees of freedom i.e., those of the exfoliated layer with itself and the environment. The other type of energy barrier is attributed to the interactions of the non perturbed degrees of freedom, i.e., those of the solvent (water and ions). These interactions cannot be smoothed by adding intermediate states or with additional sampling. In this respect, the rearrangement of the solvent molecules can increase the system's entropy affecting the free energy calculations.
Apart from the solvent phase perturbations, discrepancies of the free energy are often attributed to the statistical variance of the molecular simulations. We can see in Table 2 that the variances are small, however, they do not reflect the level of accuracy of the estimates. The error of a free energy computed by molecular simulation, ∆F sim , in respect to that of an experiment, ∆F exp , is given by This equation involves only the average volume V λ=0 of the fully coupled system, and average volume V λ=1 of the box of pure solvent (or of the solvent and the reference layer) containing the same number of solvent molecules as the coupled system [32,56]. The term with the ratio of the volumes, i.e., the correction, for the single and bilayer decoupling simulations, is tabulated in Table 2. In each case, the correction is between 10% and 40% of the reported variance. Because of this magnitude we can safely neglect this term in our calculations.

Conclusions
The availability of high-yield exfoliation methods to produce isolated, unoxidized and defect-free graphene has increased the interest in fundamental studies on the exfoliation processes by molecular simulations. We compute the liquid-phase exfoliation free energy of three graphene layers with different sizes using two simulation protocols, namely decoupling molecular dynamics and umbrella sampling. Using the two simulations, we estimate very similar free energies regardless of the size of the layers. This reflects the thermodynamic consistency of the methodologies and confirms that we designed the corresponding exfoliation paths to be equally reversible. The modeled systems, containing graphenes and water and ions, are reasonably simple, to be able to obtain good statistical sampling, allowing us to realize a lower bound on the amount of the sampling necessary to simulate more complex structures like functionalized graphenes and other 2D layered materials and metal organic nanosheets.