Molecular Dynamics Study on the Density Fluctuation of Supercritical Water

Kazuyoshi UEDA, Taro KOMAI, Isseki YU and Haruo NAKAYAMA


Return

1 Introduction

With the growing concern for the environment, much attention has been focused on the supercritical fluid as an environmentally benign solvent. In order to use a supercritical fluid effectively, understanding of the physico-chemical properties of it is important. Especially, the information about the cluster size, their distribution and the inhomogeneity are important because they relate to the solubility and the reactivity in chemical reaction [1].
Recently, Nishikawa measured the density fluctuation of CO2 and CF3H using the small angle x-ray scattering method and found the ridge of the density fluctuation exists along the extension of the coexistence curve of gas and liquid [2, 3]. They also measured the density fluctuation of supercritical water and obtained its maxima at the lower density point than the critical isochore [4].
On the other hand, many simulational works have been done especially focused on the number of hydrogen bonds of supercritical water and their orientation [5 - 8]. For cluster size and its distribution, Yoshii and Okazaki performed calculations by using large number of Lennard-Jones fluid with NVT ensemble [9]. They also investigated the microscopic properties of water in the sub- and supercritical states with polarizable water model [10]. However, the studies on the cluster size and its distribution are still limited in number compared to the works on the hydrogen bond of water. Especially, the density fluctuation of water at the supercritical state has not been investigated yet.
In this work, we performed a molecular dynamics simulation of water on a wide range of sub- and supercritical states and investigated the behavior of density fluctuation directly. As a first step of this study, we applied a simple "bin method" for the analysis of the density fluctuation. The result showed that the density fluctuation was observed in a large area from the critical point to high pressure and temperature regions and it confirmed the existence of the ridge of the density fluctuation which was found by Nishikawa et al [4].

2 Simulation procedures

All the calculations were performed using the CHARMM25 program [11]. The system consisted of a cubic cell containing 512 of TIP3P water molecules [12]. Cubic periodic boundary conditions were imposed in all directions. The calculations were performed using NPT ensemble using Nose-Hoover's temperature controlling method [13] and Andersen's pressure controlling method [14]. Electrostatic and van der Waals interactions were taped to zero by applying a switching function between 10.0 and 12.0 A [15]. Equation of motion was integrated by leap-flog method with a time step of 1 fs. A 70ps production run was performed after the equilibration at each temperature and pressure state.
Time averaged densities were evaluated at 144 thermodynamic states at the pressure region between 5 and 100 MPa and temperature region between 300 and 1000K, and the isochore curves were calculated from these values. Density fluctuations were analyzed along the critical isochore line and some isobaric lines.
For the analysis of the density fluctuation along the critical isochore, some (P, T) points were selected from the critical isochore curve obtained in the above calculation and other production runs were performed on these thermodynamics states. Similarly, temperature dependence of the density fluctuation along the isobaric lines at the pressures of 25, 30, 40 and 50 MPa were analyzed from the further simulation runs produced at these thermodynamic points.

3 Analytical method

The density fluctuation was evaluated using a method named here as "bin analysis". This method divides the cubic box filled with water into little cubes (bin) with equal volumes, and the number of the water molecules in each bin was calculated at each time step. In this work, each side of the box was divided equally by five and made 125 bins with equal volumes. Although the method is simple, this enables one to obtain the direct information of the density fluctuation of water.
Moreover, if we change the number of division, that is, the size of bin, it may give information about the size of cluster formed among water molecules. However, the size of water box used in this work is not large enough to examine the effect of the size of bin. Therefore, we only analyzed the case of division by 5 as a first step check for this method. The analyses were separately performed, one was along the critical isochore curve with density 0.3 g.cm-3 which was selected based on the critical density value of water 0.32 g.cm-3 and the other was along the isobaric lines at pressures of 25, 30, 40 and 50 MPa.

4 Results and Discussion

Time averaged density of water
Molecular dynamics simulation of water has been performed at various temperature and pressure points and the time averaged densities at these points were obtained. The density values were plotted on the T-P plane as isochore contour curves as shown in Figure 1. Above the value of 0.6 g.cm-3, the density did not show any significant change with increasing pressure under constant temperature. This indicates that the water behaves as liquid in here. On the other hand, it can be seen that the isochore lines under the density values of 0.6 g.cm-3 converge at the temperature range around 550 K to 650 K according to the decrease of the pressure. In these areas, the density changes drastically according to a small change of temperature and pressure. This indicates that the phase transition between the liquid and the gas occurred and the gas-liquid coexistence curve could be reproduced by this simulation.


Figure 1. Isochore lines of water at sub- and supercritical regions obtained by molecular dynamics simulation. The density values are shown in the figure whose dimensions are all g.cm-3. Experimental value of the critical point is shown as small circle.

The convergence of the isochore lines disappeared above the pressure of 20 MPa and the temperature of 650 K. That is, the gas-liquid coexistence curve disappeared at this point. This indicates that the critical point exists at this point in our simulation. Although it is difficult to specify the critical point exactly by the simulation, the experimentally obtained value of the critical point of water ( T = 647 K, P = 22 MPa, r = 0.32 g.cm-3 ) also exists at the same position where the convergence of the density lines disappeared. Though TIP3P is a simple model of the water, it can reproduce well the critical point. Guissani and Guillot obtained the critical point by using another water model of SPCE [5]. They also found that the critical value of their simulation corresponded well with the experimental ones. These two models are both parameterized for the water at ambient condition, but they can reproduce the critical phenomena with good accuracy for our purpose. Therefore, we next analyzed the density fluctuation from sub- to supercritical regions.
Analysis of the density fluctuation
In order to investigate the density fluctuation along the critical isochore curve, some of the state points (T, P) with the average density of 0.3 g.cm-3 were selected from Figure 1. The selected state points were shown in Table 1. The density fluctuations of these points were analyzed using bin analysis. In this analysis, the number of water molecules in each bin was calculated at each time step. Then, the time averaged percentage of the number of bins which include different numbers of water molecules were calculated. These values were plotted as ordinate and the number of water molecules in a bin was plotted as abscissa in Figure 2. In order to clear the figure, the results of the five representative state points out of ten were shown in the figure. As the simulation box was divided equally by 125 in this analysis, the average number of water molecules per a bin should be 4.09. Therefore, the histogram in Figure 2 should show a sharp peak at 4.09 if the density of the system is homogeneous. Figure 2 showed that the peak and the width of the histogram become lower and broader when the thermodynamic state approaches the critical point. Especially, at the state point of T=640K, P=23MPa, which is the closest calculated state to the critical point, the shape of the histogram distorted and showed a flattening of the peak. This indicates that the density fluctuation of the water is large near the critical point. With increasing of the pressure and the temperature, the shape of the histogram became sharp and the peak position approached the value of 4.09. These results qualitatively indicate that the density fluctuation is large near the critical point and gradually decreases with the increment of the temperature and the pressure.

Table 1. Trajectory averaged values of the temperature and pressure calculated at some selected points on the isochore curve, which were used for the bin analysis.
Average TemperatureAverage PressureAverage Density
( K )( MPa )( g.cm-3 )
63119.30.298
64022.80.308
65026.20.295
67031.30.307
70040.20.319
75149.70.301
79966.00.329
84980.60.315
89891.60.309
951100.60.301


Figure 2. Average percentage of the number of bins which include different number of water molecules were plotted against the number of water in a bin. The measured (T, P) points were on the critical isochore whose values are shown in the figure.

In order to investigate the density fluctuation quantitatively, the density fluctuation defined in the next equation [3] was calculated from the histograms in Figure 2.

The calculated values of the density fluctuation were plotted against temperatures in Figure 3. It can be seen that the density fluctuation gradually decreased with the increment of the temperature along the critical isochore line from the critical point. Nishikawa et al measured the density fluctuation of supercritical water near the critical point ( T = 660 ~ 690K and P = 22.5 ~ 38.8MPa ) using the small angle x-ray scattering method [4]. In their measurement, density fluctuation took a maximum value of 31.3 at T = 659.6K, P = 25.22 MPa and r = 0.2914 g.cm-3, and the value rapidly decreased to 1.17 at their most remote state from the critical point at T=687.4K, P=23.69, r=0.1244 g.cm-3. Their value is ten times larger than ours near the critical point, though both took similar values at the peripheral points. These differences should attribute to the number of water and the size of bin used in this study. The occurrence of the large cluster and the density fluctuation near the critical point could not be represented by the size of our 512 water molecules. Nishikawa also obtained the value of 16A for the correlation length near the critical point of water [4]. The estimated size of our bin is 7 ~ 8A in average near the critical point in our simulation. In order to investigate the fluctuation near the critical point, the simulation box with more than 3000 water molecules at least should be considered. Although the quantitative analysis of the density fluctuation near the critical point can not be achieved in this study, the rest of the wide supercritical region can be represented well by this bin analysis. Further studies with a larger scale simulation will give us a more accurate picture around the critical point.


Figure 3. Temperature dependence of the density fluctuations along the critical isochore line.

The temperature dependences of the density fluctuation along the isobaric lines were shown in Figure 4. At the pressure of 25 MPa, density fluctuation began to rapidly increase near the critical temperature, and reached a maximum value at 700K, then gradually decreased again. As isobaric lines were taken at successively higher pressures, the peak of the density fluctuation shifted toward the high temperature gradually. The temperature and the pressure values at the peak point along with each isobaric line were listed in Table 2 with its peak value of the density fluctuation.


Figure 4. Temperature dependence of the density fluctuations along the isobaric curves. The pressures are 25MPa ( ), 30MPa ( ), 40MPa ( ), and 50MPa ( ) , respectively.

Table 2. The peak values of the density fluctuation and their (T, P) value along the isobaric lines which were obtained from Figure 4.
Peak TemperaturePeak PressurePeak Density Fluctuation
( K )( MPa )
70025.01.45
75030.01.48
80040.01.32
85050.01.19

By combining the data along the isochore and isobaric lines, the density fluctuations were represented as a contour map and superimposed on the isochore lines on the P-T phase diagram in Figure 5. The contour line shows that the highest point of the density fluctuation occurs in a little high temperature region from the critical point. However, the highest density fluctuation should be at the critical point. As the simulation of this work does not have a sufficient accuracy for the fluctuation analysis at the critical point, the highest area of the fluctuation should include the critical point and the area may spread from critical point to a little high temperature region.


Figure 5. Contour plot of the density fluctuation ( - ) of supercritical water superimposed on the isochore lines ( ... ). The peak points of the density fluctuation along the isobaric lines shown in Table 2 were plotted as ( × ). Experimental value of the critical point is shown as small circle.

Compared to the isochore lines, the contour lines of density fluctuation below the value of 0.8 behave almost parallel to the isochore lines of the average density above 0.4 g.cm-3. However, the contour lines of density fluctuation above the value of 1.0 are not parallel to the isochore lines, and instead, it forms a maximum near the critical point. The threshold line locates around the isochore line with critical density of 0.3 g.cm-3. From the density fluctuation point of view, this threshold line may be regarded as the boundary of the sub- and supercritical regions.
The maxima points of the density fluctuation along the isobaric lines shown in Table 2 are also plotted in the same figure. The connected line of these peak values of the density fluctuation, that is, the ridge of the density fluctuation, roughly lies along the isochore line at density of 0.2 g.cm-3. This indicates that the ridge of the density fluctuation deviates to a little lower density side from the critical isochore line of 0.3 g.cm-3. This deviation becomes large with the increase of the pressure. Nishikawa et al [4] showed that the ridge of the density fluctuation exists and deviates to a little lower density region than the critical isochore line. The result of our simulation is in accord with their results. However, the ( P, T ) region investigated by their experiment was limited because of the experimental difficulty. They investigated the temperature ranging from 660K to 690K and pressure ranging from 22.5 MPa to 38.8 MPa. Compared to the experimental method, simulation has an advantage in that it can investigate at the extreme conditions. The results obtained in this simulation could predict that the ridge of the density fluctuation spreads to a wide area of the supercritical field. Futhermore, the ridge line deviates from the critical isochore line and lies along the isochore line at 0.2 g.cm-3.
In order to investigate the water structure around the ridge line of the density fluctuation, the coordination number of water in the first hydration shell was calculated. The coordination number, <N>, was evaluated by integrating the oxygen-oxygen radial distribution function until it takes the first minimum at 4.0 A. The results were shown in Figure 6(a). It can be seen that the coordination number decreases with increasing temperature along each isobaric line. Along the isobaric line at 25MPa, the coordination number rapidly decreases at around the critical temperature and reaches a constant value at high temperature. With the elevation of the pressure, the coordination number decreases monotonically. As the bulk density changes according to the elevation of the temperature along the isobaric line, the coordination number should be affected by the bulk density. Therefore, in order to analyze the coordination behavior in more detail, the ratio of <N>/r was plotted in Figure 6(b). It can be seen that all the values of ratio <N>/r take almost constant value below the critical temperature. This means that the coordination number changes proportionally to the bulk density. However, the ratio <N>/r increases above the critical temperature in all cases. Especially, the increment is largest at 25 MPa among the isobaric lines. This could show that the preferential aggregation of water molecules occurred above the critical temperature. This behavior is in accord with that reported by Yoshii et al [10]. Although the curve greatly waves at 25MPa due to the large fluctuation, the peak exists in all cases of the isobaric lines and their values shift toward the higher temperature with the increment of the pressure. This behavior is similar to that of density fluctuation, which was shown in Figure 4. These results may suggest that the preferential aggregation of water becomes maximum along the ridge line of the density fluctuation.


Figure 6. Temperature dependence of the coordination number of water, <N>, in the first hydration shell within the radius of 4.0 A along the isobaric curves (a). The ratio of the coordination number of water and the bulk density, <N>/r, was also plotted in (b). The pressures are 25MPa ( ), 30MPa ( ), 40MPa ( ), and 50MPa ( ) , respectively.

5 Summary

Molecular dynamics simulation of water in the supercritical region was performed to investigate the density fluctuation. The density fluctuation along the critical isochore line and some isobaric lines was analyzed using bin analysis, and the results were plotted on the T-P phase diagram. The maxima of the density fluctuation exists from the critical point to a little high temperature region. The ridge of the density fluctuation shifts to the lower density from the critical isochore line according to the increment of the temperature and pressure. As the number of water molecules was small in this study, the size of bin for the analysis of density fluctuation was also limited. In order to evaluate the behavior near the critical state, simulation with a large number of molecules, where the size of bin is larger than that of cluster, will be needed. Nevertheless, the behavior of the density fluctuation for the wide area of sub- and supercritical water except near the critical point can be obtained in this work. Study with different size of bin would also give information on the size of cluster of water. Although this is a first step of the bin analysis, the result shows that this is an effective way to analyze the density fluctuation of supercritical water.

References

[ 1] Okitsugu Kajimoto, Chem. Rev., 99, 355 (1999).
[ 2] Keiko Nishikawa, Ibuki Tanaka, Chem. Phys. Lett., 29, 149 (1995).
[ 3] Keiko Nishikawa, Takeshi Morita, J. Supercrit. Fluids, 13, 143 (1998).
[ 4] Takeshi Morita, Kohei Kusano, Hiroto Ochiai, Ken-ichi Saitow, Keiko Nishikawa, J. Chem. Phys., 112, 4203 (2000).
[ 5] Yves Guissani, Bertrand Guillot, J. Chem. Phys., 98, 8221 (1993).
[ 6] Ariel A. Chialvo, Peter T. Cummings, J. Chem. Phys., 101, 4466 (1994).
[ 7] Tahmid I. Mizan, Phillip E. Savage, Robert M. Ziff, J. Phys. Chem., 100, 403 (1996).
[ 8] Chee Chin Liew, Hiroshi Inomata, Kunio Arai, Shozaburo Saito, J. Supercrit. Fluids, 13, 83 (1998).
[ 9] Noriyuki Yoshii, Susumu Okazaki, J. Chem. Phys., 107, 2020 (1997).
[10] Noriyuki Yoshii, Hiromi Yoshie, Shinichi Miura, Susumu Okazaki, J. Chem. Phys., 109, 4873 (1998).
[11] Bernard R. Brooks, Robert E. Bruccoleri, Barry D. Olafson, David J. States, S. Swaminathan, Martin Karplus, J. Comput. Chem., 4, 187 (1983).
[12] William L. Jorgensen, Jayaraman Chandrasekhar, Jeffry D. Madura, Roger W. Impey, Michael L. Klein, J. Chem. Phys., 79, 926 (1983).
[13] Shuichi Nose, M. L. Klein, Mol. Phys., 50, 1055 (1983).
[14] Hans C. Andersen, J. Chem. Phys., 72, 2384 (1980).
[15] K. Tasaki, S. McDonald, J. W. Brady, J. Comp. Chem., 14, 278 (1993).


Return