Pressure-Induced Phase Transitions of Proton-Ordered Ices at the Low Temperature Limit
Yosuke KATAOKA and Yuri YAMADA
There are many polymorphs of ice under high pressure, in some of which the protons are ordered at low temperatures . Monte Carlo simulations have been used to obtain a phase diagram of water above 100 K, where the hydrogen atoms are disordered . In this paper, the pressure induced phase transitions in proton-ordered ices were studied, because the ordered structures are the basis of studies of the phase transitions in solids, although some disordered ices are frequently observed in the laboratory.
The TIP4P model proposed by Jorgensen et al. was assumed . This is a rigid model and is used in many molecular dynamics simulations (MD ) of liquid and solid water systems. The internal energy and volume were obtained by NTp molecular dynamics simulations at T = 1 K in our previous work . The enthalpy H was used instead of the Gibbs free energy G to estimate the phase transition pressure of proton-ordered ices at low temperatures . Although qualitative agreements were obtained in that study, ice IV appeared between ice III and VI .
In this paper, the N-dependence of the internal energy at T = 1 K was studied by NTV MD simulations at the density where the internal energy has a minimum. The limit of the infinite number is obtained for the internal energy at T = 1 K. The temperature dependence is calculated from the enthalpy, and the low temperature limit is estimated to obtain the Gibbs free energy at low temperatures.
2 Model and Calculation
2. 1 Potential Function
The pair potential energy in the TIP4P model has the following form :
The quantities e and s are the potential constants, which have dimensions of energy and length, respectively, and their values are given in Table 1. The lowest energy configuration of one water pair is given in Ref. . The negative charge -1.04 e sits on the bisection of the HOH angle in a water molecule, where e is the unit charge of protons. The intermolecular distance is 2.7×10-10 m in the lowest energy configuration. The lowest pair energy is given as:
Table 1. Potential parameters , where e is the charge of a proton.
2. 2 How to Obtain the Proton-Ordered Structure in Ices
A periodic boundary condition is assumed . The Ewald method was adopted for calculation of the coulombic interaction . The time step used was 0.1 fs and the MD length was 100,000 steps. The initial configuration of the small cell is prepared to be consistent with the observed oxygen positions in ice. The orientation has a uniform order, as shown in Figure 1.
Figure 1. Initial configuration of Ice V used in the NTV molecular dynamics simulation at T = 1 K, N = 28.
More stable ordered structures at T = 1 K were obtained by NTV molecular dynamics simulation, where the cut off distance for short range interaction was 1.4×10-9 m . This is longer than the conventional half-length of the cell. The number of molecules in the small cells used to obtain the ordered structure N0 is given in Table 2. The Materials Explorer (v3 and v4) program was used . An example of the relaxation of structure to the ordered structure is shown in Figures 1 - 3 for the case of ice V. Although proton-ordered ice V is sometimes referred to as XIII, the name ice V is used in this paper even for ordered ice V. The relaxed structure, which has the proton-ordered structure, is shown in Figure 2. This structure is compatible with that reported [8, 9]. The relaxation process is shown in Figure 3. The average potential energy per molecule is -8.722×10-20 J. This value indicates that there are 4 hydrogen bonds per molecule in this structure (see Eq. 2).
Table 2. Number of molecules in the cell of each ice type.
N0 is the number of molecules in the cell used to determine the ordered structure. N1 is the number of molecules in the cell used in the previous work . N2 is the number of molecules in the cell used in this study.
Figure 2. Relaxed configuration of Ice V obtained by the NTV molecular dynamics simulation at T = 1 K, N = 28.
Figure 3. Internal energy U, volume V, pressure p and temperature T vs. time for Ice V. An example of the relaxation of the internal energy obtained by the NTV molecular dynamics simulation at T =1 K, N = 28.
2. 3 Proton-Ordered Structure of Ices in This Work
The explicit structures used in this work are shown in Figures 4 - - 11. Figure 11 shows the layered structure found in the present calculation, which corresponds to graphite in the carbon system.
Figure 4. Structure of ice Ih  used in the present work.
Figure 5. Structure of ice Ic  used in the present work.
Figure 6. Structure of ice II  used in the present work.
Figure 7. Structure of ice III  used in the present work.
Figure 8. Structure of ice IV  used in the present work.
Figure 9. Structure of ice VI  used in the present work.
Figure 10. Structure of ice VIII  used in the present work.
Figure 11. Structure of ice in the layer structure found in the present work. This structure corresponds to the graphite in the carbon system.
2. 4 Methods used to Obtain Thermodynamic Properties
NTp MD simulations at T = 1 K were performed to obtain the average volume V and internal energy U at a given pressure p, using a time step of 0.1 fs and a standard MD length of 50,000 steps. The velocity scaling method was used to control temperature. The number of molecules in the cell is given in Table 2. The cut off distance for short range interaction was half of the cell and the Ewald method was used to determine coulombic interaction.
The enthalpy H is calculated in the standard way:
2. 5 Infinite Number Limit
The N-dependence of the internal energy was studied by NTV MD simulations at the density where the internal energy has a minimum. An example is shown for the case of ice Ih in Figure 12. The infinite number limit is estimated by linear fitting of the U-1/N plot, as shown in Figure 12.
Although the linearity is not good around the large number cases, the relative error is improved as follows: the relative error of the U/N at N = 324  is 0.001 and that is 0.0003 at N = 1500 in the present work. The internal energy of ice Ih calculated using an MD cell with number N = 1500 is corrected in this way under other pressures, although this is an assumption. The number of molecules in the basic cell N2 for the other ices is shown in Table 2.
Figure 12. Internal energy per molecule U/N vs. 1/N plot at T = 1 K, ice Ih. The infinite number limit is estimated by linear fitting of the U/N-1/N plot. The obtained equation is shown in the Eq. (4).
2. 6 Low Temperature Limit
NTp simulations at T = 2, 3, 4 and 5 K were also conducted to obtain the temperature dependence of the enthalpy H under pressure p = 0.1 MPa (=1 atm). An example is shown in Figure 13. The enthalpy per molecule H/N is plotted as a function of temperature T. The simulation result was fitted with a linear function of temperature T. The tangent of the linear function is the heat capacity per molecule Cp/N.
Figure 13. Enthalpy per molecule H/N vs. temperature T plot at p = 0.1 MPa, ice Ih. The heat capacity per molecule is estimated by linear fitting of the H/N-T plot. The obtained equation is shown in the Eq. (5).
The heat capacities of various types of ices are compared in Figure 14. This figure shows that the heat capacities of these ices are 6R within the calculation error, where R is the gas constant. From this figure, the low temperature limit of the enthalpy H can be obtained by the linear function of temperature T with the tangent 6R. This limit is the Gibbs free energy G at T = 0 K. This extrapolation of H was applied to T = 0 K under other pressures using the same tangent 6R. Although this is an assumption, the heat capacities of the ordered ices are reasonably expected to be 6R. This value is expected in the classical harmonic oscillator in the translational and the rotational degrees of freedom. The real system is not always the case because of the finite temperatures.
Figure 14. Heat capacity at constant pressure per molecule Cp/N for various ice types at p = 0.1 MPa. As for the meaning of the value 6R, see the text.
3 Results of Calculations
The internal energy per molecule U/N of ice Ih, Ic and others are shown as a function of pressure p at T = 1 K in Figure 15. The curve of ice Ic is very close to that of ice Ih, and ice Ih is stable even under negative pressure. The internal energies per molecule U/N of other ices are compared under higher pressures at T = 1 K in Figure 16, in which the curves of ices VI, VIII, IV and the layer structure are found under very high pressures.
Figure 15. Internal energy per molecule U/N vs. pressure p plot in the low temperature limit under low pressures obtained by the NTp molecular dynamics simulations.
Figure 16. Internal energy per molecule U/N vs. pressure p plot in the low temperature limit for several types of ice, obtained by the NTp molecular dynamics simulations.
The volume per molecule V/N is plotted against pressure p at T = 1 K in Figure 17, where the volume is relatively small in ice VIII and the layer structure. The volume of the layer structure is asymptotic to that of the high pressure ice VIII. One of the important features in Figure 17 is that the volume of ice VI is found in a very wide pressure range, and its value is close to that of ice VIII. In addition, ices II, III, IV and V appear in the intermediate region between the low pressure ice Ih and the high pressure ice VIII.
Figure 17. Volume per molecule V/N vs. pressure p at T = 1 K obtained by the NTp molecular dynamics simulations.
Figure 18 shows the Gibbs free energy per molecule G/N as a function of pressure p in the low pressure region. In this figure, G/N- v0.p is plotted for clarity instead of G/N, because the G/N curves overlap. The typical volume v0 is defined as follows:
This figure shows that ice Ih is most stable around p = 0. Ice Ih changes to ice II at a pressure p of 0.1 GPa. The next phase transition to ice VI then occurs at p = 1 GPa.
The high pressure part of the G/N vs. p plot is shown in Figure 19, which describes the phase transition from ice VI to ice VIII at p = 6 GPa at the low temperature limit. The internal energy of ice VIII is not the lowest in Figure 16. The volume of ice VIII is smallest in Figure 17, which is the reason why ice VIII is most stable under high pressure of p >6 GPa.
Figure 18. Gibbs free energy per molecule G/N vs. pressure p at the low temperature limit under low pressures.
Figure 19. Gibbs free energy per molecule G/N vs. pressure p at the low temperature limit.
Table 3. Phase transition pressures at the low temperature limit estimated by the Figures 18, 19.
p(exp) is the value extrapolated from the experimental results .
Phase transition pressures were obtained and are given Table 3 in comparison with the experimental results. The comparison is at least qualitatively satisfactory. The temperature effect on the phase transition pressures should be obtained by the molecular simulations at higher temperatures, which is a future study. The quantum effects like zero point motion and the tunneling of the proton are also outside of the present classical work.
The authors thank the Research Center for Computing and Multimedia Studies, Hosei University.
[ 1] B. Kamb, Ice Polymorphism and the Structure of Water, Structural Chemistry and Molecular Biology, ed. by A. Rich and N. Davidson, Freeman and Co., San Francisco (1968), pp. 507-542.
There are many reports on ice and liquid water at:
[ 2] E. Sanz, C. Vega, J. L. F. Abascal, and L. G. MacDowell, Phase Diagram of Water from Computer Simulation, Phys. Rev. Letters, 92, 255701-1-4 (2004).
[ 3] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impley, and M. L. Klein, J. Chem. Phys., 79, 926 (1983).
[ 4] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford Univ. Press, Oxford (1987).
[ 5] Y. Kataoka and Y. Yamada, Phase Transitions in Proton Ordered Ice under Pressure by Molecular Dynamics Simulations, Tech. Rep. Res. Cent. Comp. Multi. Med. Stud., 21 (2008).
[ 6] V. Buch, R. Martonak and M. Parrinello, A New Molecular-Dynamics based Approach for Molecular Crystal Structure Search, J. Chem. Phys., 123, 051108-1-4 (2005).
[ 7] http://venus.netlaboratory.com/me4/
[ 8] B. Kamb, A. Prakash and C. Knobler, Structure of ice V, Acta Crystallogr., 22, 706 (1967).
[ 9] C. G. Salzmann, P. G. Radaelli, A. Hallbrucker, E. Mayer and J. L. Finney, The preparation and structure of hydrogen ordered phases of ice, Science, 31, 1758 (2006).