Graphene exfoliation in ionic liquids: unified methodology

Abstract

Exfoliation constitutes a promising and straightforward technique to obtain a high-quality product (graphene, GRA) from an affordable source (graphite). Success of the exfoliation is blocked by inefficient solvents, which must prevent backward aggregation of GRA. The search for a solvent with optimal physicochemical properties is currently underway. This work applies a systematic methodology to rate solvents with respect to their performance in GRA exfoliation. The established methodology was applied to 1-ethyl-3-methylimidazolium tetrafluoroborate, [EMIM][BF₄], an ionic liquid (IL). Having derived thermodynamic descriptors during a real-time exfoliation, we conclude that the performance of [BMIM][BF₄] is mediocre. Nevertheless, [EMIM][BF₄] is a suitable starting point to engineer more advanced solvents. According to binding energy decomposition, the imidazole ring drives the GRA–[EMIM][BF₄] interactions, while the impact of [BF₄] is significantly smaller. We discuss possible pathways to improve the performance of ILs for GRA exfoliation.

---

Introduction

Graphene (GRA) possesses a marvelous set of mechanical, electrical, and thermal properties. These properties allow us to devise a wide variety of potential applications, including electronic and optoelectronic devices, chemical sensors, nanocomposites, and energy storage.^1–5 It is important to develop an efficient and affordable approach to produce GRA in industrial quantities to foster applications.⁶

Lots of different techniques to produce graphene sheets were introduced during the last decade.^6–8 Direct exfoliation of GRA from naturally occurring graphite flakes in the liquid phase showed its superiority for both scaling-up and easier covalent functionalization of GRA.^9–11 Functionalization of GRA is of drastic important for many applications, since it modulates electronic properties.^12,13 Exfoliation may potentially produce GRA, which is free of defects or oxygen containing groups. Oxidation of graphite to foster exfoliation is a large parallel research direction. Its major drawback is inability to remove all oxygen from the GRA surface after an exfoliation took place. In fact, understanding exfoliation is important for technological applications, such as top-down approaches to electronics. Besides, it can also be used as a higher-quality alternative to the chemical vapor deposition (CVD).^9–11 CVD is currently a popular method that produces relatively high-quality GRA with a potential for a large-scale implementation.

Despite offering certain advantages as discussed above, direct liquid-phase exfoliation, still results in a low concentrated colloidal suspensions of GRA.^14–16 Therefore, alternative liquid-phase techniques are desirable. Such techniques are expected to produce concentrated and kinetically stable suspensions of the GRA sheets in some solvent. Furthermore, extraction of the exfoliated GRA must be doable. Polar solvents are unable to deal with this task due to hydrophobicity of pristine GRA. Popular organic solvents – such as toluene, benzene, cyclohexane, octane – appear also inefficient because of a size difference between the solute and the solvent molecules (unfavorable entropic contribution to the free energy of solvation).

Ionic liquids (ILs) differ from conventional organic solvents in a few aspects.^17–21 ILs were successfully used for covalent and non-covalent functionalizing carbon nanotubes (CNTs) and GRA.^22–24 Many ILs exhibit similar surface tension to the surface energy of graphite. This is believed to be an essential feature for a direct solvent-mediated exfoliation of graphite.^21,25 In addition, ILs stabilize exfoliated GRA.^6,17,26,27 Such advantages over most other solvents make ILs interesting for GRA synthesis.

Recently, several studies employing ILs as assisting solvent media in GRA exfoliation have been reported.^6,28–34 Direct exfoliation of graphite flakes assisted by 1-butyl-3-methylimidazolium bis(trifluoromethanesulfonyl)amide, was conducted by Wang and coworkers.³¹ The applied technique produced GRA suspension with a concentration of 0.95 g L⁻¹. Mariani and coworkers³² used the 1-hexyl-3-methylimidazolium hexafluorophosphate IL as a solvent medium to exfoliate graphite at mild conditions. Concentrations of GRA up to 5.33 g L⁻¹ were recorded. Note that both successful research dealt with the imidazolium family of room-temperature ILs.

IL-rich aqueous mixtures – x(H₂O) < 10% – were demonstrated to functionalize carbonaceous nanomaterials by IL moieties.³⁵ According to the same study, high-quality carbon nanoribbons can be produced from graphite. Najafabadi and Gyenge³⁶ combined efficient cationic and anionic intercalating agents to increase GRA production. This advance was achieved using a specific technique combining cathodic and anodic exfoliation.

Atomistic-level description of structural and thermodynamic changes in the IL-containing systems during exfoliation is scarce.^27,37 It is important to understand which part of IL (cation, anion, aromatic ring, hydrocarbon chains) is primarily responsible for exfoliation and how the corresponding interaction can be enhanced. The van de Waals and π-stacking interactions (GRA–EMIM) are expected for the exfoliation. These interactions must be maximized between the successful solvent and GRA. This paper employs atomistic-precision molecular dynamics (MD) simulations based on the pairwise interaction potentials to characterize interaction between GRA and the chosen IL, 1-ethyl-3-methylimidazolium tetrafluoroborate [EMIM][BF₄]. We compute the potential of mean force (free energy vs. distance) for the exfoliation reaction and discuss this result in terms of structure and thermodynamics. The established methodology can be used to rate other ILs in view of GRA exfoliation.

Computational details

We immersed the graphene sheet (3.5 × 3.5 nm²) into 1-ethyl-3-methylimidazolium tetrafluoroborate IL, [EMIM][BF₄], 500 ion pairs. Fig. 1 depicts a molecular configuration corresponding to free energy minimum at 400 K, after 20 × 10⁶ time-steps of MD simulation.

Fig. 1 GR immersed in the [EMIM][BF₄] IL. The anions are blue, the cations are yellow, GR is light blue. GR edges were passivated by hydrogen atoms. The depicted molecular configuration corresponds to an equilibrium part of MD trajectory, after 20 × 10⁶ time-steps. Only a half of liquid volume is shown for clarity.

The MD simulations were performed in the isothermal–isobaric ensemble (N, P, T). We simulated the system at a higher temperature, 400 K, as compared to what we usually do. This was done to accelerate dynamics and avoid multiple metastable states, which are abundant in many ionic liquids at room temperature. Furthermore, GR increases shear viscosity of IL. Simulation at room conditions would require immense sampling times and, therefore, excessive computational resources. An interested reader can assess an effect of the elevated temperature using an entropic free energy component for the exfoliation process discussed below. All simulations were fulfilled at 1 bar controlled by the Parrinello–Rahman barostat³⁸ with a relaxation time of 4.0 ps and compressibility constant of 4.5 × 10⁻⁵ bar⁻¹. Constant temperature was maintained by the velocity rescaling thermostat with a relaxation time of 500 fs.³⁹

GR was simulated as a non-polarizable assembly of carbon atoms linked by harmonic potentials with the following bonded parameters: distance (C–C) = 1.41 nm, angle (C–C–C) = 120°, dihedral (C–C–C–C) = 0°. The Lennard-Jones (12, 6) parameters for the sp² carbon were taken as follows: σ = 0.34 nm, ε = 0.36 kJ mol⁻¹, which is in concordance with all popular general-purpose force fields. [EMIM][BF₄] was simulated using the force field introduced by one of us.⁴⁰ This model provides correct thermodynamic and transport properties in a wide temperature range. Since carbon atoms of GRA contain zero point charges, IL–GRA binding is purely of van der Waals nature. The systems were pre-equilibrated during 1.0 ns. Equilibrium properties were determined from the 10 ns long production runs. The integration time-step of 2.0 fs was used. Coordinates and intermediate thermodynamic quantities were saved every 200 fs.

To calculate the solvation free energy, the Bennett Acceptance Ratio (BAR) method was employed.⁴¹ A gradual decoupling of the solute molecule from its equilibrium environment (solvation shell) was performed in 21 steps independent MD simulations, Δλ = 0.05, 0 < λ < 1. The BAR method is used to enhance phase space sampling by the simulated MD system. Scaling solvent–solute interactions is a useful strategy to observe rare events and to account for their contribution into solvation free energy. With every λ, the system was equilibrated using stochastic MD during 1.0 ns. The averaged was derived from the 10 ns long production trajectory part. is in a direct relation to solvation free energy. The soft-core repulsion term for the LJ interactions was applied to avoid fusion of atomic nuclei:^42–44

V_SC(r) = (1 − λ)V([ασ⁶λ^p + r⁶]^1/6)

here, V_SC(r) is the normal hard-core pair potential, α is the LJ size parameter of the atom pair, whereas λ = 0 and λ = 1 correspond to the fully coupled and uncoupled states, respectively. The parameters for the soft-core potential were α = 0.5, p = 1.0, and σ = 0.3. Friction coefficient of 1 ps⁻¹ was used in the Langevin thermostat. Decoupling of the Coulomb and van der Waals terms upon thermodynamic integration can lead the system to instability when both above terms are decoupled together.⁴⁵ To avoid this undesirable behavior, the electrostatic interactions were turned on separately after the Lennard-Jones interactions.

Potential of mean force (PMF) is a one-dimensional function. PMF conveniently describe exfoliation of the two GRA sheets in [EMIM][BF₄]. We employed the umbrella sampling technique⁴⁶ with the weighted histogram analysis method (WHAM)⁴⁷ to derive PMF in the considered MD system. This technique requires restraints to be applied at specific positions along the reaction coordinate, so that all regions are appropriately sampled. An additional MD box was constructed and equilibrated with two parallel GRA sheets. Consequently, steered MD in the z-axis was used to generate 41 multiple configurations (Fig. 2), with a step of 0.1 nm. The distances were measured between centers-of-mass of the two GRA sheets. These distances were restrained using a one-center harmonic potential with a penalty force constant of 3000 kJ mol⁻¹ nm⁻². The reaction coordinate coincided with the normal vector towards the GRA plane. Each of the 41 positions of the reaction coordinate was sampled for 20 ns, whereas the first five nanoseconds of each such simulation were regarded as an equilibration.

Fig. 2 Three representative configurations visualizing steered MD simulations. The obtained molecular configurations were subsequently used as starting points for the umbrella sampling simulations. Anions are blue, cations are yellow. Only a half of liquid volume is shown for clarity. Note the flexibility of GRA. Simulation of rigid GRA is expected to provide totally different results in this context. The force favoring exfoliation acts along the normal to the GRA plane, since the vdW attraction between layers of graphite has to be overcome.

All reported MD simulations were performed using the GROMACS 5 molecular simulation engine.^43,48

Results and discussion

The imidazolium ring is expected orient flat near the carbonaceous surface thanks to a cation–π interaction. This interaction is a driving force coupling particles with positive charges and molecules with π-electrons.⁴⁹ Radial distribution function (RDF) derived for the centers-of-mass of GRA and the imidazole ring of [EMIM]⁺ exhibits a peak at 0.34 nm (Fig. 3) that corresponds precisely to the most favorable position cation–π interaction. The distance of 0.34 nm is approximately equal to sum of the van der Waals radii of carbon (belonging to GRA) and nitrogen (belonging to EMIM⁺). No second peak exists in this RDF. Thus, the second solvation layer of GRA is absent, confirming a short-range nature of this interaction. To understand other interactions in this system, Fig. 3 depicts RDFs between the GRA surface and selected interaction sites of [EMIM][BF₄]. Correlations of local density are observed in these cases, which however differ by their intensity. For instance, the first peaks corresponding to GRA–F (0.32 nm) and GRA–N (0.36 nm) interactions are reasonably sharp. The GRA–F interaction exhibits two peaks due to rigidity of the anion, i.e. the second peak is created by an opposite fluorine atom of the same particle.

Fig. 3 (Left) Radial distribution function (RDF) derived between GRA and the imidazole ring of the cation. (Right) RDF derived between GRA and selected atoms of [EMIM][BF₄].

Isolated GRA in IL is immobile exhibiting limited flexibility, according to Fig. 4. The root-mean-squared deviation (RMSD) was computed over an entire trajectory length, 100 ns. RMSD rises to 0.3 nm within ca. 4 ns and fluctuates around this value until the end of the MD simulation. Fluctuations appear significant, 0.1 nm < RMSD < 0.4 nm. Given the size of the simulated GRA model, the recorded RMSD is marginal. The intramolecular carbon–carbon–carbon–carbon potential plays an important role in reproducing natural dynamics of GRA over solvation and thermal motion.

Fig. 4 Root mean-squared (RMSD) and mean-squared deviations (MSD) of GRA as a function of the simulated time. The depicted MSD covers the sub-diffusion region of GRA.

Absolute translational mobility of GRA in IL is small. First, the size of GRA significantly exceeds the size of [EMIM][BF₄]. Second, [EMIM][BF₄] is a viscous IL (largely due to the cation–anion hydrogen bonding), whose genuine mobility is inferior to that of almost all molecular solvents. Third, the electrostatically maintained network of cations and anion, despite being to certain extent mobile, decreases a mean free diffusion path of any solute. Fig. 4 depicts the sub-diffusion part of the trajectory for GRA, which is observed up to more than 1000 ps. At longer times, the self-diffusion of GRA joins linear regime with a tiny slope.

The pairwise energy analysis in terms of decomposed energies between GRA and selected interaction sites of [EMIM][BF₄] is provided in Fig. 5. GRA is simulated without partial electrostatic charges with the dispersion energy represented by means of the 12-6 Lennard-Jones potential. Thus, Coulombic component of the interaction energy, which is obviously present in the cation–anion, cation–cation, and anion–anion interactions, is omitted in the GRA–IL interaction. Most of the sites feature binding energies between −150 and −250 kJ per mole of GRA. The exception is FB whose contribution amounts to −75 kJ mol⁻¹. By analyzing interactions between specific atom groups of [EMIM][BF₄] and GRA, we observed that the anion–GRA interaction is only −604 kJ mol⁻¹, while the interaction with the [EMIM]⁺ cation is much higher, −1902 kJ mol⁻¹. It is responsible for the largest fraction of the total solvation energy of GRA. The decomposition of the GRA–cation binding per intramolecular group (imidazole ring, CH₃ tail and C₂H₅ tail) confirms a dominant contribution of the imidazole ring. While anion can be further alternated to improve solubilization performance of IL, the imidazole ring must be kept unchanged. The evolution of these energies over simulated time is presented in Fig. 6.

Fig. 5 Decomposition of pairwise interaction energy between the solute and the solvent. (Left) per interaction site; (right) per group. BF – boron, FB – fluorine, CA – carbon in the imidazole ring, ME – carbon of the methyl tail, CH₂ – carbon of the methylene group of the ethyl tail, ET – carbon of the CH₃ group of the ethyl chain.

Fig. 6 The [EMIM][BF₄]–GRA interaction energy as a function of time.

Analysis of pairwise interactions reveals important aspects of solution energetics. However, it does not fully describe solvation thermodynamics, i.e. entropic fraction of the free energy. Therefore, we calculated solvation free energy of GRA in [EMIM][BF₄]. The considered thermodynamic process is transfer of GRA from vacuum to IL. This energy was found to be −969 kJ mol⁻¹ corresponding to −77.3 kJ per squared nanometer of the carbonaceous surface. Furthermore, exfoliation was simulated in real-time to record its cost step-by-step. Two GRAs were equilibrated parallel to one another in [EMIM][BF₄]. Next, they were gradually separated (Fig. 2) employing an external force (the steered molecular dynamics simulation). Each window was sampled extensively as described in the methodology to produce PMF (Fig. 7). The exfoliation free energy was defined as a difference between the plateau region (see 2.38 nm) and the minimum energy (see 0.35 nm), providing −70.1 kJ nm⁻². Note that PMF constitutes a smooth curve without local minima, which might have been related to indirect interactions mediated by the IL. Decomposition of free energy into enthalpic and entropic components shows that the enthalpic component prevails. It must be kept in mind that normalization of entropy per unit of surface area is not completely physical. Larger surfaces may induce larger entropy depending on the size of solvent particles and other structural considerations (reorganization of solution microstructure due to incorporation of the large solute). Note that free energy plateau is located far beyond 1.2 nm, i.e. where the GRA–[EMIM][BF₄] van der Waals binding occurs directly. The perturbation of the IL microstructure following exfoliation remains significant until the carbonaceous sheets separate by over 2 nm in the lateral direction. This distance may appear even larger in more ordered ionic liquids or at lower temperatures. ΔG = −70.1 kJ nm⁻² is a significant free energy cost, which motivates search for improved solvents.

Fig. 7 The exfoliation potential of mean force (black line), enthalpy (blue line), entropic contribution (red line), ΔG = ΔH − TΔS.

Conclusions

In this work, we devised a unified methodology to compare ILs as solvents for exfoliation of GRA. Exfoliation free energy is a fundamental thermodynamic descriptor of the coveted process. Exfoliation free energy vs. distance (potential of mean force) provides detailed information regarding an energetic cost of an entire process. Thus, all available ILs can be readily rated based on these parameters: the smaller energy is consumed for exfoliation, the more suitable solvent the corresponding IL is. The proposed descriptor implicitly accounts for all meaningful factors, e.g. cation–anion interaction strength, thermal motion effect, affinity to GRA, etc. As any other thermodynamics-based quantity, it does not provide any information regarding kinetics and kinetically forbidden events. Exfoliation free energy can be simply related to the exfoliation equilibrium constant, K_exf = exp(−G_exf/RT), in the ideal-behavior limit. K_exf, in turn, allows to compute GRA concentration after exfoliation.

Spontaneous exfoliation of GRA will occur if PMF decreases as distance between centers-of-mass of the neighboring sheets increases. This case, however, is unlikely possible due to the entropic factor, which prohibits distribution of large molecules within small molecules. While interactions between GRA and nonpolar organic solvents are similar with the interactions between two GRA sheets, per unit of volume, solubilization of GRA in these solvents is poor. This somewhat counter-intuitive phenomenon can be understood as a complication for smaller solvent particles to create voids, which could accommodate a large solute particle. Solvents must be used in conjunction with other established techniques of colloid chemistry, such as centrifugation and sonication. It is a particularly important role of the solvent to prevent backward aggregation of GRA.

Simulations of GRA exfoliation in [EMIM][BF₄] IL indicate that the energy investment of 70 kJ nm⁻² is necessary to separate one GRA sheet from another. In the presence of more GRA sheets, the needed energy will be even higher. A possible way to decrease this value is to decrease the cation–anion binding, e.g. via spatial separation of charged moieties. A perspective of the IL–GRA binding energy increase is unclear, as long as only van der Waals interactions can be relied on. Halogenation of the cation may be of certain interest for chasing this purpose.

Acknowledgements

V. V. C. is CAPES fellow. E. E. F. is funded through CNPq and FAPESP.

References

1. V. Palermo, Chem. Commun., 2013, 49, 2848–2857.
2. K. S. Novoselov, V. I. Fal'ko, L. Colombo, P. R. Gellert, M. G. Schwab and K. Kim, Nature, 2012, 490, 192–200.
3. O. V. Prezhdo, Surf. Sci., 2011, 605, 1607–1610.
4. P. A. Denis, J. Phys. Chem. C, 2013, 117, 3895–3902.
5. P. A. Denis and F. Iribarne, J. Mater. Chem., 2012, 22, 5470.
6. S. Ravula, S. N. Baker, G. Kamath and G. A. Baker, Nanoscale, 2015, 7, 4338–4353.
7. U. Maitra, H. S. Matte, P. Kumar and C. N. Rao, Chimia, 2012, 66, 941–948.
8. P. Kumar, RSC Adv., 2013, 3, 11987.
9. X. Cui, C. Zhang, R. Hao and Y. Hou, Nanoscale, 2011, 3, 2118–2126.
10. Y. Hernandez, V. Nicolosi, M. Lotya, F. M. Blighe, Z. Sun, S. De, I. T. McGovern, B. Holland, M. Byrne, Y. K. Gun'Ko, J. J. Boland, P. Niraj, G. Duesberg, S. Krishnamurthy, R. Goodhue, J. Hutchison, V. Scardaci, A. C. Ferrari and J. N. Coleman, Nat. Nanotechnol., 2008, 3, 563–568.
11. J. N. Coleman, Acc. Chem. Res., 2013, 46, 14–22.
12. G. Colherinhas, E. E. Fileti and V. V. Chaban, Phys. Chem. Chem. Phys., 2015, 17, 17413–17420.
13. G. Colherinhas, E. E. Fileti and V. V. Chaban, J. Phys. Chem. Lett., 2015, 6, 302–307.
14. S. Stankovich, D. A. Dikin, R. D. Piner, K. A. Kohlhaas, A. Kleinhammes, Y. Jia, Y. Wu, S. T. Nguyen and R. S. Ruoff, Carbon, 2007, 45, 1558–1565.
15. J. R. Lomeda, C. D. Doyle, D. V. Kosynkin, W. F. Hwang and J. M. Tour, J. Am. Chem. Soc., 2008, 130, 16201–16206.
16. V. C. Tung, M. J. Allen, Y. Yang and R. B. Kaner, Nat. Nanotechnol., 2009, 4, 25–29.
17. R. Hayes, G. G. Warr and R. Atkin, Chem. Rev., 2015, 115, 6357–6426.
18. V. V. Chaban and O. V. Prezhdo, J. Phys. Chem. Lett., 2013, 4, 1423–1431.
19. N. S. M. Vieira, P. M. Reis, K. Shimizu, O. A. Cortes, I. M. Marrucho, J. M. M. Araújo, J. M. S. S. Esperança, J. N. C. Lopes, A. B. Pereiro and L. P. N. Rebelo, RSC Adv., 2015, 5, 65337–65350.
20. K. Shimizu, C. E. Bernardes, A. Triolo and J. N. Canongia Lopes, Phys. Chem. Chem. Phys., 2013, 15, 16256–16262.
21. M. Tariq, M. G. Freire, B. Saramago, J. A. Coutinho, J. N. Lopes and L. P. Rebelo, Chem. Soc. Rev., 2012, 41, 829–868.
22. H. Yang, C. Shan, F. Li, D. Han, Q. Zhang and L. Niu, Chem. Commun., 2009, 3880–3882,  10.1039/b905085j.
23. E. Vazquez, F. Giacalone and M. Prato, Chem. Soc. Rev., 2014, 43, 58–69.
24. J. Texter, Curr. Opin. Colloid Interface Sci., 2014, 19, 163–174.
25. A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao and C. N. Lau, Nano Lett., 2008, 8, 902–907.
26. A. S. Pensado, F. Malberg, M. F. C. Gomes, A. A. H. Pádua, J. Fernández and B. Kirchner, RSC Adv., 2014, 4, 18017.
27. Y. Zhao and Z. Hu, J. Phys. Chem. B, 2013, 117, 10540–10547.
28. M. Tunckol, J. Durand and P. Serp, Carbon, 2012, 50, 4303–4334.
29. A. T. Najafabadi and E. Gyenge, Carbon, 2014, 71, 58–69.
30. X. Zhou, T. Wu, K. Ding, B. Hu, M. Hou and B. Han, Chem. Commun., 2010, 46, 386–388.
31. X. Wang, P. F. Fulvio, G. A. Baker, G. M. Veith, R. R. Unocic, S. M. Mahurin, M. Chi and S. Dai, Chem. Commun., 2010, 46, 4487–4489.
32. D. Nuvoli, L. Valentini, V. Alzari, S. Scognamillo, S. B. Bon, M. Piccinini, J. Illescas and A. Mariani, J. Mater. Chem., 2011, 21, 3428–3431.
33. N. Liu, F. Luo, H. Wu, Y. Liu, C. Zhang and J. Chen, Adv. Funct. Mater., 2008, 18, 1518–1525.
34. T. Kim, H. Lee, J. Kim and K. S. Suh, ACS Nano, 2010, 4, 1612–1618.
35. J. Lu, J. X. Yang, J. Wang, A. Lim, S. Wang and K. P. Loh, ACS Nano, 2009, 3, 2367–2375.
36. A. Taheri Najafabadi and E. Gyenge, Carbon, 2015, 84, 449–459.
37. G. Kamath and G. A. Baker, Phys. Chem. Chem. Phys., 2012, 14, 7929–7933.
38. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182.
39. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
40. V. V. Chaban, I. V. Voroshylova and O. N. Kalugin, Phys. Chem. Chem. Phys., 2011, 13, 7910–7920.
41. C. H. Bennett, J. Comput. Phys., 1976, 22, 245–268.
42. T. C. Beutler, A. E. Mark, R. C. van Schaik, P. R. Greber and W. F. van Gunsteren, Chem. Phys. Lett., 1994, 222(6), 529–539.
43. E. Lindahl, B. Hess and D. van Der Spoel, J. Mol. Model., 2001, 7, 306–317.
44. H. J. C. Berendsen, D. van der Spoel and R. van Drunen, Comput. Phys. Commun., 1995, 91, 43–56.
45. M. R. Shirts, J. W. Pitera, W. C. Swope and V. S. Pande, J. Chem. Phys., 2003, 119, 5740.
46. G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187.
47. S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman and J. M. Rosenberg, J. Comput. Chem., 1992, 13, 1011.
48. H. J. C. Berendsen, D. Vanderspoel and R. Vandrunen, Comput. Phys. Commun., 1995, 91, 43–56.
49. V. V. Chaban, C. Maciel and E. E. Fileti, J. Solution Chem., 2014, 43, 1019–1031.

---