Icephobicity of Functionalized Graphene Surfaces
Manipulating the ice nucleation ability of liquid water by solid surface is of fundamental importance, especially in the design of icephobic surfaces. In this paper, the icephobicity of graphene surfaces functionalized by sodium ions, chloride ions, or methane molecules is investigated using molecular dynamics simulations. The icephobicity of the surface is evaluated by the freezing temperature. The freezing temperature on surface functionalized by methane molecules decreases at first and then increases as a function of the number groups, while the freezing temperature increases monotonically as a function of the number groups upon surfaces functionalized by sodium ions or chloride ions. The difference can be partially explained by the potential morphologies near the surfaces. Additionally, the validity of indicating the ice nucleation ability of water molecules using the number of six rings in the system is examined. Current study shows that the ice nucleation upon functionalized surfaces is inhibited when compared with smooth graphene substrate, which proves the feasibility of changing the icephobicity of the surfaces by functionalizing with certain ions or molecules.
Figuring out the mechanism of ice nucleation on foreign surfaces is fundamentally important in many fields. A highly valuable application is the rational design of icephobic surfaces. In recent years, considerable research attention has been focused on the development of superhydrophobic surfaces which are expected to exhibit icephobicity because of the materials’ water repellent properties [1–4]. However, without detailed knowledge of the dependence of the icephobicity on surface morphology and surface physicochemical property, reliable and durable icephobic surfaces have not been developed yet. Alternatively, there is another possibility of developing icephobic surfaces by imitating the icephobic ability of antifreezing proteins (AFPs) [5, 6]. Experiments show that the antifreezing ability of the AFPs is attributed to some specific functional groups , which suggests that one might be able to obtain icephobic surfaces through modifying surfaces with some functional groups without the function of protein folding.
To uncover what kinds of groups determine the anti-icing ability of AFPs, studies in molecular level are needed. However, all simulation studies to date are focused on the influence of AFPs on ice crystal growth [7–12], while studies examining the influence of functional groups on ice nucleation ability are absent. One of the reasons is the extensive computational cost when trying to produce a nucleation process using common water models with charged atoms, such as SPC/E model  and TIP4P model .
To address this limitation, scholars have created many coarse-grained water models [15–18]. The monoatomic water (mW) model established by Molinero and Moore has attracted massive attention . The mW model does not only reproduce the properties of liquid water, such as surface tension, maximum density temperature, and accurate melting point but also describe the structure of water molecules in liquid or crystal state very well [18–20]. The mW model has been used in studying the ice nucleation of bulk water , undercooled drops , solutions , and icing on structured surface  and on modified graphene sheets . Meanwhile, reliable interactions between mW water and Na+ , , CH4 , and LiCl  enable the quantitative examination of ice nucleation behavior on functionalized surfaces.
In this work, we perform a series of ice nucleation simulations on graphene surfaces functionalized with Na+, , and CH4 groups. The groups are arranged randomly on a graphene sheet. Results show that ice nucleation is suppressed on the functionalized surfaces, and surface wettability does not appear as a reliable predictor of the icephobicity. The ice nucleation ability of the supercooled water could be determined by the ability to form two-dimensional ice nuclei near the surface. These findings offer new insights into the influence of functional groups on ice nucleation and facilitate the rational design of icephobic surfaces.
2. Simulation Method
2.1. Simulation Domain
As shown in Figure 1, water molecules (5072 mW) are placed on a smooth or modified graphene sheet. The number of carbon atoms in the sheet is 1008. To model functionalized graphene surface, 126, 168, 252, or 504 carbon atoms in the sheet are randomly replaced by sodium ions, chloride ions, or methane molecules, respectively. The positions of the replaced carbon atoms are maintained after replacement with different ions or molecules. This method has been successfully used to investigate the influence of hydroxyl-like groups on ice nucleation . The sodium and chloride ions are not able to form stable tetrahedral structure with the water molecules  and will be dissolved out from the solution during ice formation , while the methane molecule is a kind of hydrophobic molecule . The surfaces modified by sodium ions, chloride ions, or methane molecules are expected not to benefit ice nucleation. In organic chemistry, functional groups determine the characteristic chemical reactions of corresponding molecules. In our study, the difference in ice nucleation ability of modified and smooth graphene is due to the ions or methane molecules. We refer to ions and molecules as functional groups for convenience in the following discussion.
The dimensions of the simulation box are 5.518, 3.105, and 4.30 nm in the -, -, and -axis directions, respectively. NVT ensemble is used in the simulation and periodic boundary conditions are applied in and directions only. The interaction potential between water molecules and sodium or chloride ions is extracted from previous study , whereas the potential between water and methane molecules is reported elsewhere . The interaction between water molecules and carbon atoms in the graphene is described by the simple 12-6 Lennard-Jones (LJ) potential  , where is the distance between the carbon atoms and water molecules, nm is the distance parameter, and kJ/mol is the depth of potential well. The cut-off distance of the LJ potential is 1.0 nm. The interactions between carbon atoms are ignored. The surface atoms are fixed during the simulation.
In the simulations, the system is thermalized at 300 K for 1 ns. Then, the system is cooled at a constant cooling rate of 1 K/ns until crystallization is achieved. Moore et al.  suggested that the highest cooling rate at which ice nucleation can be observed is about 2 K/ns. The Nose-Hoover thermostat is applied to control the temperature and the velocity-Verlet algorithm is used to update the position and velocity of the molecules. To monitor ice nucleation, configurations are stored every 200 steps with the time step of 10 fs. Four simulations are conducted for each kind of surface. All simulations are carried out by using the LAMMPS package .
2.2. Determination of Freezing Temperature
To monitor the ice nucleation, we require an effective identification parameter which is not sensitive to the crystal structure and orientation for distinguishing between the liquid and solid atoms. Such parameter must also be effective in recognizing the ice molecules near the surface. Based on these considerations, we choose the bond order parameter proposed by Steinhardt et al.  and developed by ten Wolde et al.  who defined a local orientation parameter for each particle where and is the spherical harmonics and is the unit vector connecting water molecule and one of its nearest neighboring atoms . The neighboring atoms of are those within a cut-off distance nm of . is the number of the nearest neighbors of atom . Herein, is selected to identify the structures. Then, if product between atoms and is greater than 0.5, they are considered to be solid connected, where is the complex conjugate of . If one atom involves at least solid-connected neighboring atoms, the atom is a solid atom. To recognize surface solid atoms effectively, is chosen to be 7. Two connected solid atoms belong to the same crystal cluster. The evolution of the maximum-sized solid cluster is monitored and the temperature when the size of the maximum cluster changed abruptly is recognized as the freezing temperature [19, 24]. The homogeneous freezing temperature for bulk water is also determined using the same procedure.
2.3. Hydrophobicity of Different Surfaces
The hydrophobicity of the surfaces can be characterized by the contact angles of water drops on these surfaces. To determine the contact angle of the water droplet on different surfaces, we first obtain the density field from simulation trajectories at 300 K in a cylindrical coordinate system by cutting plane into square bins (width 0.1 nm) and delivering an average density using 1000 sequential configurations. Then the droplet periphery is defined by the counter line with a density of about 0.5 g/cm3 in the density field. As shown in Figure 2, the periphery is fitted by circle function to determine the contact angle. The position of the surface is at nm. To eliminate the influence of density fluctuation near the surface, only the data with nm are accounted in the fitting . The fitted lines are extended to intersect with polar axis ( nm). As shown in Figure 2, the angle between the polar axis and the tangent line of the circle in the intersection point is defined as contact angle. As shown in Figure 2, the contact angle of the water drop on the smooth graphene (°) is the largest, while those on surfaces modified by sodium ions and chloride ions are the smallest. The degree of contact angle on the smooth graphene surface agrees well with that of quantum simulation result (87°) and slightly smaller than that on graphene surface in experiments (92°) [34, 35]. The contact angles on the methane-modified surface decrease along with the number of methane molecules. The contact angles for the surfaces not shown in the figure are zero.
3. Results and Discussion
The freezing temperatures are shown in Figure 3(a). The homogeneous freezing temperature and the freezing temperature on smooth graphene are K and K, respectively, which are consistent with those in previous studies [19, 24]. As shown in the figure, all freezing temperatures on the solid surfaces are higher than that in homogeneous condition, suggesting that ice nucleation is promoted in heterogeneous conditions. Although the contact angle is the largest on the smooth graphene, the freezing temperature is higher than that on the modified surfaces. Compared with the icing on smooth graphene, ice nucleation on functionalized surfaces is inhibited, thereby demonstrating that icephobic surfaces could be obtained by treating surfaces with specific ions or molecules. The icephobicity of a surface can be characterized by the retardation of ice formation or the strength of ice adhesion, whereas the hydrophobicity of a surface is mainly determined by the static contact angle. Although many studies suggest that the freezing of water droplets is inhibited on hydrophobic surfaces [36, 37], results herein show that the hydrophobicity of a surface is not a reliable indicator of its icephobicity. Li et al. also found that the nucleation rate of water droplets on hydrophobic surfaces is higher than that on hydrophilic surfaces in experiments . The freezing temperature on the surface modified with methane molecules decreases first and then increases with the number of groups, while on the sodium or chloride ion-functionalized surfaces it increases along with the number of ions. To check the dependence of the freezing temperature on the distribution of the replaced positions of the carbon atoms, we construct three different distributions for the same number of groups. The results for the surfaces modified with sodium ions are shown in Figure 3(b). As presented, the distribution of the replaced positions does not change the freezing temperature obviously.
To discern the influence of functional groups on freezing, we analyze the growth of ice nuclei on different surfaces. As shown in Figure 4, during the ice nucleation process, a layer of two-dimensional (2D) ice molecules emerges first and then a three-dimensional (3D) ice nucleus grows steadily from the top of it. Ice nucleation on all surfaces follows this kind of pattern which is called the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) model . This kind of nucleation phenomenon deviates from the description of classical nucleation theory and is also observed in experiments  and simulations .
The ice nucleation ability near the surface likely depends on the ability to form 2D ice molecules on the surface. To verify this conjecture, we compare the ability to form dimensional ice molecules among different surfaces under different temperatures. The systems are equilibrated at 300 K for 2 ns and then simply quenched to 205, 210, 215, and 220 K for 2 ns, respectively. The number of ice-like molecules in the layer where 2D ice molecules form is accounted for. The remaining settings of the simulation are the same as those described in Section 2. Figure 5 shows that the fraction of ice-like molecules increases when the temperature decreases on the same surface. The ability to form ice-like molecules near the surface is higher at lower temperature. Furthermore, under the same temperature, the relative magnitudes of the fractions of ice-like molecules on different surfaces are consistent with the trend of freezing temperature. That is, the freezing temperature is higher as the fraction of ice-like molecules is lower.
Seeley and Seidler suggested that the icing ability of water droplets on alkanes with carbon atoms ranging from 25 to 28 depends on the 2D ice molecules near the surface . The ice nucleation rate in the drops is determined by the diffusion ability of water molecules near the surface. By freezing water drops on a silicon surface and surface modified by 1H,1H,2H,2H-perfluorodecyl-triethoxysilane (FAS-17), Li et al. also pointed out that the ice nucleation rate near the surfaces at low temperature depends on the diffusion of water molecules near the surface . The nucleation rate defined as the number of nucleation events occurring in unit time and unit volume can be expressed as , where is the kinetic prefactor, is the nucleation barrier, is the Boltzmann constant, and is the freezing temperature. The kinetic prefactor represents the diffusion of water molecules near the ice nuclei and is proportional to the inverse of the viscosity of water molecules . The apparent viscosity of water molecules in the layer where 2D ice-like molecules form can be calculated bywhere is the volume of water molecules in the occupied layer and is the pressure tensor on each water molecule, with and denoting , , or . The apparent viscosities at 240 K along with the number of functional groups are shown in Figure 6.
Figure 6 implies that the apparent viscosity of water molecules near smooth graphene surface (unfunctionalized) is smaller than that of bulk water because of the smaller interaction strength between water molecules and carbon atoms. However, the freezing temperature on the graphene surface is much higher than that of bulk water. This finding suggests that the nucleation barrier is reduced apparently in the heterogeneous ice nucleation on the graphene sheet. The viscosities near the surfaces functionalized with sodium and chloride ions increase along with the number functional groups and are about 2~3 orders higher than that on the smooth graphene sheet, which can partially explain the lower freezing temperature on the former surfaces. Forming tetrahedral structures with water molecules is also more difficult with sodium ions than with chloride ions. This behavior may lead to the lower freezing temperatures on the surface modified by sodium ions than those by chloride ions. However, the formation of 2D ice-like water molecules is not only controlled by the viscosity of water molecules in the vicinity of the surface. For instance, the freezing temperatures on surfaces modified by 504 methane molecules and sodium ions are fairly close but the viscosities near the corresponding surfaces differ considerably.
The freezing temperature on the methane-modified surfaces decreases initially () but increases while the number of methane groups increases further. The potential morphology near the surface is another factor that may affect the ice nucleation behavior of water molecules near the modified surfaces. As mentioned above, the freezing temperature is mainly affected by the formation of 2D ice molecules. When the number of methane groups is small, the morphology of the potential surface is controlled by the carbon atoms. The potential morphology becomes rougher as the number of methane groups increases. It disturbs the formation of 2D ice molecules and lowers the freezing temperature accordingly. If the number of methane groups becomes sufficiently large, the potential morphology will be controlled by the methane groups. The formation of the 2D ice becomes easier because of the smoother potential morphology. The freezing temperature will increase again with increased number of methane groups. The interaction strength between water molecules and methane groups is slightly stronger than that between water molecules and carbon atoms; therefore, the critical number of methane groups will be smaller than half of the number of carbon atoms in smooth graphene. On the surfaces modified with sodium or chloride, the potential morphology is constantly controlled by the sodium or chloride ions due to their quite strong interaction with the water molecules. Thus, the freezing temperature should increase along with the number of sodium or chloride ions.
Numerous previous studies tackled which kind of parameter is suitable for predicting the icephobicity of a surface. The proposed parameters include the tetrahedral parameter , the number of hydrogen bonds , and the number of six rings that connects water molecules in the system. Among these factors, the number of six rings have attracted the greatest attention [45, 46]. Scholars believe that a system is easier to freeze as the number of six rings increases. However, our simulations show that the number of six rings is not a reliable indicator of the ice nucleation ability of the surface. As shown in Figure 7, the distribution of six rings on the sodium-modified surface with is apparently stronger than that with . However, the freezing temperature on the surface with is a little higher than that with , suggesting that the ice nucleation more easily occurs on the surface with 168 sodium ions. Since the behavior of ice nucleation near the surface is affected by many factors, such as viscosity and potential morphology as mentioned above, it might be impossible to indicate the ability of ice nucleation upon the surface by one simple parameter only.
In this paper, the icephobicity of graphene surfaces functionalized with sodium ions, chloride ions, or methane molecules is investigated by using molecular dynamics simulations. The results reveal that the ice nucleation ability of the functionalized surfaces is weakened compared with that of the smooth graphene surface, depending on the number functional groups and type of functional group. These findings prove the feasibility of enhancing the icephobicity of surfaces by functionalizing with certain ions or molecules. Under simulation conditions, the ice nucleation process follows the KTHNY model which is characterized by an initial formation of a layer of 2D ice nucleus formed near the surface and subsequent steady growth of a 3D ice nucleus on top of the former structure. The ice nucleation ability of the surfaces is determined by the ability to form 2D ice layer. The ability of forming 2D ice nuclei is controlled by many factors, such as the apparent viscosity of the water molecules and the potential morphology near the solid-water interface. In this study, we also examine the validity of the use of the number of six rings in a system to indicate the ice nucleation ability of water molecules. The number of six rings fails to predict the icephobicity. It is difficult to indicate the ice nucleation ability of water molecules by one single parameter.
The authors declare that they have no competing interests.
This work is financially supported by the National Basic Research Program of China (no. 2015CB755801) and the National Natural Science Foundation of China (no. 51321002). The computations are carried out at the Tsinghua National Laboratory for Information Science and Technology, China.
H. Saito, K.-I. Takai, and G. Yamauchi, “A Study on ice adhesiveness to water-repellent coating,” Materials Science Research International, vol. 3, no. 3, pp. 185–189, 1997.View at: Google Scholar
J. B. Lu, Y. Q. Qiu, R. Baron, and V. Molinero, “Coarse-graining of TIP4P/2005, TIP4P-Ew, SPC/E, and TIP3P to monatomic anisotropic water models using relative entropy minimization,” Journal of Chemical Theory and Computation, vol. 10, no. 9, pp. 4104–4120, 2014.View at: Publisher Site | Google Scholar
P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, “Bond-orientational order in liquids and glasses,” Physical Review B, vol. 28, no. 2, pp. 784–805, 1983.View at: Google Scholar
A.-K. Löhmann, T. Henze, and T. Thurn-Albrecht, “Direct observation of prefreezing at the interface melt-solid in polymer crystallization,” Proceedings of the National Academy of Sciences of the United States of America, vol. 111, no. 49, pp. 17368–17372, 2014.View at: Publisher Site | Google Scholar
L. H. Seeley and G. T. Seidler, “Two-dimensional nucleation of ice from supercooled water,” Physical Review Letters, vol. 87, no. 5, Article ID 055702, 2001.View at: Google Scholar