Comparative study of the ReaxFF and potential models with density functional theory for simulating hexagonal ice

Computational Materials Science 177 (2020) 109546

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Comparative study of the ReaxFF and potential models with density functional theory for simulating hexagonal ice

T

⁎

Chunyang Wanga, Yanzhuo Xuea, Chaoying Wangb, , Duanfeng Hana a b

College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China College of Aerospace and Civil Engineering, Harbin Engineering University, Harbin 150001, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Ice Ih Molecular dynamics Water model TIP4P/ice Density functional theory Surface energy

This paper compares ReaxFF, SPC/E, TIP4P/2005 and TIP4P/Ice with Density Functional Theory (DFT) to identify which one should be used to model both the microstructures and mechanical properties of large size hexagonal Ice (Ih) over a long time scale (~10−9–10−6 s). According to the model performances with regards to determining the microstructure, surface energy and fracture toughness, the limitations and validities of ReaxFF and the three other potential models are discussed. It is found that hydrogen bond density determines the surface energy on the same surface. TIP4P/Ice has the smallest error and the ReaxFF force field has the largest error, up to 27.21% for the aforementioned quantities. If only modeling the microstructure of ice Ih, TIP4P/2005 and SPC/E are the most suitable water models. Moreover, the secondary prism face {112¯0} , with the lowest hydrogen bond density and fracture toughness, is the first cleavage plane of ice Ih. The primary prism face {1¯100} is the second cleavage plane, as shown by the surface energy and fracture toughness results of DFT and molecular dynamics (MD).

1. Introduction Ice, a compound with a diversified structure formed by simple molecules, has a rich material form. Nineteen kinds of ice phase have been discovered, such as the new ice XVIII [1] discovered in 2019 and ice XVII [2] in 2016, along with amorphous ice and glassy water [3–5]. Hexagonal ice (Ih) [6], the naturally occurring form of ice on Earth, can be seen in daily life, as it manifests as snowflakes, frost, and sea ice. Exploring the physical and mechanical properties of ice Ih is necessary for daily life, aviation, aerospace, navigation, polar science, new energy and other fields [7]. At present, ice has been investigated using several different methods, such as density functional theory (DFT) [8], molecular dynamics (MD) [9,10], the discrete-element method (DEM) [11], and peridynamics (PD) [12,13]. Since the 20th century, experimental methods have been widely used in the study of mechanical properties, such as ice microstructure [14–16], Poisson's ratio [17], the elastic modulus [18,19], and fracture toughness [20,21]. The pioneering work in brittle fractures can be traced to the classic theory of Griffith [22]. In 1969, Barker et al. first made a numerical simulation of liquid water [23]. After they used DFT to study liquid water [24], ice research became more prominent [25–28], and the experimental results of studies of ice microstructure were verified by the first-principles method [29,30]. Pan et al. [31,32] ⁎

and Sun et al. [33] studied the surface properties by DFT and reduced the surface energy range to 4.3–18.6 meV/Å2 by empirical and semiempirical methods [34–37]. It is extremely difficult to study water or ice using experiments on the nanometer scale [38], and the first-principles method cannot address large-scale or long-time-scale mechanical processes, such as stretching and fracture, due to the large atomic system and long-time-scale calculations, although it has great accuracy in its domain, without requiring empirical functions. A possible alternative for these cases may be applying an MD-based empirical potential function [39–41], which is less computationally expensive. MD simulations have been applied in the research of many new materials, such as graphene [42], single crystal silicon [43], carbon nanotube [44], silicon nanowires [45,46], and metals [47–50]. Many studies have focused on ice (or water) using MD [10,51–58], especially based on the SPC [56,59], TIP4P [60–63] and ReaxFF frameworks [57,64]. Qin et al. [57] investigated the effect of carbon dioxide on ice cracking by ReaxFF, a reactive force field for hydrocarbons. However, because this force field is mainly used in the calculation of aqueous systems [64], its accuracy in describing the properties of ice remains unverified. Vrbka et al. [56] modeled the process of brine rejection from freezing salt solutions and proposed a microscopic mechanism of brine rejection using the SPC/E water model. Abascal et al. proposed two general water models, TIP4P/2005 [60] and TIP4P/Ice [62], which

Corresponding author. E-mail address: [email protected] (C. Wang).

https://doi.org/10.1016/j.commatsci.2020.109546 Received 5 December 2019; Received in revised form 16 January 2020; Accepted 18 January 2020 0927-0256/ © 2020 Elsevier B.V. All rights reserved.

Computational Materials Science 177 (2020) 109546

C. Wang, et al.

were designed to cope with the condensed phases of water; their thermodynamic properties, phase diagrams, and properties at melting and vaporization were calculated and discussed, but their suitability to determine the microstructures and mechanical properties of Ih are uncertain. Moreira et al. [63] studied pre-melting phenomena in ice Ih, and Dong et al. [55] successfully imitated the icing process of water droplets hitting a cold wall surface using the TIP4P/Ice water model. The accuracy of ReaxFF, SPC/E, TIP4P/2005 and TIP4P/Ice in determining microstructures and mechanical properties has not yet been discussed, however, although these models are widely used in MD simulations for ice or water, as described above. To our knowledge [65,66], little work has been done for determining the best empirical potential or model for ice Ih. In this study, the microstructure properties, surface energy, fracture toughness and fracture cleavage of ice Ih are calculated using ReaxFF, SPC/E, TIP4P/2005, TIP4P/Ice and DFT. The limitations and advantages of ReaxFF and the other three potential models are discussed by comparing the above calculation results, in order to identify a model that is suitable for predicting both the microstructure and mechanical properties of ice on both a large scale (~10−9 m) and long time scale (~10−9–10−6 s) [67].

Fig. 2. Convergence verification regarding the cutoff and k-point. (a) Total energy as a function of the cutoff; (b) total energy as a function of the k-point.

In DFT simulations, the program CASTEP [68–70] and its implementation of the PW91 [71] generalized gradient approximation (GGA) function was used to perform geometric optimization and calculate the surface energies. In addition, the convergence of the relative energy to the cutoff, k-point, model size, and thickness of the cleave surface to the result, as well as the total energy or surface energy, were computed for ice Ih. It is worth noting that other untested guidelines were set to the most stringent and precise values in convergence tests. Although the energy might be expected to be different for each functional method, we observed that different guidelines provided very similar estimates of energy. As shown in Fig. 2(a), for the unit cells examined, the total energy had a nearly constant value with an increasing cut-off value in the range from 300 eV to 900 eV. The total energy approached a constant well before the cut-off value of 800 eV; thus, it was chosen to save computing time. In Fig. 2(b), the energy was constant when the k-point was set to 2 × 2 × 2. Using the cut-off of 800 eV and a k-point of 2 × 2 × 2, the number of bilayers and vacuum layer thicknesses were investigated by convergence tests. The surface energy results were 12.47 meV/Å2 and 12.28 meV/Å2, respectively, when the number of bilayers [31,32] were 6 and 12 (in the Y direction). Slabs that are six bilayers thick were expected to obtain an accurate result, because they had a surface energy difference of 0.19 meV/Å2. Up to 15bilayer-thick slabs were used in other convergence tests, and the surface energy was insensitive to the number of bilayers used; the results were stable at 2–3 bilayer thick slabs [31]. Therefore, 12 bilayers were selected in this paper to ensure reliability. In the vacuum layer thickness test, the thickness was set to 10 Å and 20 Å. The surface energy had a small variation range (12.29 meV/Å2 and 12.28 meV/Å2). The vacuum layer was selected as 20 Å in DFT, considering the successful application of a 20 Å vacuum layer thickness in a previous study [32]. The PBE

2. Computational methods In this work, the initial structure of the ice microstructure is a lowdensity structure, in which two hydrogen atoms are surrounded by oxygen atoms and the hydrogen is regularly arranged. An oxygen atom is adjacent to four surrounding oxygen atoms, and each oxygen atom is located at the center of a tetrahedron surrounded by four other oxygen atoms. Another important feature is that the oxygen atoms are concentrated on the basal face, and the direction perpendicular to the basal face is known as the c-axis. In order to employ the periodic boundary conditions (PBC) in the three directions of the DFT and MD simulations, the orthogonal models shown in Fig. 1 are used. The Z axis is the c-axis direction [0001] in the crystal system, and the X axis and Y axis are the outer normal direction of the secondary prism face {112¯0} and the primary prism face {1¯100} , respectively. Six representative surfaces (Fig. 1) were selected to determine the surface (cleavage plane) prone to cracking, according to the principle of equivalent crystal planes and hydrogen bond distribution on each surface. Specifically, there are two equivalent crystal planes with different topological structures along the direction [112¯0], one due to the so-called intermolecular interaction and one due to hydrogen bonding. The secondary prism face {112¯0} , with a hydrogen bond density of 0.1277 Å−2, is named the X1 surface, and the face with 0.0639 Å−2 is named the X2 surface. These values are also shown in Table 2. The two equivalent crystal planes along the direction [0001] are named the Z1 and Z2 surfaces by the same principle. Along the direction [1¯100], although the surfaces have the same hydrogen bond density of 0.0738 Å−2, they have different topological structures. Thus, the primary prism faces {1¯100} are identified as the Y1 and Y2 surfaces.

Fig. 1. Six representative surfaces that may break. (a) Surfaces X1 and X2; (b) Surfaces Y1 and Y2; (c) Surfaces Z1 and Z2. The solid black line in the three panels corresponds to the secondary prism face, primary prism face, and basal face, respectively. 2

Computational Materials Science 177 (2020) 109546

C. Wang, et al.

DFT molecular angle θHintra - O - H = 107.130 was close to the DFT results of Kuo et al. (106.34 ± 0.36° [29]) and the experimental results of Kuhs et al. (107 ± 1° [15]). The DFT molecular bond and angle results in this work are approximately equal to the existing results of other scholars. The difference was less than 0.28% between the angle formed inter by the nearest three oxygen atoms along the [0001] θO[0001] = 109.633° inter and perpendicular to [0001] θO{0001} = 109.308°, with other DFT results inter inter (θO[0001] =109.4° and θO[0001] = 109.3°) [30] and experimental values inter inter (θO[0001]=109.33° and θO[0001] = 109.61°) [14]. In terms of the distance between the two nearest oxygen atoms, our DFT result, lOinter - O = 2.680 Å, was similar to the DFT result of 2.751 ± 0.007 Å [29] by Kuo et al. and the experimental results of 2.750 Å [14] by Kuhs et al. These slight differences may be induced by differences in the temperature conditions and initial structure (idealized model in DFT and experimental model in the experiment). After a comprehensive comparison between the experiment, our results and other DFT outcomes, we found that the distance errors between the two nearest oxygen atoms lOinter - O was less than 0.071 Å, and the errors of the other four microstructure parameters were not more than 1.10%. Ice is a brittle material [11,72], and surface energy is a one of its most important quantities [32]. Using Eq. (1), the surface energy per unit area γ is calculated and the results are shown in Table 2:

Fig. 2. (continued)

functional scheme provided a similar estimated energy when compared with the PW91 functional (see Section 3.1 for details). Apart from the formula employed above, the quality of convergence tolerance was set as “fine” (a quality option in CASTEP); that is, the optimization convergence criterion was that the energy difference and displacements between cycles were less than 10−5 eV/atom and 0.001 Å, and the force and stress were less than 0.03 eV/Å and 0.05 GPa. In addition, the order of magnitude of energy and self-consistent field (SCF) tolerance were set as 10−6 eV/atoms. In order to better describe dispersion relations and the hydrogen bonding between water molecules, and to obtain more accurate results, the Grimme method for DFT-D correction was employed. MD simulations were performed using the Large-scale AtomicMolecular Massively Parallel Simulator (LAMMPS) software package [57]. ReaxFF, SPC/E, TIP4P/2005, and TIP4P/Ice were selected for comparisons with each other and the DFT results by first performing geometric optimization and calculating the surface energies. The MD model used the orthogonal structure, while its three-sided vector and the researched surfaces were the same as the DFT simulation. The surface was constructed with a large vacuum layer thickness of 40 Å. In order to relax the ice structures and obtain stable structures, MD simulations were conducted in the NPT ensemble with 2 fs timesteps at 1.0 K and atmospheric pressure. Three-dimensional periodic boundary conditions were applied in all directions of the simulation box. A NoseHoover thermostat and barostat were employed to control the temperature and pressure with damping times of 1 ps. The relaxation was performed with a simulation time of 2 ns to ensure that the convergence of the energy and volume was achieved in order to obtain a stable system.

γ=

Esurface − E0 2A

(1)

where E0 is the total energy of the supercell, Esurface is the total energy of the supercell with new surfaces, and A is the surface area of one side, which depends on the considered face. The results in Table 2 indicate that the X1 surface energy, 15.98 meV/Å2, was approximately 3.85 meV/Å2 larger than the X2 surface energy, 12.13 meV/Å2, because of the higher hydrogen bond density on the X1 surface than the X2 surface. Therefore, if there is a crack in the secondary prism face, it would more easily form on the X2 surface. The two surfaces of the primary prism face (Y face) were completely equal in hydrogen bond density, and the Y1 surface energy, 12.32 meV/Å2, was in excellent agreement with the Y2 surface energy, 12.20 meV/Å2. The main axis directions were the close packed direction inducing the largest surface energy of the Z surfaces in all three directions. The surface energy of the surface Z1 is larger than Z2 by 11.48 meV/Å2, which was induced by the higher hydrogen bond density on Z1 than Z2. Our X and Y surface energy calculations using DFT agreed with other calculations (4.3–18.6 meV/Å2) by empirical and semi-empirical methods. Pan et al. and Sun et al. performed DFT calculations that produced similar values (12.0–18.2 meV/Å2). The X2, Y1 and Y2 surfaces are the most likely cleavage surfaces, where the crack expands in the six selected representative surfaces because the energy is the smallest. The accuracy of the DFT calculation results are related to the selected functional form. Therefore, in order to study the effect of this functional on the surface energy, three surfaces were randomly selected to calculate the surface energy using a PBE function. The surface energy values given in Table 2 show that the two predicted results with PW91 and PBE functions were relatively close. The differences were 0.24 meV/Å2, 1.01 meV/Å2 and 0.86 meV/Å2, respectively, which suggest that the selection of the functional has little effect on the accuracy of the results. Thus, the results using function PW91 will be used as a reference.

3. Results and discussion 3.1. DFT studies Ice Ih exhibited quasi-regular tetrahedral microstructural features in the DFT simulation, as shown in Fig. 3. The five oxygen atoms were located at the center of the tetrahedron and the four vertices. There were two hydrogen atoms in the vicinity of each oxygen atom and the hydrogen atoms were located on an approximate line between the two nearest neighboring oxygen atoms. Each water molecule formed a hydrogen bond with the nearest four water molecules. The average molecular geometry of a water molecule is summarized in Table 1. The calculated H–O bond length in the water moleculeslHintra - O = 1.010 Å agreed well with the experimental value of 1.005 Å [14] and other DFT calculations 0.9990 ± 0.0008 Å [29]. Our

3.2. MD studies The stable structure and surface energy of ice Ih were studied by MD, and the microstructure results are summarized in Table 3. The three water models, SPC/E, TIP4P/2005 and TIP4P/Ice, shook the H–O bond and the H-O-H angle in the same water molecule [38] and considered the interaction between molecules [73]. The H–O bond length in the same molecule was 1.003 Å by ReaxFF, whose difference with 3

Computational Materials Science 177 (2020) 109546

C. Wang, et al.

Fig. 3. Quasi-tetrahedral microstructure of ice Ih. The hydrogen bond is shown with a dashed line, red for oxygen, white for hydrogen. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

DFT (1.010 Å) and the experimental value (1.005 Å) [14] was less than 0.01 Å. However, the error in the H–O–H bond angle by ReaxFF was less than 1.42°, as compared with 107.13° for DFT and 107 ± 1° [15] for the experimental value. Regarding structural properties between inter molecules, θO[0001] , the angle formed by the three nearest oxygens along [0001], a TIP4P/2005 result of 108.426° was 1.1% less than our DFT result of 109.633° and the experimental value of 109.33° [14]. The water model TIP4P/Ice value of 107.469° and the SPC/E result of 107.658° were similar; differences with DFT and the experimental vainter lues were not more than 1.98%, and the θO[0001] for ReaxFF had the inter inter worst accuracy. Interestingly, θO{0001} had the same trend with θO[0001] , inter where TIP4P/2005 θO{0001} 110.435° was the closest to DFT 109.308° and the experimental value 109.61° [14]. 111.195° for TIP4P/Ice and 110.999° for the SPC/E water model were acceptable, with differences less than 1.89°. ReaxFF did not have good simulation accuracy. ReaxFF and TIP4P/Ice in MD mirrored the experimental value of the O–O bond length lOinter − O accurately. Surface energies calculated by DFT and MD are shown in Fig. 4. The surface energies calculated by MD approached the DFT values in the following order (with decreasing error): ReaxFF, SPC/E, TIP4P/2005, and TIP4P/Ice. We mainly discuss the TIP4P/2005 and TIP4P/Ice water models because of the large errors (up to 27.21%) in ReaxFF and SPC/E. For the secondary prism face X2, a surface energy of 11.58 meV/Å2 was obtained by TIP4P/Ice, which was close to the DFT value of 12.13 meV/ Å2. Therefore, the TIP4P/Ice model can calculate the secondary prism face X2 accurately. However, the difference of TIP4P/2005, which was closely followed by TIP4P/Ice, was up to 15.25% when compared with DFT. The primary prism faces, Y1 and Y2, revealed the same trend as the secondary prism face X2. The TIP4P/Ice result of 12.63 meV/Å2 on

Table 2 Surface energies of ice Ih obtained by DFT. Indices of the crystal face Surface name

{112¯0} X1

X2

{1¯100} Y1

Y2

{0001} Z1

Z2

γ (PW91) (meV/Å ) γ (PBE) (meV/Å2) Hydrogen bond density (Å−2)

15.98 − 0.1277

12.13 12.37 0.0639

12.32 13.33 0.0738

12.20 − 0.0738

51.96 − 0.1814

40.48 41.34 0.0605

2

Table 3 Comparison of the averaged molecular geometry of ice Ih by MD.

DFT Experimental TIP4P/Ice SPC/E TIP4P/2005 ReaxFF

Intramolecular

Intermolecular

lHintra - O (Å)

θHintra - O - H (°)

inter (°) θO[0001]

inter θO{0001} (°)

lOinter - O (Å)

1.010 1.005 − − − 1.003

107.130 107 ± 1 − − − 105.713

109.633 109.33 107.469 107.658 108.426 106.922

109.308 109.61 111.195 110.999 110.435 111.433

2.680 2.750 2.752 2.714 2.733 2.748

Y1 was approximately equal to the DFT value of 12.32 meV/Å2. Although the TIP4P/2005 results of 11.46 meV/Å2 were satisfactory, the difference was much larger than that of TIP4P/Ice. For the Y2 surface, the TIP4P/Ice and TIP4P/2005 calculations were similar to the DFT value, and their differences were 0.55 meV/Å2 and 0.50 meV/Å2, respectively. The results of ReaxFF and SPC/E differed from the DFT

Table 1 intra Summary of the geometric parameters for simulated ice Ih by DFT. lHintra - O and θH - O - H are the H-O bond length and the H-O-H angle in the water molecule, respectively; inter inter and θO[0001] and θO {0001} are the angles parallel and perpendicular to [0 0 0 1] formed by the three nearest oxygen atoms, where lOinter - O represents the length of the two nearest oxygen atoms. Intramolecular Parameters

Our DFT Other DFT Experimental

Intermolecular Parameters

lHintra - O (Å)

θHintra - O - H (°)

inter (°) θO[0001]

inter θO{0001} (°)

lOinter - O (Å)

1.010 0.9990 ± 0.0008 1.005

107.130 106.34 ± 0.36 107 ± 1

109.633 109.4 109.33

109.308 109.3 109.61

2.680 2.751 ± 0.007 2.750

4

Computational Materials Science 177 (2020) 109546

C. Wang, et al.

Fig. 4. Comparison of the surface energies obtained by MD and DFT.

result significantly (difference was above 11.9%). In addition, the surface energies of the different surfaces (X2, Y1 and Y2) calculated by the same water model had the same trend. This trend shows that the surface energy of the secondary prism face X2 was the lowest of all the surfaces. For the surface energies of the primary prism faces, Y1 and Y2, the calculations were close, which was mainly due to equal hydrogen bond density on Y1 and Y2.

Fig. 5. Fracture toughness KIc of each surface.

the different methods between the experiment and DFT. As the only water model that can adequately reproduce the primary and secondary prism surfaces at the same time, the TIP4P/Ice calculation in fracture toughness agreed well with DFT, with a difference less than 1.52 kPa·m1/2. TIP4P/2005 had poorer accuracy than TIP4P/Ice; the differences were about 7.94%, 3.55%, and 2.07%, respectively. The applicability of the SPC/E and ReaxFF water models for ice Ih is not optimistic, as the errors were the largest of the four models. Finally, the fracture toughness of the secondary prism surface X2 {112¯0} was the lowest. It had the worst ability to resist the extension of cracks and was the first cleavage plane for crack propagation in ice Ih. The primary prism surfaces, Y1 and Y2 {1¯100} , were the second cleavage planes of crack propagation in ice Ih, considering their difference in fracture toughness was only 0.49%, and that they had the same hydrogen bond density.

3.3. Comparison of the fracture toughness Based on the above research, the fracture toughness of X2, Y1 and Y2 and where they might fracture were studied. According to linear elastic fracture mechanics, the fracture toughness KIc of ice Ih can be calculated as [57]:

KIc =

2Eγ (1 − ν 2)

(2)

where E is the elastic modulus, ν is Poisson’s ratio, and γ is the surface energy per unit area. Langleben and Pounder [18] found that the elastic modulus of ice was 9–10 GPa, in the case of a low brine volume. They fit the relationship between the volume of the brine formed (vb ) and the elastic modulus as follows: E = 10 − 0.0351vb (in GPa). Therefore, the elastic modulus of brine free ice Ih used in this paper was 10 GPa, which is consistent with the experimental estimate [19] of freshwater ice at a high strain rate. Weeks and Assur [17] derived the relationship between Poisson's ratio for sea ice and temperature as ν = 0.333 + 0.06105 exp(Ti /5.48) by analyzing Lin'kov's experimental results, where Ti is the temperature of ice (°C). We confirmed that Poisson's ratio in this paper was 0.333 by the above relationship, as Poisson's ratio tended to be 0.333 when the temperature was lower than −30 °C. The fracture toughness of each surface using DFT and MD is shown in Fig. 5. The fracture toughness of the three surfaces by DFT were similar, wherein the fracture toughness of the secondary prism face X2 was the lowest (66.115 kPa·m1/2); the fracture toughness of Y2 was 66.305 kPa·m1/2. The fracture toughness of the primary prism surface Y2 was the highest at 66.631 kPa·m1/2. Petrovic [20] found that the fracture toughness of ice is a function of temperature, grain size and loading rate, and the fracture toughness is generally between 50 and 150 kPa·m1/2. Urabe et al. [21] found that the fracture toughness of ice was 64.724 kPa·m1/2 using a three-point bending test when the strain rate was 2.91 × 10−3 s−1, which is approximately 2.15% smaller than the minimum fracture toughness (X2, 66.115 kPa·m1/2) obtained using the DFT method. Our calculation results were satisfactory considering

4. Conclusion In summary, this article compared the ReaxFF, SPC/E, TIP4P/2005 and TIP4P/Ice models with DFT and discussed their limitations and advantages in both the microstructures and mechanical properties of large size hexagonal Ice (Ih) over a long time scale. Similar DFT and MD results (microstructure, surface energy) were obtained in this work as compared to experiments conducted by other researchers, indicating that the results are credible. The microstructure of ice Ih was predicted using MD and DFT, and we compared them with other DFT and experimental results. It was found that the SPC/E and TIP4P/2005 water models could predict the inter inter angle formed by the three nearest oxygen atoms, θO[0001] and θO(0001) , whose errors were less than 2% but were poor at predicting the O–O inter bond length lOinter - O . Except for lO - O , the errors in the ReaxFF force field were largest. The TIP4P/Ice potential model had a comprehensively better performance when modeling the microstructures of ice Ih, which suggests that TIP4P/Ice can best characterize the actual properties of ice Ih. In the study of surface energy, the secondary prism face (X2), with the lowest surface energy of 12.13 meV/Å2 by DFT and 11.58 meV/Å2 by TIP4P/Ice, was the most probable cleavage plane for crack growth in ice Ih. Because of equal hydrogen bond density (0.0738 Å−2), the primary prism faces (Y1 and Y2) had similar surface energy, 12.32 meV/Å2 and 12.20 meV/Å2, respectively. It was found that hydrogen bond density determines the surface energy on the same 5

Computational Materials Science 177 (2020) 109546

C. Wang, et al.

surface. The basal face Z had a high surface energy (up to 51.96 meV/ Å2) and hydrogen bond density (up to 0.1814 Å−2), which was induced by the close packed direction of the main axis. The fracture toughness of X2, Y1 and Y2 were studied using linear elastic fracture mechanics. The TIP4P/Ice results were satisfactory, and its difference was less than 1.52 kPa·m1/2 with DFT. TIP4P/2005 had poorer accuracy than TIP4P/Ice. The applicability of the SPC/E and ReaxFF water models were not optimistic, as the errors were the largest, up to 14.68%, for the four models. In general, it was found that TIP4P/ Ice was the most suitable model for determining both the microstructures and mechanical properties of ice Ih, concerning the performance in surface energy, fracture toughness and microstructures. The water models TIP4P/2005 and SPC/E are only acceptable for predicting the microstructure of ice Ih. Although different water models and force fields predicted different fracture toughness, the same trend shows that the secondary prism face, with the lowest hydrogen bond density (0.0639 Å−2) and fracture toughness (66.115 kPa·m1/2), was the first cleavage plane of ice Ih, and the primary prism face was the second cleavage plane, whose fracture toughness was 66.305 kPa·m1/2 and hydrogen bond density was 0.0738 Å−2.

2010.04.004. [9] Q. Yin, L. Hu, X. Wu, K. Xiao, C. Huang, Temperature-dependent phase transformation of ice-1h under ultrafast uniaxial compression: a molecular dynamics simulation, Comput. Mater. Sci. 162 (2019) 340–348, https://doi.org/10.1016/j. commatsci.2019.03.013. [10] S. Xiao, J. He, Z. Zhang, Nanoscale deicing by molecular dynamics simulation, Nanoscale 8 (2016) 14625–14632, https://doi.org/10.1039/c6nr02398c. [11] S.C. Di, Y.Z. Xue, X.L. Bai, Q. Wang, Effects of model size and particle size on the response of sea-ice samples created with a hexagonal-close-packing pattern in discrete-element method simulations, Particuology 36 (2018) 106–113, https://doi. org/10.1016/j.partic.2017.04.004. [12] Q. Wang, Y. Wang, Y. Zan, W. Lu, X. Bai, J. Guo, Peridynamics simulation of the fragmentation of ice cover by blast loads of an underwater explosion, J. Mar. Sci. Technol. 23 (2018) 52–66, https://doi.org/10.1007/s00773-017-0454-x. [13] M.H. Liu, Q. Wang, W. Lu, Peridynamic simulation of brittle-ice crushed by a vertical structure, Int. J. Nav. Archit. Ocean Eng. 9 (2017) 209–218, https://doi. org/10.1016/j.ijnaoe.2016.10.003. [14] W. Kuhs, M. Lehmann, The structure of the ice Ih by neutron diffraction, J. Phys. Chem. 87 (1983) 4312–4313, https://doi.org/10.1021/j100244a063. [15] W. Kuhs, M. Lehmann, Water Science Reviews, Cambridge Univ. Press, Cambridge, 1986. [16] R. Gagnon, H. Kiefte, M. Clouter, E. Whalley, Acoustic velocities and densities of polycrystalline ice Ih, II, III, V, and VI by Brillouin spectroscopy, J. Chem. Phys. 92 (1990) 1909–1914, https://doi.org/10.1063/1.458021. [17] G.W. Timco, W.F. Weeks, A review of the engineering properties of sea ice, Cold Reg. Sci. Technol. 60 (2010) 107–129, https://doi.org/10.1016/j.coldregions.2009. 10.003. [18] M.P. Langleben, E.R. Pounder, Elastic Parameters of Sea Ice, MIT, USA, 1963. [19] L.W. Gold, Engineering properties of fresh-water ice, J. Glaciol. 19 (1977) 197–212, https://doi.org/10.3189/S0022143000215608. [20] J.J. Petrovic, Mechanical properties of ice and snow, J. Mater. Sci. 38 (2003) 1–6, https://doi.org/10.1023/a:1021134128038. [21] N. Urabe, T. Iwasaki, A. Yoshitake, Fracture toughness of sea ice, Cold Reg. Sci. Technol. 3 (1980) 29–37, https://doi.org/10.1016/0165-232X(80)90004-X. [22] A. Griffith, The phenomena of flow and rupture in solids, Philos. Trans. R. Soc. Lond. Ser. A 221 (1920) 163–198. [23] J.A. Barker, R.O. Watts, Structure of water; a Monte Carlo calculation, Chem. Phys. Lett. 3 (1969) 144–145, https://doi.org/10.1016/0009-2614(69)80119-3. [24] K. Laasonen, M. Sprik, M. Parrinello, R. Car, ‘‘Ab initio’’ liquid water, J. Chem. Phys. 99 (1993) 9080–9089, https://doi.org/10.1063/1.465574. [25] J. Chen, G. Schusteritsch, C.J. Pickard, C.G. Salzmann, A. Michaelides, Two dimensional ice from first principles: structures and phase transitions, Phys. Rev. Lett. 116 (2016) 025501, , https://doi.org/10.1103/PhysRevLett. 116.025501. [26] B. Santra, J. Klimes, D. Alfe, A. Tkatchenko, B. Slater, A. Michaelides, R. Car, M. Scheffler, Hydrogen bonds and van der Waals forces in ice at ambient and high pressures, Phys. Rev. Lett. 107 (2011) 185701, , https://doi.org/10.1103/ PhysRevLett. 107.185701. [27] E. Schwegler, M. Sharma, F. Gygi, G. Galli, Melting of ice under pressure, Proc. Natl. Acad. Sci. U.S.A. 105 (2008) 14779–14783, https://doi.org/10.1073/pnas. 0808137105. [28] Y.A. Mantz, F.M. Geiger, L.T. Molina, B.L. Trout, First-principles theoretical study of molecular HCl adsorption on a hexagonal ice (0001) surface, J. Phys. Chem. A 105 (2001) 7037–7046, https://doi.org/10.1021/jp010817u. [29] J.L. Kuo, M.L. Klein, W.F. Kuhs, The effect of proton disorder on the structure of iceIh: a theoretical study, J. Chem. Phys. 123 (2005) 134505, , https://doi.org/10. 1063/1.2036971. [30] T.K. Hirsch, L. Ojamae, Quantum-chemical and force-field investigations of ice Ih: computation of proton-ordered structures and prediction of their lattice energies, J. Phys. Chem. B 108 (2004) 15856–15864, https://doi.org/10.1021/jp048434u. [31] D. Pan, L.M. Liu, G.A. Tribello, B. Slater, A. Michaelides, E. Wang, Surface energy and surface proton order of the ice Ih basal and prism surfaces, J. Phys.: Condens. Matter 22 074209 (2010), https://doi.org/10.1088/0953-8984/22/7/074209. [32] D. Pan, L.M. Liu, G.A. Tribello, B. Slater, A. Michaelides, E. Wang, Surface energy and surface proton order of ice Ih, Phys. Rev. Lett. 101 (2008) 155703, , https://doi. org/10.1103/PhysRevLett. 101.155703. [33] Z.R. Sun, P. Ding, X.U. Limei, E.G.J.S.S. Wang, Progress in the structural and physical properties of ice surface, Sci. Sin. 43 (2013) 1144–1150, https://doi.org/ 10.1360/132013-336. [34] A.V.S. de Reuck, The surface free energy of ice, Nature 179 (1957) 1119, https:// doi.org/10.1038/1791119a0. [35] C. Van Oss, R. Giese, R. Wentzek, J. Norris, E. Chuvilin, Surface tension parameters of ice obtained from contact angle data and from positive and negative particle adhesion to advancing freezing fronts, J. Adhes. Sci. Technol. 6 (1992) 503–516, https://doi.org/10.1163/156856192X00827. [36] J. Douillard, M. Henry, Calculation of surface enthalpy of solids from an ab initio electronegativity based model: case of ice, J. Colloid Interface Sci. 263 (2003) 554–561, https://doi.org/10.1016/S0021-9797(03)00093-6. [37] M. Henry, First-principles derivation of vacuum surface energies from crystal structures, Solid State Sci. 5 (2003) 1201–1205, https://doi.org/10.1016/S12932558(03)00183-3. [38] S. Meng, E.G. Wang, Water Science Theory and Experiment, Peking University Press, Beijing, 2014. [39] M.A. Carignano, Formation of stacking faults during ice growth on hexagonal and cubic substrates, J. Phys. Chem. C 111 (2007) 501–504, https://doi.org/10.1021/ jp067388q. [40] J. Gelman Constantin, M.A. Carignano, H.R. Corti, I. Szleifer, Molecular dynamics

CRediT authorship contribution statement Chunyang Wang: Conceptualization, Methodology, Software. Yanzhuo Xue: Validation, Visualization. Chaoying Wang: Supervision, Writing - review & editing. Duanfeng Han: Data curation, Writing - original draft. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments Funding: This work was supported by the Key Program of the National Natural Science Foundation of China [grant numbers 51639004], the National Key Research and Development Program of China [grant numbers 2018YFC1406000, 2017YFE0111400, 2016YFE0202700], the National Natural Science Foundation of China [grant numbers 51579054], and the Fundamental Research Funds for the Central Universities of Ministry of Education of China. References [1] M. Millot, F. Coppari, J.R. Rygg, A. Correa Barrios, S. Hamel, D.C. Swift, J.H. Eggert, Nanosecond X-ray diffraction of shock-compressed superionic water ice, Nature 569 (2019) 251–255, https://doi.org/10.1038/s41586-019-1114-6. [2] L. del Rosso, M. Celli, L. Ulivi, New porous water ice metastable at atmospheric pressure obtained by emptying a hydrogen-filled ice, Nat. Commun. 7 (2016) 13394, https://doi.org/10.1038/ncomms13394. [3] T. Matsui, M. Hirata, T. Yagasaki, M. Matsumoto, H. Tanaka, Communication: hypothetical ultralow-density ice polymorphs, J. Chem. Phys. 147 (2017) 091101, , https://doi.org/10.1063/1.4994757. [4] Y.Y. Huang, C.Q. Zhu, L. Wang, J.J. Zhao, X.C. Zeng, Prediction of a new ice clathrate with record low density: a potential candidate as ice XIX in guest-free form, Chem. Phys. Lett. 671 (2017) 186–191, https://doi.org/10.1016/j.cplett. 2017.01.035. [5] Y.Y. Huang, C.Q. Zhu, L. Wang, X.X. Cao, Y. Su, X. Jiang, S. Meng, J.J. Zhao, X.C. Zeng, A new phase diagram of water under negative pressure: the rise of the lowest-density clathrate s-III, Sci. Adv. 2 (2016) e1501010, , https://doi.org/10. 1126/sciadv.1501010. [6] K. Röttger, A. Endriss, J. Ihringer, S. Doyle, W. Kuhs, Lattice constants and thermal expansion of H2O and D2O ice Ih between 10 and 265 K, Acta Crystallogr. B 50 (1994) 644–648, https://doi.org/10.1107/S0108768111046908. [7] L. Zhou, H.J. Xu, S.K. Long, D.W. Li, Research of aircraft icing characteristics and anti-icing and de-icing technology, China Saf. Sci. J. 20 (2010) 105–110, https:// doi.org/10.16265/j.cnki.issn1003-3033.2010.06.018. [8] X. Fan, D. Bing, J. Zhang, Z. Shen, J.-L. Kuo, Predicting the hydrogen bond ordered structures of ice Ih, II, III, VI and ice VII: DFT methods with localized based set, Comput. Mater. Sci. 49 (2010) S170–S175, https://doi.org/10.1016/j.commatsci.

6

Computational Materials Science 177 (2020) 109546

C. Wang, et al.

[41] [42]

[43]

[44]

[45]

[46]

[47]

[48]

[49]

[50]

[51]

[52]

[53]

[54]

[55]

[56]

[57] Z. Qin, M.J. Buehler, Carbon dioxide enhances fragility of ice crystals, J. Phys. D Appl. Phys. 45 (2012), https://doi.org/10.1088/0022-3727/45/44/445302. [58] M.S. Korlie, 3D simulation of cracks and fractures in a molecular solid under stress and compression, Comput. Math. Appl. 54 (2007) 638–650, https://doi.org/10. 1016/j.camwa.2005.09.008. [59] H.J. Berendsen, J.P. Postma, W.F. van Gunsteren, J. Hermans, Interaction Models for Water in Relation to Protein Hydration, in: B. Pullman (Ed.), Intermolecular Forces, Springer, Dordrecht, 1981, pp. 331–342. [60] J.L.F. Abascal, C. Vega, A general purpose model for the condensed phases of water: TIP4P/2005, J. Chem. Phys. 123 (2005) 234505, , https://doi.org/10.1063/1. 2121687. [61] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys. 79 (1983) 926–935, https://doi.org/10.1063/1.445869. [62] J.L.F. Abascal, E. Sanz, R.G. Fernandez, C. Vega, A potential model for the study of ices and amorphous water: TIP4P/Ice, J. Chem. Phys. 122 (2005) 234511, , https:// doi.org/10.1063/1.1931662. [63] P.A. Franco Pinheiro Moreira, R.G. de Aguiar Veiga, I.D.A. Ribeiro, R. Freitas, J. Helfferich, M. de Koning, Anomalous diffusion of water molecules at grain boundaries in ice I-h, PCCP 20 (2018) 13944–13951, https://doi.org/10.1039/ c8cp00933c. [64] A.C.T. van Duin, S. Dasgupta, F. Lorant, W.A. Goddard, ReaxFF: a reactive force field for hydrocarbons, J. Phys. Chem. A 105 (2001) 9396–9409, https://doi.org/ 10.1021/jp004368u. [65] Pedro Augusto Franco Pinheiro Moreira, Roberto Gomes de Aguiar Veiga, Maurice de Koning, Elastic constants of ice I h as described by semi-empirical water models, J. Chem. Phys. 150 (4) (2019) 044503, https://doi.org/10.1063/1. 5082743. [66] P.A. Santos-Flórez, C.J. Ruestes, M. de Koning, Uniaxial-deformation behavior of ice I h as described by the TIP4P/Ice and mW water models, J. Chem. Phys. 149 (2018) 164711, , https://doi.org/10.1063/1.5048517. [67] A. Gooneie, S. Schuschnigg, C. Holzer, A review of multiscale computational methods in polymeric materials, Polymers 9 (2017), https://doi.org/10.3390/ polym9010016. [68] V. Milman, B. Winkler, J. White, C. Pickard, M. Payne, E. Akhmatskaya, R. Nobes, Electronic structure, properties, and phase stability of inorganic crystals: a pseudopotential plane-wave study, Int. J. Quantum Chem. 77 (2000) 895–910, https:// doi.org/10.1002/(SICI)1097-461X(2000) 77:5<895::AID-QUA10>3.0.CO;2-C. [69] S.J. Clark, M.D. Segall, C.J. Pickard, P.J. Hasnip, M.J. Probert, K. Refson, M.C. Payne, First principles methods using CASTEP, Z. Kristallogr. 220 (2005) 567–570, https://doi.org/10.1524/zkri.220.5.567.65075. [70] K. Zhang, P. Zhang, Z.R. Wang, X.L. Zhu, Y.B. Lu, C.B. Guan, Y.H. Li, DFT simulations of the vibrational spectrum and hydrogen bonds of ice XIV, Molecules 23 (2018), https://doi.org/10.3390/molecules23071781. [71] D. Liu, J. Deng, Y. Jin, Interactions of water molecule with HfB2 and TaB2 (0001) surfaces: a first-principles investigation, Comput. Mater. Sci. 97 (2015) 116–120, https://doi.org/10.1016/j.commatsci.2014.10.006. [72] Paul Duval, Creep and Fracture of Ice, Cambridge University Press, 2009. [73] B. Guillot, A reappraisal of what we have learnt during three decades of computer simulations on water, J. Mol. Liq. 101 (2002) 219–260, https://doi.org/10.1016/ s0167-7322(02)00094-6.

simulation of ice indentation by model atomic force microscopy tips, J. Phys. Chem. C 119 (2015) 27118–27124, https://doi.org/10.1021/acs.jpcc.5b10230. R.G. Pereyra, M.A. Carignano, Ice nanocolumns: a molecular dynamics study, J. Phys. Chem. C 113 (2009) 12699–12705, https://doi.org/10.1021/jp903404n. Y. Liu, J. Liu, S. Yue, J. Zhao, B. Ouyang, Y. Jing, Atomistic simulations on the tensile deformation behaviors of three-dimensional graphene, Phys. Status Solidi B 255 (2018), https://doi.org/10.1002/pssb.201700680. Q. Liu, S. Shen, Interaction of voids and nano-ductility in single crystal silicon, Comput. Mater. Sci. 67 (2013) 123–132, https://doi.org/10.1016/j.commatsci. 2012.08.039. Y. Jing, N.R. Aluru, Size effect on brittle and ductile fracture of two-dimensional interlinked carbon nanotube network, Phys. B 520 (2017) 82–88, https://doi.org/ 10.1016/j.physb.2017.06.026. Y. Jing, C. Zhang, Y. Liu, L. Guo, Q. Meng, Mechanical properties of kinked silicon nanowires, Phys. B 462 (2015) 59–63, https://doi.org/10.1016/j.physb.2015.01. 018. Q. Liu, L. Wang, S. Shen, Effect of surface roughness on elastic limit of silicon nanowires, Comput. Mater. Sci. 101 (2015) 267–274, https://doi.org/10.1016/j. commatsci.2015.02.009. W. Yu, S. Shen, Q. Liu, Energetics of point defect interacting with bi-crystal Σ3 copper grain boundaries, Comput. Mater. Sci. 118 (2016) 47–55, https://doi.org/ 10.1016/j.commatsci.2016.02.038. G. Quan, Y. Mao, G. Li, W. Lv, Y. Wang, J. Zhou, A characterization for the dynamic recrystallization kinetics of as-extruded 7075 aluminum alloy based on true stress–strain curves, Comput. Mater. Sci. 55 (2012) 65–72, https://doi.org/10. 1016/j.commatsci.2011.11.031. G. Quan, Y. Tong, G. Luo, J. Zhou, A characterization for the flow behavior of 42CrMo steel, Comput. Mater. Sci. 50 (2010) 167–171, https://doi.org/10.1016/j. commatsci.2010.07.021. G. Quan, G. Luo, J. Liang, D. Wu, A. Mao, Q. Liu, Modelling for the dynamic recrystallization evolution of Ti–6Al–4V alloy in two-phase temperature range and a wide strain rate range, Comput. Mater. Sci. 97 (2015) 136–147, https://doi.org/10. 1016/j.commatsci.2014.10.009. N.J. English, Structural properties of liquid water and ice Ih from ab-initio molecular dynamics with a non-local correlation functional, Energies 8 (2015) 9383–9391, https://doi.org/10.3390/en8099383. J.A. Hayward, A.D.J. Haymet, The ice/water interface: molecular dynamics simulations of the basal, prism, {2021}, and {2110} interfaces of ice Ih, J. Chem. Phys. 114 (2001) 3713–3726, https://doi.org/10.1063/1.1333680. T. Ikeda-Fukazawa, S. Horikawa, T. Hondoh, K. Kawamura, Molecular dynamics studies of molecular diffusion in ice Ih, J. Chem. Phys. 117 (2002) 3886–3896, https://doi.org/10.1063/1.1495844. H.B. Hu, Q. He, S.X. Yu, Z.Z. Zhang, D. Song, Freezing behavior of droplet impacting on cold surfaces, Acta Phys. Sin. 65 (2016) 104703, , https://doi.org/10. 7498/aps.65.104703. Q.Q. Dong, H.B. Hu, S.Q. Chen, Q. He, L.Y. Bao, Molecular dynamics simulation of freezing process of water droplets impinging on cold surface, Acta Phys. Sin. 67 (2018) 054702, , https://doi.org/10.7498/aps.67.20172174. L. Vrbka, P. Jungwirth, Brine rejection from freezing salt solutions: a molecular dynamics study, Phys. Rev. Lett. 95 (2005) 148501, , https://doi.org/10.1103/ PhysRevLett. 95.148501.

7