Simulation Details and Additional Simulation Results
Stabilization of Graphene Sheets by a Structured Benzene / Hexafluorobenzene Mixed Solvent Andrew J. Oyer,† Jan-Michael Y. Carrillo,†,‡ Chetan C. Hire,† Hannes C. Schniepp,∫ Alexandru D. Asandei,†,¶ Andrey V. Dobrynin,†,‡ and Douglas H. Adamson†,¶,* †
Institute of Materials Science, Polymer Program, University of Connecticut, Storrs, CT 06269, ‡Physics Department, University of Connecticut, Storrs, CT, 06269, ¶Chemistry Department, University of Connecticut, Storrs, CT 06259, ∫Department of Applied Science, The College of William and Mary, Williamsburg, VA, 23185
Simulation Details and Additional Simulation Results We have performed molecular dynamics simulations of the adsorption of water, benzene, hexafluorobenzene (C6F6), xylenes, and heptane on graphene sheets. The Generalized Amber Force Field (GAFF)(1) parameters were used for atomistic models of solvents and graphene while water was simulated by using a modified TIP3P model optimized for simulations with Ewald summation.(2) The partial charge distributions for the solvents, except water, were obtained by performing abinitio calculations using the Gaussian 09 (G09) simulation package(3) with 6-31G(d) basis set and B3LYP DFT method. The total potential energy of the system consisted of the bonded, bond angle, dihedral angle, improper angle and non-bonded interaction potentials. The interaction parameters for the van der Waals potential between heterogeneous atomic pairs were calculated as the geometric mean of the interaction parameters for each atom. The default AMBER force field weighing coefficients for pair-wise energy and force contributions were used to account for contributions from the van der Waals and electrostatic interactions. U TOTAL
K r r
2
r
BONDS
eq
K
2
eq
ANGLES
Aij Bij qi q j (SI.1) Vn 2 1 cos n K 12 6 eq Rij Rij DIHEDRALS 2 IMPROPER i j Rij
The simulation box was built using Chem3D,(4) G09, Antechamber(5) and AMBER2LAMMPS python script that is included with LAMMPS.(6) The G09 input file for the solvent molecule, e.g. benzene, was built in Chem3D, then G09 calculations
SI 1
with geometry optimization using AM1 semi-empirical method were performed. The Gaussian output from the calculation was used as an input for Antechamber to determine charges, atom type, bond type, angle and dihedral type assignments.
Figure SI 1 Partial charge distributions used in simulations: water (a), benzene (b), xylene (c), heptane (d) and C6F6 (e). All charges except for water were obtained by using Mulliken population analysis from ab initio calculations with 6-31G(d) basis set and B3LYP DFT method. The water charges were obtained from Price and Brooks.(2)
SI 2
Figure SI 2 Initial system configuration of equimolar mixture of benzene and C 6F6 molecules. The graphene sheet is located at z=0. Hydrogen atoms are shown as green beads, fluorine atoms as yellow beads, aromatic carbon atoms as grey beads, and carbon atoms in the graphene sheet are colored in black.
The AMBER topology file was created by using LEAP that was included in the Antechamber package. The AMBER topology file was converted into a LAMMPS data file using the python script AMBER2LAMMPS. Using output of the AMBER2LAMMPS script as a template, the solvent molecule was replicated and distributed in the simulation box using in-house code. The graphene sheet was modeled as a neutral xy-periodic “macromolecule” using the GAFF definition of the aromatic carbon for the van der Waals interaction parameters. The partial charges of the solvents, except water, were obtained with the Mulliken population analysis from ab initio calculations using G09 with 6-31G(d) basis set and B3LYP DFT method. The results are summarized in Figure SI 1.
SI 3
The NPT ensemble simulations were performed using GPU accelerated LAMMPS code.(7) The equations of motion were integrated by using the velocity Verlet algorithm with a time step of 1.0 fs. The system was periodic in the x, y and z directions. The standard PPPM(8) method with an accuracy of 1.0 × 10-5, and the near-field cutoff set to 10.0 Å, was used to account for contributions from the longrange electrostatic interactions. The graphene sheet was placed at z = 0 Å and spanned the xy-plane. The coordinates of the carbon atoms forming a graphene sheet were fixed during the entire simulation run. Solvent molecules were distributed over the volume of the simulations box. Figure SI 2 shows a snapshot of the initial system configuration. The simulation box sizes and the number of atoms in a system are given in Table SI 1. The system was equilibrated for 10 ns to achieve the equilibrium box volume, with an average system pressure of 1 atm and a temperature of 300 K. A Nose-Hoover thermostat and barostat with relaxation times of 0.1 ps and 1.0 ps respectively were used to maintain the temperature and pressure in the system. The Nose-Hoover barostat was applied along the z-direction only because the size of the simulation box in the x and y-directions was fixed by the size of the graphene sheet. The data collection was performed during the final 3 ns of the simulation run (production run). Figure SI 3 shows the evolution of the simulation box size along the z-axis, Lz, and system density during the entire simulation run.
SI 4
Figure SI 3 Time dependence of the box size along the z-axis and system density in C6F6/ graphene simulation.
Figure SI 4 Density distribution along z-axis in C6F6/ C6H6 mixtures with different mole fraction of C6F6.
SI 5
Table SI 1 nCarbon
nMolecules
nAtoms
Graphene
Solvent
Total
Water 116.5 115.3 142.1 0.2
5376
67200
206976
Benzene 116.5 115.3 146.1 0.3
5376
12096
150528
Xylenes 116.5 115.3 200.2 0.3
5376
12096
223104
Heptane 116.5 115.3 170.4 0.3
5376
8736
206304
C6F6 116.5 115.3 173.4 0.2
5376
12096
150528
C6F6/C6H6 (3:1) 116.5 115.3
166.7 0.2
5376
12096
150528
C6F6/C6H6 (1:1) 116.5 115.3
158.7 0.3
5376
12096
150528
C6F6/C6H6 (1:3) 116.5 115.3 154.0 0.3
5376
12096
150528
Solvent
Lx(Å)
Ly(Å)
(Å)
To elucidate the affinity between different solvent mixtures and graphene, we have calculated the density distribution function of solvent molecules along the z-direction. Figure SI 4 shows the density distribution in the different molar mixtures of C6F6 and benzene. This density distribution was obtained by binning atom positions along the z-axis with a bin size of =0.1 Å. As in the case with other solvents (see Figure 1 in main text) the 1:1 mixture of C6F6/C6H6 has the highest degree of ordering among all C6F6/ C6H6 mixtures. This degree of ordering is due to the orientation of molecules parallel to the graphene sheet plane. This is clearly seen in Figure SI 5 showing the distribution of the orientational order parameter:
S ( z)
3 cos 2 ( z ) 1 2 SI 6
(SI 2)
where z is an angle between the z-axis and normal to the carbon ring of the aromatic solvents, and brackets correspond to the ensemble over all molecule orientations at distance z from the surface during the production run. The positive values of the order parameter S(z) correspond to the parallel alignment of molecules with the graphene sheet while the negative values of the order parameter correspond to perpendicular orientations of molecules with respect to the graphene sheet. The high degree of parallel alignment indicates columnar stacking of the C6F6 and benzene molecules.
Figure SI 5 Distribution of the orientational order parameter, S(z), along z-axis for different mole fractions of C6F6 in C6F6/ C6H6 solutions.
SI 7
Figure SI 6 Charge distribution, Q(z), and excess density distribution, (z)-b in 1:1 mixture of benzene and C6F6, where b is a bulk density.
It is interesting to point out that there is a correlation between density excess Figure in the adsorbed layer and local charge distribution. This is illustrated in Figure SI 6. The excess of negative charge correlates with an excess of the local density. The main contribution to this excess comes from fluorine atoms. Hydrogen atoms, with the smallest van der Waals radius, generate the positive charge excess observed at the graphene surface. The shoulders and double peaks seen in the charge distribution between peaks in the density distribution correspond to positively charged carbon atoms in C6F6 and positively charged hydrogen atoms in benzene.
SI 8
Figure SI 7 Snapshots of the first layer structure relative to the graphene sheet for water (left), heptane (center), and xylene (right).
Figure SI 7 shows snapshots of the structure of first adsorbed layer of molecules at the graphene surface. For these pictures we used only atoms which are located within 6 Å from the substrate. As one can see from these pictures, the number of molecules within these layers decreases with increasing molecular size. The highest number of molecules is observed for water. In addition, heptane and xylene molecules show weaker orientational alignment with respect to the graphene surface as compared with that seen in C6F6, C6F6/ C6H6, and C6H6 solvents as shown in Figure 2 in the main text. Quantitative information about the affinity of different solvents to graphene can be obtained from solvent surface excess: 2
Lz / 2
( z ) dz b
(SI.3)
d
where b is the bulk density and d=2.5 Å. The factor of two accounts for two graphene surfaces exposed to a solvent. The integration of equation SI.3 was performed numerically with the integration step =0.1 Å. Table SI 2 summarizes our results for solvent surface excess.
SI 9
Table SI 2 Solvent surface excess
Solvent
b
b g/cm
Water
1.04
0.67
0.65
Benzene
0.82
2.40
2.94
Xylenes
0.80
2.65
3.30
Heptane
0.64
2.43
3.78
C6F6
1.64
2.46
1.5
C6F6/C6H6 (1:1)
1.26
3.98
3.15
C6F6/C6H6 (3:1)
1.46
2.76
1.89
C6F6/C6H6 (1:3)
1.04
3.39
3.27
The largest excess mass per unit area is observed for the 1:1 mixture of C6F6/ C6H6. Note also that this value is larger than an average solvent surface excess of pure benzene and pure C6F6. This can be a macroscopic manifestation of the solvent mixture structuring at the graphene surface. The lowest surface excess was observed for water. The normalized quantity of the surface excess has the largest value for the 1:3 mixture of C6F6/ C6H6. This is a reflection of two tendencies: the growth of the density of the solvent mixture and the growth of the solvent surface excess. However, even for the normalized quantity, the value of the solvent surface excess for the 1:1 C6F6/ C6H6 mixture is larger than the average quantity for pure components.
SI 10
Figure SI 8 Generation eight (G8) coronene-like molecule C384H48. Carbon atoms are shown in black and hydrogen atoms are colored in green.
To show the superior graphene solubility of the 1:1 C6F6/ C6H6 mixture in comparison with pure solvents, we used the Weighted Histogram Analysis Method to calculate the potential of mean force between a graphene sheet and a graphene flake in different solvents.(9) The graphene flake was modeled by a G8 coronenelike molecule consisting of eight generations of carbon rings and terminated by hydrogen (see Figure SI 8). The partial charges of the coronene-like molecule were obtained with the Mulliken population analysis from ab initio calculations using G09 with 6-31G(d) basis set and B3LYP DFT method without geometry optimization. The C-C bond is set to 1.387 Å while the C-H bond is set to 1.087 Å. This is consistent with the bond lengths definition for aromatic carbons and hydrogens in GAFF.
SI 11
Table SI 3 – Systems used in PMF simulations nCarbon
nMolecules
nAtoms
System
Lx(Å)
Ly(Å)
Lz (Å)
Graphene
Solvent
Total
C6H6
58.3
57.7
75.7
1344
1512
19920
C6F6/C6H6 (1:1)
58.3
57.7
82.3
1344
1512
19920
C6F6
58.3
57.7
90.5
1344
1512
19920
In these simulations, the graphene sheet was placed parallel to the xy-plane at z=0 Å and spanned the xy-plane. The positions of the carbon atoms forming a graphene sheet were fixed during the entire simulation run. Initially a graphene flake was placed at z= 18.254 Å followed by a solvent being added to the simulation box. The initial size of the simulation box along the z direction was equal to Lz= 159.9 Å. The periodic boundary conditions were used in the x, y, and z directions. The number of atoms used in these simulations is given in Table SI 3. The system was equilibrated for 3 ns to achieve the equilibrium box volume, the average system pressure (1 atm) and the temperature (300 K). A Nose-Hoover thermostat and barostat with relaxation times 0.1 ps and 1.0 ps respectively were used to maintain temperature and pressure in the system. The Nose-Hoover barostat was applied along the z-direction only because the size of the simulation box in the x and ydirections was fixed by the size of the graphene sheet. This was followed by simulation runs where we calculated the potential of the Mean Force by using WHAM.(9) These simulations were performed at constant temperature and volume (system sizes are listed in Table SI 3). The constant
SI 12
temperature was maintained by coupling a system to the Nose-Hoover thermostat with relaxation time 0.1 ps. In these simulations the z-coordinate of the center of mass of the graphene flake with coordinates (xcm, ycm, zcm) was tethered to (0, 0, z*) by harmonic springs:
Uspring =
K1 2 K ( zcm - z *) + 2 ( xcm2 + ycm2 ) 2 2
(SI.4)
where the values of the spring constants are K1= 250 Kcal/mole/Å2 and K2 = 2000 Kcal/mole/Å2. During these simulation runs we varied the location of the tethering point z* with the increment z*=0.1 Å, until z*= 3.354 Å was reached. For each location of the tethered point, the system was equilibrated for 0.1 ns. The equilibration step was followed by the production run with a duration of 0.3 ns. During this step, we obtained the distribution of the center of mass locations of the graphene flake for the potential of mean force WHAM calculations. The simulation results are shown in Figure SI 9. The potential of the mean force calculations do not give the absolute value of the system Helmholtz free energy but provide information about the change in the system Helmholtz free energy as the graphene flake approaches the larger graphene sheet. The largest positive change in the Helmholtz free energy is observed for the 1:1 C6F6/C6H6 mixture. The positive change indicates that the mixture is a good solvent for graphene. In good solvents, a solvated state, with both sides of graphene covered with solvent, has a lower free energy than the aggregated (stacked) state, when the graphene sheets come into contact. Similar trend is observed for pure hexafluorobenzene. However, these values of F are smaller than in the mixed
SI 13
solvent, with the maximum values of F being achieved at a distance z about 5.2 Å. At such small separations, the graphene flake experiences large bending deformations, required for the expulsion of the last layer of the adsorbed molecules (see Figure SI 10). Note that in the 1:1 C6F6/C6H6 mixture, changes in F occur in a step-like fashion, reflecting the removal of well-organized layers of solvent molecules (see Figure SI 5). Steps in F are also seen for pure hexafluorobenzene, but they are less pronounced, and appear at smaller separations.
Figure SI 9 Dependence of the difference in the Helmholtz free energy F on the distance between graphene sheet and graphene flake. Letters correspond to snapshots of graphene sheet and flake configurations shown in Figure SI 10.
The opposite trend is observed for benzene. In this case, F first decreases with decreasing the distance between graphene sheet and flake. Thus, expelling benzene
SI 14
from the gallery between the sheet and flake is a thermodynamically favorable process. The value of F begins to increase at z~7 Å. This increase is due to the bending of the flake, which occurs prior the expulsion of the last layer of the benzene molecules located between the graphene sheet and graphene flake (see Figure SI10). Finally, at z
Institute of Materials Science, Polymer Program, University of Connecticut, Storrs, CT 06269, ‡Physics Department, University of Connecticut, Storrs, CT, 06269, ¶Chemistry Department, University of Connecticut, Storrs, CT 06259, ∫Department of Applied Science, The College of William and Mary, Williamsburg, VA, 23185
Simulation Details and Additional Simulation Results We have performed molecular dynamics simulations of the adsorption of water, benzene, hexafluorobenzene (C6F6), xylenes, and heptane on graphene sheets. The Generalized Amber Force Field (GAFF)(1) parameters were used for atomistic models of solvents and graphene while water was simulated by using a modified TIP3P model optimized for simulations with Ewald summation.(2) The partial charge distributions for the solvents, except water, were obtained by performing abinitio calculations using the Gaussian 09 (G09) simulation package(3) with 6-31G(d) basis set and B3LYP DFT method. The total potential energy of the system consisted of the bonded, bond angle, dihedral angle, improper angle and non-bonded interaction potentials. The interaction parameters for the van der Waals potential between heterogeneous atomic pairs were calculated as the geometric mean of the interaction parameters for each atom. The default AMBER force field weighing coefficients for pair-wise energy and force contributions were used to account for contributions from the van der Waals and electrostatic interactions. U TOTAL
K r r
2
r
BONDS
eq
K
2
eq
ANGLES
Aij Bij qi q j (SI.1) Vn 2 1 cos n K 12 6 eq Rij Rij DIHEDRALS 2 IMPROPER i j Rij
The simulation box was built using Chem3D,(4) G09, Antechamber(5) and AMBER2LAMMPS python script that is included with LAMMPS.(6) The G09 input file for the solvent molecule, e.g. benzene, was built in Chem3D, then G09 calculations
SI 1
with geometry optimization using AM1 semi-empirical method were performed. The Gaussian output from the calculation was used as an input for Antechamber to determine charges, atom type, bond type, angle and dihedral type assignments.
Figure SI 1 Partial charge distributions used in simulations: water (a), benzene (b), xylene (c), heptane (d) and C6F6 (e). All charges except for water were obtained by using Mulliken population analysis from ab initio calculations with 6-31G(d) basis set and B3LYP DFT method. The water charges were obtained from Price and Brooks.(2)
SI 2
Figure SI 2 Initial system configuration of equimolar mixture of benzene and C 6F6 molecules. The graphene sheet is located at z=0. Hydrogen atoms are shown as green beads, fluorine atoms as yellow beads, aromatic carbon atoms as grey beads, and carbon atoms in the graphene sheet are colored in black.
The AMBER topology file was created by using LEAP that was included in the Antechamber package. The AMBER topology file was converted into a LAMMPS data file using the python script AMBER2LAMMPS. Using output of the AMBER2LAMMPS script as a template, the solvent molecule was replicated and distributed in the simulation box using in-house code. The graphene sheet was modeled as a neutral xy-periodic “macromolecule” using the GAFF definition of the aromatic carbon for the van der Waals interaction parameters. The partial charges of the solvents, except water, were obtained with the Mulliken population analysis from ab initio calculations using G09 with 6-31G(d) basis set and B3LYP DFT method. The results are summarized in Figure SI 1.
SI 3
The NPT ensemble simulations were performed using GPU accelerated LAMMPS code.(7) The equations of motion were integrated by using the velocity Verlet algorithm with a time step of 1.0 fs. The system was periodic in the x, y and z directions. The standard PPPM(8) method with an accuracy of 1.0 × 10-5, and the near-field cutoff set to 10.0 Å, was used to account for contributions from the longrange electrostatic interactions. The graphene sheet was placed at z = 0 Å and spanned the xy-plane. The coordinates of the carbon atoms forming a graphene sheet were fixed during the entire simulation run. Solvent molecules were distributed over the volume of the simulations box. Figure SI 2 shows a snapshot of the initial system configuration. The simulation box sizes and the number of atoms in a system are given in Table SI 1. The system was equilibrated for 10 ns to achieve the equilibrium box volume, with an average system pressure of 1 atm and a temperature of 300 K. A Nose-Hoover thermostat and barostat with relaxation times of 0.1 ps and 1.0 ps respectively were used to maintain the temperature and pressure in the system. The Nose-Hoover barostat was applied along the z-direction only because the size of the simulation box in the x and y-directions was fixed by the size of the graphene sheet. The data collection was performed during the final 3 ns of the simulation run (production run). Figure SI 3 shows the evolution of the simulation box size along the z-axis, Lz, and system density during the entire simulation run.
SI 4
Figure SI 3 Time dependence of the box size along the z-axis and system density in C6F6/ graphene simulation.
Figure SI 4 Density distribution along z-axis in C6F6/ C6H6 mixtures with different mole fraction of C6F6.
SI 5
Table SI 1 nCarbon
nMolecules
nAtoms
Graphene
Solvent
Total
Water 116.5 115.3 142.1 0.2
5376
67200
206976
Benzene 116.5 115.3 146.1 0.3
5376
12096
150528
Xylenes 116.5 115.3 200.2 0.3
5376
12096
223104
Heptane 116.5 115.3 170.4 0.3
5376
8736
206304
C6F6 116.5 115.3 173.4 0.2
5376
12096
150528
C6F6/C6H6 (3:1) 116.5 115.3
166.7 0.2
5376
12096
150528
C6F6/C6H6 (1:1) 116.5 115.3
158.7 0.3
5376
12096
150528
C6F6/C6H6 (1:3) 116.5 115.3 154.0 0.3
5376
12096
150528
Solvent
Lx(Å)
Ly(Å)
(Å)
To elucidate the affinity between different solvent mixtures and graphene, we have calculated the density distribution function of solvent molecules along the z-direction. Figure SI 4 shows the density distribution in the different molar mixtures of C6F6 and benzene. This density distribution was obtained by binning atom positions along the z-axis with a bin size of =0.1 Å. As in the case with other solvents (see Figure 1 in main text) the 1:1 mixture of C6F6/C6H6 has the highest degree of ordering among all C6F6/ C6H6 mixtures. This degree of ordering is due to the orientation of molecules parallel to the graphene sheet plane. This is clearly seen in Figure SI 5 showing the distribution of the orientational order parameter:
S ( z)
3 cos 2 ( z ) 1 2 SI 6
(SI 2)
where z is an angle between the z-axis and normal to the carbon ring of the aromatic solvents, and brackets correspond to the ensemble over all molecule orientations at distance z from the surface during the production run. The positive values of the order parameter S(z) correspond to the parallel alignment of molecules with the graphene sheet while the negative values of the order parameter correspond to perpendicular orientations of molecules with respect to the graphene sheet. The high degree of parallel alignment indicates columnar stacking of the C6F6 and benzene molecules.
Figure SI 5 Distribution of the orientational order parameter, S(z), along z-axis for different mole fractions of C6F6 in C6F6/ C6H6 solutions.
SI 7
Figure SI 6 Charge distribution, Q(z), and excess density distribution, (z)-b in 1:1 mixture of benzene and C6F6, where b is a bulk density.
It is interesting to point out that there is a correlation between density excess Figure in the adsorbed layer and local charge distribution. This is illustrated in Figure SI 6. The excess of negative charge correlates with an excess of the local density. The main contribution to this excess comes from fluorine atoms. Hydrogen atoms, with the smallest van der Waals radius, generate the positive charge excess observed at the graphene surface. The shoulders and double peaks seen in the charge distribution between peaks in the density distribution correspond to positively charged carbon atoms in C6F6 and positively charged hydrogen atoms in benzene.
SI 8
Figure SI 7 Snapshots of the first layer structure relative to the graphene sheet for water (left), heptane (center), and xylene (right).
Figure SI 7 shows snapshots of the structure of first adsorbed layer of molecules at the graphene surface. For these pictures we used only atoms which are located within 6 Å from the substrate. As one can see from these pictures, the number of molecules within these layers decreases with increasing molecular size. The highest number of molecules is observed for water. In addition, heptane and xylene molecules show weaker orientational alignment with respect to the graphene surface as compared with that seen in C6F6, C6F6/ C6H6, and C6H6 solvents as shown in Figure 2 in the main text. Quantitative information about the affinity of different solvents to graphene can be obtained from solvent surface excess: 2
Lz / 2
( z ) dz b
(SI.3)
d
where b is the bulk density and d=2.5 Å. The factor of two accounts for two graphene surfaces exposed to a solvent. The integration of equation SI.3 was performed numerically with the integration step =0.1 Å. Table SI 2 summarizes our results for solvent surface excess.
SI 9
Table SI 2 Solvent surface excess
Solvent
b
b g/cm
Water
1.04
0.67
0.65
Benzene
0.82
2.40
2.94
Xylenes
0.80
2.65
3.30
Heptane
0.64
2.43
3.78
C6F6
1.64
2.46
1.5
C6F6/C6H6 (1:1)
1.26
3.98
3.15
C6F6/C6H6 (3:1)
1.46
2.76
1.89
C6F6/C6H6 (1:3)
1.04
3.39
3.27
The largest excess mass per unit area is observed for the 1:1 mixture of C6F6/ C6H6. Note also that this value is larger than an average solvent surface excess of pure benzene and pure C6F6. This can be a macroscopic manifestation of the solvent mixture structuring at the graphene surface. The lowest surface excess was observed for water. The normalized quantity of the surface excess has the largest value for the 1:3 mixture of C6F6/ C6H6. This is a reflection of two tendencies: the growth of the density of the solvent mixture and the growth of the solvent surface excess. However, even for the normalized quantity, the value of the solvent surface excess for the 1:1 C6F6/ C6H6 mixture is larger than the average quantity for pure components.
SI 10
Figure SI 8 Generation eight (G8) coronene-like molecule C384H48. Carbon atoms are shown in black and hydrogen atoms are colored in green.
To show the superior graphene solubility of the 1:1 C6F6/ C6H6 mixture in comparison with pure solvents, we used the Weighted Histogram Analysis Method to calculate the potential of mean force between a graphene sheet and a graphene flake in different solvents.(9) The graphene flake was modeled by a G8 coronenelike molecule consisting of eight generations of carbon rings and terminated by hydrogen (see Figure SI 8). The partial charges of the coronene-like molecule were obtained with the Mulliken population analysis from ab initio calculations using G09 with 6-31G(d) basis set and B3LYP DFT method without geometry optimization. The C-C bond is set to 1.387 Å while the C-H bond is set to 1.087 Å. This is consistent with the bond lengths definition for aromatic carbons and hydrogens in GAFF.
SI 11
Table SI 3 – Systems used in PMF simulations nCarbon
nMolecules
nAtoms
System
Lx(Å)
Ly(Å)
Lz (Å)
Graphene
Solvent
Total
C6H6
58.3
57.7
75.7
1344
1512
19920
C6F6/C6H6 (1:1)
58.3
57.7
82.3
1344
1512
19920
C6F6
58.3
57.7
90.5
1344
1512
19920
In these simulations, the graphene sheet was placed parallel to the xy-plane at z=0 Å and spanned the xy-plane. The positions of the carbon atoms forming a graphene sheet were fixed during the entire simulation run. Initially a graphene flake was placed at z= 18.254 Å followed by a solvent being added to the simulation box. The initial size of the simulation box along the z direction was equal to Lz= 159.9 Å. The periodic boundary conditions were used in the x, y, and z directions. The number of atoms used in these simulations is given in Table SI 3. The system was equilibrated for 3 ns to achieve the equilibrium box volume, the average system pressure (1 atm) and the temperature (300 K). A Nose-Hoover thermostat and barostat with relaxation times 0.1 ps and 1.0 ps respectively were used to maintain temperature and pressure in the system. The Nose-Hoover barostat was applied along the z-direction only because the size of the simulation box in the x and ydirections was fixed by the size of the graphene sheet. This was followed by simulation runs where we calculated the potential of the Mean Force by using WHAM.(9) These simulations were performed at constant temperature and volume (system sizes are listed in Table SI 3). The constant
SI 12
temperature was maintained by coupling a system to the Nose-Hoover thermostat with relaxation time 0.1 ps. In these simulations the z-coordinate of the center of mass of the graphene flake with coordinates (xcm, ycm, zcm) was tethered to (0, 0, z*) by harmonic springs:
Uspring =
K1 2 K ( zcm - z *) + 2 ( xcm2 + ycm2 ) 2 2
(SI.4)
where the values of the spring constants are K1= 250 Kcal/mole/Å2 and K2 = 2000 Kcal/mole/Å2. During these simulation runs we varied the location of the tethering point z* with the increment z*=0.1 Å, until z*= 3.354 Å was reached. For each location of the tethered point, the system was equilibrated for 0.1 ns. The equilibration step was followed by the production run with a duration of 0.3 ns. During this step, we obtained the distribution of the center of mass locations of the graphene flake for the potential of mean force WHAM calculations. The simulation results are shown in Figure SI 9. The potential of the mean force calculations do not give the absolute value of the system Helmholtz free energy but provide information about the change in the system Helmholtz free energy as the graphene flake approaches the larger graphene sheet. The largest positive change in the Helmholtz free energy is observed for the 1:1 C6F6/C6H6 mixture. The positive change indicates that the mixture is a good solvent for graphene. In good solvents, a solvated state, with both sides of graphene covered with solvent, has a lower free energy than the aggregated (stacked) state, when the graphene sheets come into contact. Similar trend is observed for pure hexafluorobenzene. However, these values of F are smaller than in the mixed
SI 13
solvent, with the maximum values of F being achieved at a distance z about 5.2 Å. At such small separations, the graphene flake experiences large bending deformations, required for the expulsion of the last layer of the adsorbed molecules (see Figure SI 10). Note that in the 1:1 C6F6/C6H6 mixture, changes in F occur in a step-like fashion, reflecting the removal of well-organized layers of solvent molecules (see Figure SI 5). Steps in F are also seen for pure hexafluorobenzene, but they are less pronounced, and appear at smaller separations.
Figure SI 9 Dependence of the difference in the Helmholtz free energy F on the distance between graphene sheet and graphene flake. Letters correspond to snapshots of graphene sheet and flake configurations shown in Figure SI 10.
The opposite trend is observed for benzene. In this case, F first decreases with decreasing the distance between graphene sheet and flake. Thus, expelling benzene
SI 14
from the gallery between the sheet and flake is a thermodynamically favorable process. The value of F begins to increase at z~7 Å. This increase is due to the bending of the flake, which occurs prior the expulsion of the last layer of the benzene molecules located between the graphene sheet and graphene flake (see Figure SI10). Finally, at z
