

Table 1. Setting conditions required for MD method.
| Numerical Integration | Gear Method |
| Temperature Control | Velocity Scaling |
| Ensemble | NTP |
| Initial Density | 0.00142g/cm3 |
| MD Steps | 100000 |
| Time Step | 1 fs |
| 3D Periodic Boundary Condition | |
| Cubic Cell |
In addition, the software and environment of the personal computer used for MD simulations are shown in the following Table 2.
Table 2. PC spec. for MD simulations.
| Software | WinMASPHYC-Pro1.0 |
| OS | Windows98SE |
| CPU | PentiumIII600MHz |
| MEMORY | SDRAM256MB |
The thermodynamic quantities per one Ar molecule were calculated in the following procedure. First, the thermodynamic quantities of Ar and He mixed system were calculated. Next, the thermodynamic quantities of the pure He system were calculated. And we subtracted the latter from the former. Furthermore, the result was divided by number of Ar molecules (see eq.3):

The reliability of the thermodynamic quantities per one Ar molecule solved by using the above procedure is explained in the following section. The thermodynamic quantities were shown separately by number of career gas in Figures 1 - 3.

Figure 1. Ue values per one Ar molecule varying with number of career gas, as a function of temperature. (N=8)

Figure 2. PeV values per one Ar molecule varying with number of career gas, as a function of temperature. (N=8)

Figure 3. V values per one Ar molecule varying with number of career gas, as a function of temperature. (N=8)
Hereinafter how to calculate thermodynamic quantities of each figure are explained briefly. Concerning the values Ue and PeV already shown in eq. (1), the ideal gas energy values were subtracted from the total energy values, because calculation of the interaction terms generated by Ar is the purpose of this research.
For comparison purposes, Figures 1 - 3 include the results obtained by NTV ensemble. In the NTV ensemble, the career gas and the periodic boundary condition are disabled. The velocity scaling method is used for a temperature control. On the result of NTV ensemble, Ue value increases around T = 0.48 and 0.72 e/kB when the structure of Ar microcluster changed. However decomposition of Ar microcluster could not be observed within the range of set-up temperature on NTV ensemble.
Next, in NTP ensemble, periodic boundary condition and disposition of career gas are enabled. In this setting, Ue values per one Ar molecule increase drastically from the point at which the temperature exceeds T = 0.3 e/kB and such values finally reach Ue = 0 e. This result means that the microcluster has collapsed and evaporated with first order phase transition. However the number of career gas influences Ue value. The Ue curves about the different numbers of career gas show relevance between increase of the number of career gas and drop in the phase transition temperature. Although such number increase does not unlimitedly bring the phase transition temperature down it seems that such temperature converges into the real phase transition temperature. Such tendency is explained by the following observations. In the case of 50 He molecules, the phase transition temperature was around T = 0.43 e/kB due to obvious lack of career gas number. However in the cases of 100 and 200 He molecules, their respective phase transition temperatures reached T = 0.34 and 0.32 e/kB, and both temperatures showed little difference. At the same time, in the case of 200 He molecules their phase transition temperatures ranged shortly from T = 0.31 to 0.32 e/kB whereas on the contrary in case of 100 He molecules the temperatures ranged widely from T = 0.34 to 0.45 e/kB. It is considered that lack of energy is the cause. Ar microcluster became unstable due to high energy generated, so it released one molecule and became stable again. At the same time the energy that microcluster received from career gas 100 He molecules became insufficient to maintain its phase transition, and such transition terminated at such temperature. From these results, the large number of career gas brings the phase transition temperature down to its accurate temperature, and shortens the phase transition temperature range and the required time. Figure 2 shows almost constant PeV value notwithstanding change of the number of career gas molecules. According to Figure 3, the assembly volume grows in proportion to the number of career gas. However, the pure Ar system volume was not influenced by the number of career gas molecules.

Figure 4. Ue values per one Ar molecule varying with pressure intensity, as a function of temperature. (N=8)

Figure 5. PeVvalues per one Ar molecule varying with pressure intensity, as a function of temperature. (N=8)

Figure 6. Vvalues per one Ar molecule varying with pressure intensity, as a function of temperature. (N=8)
Figures 4 - 6 show the thermodynamic quantities separated by change of pressure. In the figures, we used the system that contains 200 He molecules because the most accurate thermodynamic quantities were obtained by using 200 He molecules as career gas from Figures 1 - 3.
Each thermodynamic quantity was calculated under pressures of 0.1, 1.0 and 10 MPa. According to Figure 4, at 0.1 MPa level, high energy state has been generated from the beginning and movement of Ar molecules is perceived at around T = 0.1 e/kB. In addition, such energy has been consistently unstable due to low pressure and high amplitude of expansion and contraction of the cell. And the starting point of phase transition has been unascertained because of gradual increase of Ue value until reaching zero. On the contrary, at the 10 MPa level, high pressure prevented the cell from expanding sufficiently and both Ar microcluster and career gas were forced to fit themselves in the small cell despite the temperature rise. On the preceding condition, even if such Ar had been liberated from Ar microcluster, it collided with career gas and reattached itself to microcluster immediately after such liberation. As a result, the phase transition temperature was considerably high. In accordance with Figures 5, 6, 0.1 and 10 MPa curves are unstable at PeV value and particularly such instability is outstanding after the phase transition compared to the stability of the pressure 1 MPa curve.
According to Figures 1 - 6, it is ascertained that the 8 Ar molecules microcluster selected for this simulation requires the career gas number to be 200 molecules, pressure to be 1 MPa and initial density of cell to be one hundredth of liquid density 0.0142 g/cm3.
Table 3. Setting conditions of MD method and MC method.
| Calculation Technique | MD Method | MC Method |
|---|---|---|
| Numerical "Integration" | Gear Method | Metropolis Method |
| Ensemble | NTP | NTV |
| Used Molecule | Ar(8)+He(200) | Ar(8) |
| Initial Volume | 3318.03 s3 | 3795.36 s3 |
| Pressure | 1MPa | -- |
| Periodic Boundary Condition | Enabled | Enabled |
| Time Step | 1fs | -- |
Table 4. PC spec. for MC simulations.
| Development language | Compaq Visual Fortran6.5 |
| OS | Windows98SE |
| CPU | PentiumIII600MHz |
| MEMORY | SDRAM256MB |
The development language and spec. of the personal computer used for MC simulations are shown in the Table 4. The personal computer is the same as that used for calculation of the MD method with the section 2. The program for the MC method was developed by us.
In the case of MC method, career gas was not disposed on the contrary to the MD method. Sufficient number of steps in MC simulations leads to accurate results, so it is expected that this calculation gives us thermodynamic quantities of pure Ar. In addition, we used the volume measured immediately after completion of phase transition for using MD method in MC simulations with molecular extension because the volume is constant according to NTV ensemble. Each thermodynamic quantity calculated under such setting conditions using MC method compared to MD method is shown in Figures 7 - 14 below.
In Figure 7, average internal energy U obtained by using MD method in two systems; Ar and career gas mixture system and pure career gas system is shown. In Figure 8, the product of pressure and volume PV is shown as Figure 7. On one hand, the mixed system curves change with phase transition, on the other hand, the pure career gas system curves are linear. As referring to eq. (3), we obtain thermodynamic quantity per one Ar molecule by subtracting pure career gas system value from mixed system value described in these figures and dividing it by Ar molecules number 8. It is not easy to determine whether the bold procedure using subtraction gives an accurate calculation result of Ar value extracted from mixed phase. However, as the thermodynamic quantity curve of pure career gas phase is the linear function of temperature, so we may say that such result obtained by subtraction is the pure thermodynamic quantity of Ar.

Figure 7. Assembly value and career gas value U obtained by MD method, as a function of temperature.

Figure 8. Assembly value and career gas value PV obtained by MD method, as a function of temperature.
Figures 9, 10 compare the results of Ue and PeV in MD method obtained based on Figures 7, 8 with such results obtained by using MC method. As already stated in section 2, U and PV are constituted by two elements, i.e. ideal gas terms and pure interaction terms, see eq. (4).

Upon comparison of MD method and MC method, interaction energy plays an important role. Therefore, the thermodynamic quantities were compared as interaction terms in Figures 9 - 14. The foregoing considered, when we refer to Figures 9, 10, both curves show similar tendency except that PeV fluctuation on the MD method curve is relatively greater than MC method.

Figure 9. Comparison of Ue values per one Ar molecule using MD and MC method. (N=8)

Figure 10. Comparison of PeV values per one Ar molecule using MD and MC method. (N=8)
In Figure 11, the interaction term of enthalpy He using MD method in accordance with eq. (5) and fitting curve obtained from He are shown.

Figure 11. He value per one Ar molecule using MD method and fittint curve. (N=8)

The fitting below is necessary to calculate the interaction term of Gibbs free energy Ge, because the calculation procedure to obtain Ge uses fitting function. In addition, Figure 11 does not include the MC method curve because the interaction term of Helmholtz free energy Ae making pair with Ge has been calculated for MC method, therefore He is not required to obtain Ae.
In Figure 12, we obtained two kinds of interaction terms of heat capacities Cpe and Cve. The heat capacity at constant pressure Cpe is calculated by MD method and another heat capacity at constant volume Cve is calculated by MC method. Each heat capacity was calculated by partial differential in the following eq. (6).

Figure 12. Comparison of heat capacities per one Ar molecule using MD method and MC method.

Although their peaks are not in the same position, such difference is not problematic since they correspond to each phase transition temperature variation.
Interaction term of entropy Se by usage of numerical integration described in eq. (7) below is put in Figure 13:


Figure 13. Comparison of Se values per one Ar molecule using MD method and MC method. (N=8)
In eq. (7), the temperature interval DT is constant and it shall be set in accordance with the equation DT=T1. Accordingly, the first term shall be as (Ce1/T1)×T1=Ce1, and Se curve and Ce curve shall have the same starting point.
MD method and MC method were compared using relative value [5] defined by Se1= Ce1. The MD method has a relatively higher Se value than the MC method because their values are based on difference of heat capacity Cpe and Cve. As the MD method takes into consideration the energy contribution from PeV, the Se value has become higher than that obtained by using MC method.
In Figure 14, comparison of free energy interaction terms obtained from MD method and MC method is indicated.

Figure 14. Comparison of Free Energies per one Ar molecule using MD method and MC method. (N=8)

As entropy is used also in eq. (8), the relative values of free energy were compared in the same manner described in the preceding paragraph. There is a difference between constant pressure and constant volume conditions in these two energies. It was able to confirm that phase transition from microcluster to vapor is a voluntary transition because both free energy values show similar curves whose slopes incline in a more negative direction after phase transition. Consequently, we shall conclude that the simulations results of MD method and MC method correspond to each other.