Phase Transition of Minute Argon Microcluster in Molecular Dynamics Simulations

Yosuke UEDA, Yuri YAMADA and Yosuke KATAOKA


1 Introduction

The pure substances existing around us generally stay in three stable states such as solid, liquid or vapor phase under variable temperature and pressure. Transfer to different phase due to changed intensive variables is referred to as phase transition. The phase transition plays an important role upon analyzing properties of the substance. And the analyses on the microscopic level teach facts that are cooperative phenomena produced by interaction of consisting elements such as atom and molecule. In order to examine the phase transition phenomenon at the microscopic level, we analyzed thermodynamic quantities of microcluster [1] changing as functions of temperature.
The microcluster simulations using MD method [2] generally require disposition of target microcluster in a sufficiently large size cell. And the cell is applied to free boundary condition. However, in cases where the size of such target microcluster is small and molecules in the microcluster are monatomic, the energy source is limited only to intermolecular interaction and kinetic energy. Under such condition, low temperature preset for the simulations may prevent the microcluster from obtaining sufficient energy and such microcluster may remain unchanged. To avoid the result that such microcluster is not able to obtain sufficient energy from outside, the model [3] disposing numerous career gases in addition to molecules which constitute the microcluster in the cell is often used. And the cell is applied to periodic boundary condition. By this method, when career gas particles collide with microcluster, the microcluster can obtain energy from the exterior. As a result, the phase transition phenomenon is promoted. However, this method may cause the behavior of microcluster with preset pressure and the career gas number. And accurate thermodynamic quantities of microcluster may not be calculated by existence of career gas. Accordingly, since the phase transition phenomenon is intentionally promoted using career gas, the calculation result is not necessarily right.
Therefore in this paper, argon microcluster constituted monatomic was taken up for an example, MD simulations using career gas were calculated, and the thermodynamic quantities of pure Ar microcluster on phase transition were obtained. In regard to career gas, we used helium molecule belonging to rare gases. However, Lennard-Jones parameters of career gas were Ar ones for the purpose of simplifying the model. Specifically, interaction terms of internal energy, enthalpy, heat capacity, entropy and free energy were obtained. We also performed MC simulations [4] without career gas and obtained the quantities like the MD method, in order to compare the results from two methods and evaluate the suitability of MD method.

2 Molecular Dynamics Simulations of Argon Microcluster Phase Transition

MD simulations were carried out under various career gas numbers and pressures. We observed the evaporation process from cohesive phase microcluster. In this process, average values of sums of intermolecular potential Ue, volume V and interaction term of pressure Pe by multiplying volume PeV were calculated. The variables such as Ue and PeV were calculated in eq. (1) below and obtained respectively by deducting ideal gas energy values Ui and PiV from the total energy values U and PV.

With respect to temperature, we preset it at 1 K (= 0.00834 e/kB, see below) and raised the temperature by 1 K serration after sufficient stabilization of 8-molecules microcluster constituted by Ar molecules until reaching 100 K during this simulations period. However, when Ue of the microcluster did not reach an equilibrium state in one period of calculation, the next calculation was continued at the same temperature. Therefore, around at phase transition temperature, the number of times of calculations at the same temperature increases inevitably. In this simulation, we think obtaining a statistics average is important, because this simulations result will be compared with the MC simulations result later. Therefore, the factor of time is seldom taken into consideration. In addition , in relation to the initial setting of molecules in the cell, 8 Ar molecules were placed collectively in the center of the cell in order to facilitate microcluster constitution and He as career gas was disposed around them. Concerning intermolecular potential, the eq. (2) of Lennard-Jones potential was applied thereto because both Ar and He are monatomic molecules having the property of similar spherical structure.

In eq. (2), r means distance between two molecules, e and s are both Lennard-Jones parameters when e means well's depth of potential curve and s means intermolecular distance from which potential energy decreases to zero. In this case, the parameters e =1.65×10-21 J and s =3.405×10-10 m were used as suitable for Ar. Furthermore, the parameters e and s were used as units of energy and length respectively. However, in connection with intermolecular potentials between Ar-He or He-He, attractive force term was deleted and only repulsive force term of Lennard-Jones potential was used in order to prevent clusterization of career gas molecules. Table 1 shows setting conditions in relation to the other elements.

Table 1. Setting conditions required for MD method.
Numerical IntegrationGear Method
Temperature ControlVelocity Scaling
Initial Density0.00142g/cm3
MD Steps100000
Time Step1 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.

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.

3 Comparison of Molecular Dynamics Method and Monte Carlo Method

Repeating simulations employing MD method provided in section 2 enables one to identify the most suitable calculation condition for 8 Ar molecules microcluster and we obtained the thermodynamic quantities under such condition. However, it is not ascertained whether thermodynamic quantities of the pure Ar system can be extracted correctly from the Ar and career gas mixed system. The MC simulations were calculated in order to evaluate from a different perspective. If the same results of thermodynamic quantities are obtained from the two different methods, such results shall be reliable. The comparison of calculation conditions of MD method and MC method is shown in Table 3 below:

Table 3. Setting conditions of MD method and MC method.
Calculation TechniqueMD MethodMC Method
Numerical "Integration"Gear MethodMetropolis Method
Used MoleculeAr(8)+He(200)Ar(8)
Initial Volume3318.03 s33795.36 s3
Periodic Boundary ConditionEnabledEnabled
Time Step1fs--

Table 4. PC spec. for MC simulations.
Development languageCompaq Visual Fortran6.5

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.

4 Conclusion

We affirmed that usage of career gas could shorten phase transition calculation of microcluster without intramolecular interaction. In addition, we could extract the thermodynamic quantity of pure Ar correctly from the Ar and career gas mixed system. The results calculated by MD simulations are approximately equivalent to the results using MC method.


[ 1] Koa Kajimoto, Chemistry of Cluster, Baihukan, Tokyo (1992).
[ 2] WinMASPHYC-Pro1.0 user's guide, Fujitsu (2000).
[ 3] Kenji Yasuoka, Mitsuhiro Matsumoto, J. Chem. Phys., 109, 8451-8462 (1998).
[ 4] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford (1987).
[ 5] J.M. Haile, Molecular Dynamics Simulation -Elementary Methods-, Wiley-Interscience, New York (1997).