Investigation of graphene oxide nanosheets dispersion in water based on solubility parameters: a molecular dynamics simulation study

Abstract

In this work, the dispersion behavior of graphene oxide (GO) nanosheets in water was studied by a series of MD simulations for water, GO, and the water/GO mixed systems. The simulation results showed that the water/GO system has a well-ordered structure with strong H-bond interactions. The initial value (4.0 Å) of the interlayer distance between two GO sheets at the beginning of the simulation clearly increased (∼7.2 Å) in the presence of water as compared to its value after simulation in the absence of water (5.7 Å). The solubility parameter of GO was calculated as a function of both temperature and number of layers. The solubility parameter of GO at 300 K reached a plateau at 44.9 MPa^1/2 when the number of layers was five. Hence, water is a good solvent for dispersing GO since their solubility parameters are close to each other especially at temperatures close to the freezing point of water. The strong H-bond between water and oxygen-functional groups of GO makes the enthalpy favorable to form stable GO dispersions in water.

---

Introduction

Graphene oxide (GO), an oxygen-rich carbonaceous layered material, is an atomically thin sheet of graphite which is increasingly attracting chemists due to its properties. Most of the outstanding properties of GO arise from its hybrid electronic structure.^1,2

Many industrial and scientific applications benefit from the outstanding chemical, thermal, and electrical properties of GO, making its colloidal suspensions stable in various solvents, especially in water, very important.^1–5 A number of studies have reported the dispersion of GO in polar solvents.^3–5 The choice of an appropriate solvent for the dispersion of GO was based primarily on trial-and-error experiments.

Cohesive energy density is the net effect of all the interatomic/molecular interactions including van der Waals interactions, covalent bonds, ionic bonds, hydrogen bonds, and electrostatic interactions. An understanding of cohesive energy density is very important for predicting the stability of a colloid. The cohesive energy density of a material can be quantified in a number of ways. The most common approach is to use the solubility parameter, δ. The solubility parameter is one of the most useful concepts in the physical chemistry and thermodynamics of solutions and colloids. There is much interest in utilizing solubility parameters for rationally designing novel processes such as supercritical extraction, coatings, or new materials such as drugs and polymer alloys.⁶

In 1936, Hildebrand proposed a simple definition for a “solubility parameter” which provides a systematic description of the miscibility behavior of solvents.⁷ The solubility parameter is defined as the square root of the cohesive energy density (CED), i.e. the amount of energy needed to completely remove unit volume of molecules from their neighbors to infinite separation:⁸

where ΔE_v, ΔH_v, and V_m are internal energy change of vaporization, enthalpy of vaporization, and molar volume of the liquid at the temperature of vaporization. R and T are gas constant and absolute temperature, respectively.

The dispersion behavior of GO at different organic solvents were investigated by Paredes et al.³ Graphite oxide was dispersed in N,N-dimethyl formamide, tetrahydrofuran, and ethylene glycol. They mentioned that in all of these solvents, full exfoliation of the graphite oxide into individual and single-layer GO sheets was achieved by sonication. Konios et al.⁵ investigated the dispersion behavior of GO and reduced graphene oxide (rGO) in eighteen solvents by estimating the solubility parameter using UV-Vis spectroscopy. Hadadian et al.⁹ investigated the stability of GO–distilled water dispersions by a number of experiments. Although these studies provide useful preliminary information, they do not permit molecular-level design of new solvents that are capable of efficiently dispersing GO. Hence, developing a molecular-level understanding of GO–solvent interactions represents a very important step toward optimizing the design of stable GO dispersions. Molecular dynamics (MD) simulations provide a promising computational tool for elucidating the nature of GO–solvent interactions at the molecular level. Simulations can provide important insights into the structural, electronic, and chemical properties of GO.

The aim of the paper is theoretical study of GO dispersion in water using MD simulation, calculating GO and water's solubility parameters by changing the effective parameter “temperature”, in order to find the optimum conditions for having stable GO's suspension in water. To the best of our knowledge, up to now, no research has reported for calculating the solubility parameter of GO using MD simulation. Although there are several scattered reports on different aspects of MD simulation for GO's dispersion in water, none of these reports directly address the solubility parameter. For example, Bagri et al.¹⁰ used MD simulation to study the atomistic structure of progressively reduced graphene oxide. Medhekar et al.¹¹ elucidated the atomic level structure and mechanical properties of GO paper based on MD simulations using the reactive force field (ReaxFF). Shih et al.¹² carried out a series of comparative experimental and MD simulation studies to understand the fundamentals of the surface activity and colloidal stability of GO aqueous solutions at different pH values. Dimiev et al.¹³ proposed an unconventional view of GO chemistry and develop the corresponding “dynamic structural model”. Nicolaï et al.¹⁴ used the quantum mechanical calculations to develop a full set of force field parameters in order to perform MD simulations to understand and optimize the molecular storage properties inside graphene oxide frameworks. Zokaie and Foroutan¹⁵ compared the structural and dynamical properties of water confined between two GO sheets through MD simulations. They showed that the structure and dynamics of water near the GO surfaces changes under confinement conditions. However, none of the recent MD simulation studies have been calculated the solubility parameter of GO.

Computational methods

The simulations were carried out using LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) code provided by Sandia National Laboratories.¹⁶ In this work, the MD simulations of bulk water, GO, and mixed system containing GO in water (water/GO) were performed separately.

Water

The TIP4P potential model¹⁷ was used for water molecules. TIP4P is a rigid four site model which consists of three fixed point charges and one Lennard-Jones center. Vega et al.¹⁸ showed that TIP4P is a suitable model for deriving the equations of state and phase diagram of water. Zokai and Foroutan¹⁵ successfully applied TIP4P model to investigate the behavior of ordered water molecules confined between GO surfaces.

The simulations of water were run at different temperatures for 1000 water molecules with the time step of 0.1 fs. The periodic boundary conditions were applied for all three directions. All water simulations were initially run for 0.1 ns in NPT ensemble (Nose–Hoover thermostat and barostat) at 1 bar and desired temperature to reach the equilibrium and then switched to NVE ensemble for 0.1 ns. The NVE ensemble was selected for production and averaging step because the energy fluctuations are minimum in this ensemble.

Graphene oxide

The geometric structure of the single sheet of GO was extracted from the experimental data.¹⁹ According to the experiments, the structure of the single layer of GO was considered as 3 × 3 nm² containing 15 epoxy, 20 quinone, 3 carboxyl, and 140 hydroxyl groups. The oxygen functional groups were distributed on both sides of GO sheet so that the epoxy and hydroxyl groups on above and below graphene layer and the carboxylic and quinone groups at the edges and in the positions of the lattice defects. The single sheet of GO contains 412 carbon atoms and 181 oxygen atoms. Therefore, the ratio of C/O is about 2.3 which is in accord with the experimental data.¹⁹ Also, two structural defects in each GO sheet, one single and other di-vacancy were included. The initial interlayer distance of GO was chosen to be 4.0 Å for multilayer GO. Fig. 1 shows the scheme of monolayer GO containing structural defects.

Fig. 1 The schematic representation of the optimized single layer of GO. The positions of defects were shown with yellow circles. The gray, white, and red circles are carbon, hydrogen, and oxygen atoms, respectively.

The REBO potential²⁰ was used for carbon and hydrogen atoms. The bond stretching was represented by Morse potential (for C–O, O–H, and C=O bonds). The parameters for the Morse potential were taken from elsewhere.²¹ The harmonic cosine form for bending potential was used.²¹ The dihedrals were described by cosine form.²² LJ (12, 6) potential was used for describing the vdW interactions. The cutoff distance for non-bonded interactions was considered to be 15 Å. The LJ parameters for H, C, and O atoms are listed in Table 1.

Table 1 LJ parameters for C, H, and O atoms²³

Atom	σ (Å)	ε (kcal mol^−1)
C	3.80	0.080
H	2.60	0.008
O	3.60	0.150

---

The columbic interactions were also considered using Coulomb law with the cutoff radius of 15 Å. The long range columbic interactions beyond the Coulomb cutoff distance were also calculated using the particle–particle particle-mesh (PPPM) method.²⁴

To calculate the partial atomic charges, first the initial structure of GO was optimized in B3LYP/6-31G(d) level of theory. Then, the natural bond orbital (NBO) was employed to calculate the partial atomic charges using density functional theory at the same level of theory. The mean partial atomic charges for oxygen and hydrogen atoms of hydroxyl groups were calculated as −0.46 and 0.27, respectively. The mean partial atomic change for oxygen atoms of epoxy groups was calculated as −0.33. These values were 0.35 and −0.30 for carbon and oxygen atoms of quinone groups, respectively. Also, the mean atomic charges calculated for carbon and oxygen atoms due to carbonyl group and oxygen and hydrogen atoms of hydroxyl group of carboxyl groups were 0.61 and −0.40, −0.43 and 0.39, respectively. The partial atomic charges calculated using B3LYP/6-31G(d) level of theory are a preliminary and raw guess for atomic charges. Hence, we used the charge equilibration method,²⁵ as implemented in LAMMPS, to minimize the electrostatic energy of the system by adjusting the partial charge on individual atoms based on interactions with their neighbors. All the quantum mechanics calculations were performed by Gaussian 09 software program package.²⁶ After that, the simulations of GO were carried out in NPT ensemble with the total simulation time of 40 ns (time step of 1 fs). Three MD simulations for water, GO, and mixed water/GO systems were performed.

Results and discussions

Water

In order to evaluate and understand the effect of the presence of GO nanosheets on the surrounding water molecules in the mixed water/GO system, we need to compare some structural and dynamical properties of pure water with those of water in the mixed water/GO system. These properties are radial distribution functions (RDFs), g(r), and mean square displacements (MSDs). The radial distribution function, g(r), (RDF) is an example of a pair correlation function which describes how on average the particles correlations in a substance decay with increasing separation.²⁷ The three site–site radial distribution functions, g_OO(r), g_OH(r), and g_HH(r) are commonly used when the structure of the liquid water is studied.

The oxygen–oxygen, g_OO(r), oxygen–hydrogen, g_OH(r), and hydrogen–hydrogen, g_HH(r), RDFs were derived using MD simulation with considering TIP4P potential model at 300 K and 1 bar. g_OO(r), g_OH(r), and g_HH(r) are given in Fig. 2a–c, respectively. As Fig. 2 shows, there is a good agreement between our work and previous experimental and theoretical studies.^28–31

Fig. 2 (a) Oxygen–oxygen, (b) oxygen–hydrogen, (c) hydrogen–hydrogen radial distribution functions calculated for the TIP4P force field of liquid water at 300 K and 1 bar and compared with experimental and theoretical values.

The mean square displacements (MSDs) of liquid water at different temperatures were computed using a series of MD simulations. The curve of MSD versus time calculated from averaging over for 1000 water molecules and the 0.1 ns trajectories in the temperature range from 280 to 360 K (Fig. 3). Fig. 3 shows that as time elapses, the MSD grows linearly. Also, the MSD increases with raising temperature since the mobility of molecules increases.

Fig. 3 The MSD plot for liquid water molecules at different temperatures and 1 bar. The curves were calculated using the MD simulation with considering TIP4P potential model for liquid water.

Table 2 gives the calculated internal energy change of vaporization (ΔE_vap), molar volume (V_m), and solubility parameter of water (δ_water). The solubility parameters were derived using eqn (1). As Table 2 shows there is a good agreement between the calculated and experimental values of solubility parameter for liquid water at different temperatures. Table 2 also indicates that as temperature rises, the solubility parameter decreases. With increasing the temperature, the interactions between the molecules in the system weaken as a result of an increase in kinetic energies of the molecules. Therefore, the amount of energy required to overcome the interactions and vaporize the molecules decreases, and the solubility parameter decreases.

Table 2 Calculated molar volumes, internal energy change of vaporization, and solubility parameters of water at different temperatures and 1 bar

T (K)	Vm^a (cm^3 mol^−1)	ΔEvap^a (kJ mol^−1)	δwater (MPa^1/2)		% Error (δwater)
			Calc.	Exp.^b
a Calculated using MD simulation.b Experimental data was calculated based on enthalpies of vaporization extracted from ref. 32.c Ref. 6.
280	17.908	40.368	47.480	48.591	2.3
300	18.124	39.106	46.454	47.994	3.2
				47.820^c	2.8
320	18.160	37.872	45.663	47.391	3.6
340	18.630	36.564	44.304	46.779	5.3
360	18.963	35.271	43.128	46.159	6.5

---

Graphene oxide

The structural properties including RDFs and potential of mean force and thermodynamic properties such as solubility parameter of GO were calculated using MD simulations.

RDFs

Fig. 4 shows the RDF plots for oxygen of the first layer and hydrogen of the first (O₁–H₁), second (O₁–H₂), third (O₁–H₃), fourth (O₁–H₄), and fifth (O₁–H₅) layer of GO at 280 K and 1 bar. The inset plots are the corresponding potential of mean force, W(r), (PMF) per unit area. The PMF is the work needed to bring the two particles from infinite separation to a distance r in a dense system. It is related to the radial distribution function of the system by the following equation:³³

W(r) = −kT ln g(r) (2)

where k is Boltzmann constant.

Fig. 4 The RDF plots for the oxygen of the first layer and hydrogen of the (a) first (b) second (c) third (d) fourth and (e) fifth layers of GO at 280 K and 1 bar. The insets are the corresponding PMF curves per unit area.

The first peak of RDF and PMF for the O₁–H₁ appears at 1.98 Å (Fig. 4a). The negative value of PMF at this distance indicates there is a strong attraction between O₁ and H₁. The ratio of height of second peak is much less than the first one. It means there is a weak correlation between O₁–H₁ at this distance (second peak). In other words at distances greater than 2.2 Å, there is no correlation between O₁ and H₁.

Fig. 4b shows that the height of the prominent peak, first peak, is less than that of Fig. 4a. This means there is a weaker interaction between oxygen of the first layer and hydrogen of the second layer of GO compared to that of the first layer. Also, the value of PMF of O₁–H₂ (−0.365 kJ mol⁻¹ nm⁻²), is more positive than that of O₁–H₁ (−0.428 kJ mol⁻¹ nm⁻²) at the position of the first peak. As Fig. 4c to e show there is no interaction between O₁–H₃, O₁–H₄, and O₁–H₅.

Fig. 5a and b show the RDF and PMF plots for O₁–H₂ at different temperatures, respectively. Because of thermal motions, the intensity of the first peak of the RDFs decreases as temperature rises.

Fig. 5 (a) The RDF and (b) the PMF plots for O₁–H₂ of GO for different temperatures at 1 bar. The insets are a close view of (a) the first peak of RDFs (b) the well depth of the PMF curves.

Solubility parameter

The internal energy change of vaporization, ΔE_vap, of GO as functions of both temperature and number of layers was computed using MD simulations. The ΔE_vap of a substance is defined as the increase in internal energy of substance upon removal of all intermolecular interactions. The internal energy of vaporization of GO is defined as:³⁴

ΔE_vap = E_NL − NE_SL (3)

where E_NL and E_SL are the total energy of N-layer and single layer GO. The total energy of GO for different number of layers at various temperatures were computed using MD simulations (Table 3). The energy of two and three layers GO were calculated just at 300 K and are equal to −6668 and −10 022 eV, respectively.

Table 3 The energy of multilayer GO

T (K)	EGO (eV)
	N = 1	N = 4	N = 5 (bulk)
280	−3278	−13 401	−16 751
300	−3289	−13 401	−16 751
320	−3304	−13 428	−16 797
340	−3310	−13 445	−16 819
360	−3312	−13 459	−16 822

---

The volumes and internal energy of vaporizations of GO for different number of layers at 300 K was computed using MD simulations and the results are shown in Table 4. As this table shows the internal energy of vaporization rises with the increase of the number of layers because the number of interactions increases.

Table 4 The volumes and internal energy change of vaporizations of GO for different number of layers at 300 K

N	V (Å^3) × 10^3	ΔEvap (eV)
2	10.800	90
3	14.852	155
4	19.531	245
5	24.281	306

---

The average number of hydrogen bonds (H-bond) per single layer of GO can be computed by integrating g_O1–H2(r) plot up to the first minimum. Table 5 gives the average number of H-bond per single layer of GO and the number of H-bond per unit area (D_HB). As temperature increases both number of H-bond and DHB decreases because of increasing of thermal motions.

Table 5 The average number of H-bonds per a single layer of GO and the number of H-bond per unit area

T (K)	Average no. of H-bonds per single layer of GO	DHB (no. H-bond per nm^2)
280	71.6	7.99
300	71.3	7.92
320	69.5	7.72
340	68.0	7.56
360	63.0	7.00

---

Although a large number of experimental work on solubility parameter for liquids can be found in literature, for solids only a few publications reported the solubility parameter. This is probably due to the fact that experimental measurements of solubility parameters for liquids are easier to obtain than for solids. In this work, the solubility parameter of GO was calculated as a function of both temperature and number of layers using volume and ΔE_vap reported in Table 4. Fig. 6 shows the solubility parameter of GO versus the number of layers at 300 K and 1 bar. The solubility parameter of GO at 300 K reaches a plateau at 44.9 MPa^1/2 when the number of layers is five. Therefore, the five layer GO was considered as bulk state.

Fig. 6 Solubility parameter of GO as a function the number of layers at 300 K and 1 bar.

Table 6 shows the computed values of volume and solubility parameter of the bulk GO (N = 5) at different temperatures. The computed values of solubility parameters were compared to the corresponding computed data for liquid water in Table 6. This table shows both values are very close to each other especially at lower temperatures. Hence, water is a good solvent for dispersing GO especially at temperatures close to freezing point of water. This observation may be associated with the fact that hydrogen bonds between water molecules and oxygen groups of GO are stronger at lower temperatures. As mentioned before, the solubility parameter is the net effect of all interatomic/molecular interactions that here are the electrostatic (H-bond) and van der Waals. Table 6 reveals that solubility parameter of the both GO and water increases as temperature falls from 360 to 280 K but the slope of this enhancement is steeper for GO. Finally, the minimum difference between these two solubility parameters occur at temperatures close to the freezing point of water. In the other words, the interactions between water molecules and between GO layers become more similar at lower temperatures. One of the basic concepts in solution chemistry is that materials with the same interactions mix more easily with each other.

Table 6 The volumes and solubility parameters of five layers GO (bulk GO) and solubility parameters of water at different temperatures

T (K)	V (Å^3) × 10^3	δGO (MPa^1/2)	δwater (MPa^1/2)	% Error
280	24.246	48.8	47.5	2.7
300	24.281	44.9	46.4	−3.3
320	24.292	42.5	45.7	−7.5
340	24.317	41.6	44.3	−6.5
360	24.322	41.2	43.1	−4.6

---

It is worthy to mention that along with simulation studies, we investigated the dispersion of GO in water experimentally.⁹ Our experimental data are in acceptable agreement with the simulation results, i.e., the aqueous suspensions of the prepared GO were stable for more than 5 months without any sedimentation.

Water/GO mixed system

To gain insight into atomic-level structural change of GO and water in the presence of each other (water/GO system), two 3 × 3 nm² parallel GO sheets (in XZ plane) containing 30 epoxy, 40 quinone, 6 carboxyl, and 280 hydroxyl groups were constructed for simulations. Total number of carbon and oxygen atoms in two GO sheets were 824 and 362, respectively; therefore, the ratio of C/O is about 2.3. We also included two structural defects in each GO sheet, one single and other di-vacancy. The initial interlayer distance was chosen to be 4.0 Å. 1500 water molecules were randomly dispersed around GO with initial solute/solvent distance of 2.0 Å. We performed extensive MD simulations using force fields described before in Sections 2.1 and 2.2 for water and GO, respectively. The long range coulomb interactions beyond a cut of distance of 15 Å are calculated using PPPM method. The simulations were performed in NPT ensemble at 300 K and 1 bar with time step of 0.1 fs. The total simulation time was 20 ns. The snapshot of the water/GO system at 300 K and 1 bar is shown in Fig. 7. As this figure shows the interlayer distance between two GO sheets clearly increased (∼7.2 Å) in the presence of water compared to its initial value (4.0 Å) at the start of simulation and its value after simulation in the absence of water (5.7 Å). In the other words, water exfoliates the GO sheets in the water/GO system. The increase in interlayer distance between two GO sheets in the absence of water is due to the presence of oxygen-containing functional groups because they produce electrostatically charged interlamellar layer and weaken the van der Waals interactions. By introducing of water molecules, further increase in the distance between sheets of GO occurs through the formation of H-bond between the oxygen-containing functional groups and water molecules and also between water molecules trapped inside the interlamellar layer. This weakens the van der Walls interactions.³⁵

Fig. 7 Snapshot of simulation box containing bilayer GO and water molecules at 300 K and 1 bar.

To study the influence of GO in the structure of water, we first computed the oxygen–oxygen RDF(g_OO(r)) and oxygen–hydrogen RDF (g_OH(r)) of water in the presence of double layer GO and compared to that of pure water (Fig. 2a and b). The main features of Fig. 8 are:

Fig. 8 (a) g_OO(r) and (b) g_OH(r) of water in the presence and the absence of GO.

(i) The number of peaks of g_OO(r) in the water/GO system with respect to liquid water increases, meaning that the structure of water in the presence of GO is more ordered than that of pure water. The g_OO(r) plot of water in the water/GO system approaches to unity at longer distances (14 Å) compared with that of pure water (8 Å). This is due to the strong H-bond between oxygen-containing functional groups of GO and water molecules. Hence, the H-bond between GO and water has a profound effect on the structure of water.

(ii) The height of the first peak of g_OO(r) for water in the water/GO system is greater than that of pure water. It means that the number of water molecules at the nearest neighboring shell increases because of more ordered structure of water at the presence of GO.

(iii) The position of the first peak of g_OO(r) remains almost unchanged whereas the second peak shifts to longer distance in the water/GO system. The same results were obtained for water molecules confined between GO sheets.¹⁵ The shift observed for the second peak can be assigned to the disruption of hydrogen bond between water molecules in the presence of GO.¹⁵

(iv) Fig. 8b reveals an interesting fact. As this figure shows, in the bulk water, the first peak of g_OH(r) located at 1.84 Å (assigned as mean of H-bond length) was disappeared in the presence of GO sheets. Instead, some overlapped broad peaks appear in the region 2.2–3.0 Å. This finding approves that H-bond interactions between water molecules in the water/GO system were strongly disrupted and a longer (so weaker), and non-uniform H-bond forms. Medhekar et al.¹¹ shows that optimum H-bond length between water molecules in the similar water/GO system is 2.55 Å.

To get further knowledge on the structure of water/GO system, the most interesting RDFs, g_OH(r) (oxygen atom of GO and hydrogen atom of water) was computed (Fig. 9). The first peak, which indicates the hydrogen bond interactions, located at about 1.20 Å. It was shifted to shorter distances compared to H-bond distance in water (1.84 Å) (Fig. 2b) which clearly shows a strong H-bond between GO oxygen functional groups and water molecules was formed.

Fig. 9 RDF plot for oxygen atom of GO and hydrogen atom of water for water/GO system.

The MSD plots in the X, Y, and Z-directions as well as total MSD of water for water/GO system (T = 300 K and 1 bar) were also computed using MD simulations (Fig. 10). A comparison between the MSD plots of bulk water (Fig. 3) and water in water/GO system (Fig. 10) indicates that strong H-bond network between GO and water molecules makes the water to behave as a pseudo-solid material with a saturated MSD plot. In fact, the mobility of water molecules diminishes because of strong H-bonds between GO and water. The MSD values in X and Z directions are greater than that of Y direction because Y axes is perpendicular direction to the GO surface (XZ plane). The similar results were observed in ref. 17.

Fig. 10 MSD plots of water in water/GO system in X, Y, and Z directions as well as total MSD.

Conclusions

In this work, the MD simulation as a powerful tool was used to study the GO nanosheets dispersion in water. For this purpose, three different sets of MD simulations were performed, namely for water, GO, and GO in the presence of water.

The main conclusions of this work are:

(i) The internal energy change of vaporization of GO increases with the increase of layers.

(ii) Water/GO system has a well-ordered structure with strong H-bond interactions. Water molecules penetrate into the GO interlayer space and exfoliate it easily (thanks to H-bond interactions).

(iii) MD simulation results showed that water is a good solvent for dispersing GO since the values of solubility parameter of GO is close to those of water especially at temperatures close to freezing point of water. The strong H-bond between water and oxygen-functional groups of GO makes the enthalpy favorable to form stable GO dispersions in the water.

(iv) The strong hydrogen bonding between GO and water molecules is responsible to form stable colloidal suspension of GO in water. It is an enthalpy driven process. The value of solubility parameter for GO supports this tendency since it was calculated based on enthalpy of vaporizations. Hence, some sort of enthalpy driven feature of GO dispersion in water can be concluded.

Acknowledgements

Financial support from the Ferdowsi University of Mashhad (grant no. 3/31131) is gratefully acknowledged.

References

1. S. Stankovich, R. D. Piner, S. T. Nguyen and R. S. Ruoff, Carbon, 2006, 44, 3342–3347.
2. O. C. Compton and S. T. Nguyen, Small, 2010, 6, 711–723.
3. J. Paredes, S. Villar-Rodil, A. Martinez-Alonso and J. Tascon, Langmuir, 2008, 24, 10560–10564.
4. D. H. Kim, Y. S. Yun and H.-J. Jin, Curr. Appl. Phys., 2012, 12, 637–642.
5. D. Konios, M. M. Stylianakis, E. Stratakis and E. Kymakis, J. Colloid Interface Sci., 2014, 430, 108–112.
6. C. M. Hansen, Hansen Solubility Parameters: a User's Handbook, CRC press, 2012.
7. J. H. Hildebrand, Solubility of Non-electrolytes, Reinhold Publishing Corporation, 1936.
8. E. K. Goharshadi and M. Hesabi, J. Mol. Liq., 2004, 113, 125–132.
9. M. Hadadian, E. K. Goharshadi and A. Youssefi, J. Nanopart. Res., 2014, 16, 1–17.
10. A. Bagri, C. Mattevi, M. Acik, Y. J. Chabal, M. Chhowalla and V. B. Shenoy, Nat. Chem., 2010, 2, 581–587.
11. N. V. Medhekar, A. Ramasubramaniam, R. S. Ruoff and V. B. Shenoy, ACS Nano, 2010, 4, 2300–2306.
12. C.-J. Shih, S. Lin, R. Sharma, M. S. Strano and D. Blankschtein, Langmuir, 2011, 28, 235–241.
13. A. M. Dimiev, L. B. Alemany and J. M. Tour, ACS Nano, 2013, 7, 576–588.
14. A. Nicolaï, P. Zhu, B. G. Sumpter and V. Meunier, J. Chem. Theory Comput., 2013, 9, 4890–4900.
15. M. Zokaie and M. Foroutan, RSC Adv., 2015, 5, 39330–39341.
16. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
17. J. L. Abascal and C. Vega, J. Chem. Phys., 2005, 123, 234505.
18. C. Vega, C. McBride, E. Sanz and J. L. Abascal, Phys. Chem. Chem. Phys., 2005, 7, 1450–1456.
19. G. Shao, Y. Lu, F. Wu, C. Yang, F. Zeng and Q. Wu, J. Mater. Sci., 2012, 47, 4400–4409.
20. D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni and S. B. Sinnott, J. Phys.: Condens. Matter, 2002, 14, 783.
21. S. L. Mayo, B. D. Olafson and W. A. Goddard, J. Phys. Chem., 1990, 94, 8897–8909.
22. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
23. M. Freindorf and J. Gao, J. Comput. Chem., 1996, 17, 386–395.
24. B. A. Luty and W. F. van Gunsteren, J. Phys. Chem., 1996, 100, 2581–2587.
25. A. K. Rappe and W. A. Goddard III, J. Phys. Chem., 1991, 95, 3358–3363.
26. M. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci and G. A. Petersson, Gaussian 09, Wallingford, CT, 2009.
27. A. Morsali, E. K. Goharshadi, G. A. Mansoori and M. Abbaspour, Chem. Phys., 2005, 310, 11–15.
28. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
29. J. M. Sorenson, G. Hura, R. M. Glaeser and T. Head-Gordon, J. Chem. Phys., 2000, 113, 9149–9161.
30. A. Soper, J. Phys.: Condens. Matter, 2007, 19, 335206G.
31. G. Hura, D. Russo, R. M. Glaeser, T. Head-Gordon, M. Krack and M. Parrinello, Phys. Chem. Chem. Phys., 2003, 5, 1981–1991.
32. H. W. Horn, W. C. Swope, J. W. Pitera, J. D. Madura, T. J. Dick, G. L. Hura and T. Head-Gordon, J. Chem. Phys., 2004, 120, 9665–9678.
33. D. A. McQuarrie, Statistical Thermodynamics, University Science Books, California, 1973.
34. J. Gupta, C. Nunes, S. Vyas and S. Jonnalagadda, J. Phys. Chem. B, 2011, 115, 2014–2023.
35. S. Abdul Rashid, S. A. Mohd Zobir, S. Krishnan, M. Hassan and H. Lim, J. Nanopart. Res., 2015, 17, 1–11.

---

Footnote

† Electronic supplementary information (ESI) available: Table S1 (justification of validity of potential used in this work). See DOI: 10.1039/c5ra19932h

---