Confinement effects of graphene oxide nanosheets on liquid–solid phase transition of water

Abstract

In this work, the liquid–solid phase transition temperature of water confined between two graphene oxide (GO) sheets is investigated using molecular dynamics simulations. The results demonstrate that, due to the presence of functional groups of the GO sheets, at temperatures below liquid–solid phase transition temperature, the water molecules near the confining sheets are not structured like ice and remained in the liquid form. The results also reveal the confinement effects on the melting and freezing rate of water molecules. The confining conditions delay the freezing process of the water molecules compared to the bulk water molecules. These results are confirmed by the calculated total energy, the density profile of water molecules confined between two GO sheets, the number of hydrogen bonding, radial distribution function, and mean square displacement. The liquid–solid phase transition temperature of water in the presence of GO sheets was calculated as 236 K, which was 34 K less than the temperature of the bulk water.

---

1. Introduction

Experimental studies have shown that the melting and freezing points of various fluids confined within the nanostructures are different than the values of the bulk fluids. Increasing the surface to volume ratio of confined systems changes phase equilibrium temperature.^1,2 The water confined in nanometric space shows different behaviors and properties of bulk water. Many simulations and experimental studies have been done on the water confined in hydrophobic and hydrophilic environments.^3–6 The study of water confined in porous silica has shown that the phase behavior of water changes and effect of volume reduction is significant. Also, freezing point of confined water is reduced compared with bulk water.⁷ Molinero et al.⁸ used molecular dynamic simulations to study the melting temperatures and coexistence of ice and a premelted liquid layer in the nanopores with radius 2 nm and a broad range of surface-water interactions, ranging from very hydrophobic to strongly hydrophilic. Their results indicated that the melting temperatures of ice confined in the pores are strongly decreased compared with the temperature of the bulk ice, 274 K, in the mW water model. The application of molecular dynamics simulations to study the water confined between hydrophobic and hydrophilic spaces has led to great success. Heterogeneous nucleation of ice on the surface of graphite was investigated in 2014 by Molinero et al.,⁹ the results of which showed that the graphitic surfaces encouraged the heterogeneous nucleation of ice and increased the freezing temperature. Molinero et al.¹⁰ also studied heterogeneous nucleation of ice on the surface of graphene containing hydroxyl groups. These surfaces are not a suitable option for this process, because they reduce the temperature of freezing water; but, they could play an indirect role for the heterogeneous nucleation of ice because of increasing hygroscopicity.

Recently, GO sheets, as a hydrophilic substance, have received great interest due to their applications in energy storage devices, sensors, etc.¹¹ Zheng et al.¹² studied the process of ice nucleation on superhydrophilic GO nanoflakes. Various ice solids can thus be grown in ambient conditions on superhydrophilic GO. They obtained real-time imaging of the continuous phase transition amorphous ice to a transient cubic ice Ic stage, which in the end turns into the stable hexagonal ice Ih. They discovered that ice nucleation and growth can be affected by modifying the functional groups of nano GO and by intermolecular hydrogen-bonding between nano GOs. This could provide a way to control heterogeneous ice nucleation and snow crystal formation. Due to the high ability of GO to disperse in water, it is assumed to be a hydrophilic material. According to the layered structure of GO in water, it can serve as a model for studying water confined in the nanoscale between hydrophilic surfaces.¹³ In our previous work, we studied the structural and dynamical properties of water confined between two GO sheets. Our results indicated that the confined water molecules show different structural and dynamical behaviors near the sheets, as the hydrogen bonds between water molecules near the confining surfaces decreases due to bonding the water molecules to the substituents of GO. Also, the orientation of molecules in the different distances is different from the confining surfaces.¹⁴ Determining the phase transition temperature of liquid–ice water confined between two GO sheets and also the structural and behavior of liquid and ice water confined between two GO sheets are our goal in the project.

2. Computational methods

2.1. Simulation details of bulk water

In order to determine the phase transition temperature of liquid–ice water confined between two GO sheets, a reference system using coexistence between solid and liquid¹⁵ was studied. A TIP4P/Ice¹⁶ model was used for the water molecules. This model properly reproduces many properties of water. For instance, the predicted value of water density is equal to 0.993 g cm⁻³ at T – 298.15 K and P – 1 atm, which is in good agreement with the experimental data. We also set up a box of liquid water with 1024 molecules and a box of ice Ih with 1024 molecules into contact with each other for the ice/liquid water initial configuration with 2048 molecules. The typical size of the simulation box was 72 × 30 × 28 ų. We oriented the solid such that the solid face in contact was the secondary prismatic {1−210} plane. The main advantage of this choice is that this face is the fastest growing face.¹⁷ Several NPT molecular dynamic simulations at constant pressure in 1 bar were performed at different temperatures in the temperature range of 260 K and 280 K. In our simulations, the temperature was fixed with a Nosé-Hoover thermostat with a relaxation time of 1 ps. Pressure was also fixed with a Nosé-Hoover barostat and relaxation time of 10 ps. Periodic boundary conditions were applied in all directions. Simulations with time step 1 fs were conducted for 20 ns. The van der Waals forces with the cutoff of 10 Å were calculated, and Ewald summation approach was used for calculating the electrostatic interactions in the system.

2.2. Simulation details of the mixture of ice/water confined between two GO sheets

In this study, the applied GO is one-layered with 840 carbon atoms as well as OH and epoxide functional groups on both sides. Geometry optimization was performed using the hybrid B3LYP functional.¹⁸ The Gaussian 03 program was employed with the internally stored 3-21G basis set. The obtained partial charges of oxygen and hydrogen atoms of the hydroxyl group are −0.524e and 0.328e, also, for the oxygen atom of the epoxy group, it is equal to −0.4e. These values are in agreement with the partial charges of OPLS-AA reported by Jorgensen.¹⁹ The interaction between water molecules and carbon was described by the 12−6 Lennard-Jones potential, with parameters ε_C–O = 0.1889 kcal mol⁻¹ and σ_C–O = 3.19 Å at an interatomic distance r. The initial configuration of our system consisted of two parallel surfaces of GO located at the distance of 20 Å and along the X–Y plane. There were 1000 molecules of water, half of which were in the liquid form and the other half were in the form of ice Ih. The Z dimension of the box was large enough and the other two were equal to the length and height of GO. Periodic boundary conditions were applied in all respects to omit surface effects. Cut off 10 Å was applied to all non-bonded interactions. Ewald sums were used to correct long-range interactions. To solve the equations of motion, we used the Verlet algorithm. Simulations were performed at different temperatures with NPT ensemble. A Nosé-Hoover thermostat and a Nosé-Hoover barostat fixed temperature and pressure, respectively. The simulations were done at the pressure of 1 atm and, according to the anisotropy of the system, Z-direction changed independently from the x- and y-directions. After the energy minimizing of the configurations, the simulations were done at different temperatures for 50 ns. Fig. 1 shows the configuration of the system consisting of the solid/liquid mixture of water confined between GO sheets. All the MD simulations were performed using LAMMPS.²⁰

Fig. 1 Configuration of the simulated system consists of mixing of liquid and solid water confined between two GO surfaces. Green, red, and blue colors correspond to carbon, oxygen, and hydrogen atoms, respectively.

3. Results and discussions

In this Section, the simulation results for both systems, the bulk water and confined water are presented. The results of molecular dynamics simulations are used to analyze and determine the phase transition temperature for both systems using the total energy curves, the number of ice molecules, radial distribution functions, mean square displacement, coordination numbers and number of hydrogen bonding.

3.1. Total energy

Recently, an interesting approach has been reported by Vega et al. for finding the melting point of ice.²¹ They considered improving total energy (sum of kinetic and potential energies) with respect to time before reaching equilibrium of the system using the method coexistence of liquid–solid interface and calculated the melting point of several atomic models of water. For instance, the energy of the TIP4P/2005 model changed slowly with respect to time at the temperatures 248 K and 250 K. Therefore, the system was in equilibrium in this area of temperatures, and the average of these two temperatures was introduced as the phase transition temperature. The phase transition temperature of the TIP4P/Ice model was also obtained as 270 K using the same approach. Although the obtained values for the water/ice mixture in compared to the pure ice shows more deviation than experimental results, the melting and freezing processes are faster in this approach and the needed time to reach equilibrium decreases. Also, supercooling and superheating phenomena are removed.

In this work, total energy curves were used to obtain the phase transition temperature. At the temperatures above the phase transition temperature, the thermostat provides the needed energy to the system to melt ice and the temperature is lower; the thermostat takes energy out of the system; so, the water freezes. At the phase transition temperature, the total energy of the system does not change over time. In order to obtain the phase transition temperature, at the beginning, two temperatures were chosen for which complete freezing/melting was seen. Then, each of these temperatures was decreased or increased by a few degrees of celsius to find the phase transition temperature. Fig. 2 exhibits the evolution of the total energy of the bulk system at the temperatures of 262 K, 268 K, 270 K, 271 K, and 276 K. As can be observed, the energy decreases with time at 262 K and water transforms to ice crystals. While 15 ns passed, the system does not still freeze. Thus, freezing process needs more time than the melting process. At 276 K, the ice molecules starts to melt, causing the increase of the total energy with time. It can be also seen that, before reaching the equilibrium, the slope of the energy curve increases, which shows the increase of the ice melting rate as one or two layers of ice has retained. At the temperatures 268 K and 271 K, evolution of the curves is minimized; so, the mean of these two temperatures, that is 270 K, is considered the solid–liquid phase transition temperature. Fig. 3 shows the snapshots of the initial configuration and also the final configuration of the bulk system at two different temperatures.

Fig. 2 Evolution of the total energy of the bulk water at different temperatures.

Fig. 3 Snapshots of the system in (a) the initial configuration, (b) 276 K, and (c) 262 K, red and white colors correspond to oxygen and hydrogen atoms, respectively.

Fig. 4 illustrates the evolution of the total energy with respect to time for water confined between two GO sheets. Energy curves increase at the temperatures of 238 K and 239 K, indicating ice is melting. At the temperatures of 232 K and 234 K, energy decreases with respect to time. The system is in the equilibrium state between solid and liquid in the temperature area of (235–237) K and the energy curves show constant and mild fluctuations. The average of these two temperatures, 236 K, was considered the solid–liquid phase transition temperature of water confined between two GO sheets. Fig. 4 also shows that confinement conditions retard the freezing process, and more time – upper than 50 ns – is needed to freeze the water molecules at the temperatures lower than the phase transition point, although this process occurs faster in the bulk water.

Fig. 4 Evolution of the total energy of confined water between two GO sheets at different temperatures.

3.2. The number of ice molecules

The other approach to determine the melting point is to count the fraction of ice molecules or liquid water molecules as a function of time. To analyze the evolution, to obtain the mechanism of ice Ih melting process, an algorithm has been previously developed by Ba'ez and Clancy.^22,23 This algorithm determines the number of ice and water molecules of the system at each time. The tetrahedral order parameter is defined as F_i,is calculated for each water molecule i. In eqn (1), n_i is the number of water molecules in the first solvent shell of the water molecule i. j and k are indices running over these nearest neighbors and θ_jik is the angle between the oxygen atoms of j, i, and k monomers. For tetrahedral bonding, there are six angles and the cosine of θ_jik angle is close to −0.33. Thus, F_i provides a measure of deviation from perfect tetrahedral coordination (F = 0) and serves to distinguish solid-like (F < 0.4) and liquid-like (F > 0.4) water molecules. Fig. 5 shows the number of bulk water molecules in the ice phase as a function of time for different temperatures. At 270 K, the number of ice molecules, except for small fluctuations around a mean value, has almost been unchanged. The solid phase is totally melted at 276 K. Increasing the curve at the temperature of 262 K indicates the system is freezing.

Fig. 5 Number of bulk water molecules in the ice phase. The temperatures of the three upper curves are – from up to down – are 262, 268, and 270 K, respectively. The lower curve corresponds to a temperature of 276 K.

Fig. 6 depicts the change in the number of confined water molecules in the ice phase with respect to time for several temperatures. The number of confined water molecules in solid state decreases with time at 239 K which is above the ice melting temperature and this process, like the bulk water system, is faster than the freezing process at 234 K. At the temperature of 237 K which is close to the phase transition temperature, the number of ice molecules retains constant over time. Comparison of Fig. 5 and 6 also confirms that the freezing process of confined water is longer than it in the bulk water.

Fig. 6 Number of confined water molecules in the ice phase. The temperatures of the curves from up to down are 234, 237, 238 and 239 K, respectively.

3.3. Density profile

In this Section, we calculate the density values of water molecules for the bulk water in the temperature range of 262–277 K and, for the confined water, in the temperature range of 232–240 K. Liquid density decreases generally as temperature is increased at the given pressure; yet, the experimental density curve of water behaves abnormally as the water reaches the maximum density at the temperature of 4 °C. In the ice/water mixture, if the ice melts, the volume of the system will decline and the density will therefore rise; however, if water freezes the density will decrease. Fig. 7 shows the density curve with respect to temperature in the bulk and confined systems. According to these curves, the system density at the temperatures, upper and lower than the phase transition temperature shows a steady decrease. At 270 K for the bulk water and 236 K for the confined water, due to the solid–liquid equilibrium, the density of the system is between the density of the upper and lower temperatures.

Fig. 7 Density curve with respect to temperature in (a) the bulk and (b) the confined systems.

3.4. Coordination number

The area under the first peak of oxygen–oxygen radial distribution function, O–O RDF, defined as the coordination number, rises slightly upon melting of ice Ih.²⁴ In water, this number increases with temperature, yet remains not more than about 4, the coordination number in crystalline ice. The coordination number for bulk water was obtained, in agreement with experimental data and simulation calculations,¹⁴ which was equal to 4.7. The coordination number diagram is plotted in Fig. 8 for both the system of confined and bulk water. For both, this number changes interestingly at the phase transition temperature, indicating that the number of the neighbors of the central liquid water increased compared with the crystal ones and this quantity for the confined water is indeed less than that of the bulk water. The average coordination number for bulk is 4.04–4.22 in the temperature range of 262–277 K; it is also 3.95–4.1 in the temperature range of 232–240 K for the confined system. At the phase transition temperature, the coordination number equals 4.14 and 4.015 for the bulk and confined water, respectively. Besides, the data for both of the systems indicates that the number of neighbors of central confined water decreases with the bulk water one.

Fig. 8 (a) Coordination number of bulk water and (b) confined water, at different temperatures.

3.5. The Z-axis density profile

Z-density profile evaluated using the average number of atoms in a bin with thickness ΔZ in Z-direction and with width of ΔZ = 0.2 was used to evaluate density. Fig. 9 presents the density profile of the oxygen atoms of water molecule confined between two surfaces of GO at temperature of 230 K (solid line) and 240 K (dashed line). At 240 K at which the system is liquid, the density near the surfaces is more than the density in other regions. In the density profile, three different groups of water molecules can be observed. In the area between −8 Å and 8 Å, there is a straightforward line and the amount of density agrees with the bulk water density. At the distance of 8 Å to 14 Å (also the density of −8 Å to −14 Å), one peak is seen in the middle part of GO. Furthermore, one peak with maximum height in the area between 14 Å and 17.5 Å (and −14 Å and −17.5 Å) is observed beside the confining surfaces. There is no water molecule at the region of 0 A to 2.5 Å, due to the repulsion effects causing by the intersections of water molecules with the confining surfaces. At 230 K at which the system is solid, the density of water molecules along the z-axis marks by the solid line. As can be seen, the density of water molecules in the first and the second areas near the surface behaves in agreement with the temperature of 240 K, at which there is only liquid water. It can be derived that the presence of substituents on the GO surfaces affects the water molecule structure and there is no Ih structure in this area. In the area of −8 Å to 8 Å, the drastic and ordered fluctuations of the density confirm the presence of ice. Fig. 10 shows a snapshot of the solid–liquid water mixture system between GO sheets after 50 ns at 230 K. The disorder structure of Ih ice near the confining sheets can be seen in Fig. 10.

Fig. 9 Z-density distribution of the oxygen atoms of water molecule confined between two GO surfaces at temperature of 230 K (solid line) and 240 K (dashed line).

Fig. 10 Snapshot of the solid–liquid water mixture system between GO sheets after 50 ns at 230 K.

3.6. Number of hydrogen bonds

In the ice Ih, each water formed four hydrogen bonds with O⋯O distance of 2.76 Angstrom to the nearest oxygen neighbor. In ice Ih, each molecule is hydrogen bonded to four other molecules. However, in the liquid water, that number varies from molecule to molecule. To find out the structure of water molecules at different distances to the confining GO sheets, the number of hydrogen bonds of water molecules at different distances to GO sheets was calculated. These calculations at the temperature before the phase transition temperature was done to find out whether there is Ih ice or not near the GO sheets. The related data to the hydrogen bonds per molecule at the different distances to GO sheets at 232 K are presented in Fig. 11. As mentioned above, in the area close to the sheets – 14.5 Å to 17.5 Å – the number of hydrogen bonds equals 3.24 and the presence of GO substituents has disturbed the structure of Ih ice in this area. Increasing the distance from the surfaces and decreasing the presence of substituents increases the number of hydrogen bonds per water molecule and this amount gets closer to 4, which is equal to the hydrogen bonds per molecule in the ice Ih.

Fig. 11 Number of hydrogen bonds of water molecules confined between two GO sheets.

3.7. Radial distribution function (RDF)

Fig. 12a and b show the oxygen–oxygen and hydrogen–hydrogen RDFs, respectively, for the ice/water mixture at the phase transition temperature and also at temperatures upper and lower than this point. Fig. 12a shows that the RDF at the temperature of 276 K is different from the RDFs at other temperatures. By increasing the temperature from 262 K to 276 K, the height of these peaks is decreased, but about 276 K, the reduction is most severe. Also, the position of the second peak at the temperature of 276 K is different from it at other temperatures. This issue can be observed in the hydrogen–hydrogen RDFs in Fig. 12b.

Fig. 12 RDFs of bulk water. (a) O–O RDFs. (b) O–H RDFs.

This event indicates a steep decrease of the accumulation of liquid molecules near each other compared with the solid ones. Thus, the correlation of the ice molecules with the neighbors is more than that of the liquid water molecules. At the temperature of 270 K, at which the system is in equilibrium, the first peak height was between the height of the solid and liquid peaks.

Fig. 13 shows the O–O RDF for the confined water at the temperatures close to the phase transition point. As mentioned above, the height of the first peak decreases with increasing the temperature. At the lowest studied temperature, 234 K, the height of this peak is about 9.5, which is the highest. The smallest peak is observed at 239 K with a height of 7.8. The O–O RDF of the water molecules for different distances of the confining surfaces, at the temperature of 234 K, is presented in Fig. 14. According to Fig. 14, the height of these peaks in the closest layer to the GO surface shows a steep decrease with respect to the second layer.

Fig. 13 RDFs for O–O of water molecules confined between GO surfaces.

Fig. 14 RDFs for O–O of water molecules confined between GO surfaces for the nearest layer to confining surfaces (solid line), and the second layer (dashed line).

3.8. Mean square displacement

Mean square displacement (MSD) can be used to indicate the melting process of a solid. For an ordered-structure liquid, the MSD gradually increases with time; but, for a solid, it fluctuates around a mean value.²⁵ Fig. 15 shows the outcome results of the MSD of the bulk water molecules (a) and the confined water molecules (b) at different temperatures. As can be observed, the slope increases with time at 276 K for bulk water and, at 239 K, for confined water, whereas at the temperatures of 262 K and 232 K, respectively, for the bulk water and the confined water, the curves, since the system is freezing, fluctuates around a constant value. Different scales of the vertical axes in two plots mentions that the confining GO surfaces effects a lot on the displacement of the water molecules. At the temperature of 239 K, in which the ice molecules have been melted, the mobility of the confined water molecules is much lower than the mobility of the bulk water molecules at the temperature of 276 K.

Fig. 15 MSD of water molecules for (a) bulk system, and (b) confined system.

4. Conclusion

In the present work, the solid–liquid phase transition temperature of the confined water between graphene oxide nanosheets was studied using molecular dynamic simulations and compared with the bulk water results. To achieve this aim, evolution of the total energy with time and change of the number of molecules in the ice phase were calculated. Our results show the phase transition temperature of 236 K is 34 K less than the phase transition temperature of bulk water and is really close to the melting point of confined water between hydrophilic cylindrical nanopores with 2 nm radius,²⁶ which has been recently reported. These results also demonstrated that the confinement effects on the melting and freezing rate of the water molecules so that the confined water molecules need more time to freeze. The radial distribution functions showed that the correlation of the confined water molecules, like the bulk water, in the crystal state is more than this correlation in the liquid water. The structural properties of ice molecules at the temperature below the phase transition point were studied. The outcome results from considering the density of water molecules in the Z-direction, number of hydrogen bonds, and radial distribution functions showed that even at the temperature below the solid–liquid phase transition, at which the system is solid, the structure of water molecules near the GO sheets is like that of liquid water and, in this region, there is no ordered-structure ice. The results showed that the confining surfaces reduce the displacement of the water molecules compared to the bulk water.

References

1. I. Brovchenko and A. Oleinikova, Interfacial and Confined Water, Elsevier, 2008.
2. O. Petrov and I. Furo, Phys. Rev. E, 2006, 73, 11608.
3. H. K. Christenson, R. G. Horn and J. N. Israelachvili, J. Colloid Interface Sci., 1982, 88, 79.
4. C. Faivre, D. Bellet and G. Dolino, Eur. J. Phys. B, 1999, 7, 19.
5. J. C. Dore, M. Dunn, T. Hasebe and J. H. Strange, Colloids Surf., 1989, 36, 199.
6. D. W. Aksnes and L. Gjerdakbr, J. Mol. Struct., 1999, 475, 27.
7. M. Heuberger, M. Zach and N. D. Spencer, Science, 2001, 292, 905.
8. E. Moor, J. Allen and V. Molinero, J. Phys. Chem. C, 2012, 116, 7507.
9. L. Lupi, A. Hudait and V. Molinero, J. Am. Chem. Soc., 2014, 136, 3156.
10. L. Lupi and V. Molinero, J. Phys. Chem. A, 2014, 118, 7330.
11. K. Chen, F. Liu, S. Song and D. Xue, CrystEngComm, 2014, 16, 7771–7776.
12. Y. Zheng, C. Su, J. Lu and K. Ping Loh, Angew. Chem., Int. Ed., 2013, 52, 1–6.
13. A. V. Talyzin, S. M. Luzan, T. Szabo, D. Chernyshev and V. Dmitriev, Carbon, 2011, 49, 1894.
14. M. Zokaie and M. Foroutan, RSC Adv., 2015, 5, 39330.
15. R. Garcia Fernandez, J. L. F. Abascal and C. Vega, J. Chem. Phys., 2006, 124, 144506.
16. J. L. F. Abascal, E. Sanz, R. Garcia Fernandez and C. Vega, J. Chem. Phys., 2005, 122, 234511.
17. H. Nada and Y. Furukawa, J. Cryst. Growth, 2005, 283, 242.
18. A. D. Deyoung, S. W. Park, N. R. Dhumal, Y. Shim, Y. Jung and H. J. Kim, J. Phys. Chem. C, 2014, 118, 18472.
19. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 2008, 101, 065502.
20. S. Plimpton, J. Comput. Phys., 1995, 117, 1.
21. R. Garcia Fernandez, J. L. F. Abascal and C. Vega, J. Chem. Phys., 2006, 124, 144506.
22. L. A. Baez and P. Clancy, Ann. N. Y. Acad. Sci., 1994, 715, 177.
23. L. A. Baez, Ph.D thesis, Cornell University, 1996.
24. T. head gordon and M. E. Johnson, PNAS, 2006, 103, 7973–7977.
25. M. P. Allen, Computer simulation of liquids, Oxford University Press, USA, 1989.
26. E. Moor, J. Allen and V. Molinero, J. Phys. Chem. C, 2012, 116, 7507.

---