Model and Method - 液態鎵的結構、動力學及熱力學性質之研究

III.3.1 Model

Shown in Fig. III.1, the Ga interatomic pair potential at T = 323K, φ(r), has a

0.6 0.8 1 1.2 1.4 1.6

Figure III.1: Interatomic pair potential, φ(r), (solid line). The dot-dashed line is the LJ potential with the same ε and σ. The dashed line is a linear function a1(r/σ) + a2 with a1 = −115 and a2 = 109; the dotted line is the function c1exp(−c2r/σ)(r/σ)^c³ , where c¹ = 4.42 × 10⁶, c² = 16.4 and c³ = 1.5

ledge-shape repulsive core and an oscillatory part, whose long-range behavior generally follows the Friedel oscillations [58]. Two parameters of φ(r) are ε, the depth of the main attractive well, and σ, the shortest distance at φ(r) = 0; the value of ε corresponds to 47K and σ = 4.04˚A. The repulsive core of φ(r), denoted as φ0(r), is the pure repulsive potential for r smaller than r⁰, which is the position of the main attractive well, and r⁰ is about 1.07σ. The pure repulsive potential, φ0(r), can be roughly divided into three sections. As r less than 0.7σ, φ0(r) increases almost exponentially with decreasing r so

that the repulsive core in this region is extremely stiff. In the intermediate region between 0.9σ and 0.7σ, the value of φ0(r) increases roughly from 6ε to 30ε and the shape of φ0(r), becoming softer, has a ramp-like behavior with a reflection point around 0.8σ [86]. As r is larger than 0.9σ, φ⁰(r) behaves like the repulsive core of the LJ potential with the same σ and ε, and the interaction potential is less than 6ε [60].

III.3.2 Method

We have carried out a series of molecular dynamics simulations at constant NVT ensemble. In each simulation, we fix the temperature to be 6.85ε which corresponds to T = 323K, and the box size is equal to 10.2σ. All quantities given in this chapter are in units of σ, ε and the mass of Ga atom. The number of simulated particles starts at N = 3500. Then, in each simulation for a new thermodynamic state, N is reduced by 100, with the lowest N being 100; therefore, we have a series of MD generated configurations at each density. At N = 3500, the reduced system density ρ is equal to 3.305, which is close to that of liquid Ga at T = 323K and pressure of about 1 bar. At ρ = 3.305, the static and dynamic structure factor of this simulated model agree well with the experimental results at temperature close to the melting point (Tm = 303K), and dynamic anomaly of liquid Ga is reproduced [58]. In the following, the thermodynamic properties of the model fluid will be investigated with simulation and various hard-sphere perturbation theories.

Figure III.2: Radial distribution function of the model fluid at several reduced densities.

In each panel, the solid is the simulation result and the dashed line is the approximation given by Eq. III.10 with the Lado-WCA value of σHS. The dotted line indicates the σHS

value obtained by the E-MCRS theory.

III.4 Results

The radial distribution functions of these simulated systems at various reduced system density are shown in Fig. III.2. At ρ = 3.305, the main peak of g(r) is located at 0.686σ, well inside the repulsive core of φ(r). By descreasing ρ, the main peak of g(r) gets lower and shifts outwardly, and a shoulder near σ appears and grows. Moreover, around ρ = 2.0, the roles of the main peak and the shoulder exchange, because the main peak occurrs near

σ and the shoulder resides in the inner side of the repulsive core. As ρ keeps on decreasing, the in-inside shoulder finally disappears at ρ = 1.35. The new main peak near σ keeps the same shape and moves toward the position of the first attractive well with descreasing ρ. Similar results of the variation of g(r) with density are also observed for the one-scale and two-scale ramp potentials [87], introdcued by Jagla [88]. On the other hand, the approximation in Eq. III.10 for our model fluid are examined by evaluating y_HS(r) in the approximate Eq. III.11 with the Lado-WCA HS diameter for some reduced densities, and the results are also shown in Fig. III.2. For each density, Eq. III.10 is only good for g(r) inside the value of σ_HS, which is obtained by the E-MCRS theory; however, in following, the calculations need the contributions of gHS(r) inside the σHS.

w(r), giving the mean force acting on a particle in the fluid, is the other way to manifest the density effect on the variation of g(r) in our model. The definition of w(r) have been given in Section. III.2.2, and the variation of w(r) with reduced density ρ is shown in Fig. III.3. At ρ = 3.305, because of the high compactness in the fluid, w(r) has a deep attractive well, which is at the first-peak position of g(r). Furthermore, the depth of this attractive well gradually reduces and a shoulder near σ appears by decreasing ρ.

Around ρ = 1.889, the roles of this attractive well and the shoulder in w(r) switch with each other, as the similar case in Fig. III.2. Beginning at ρ = 1.889, the new attractive well in w(r) continuously moves out and the new attractive well become deep. Finally, the well is attenuated to ε.

Fig. III.4 shows the variation with density of the HS diameter, which are estimated

0.5 0.75 1 1.25 1.5 1.75 2

r / σ

-10 0 10 20 30

Potential

φ(r)

w(r) for ρ=3.305 w(r) for ρ=2.361 w(r) for ρ=1.889 w(r) for ρ=0.944 w(r) for ρ=0.094

Figure III.3: Potential of mean force w(r) of the model fluid at ρ = 3.305 (dashed line), ρ = 2.361 (short-dashed line), ρ = 1.889 (dot-dashed line), ρ = 0.944 (dot-dot-dashed line) and ρ = 0.094 (dotted line). The solid line is the interatomic pair potential φ(r).

All potentials are scaled with ε and distance is scaled with σ.

0 1 2 3

Figure III.4: The variation of effective HS diameter σHS as a function of density. σHS is estimated by the WCA (triangles), Lado-WCA (diamonds), MCRS (circles) and E-MCRS (squares) theories. The crosses stand for the main-peak position of g(r).

by the E-MCRS, MCRS, WCA, and L-WCA theories. The values of the HS diameter estimated by the four perturbation theories decrease with increasing density, and our results are consistent with the results that the effective HS diameter of a fluid is decreased with increasing the fluid density [33]. At high densites, both WCA and L-WCA σHS

values are larger than the main-peak position of g(r), which is inside the repulsive core;

however, the HS diameters become smaller than the position of main-peak of g(r) below ρ = 2.0. Not only the σHS values of the MCRS theory is limited by the peak position of g(r) but also the E-MCRS value is almost along the track of the peak position as ρ

is larger than 2. According to Fig. III.1, as the density is below ρ = 2.0, because the exponential-like behavior of average distance between particles is replaced by ramp-like behavior, the position of the main-peak of g(r) makes a jump to a larger value around ρ = 2.0. Therefore,the HS diameters of the MCRS and E-MCRS theory are apparently released from the restriction due to the main-peak position of g(r) and ascend manifestly as ρ varies from 2 to 1. The MCRS value is roughly fixed at 0.9σ as density approaches to zero, but the E-MCRS value continuously increases at very low densities and passes over the values of the WCA and Lado-WCA theories.

According to the pressure equation, the compressibility factor Z(ρ) can be calculated, and the results are shown by symbols in Fig. III.5. In order to perform numerical integration of Z(ρ), we fit Z(ρ) in Fig. III.5 with the series in Eq. III.19 truncated beyond the fouth order of ρ and treat the virial coefficients as fitting parameters. The simulated data of Z(ρ) can be separated into two parts. One is larger than ρ = 1 and the other is smaller than ρ = 1, and we give the values of the fitting parameters in the caption of Fig. III.5. The inset of Fig. III.5 shows the second Virial approximation that check the accuracy of our calculations at low densities.

In the MCRS and E-MCRS theories, the variational curves of excess Helmholtz free energy ˜A for some reduced densities are shown in Fig. III.6, and each curve has a mini-mum, which corresponds to a HS diameter. In these two theories, the approximate gHS(r), obtained by Verlet algorithm [30], is used, and yHS(r) is approximated by Eq. III.11 with

0 1 2 3

0 10 20 30 40 50

simulation data fitting function

0 0.25 0.5 0.75 1

2 4 6

Second Virial Approx.

Figure III.5: Compressibility factor, Z = βP/ρ, as a function of density. The circles are calculated by the pressure equation with the simulated g(r). The solid line is the Virial function III.19 with the Virial parameters: B2 = 0.995, B3 = 4.129, B4 = −4.969, and B5 = 4.67 for ρ > 1, and B2 = −0.323, B3 = 5.228, B4 = −0.523, and B5 = 0.055 for ρ < 1. The inset shows Z at low densities and the second virial approximation (dashed line).

0.64 0.66 0.68

Figure III.6: Variation of A^ex/N kBT on σHS at the indicated reduced density. In each panel, the solid and dashed lines are the curves for the E-MCRS and MCRS A^ex/N kBT , respectively. The dot-dashed line is the variational A^ex/N k_BT with Mon’s correction term given in Eq .III.8. At ρ = 0.472 and ρ = 0.094, the solid and dot-dashed lines are almost the same.

mum of the variational curve shifting to a larger effective HS diameter and a lower value of ˜A than those obtained by the MCRS theory. These results are expected to be true for all kinds of fluid [33]. The accuracy of the approximation in Eq. III.10 is examined again by calculating excess Helmholtz free energy. In Fig. III.6, although the corresponding variational curve of Eq. III.9 generally shifts upward relative to the one of Eq. III.16, the HS diameter σ_HS, which is estimated by Eq. III.9 is the same as that, which is es-timated by Eq. III.16, and the overall thermodynamics properties of our calcualtion are not effected. By using the fitting function Eq. III.19, we have done the integration in Eq.

III.17 and the results are indicated by the symbols in Fig. III.7. Shown in Fig. III.7, the excess Helmholtz free energies estimated by the various HS perturbation theories, includ-ing the WCA, L-WCA, MCRS and E-MCRS theries, and the comparison between these results is presented. Although the MCRS prediction slightly deviates from the simulated excess Helmholtz free energy, Mon’s correction cause the E-MCRS prediction to match the simulated data well for entire densities. Apparently, the variational approaches do a better prediction than the WCA [73, 74, 75, 76] and Lado-WCA [77] theories, which are only good for the region of ρ < 1.

Fig. III.8 shows the behaviors of simulated and approximate excess entropies ˜S. In thermodynamics, ˜S has a linear behavior at low densities because the excluded volume effect due to the repulsive core of φ(r). Intriguingly, the density curve of ˜S has a small positive slope in the intermediate range, which is roughly from ρ = 1.4 to 1.8. The anomalous behavior of excess entropy is also observed in a two-scale ramp potential [89].

0 1 2 3

0 10 20 30 40

Aex / Nk BT

simulation data E-MCRS MCRS WCA Lado-WCA

Figure III.7: The excess Helmholtz free energy ˜A = A^ex/N k_BT predicted by the E-MCRS (solid line), E-MCRS (dash line), WCA (dotted line) and Lado-WCA (dot-dashed line) theories and the result calculated by Eq .III.17 (symbols).

0 1 2 3

-5 -4 -3 -2 -1 0

Sex / Nk B

simulation data EMCRS

MCRSHS fluid wth E-MCRS σ_HS

Figure III.8: Density dependence of the excess entropy ˜S = S^ex/N kB. The solid and dashed lines are the predictions of the E-MCRS and MCRS theories, respectively. The dot-dashed line is the HS-fluid contribution in the E-MCRS theory. The symbols are the results of thermodynamic calculation.

Beyond the intermediate region, ˜S decrease with density in a slower decreasing rate, be-cause the density effect is more significant than the diameter reduced effect again. In approximation, the prediction of the MCRS theory agrees well with the result of ther-modynamic calculation at low densities, and also has smaller values of excess entropy anomaly in the intermediate range; however, the excess entropy of MCRS’s prediction are deviated from simulation data at higher densities. As the E-MCRS diameter is substi-tuted into Eq. III.23, which is the excess entropy of a HS fluid , the absolute value of the excess entropy with a manifested anomaly is larger than that estimated by the MCRS theory. It is amazing that the Mon’s correction in Eq. III.24 makes the prediction of the E-MCRS theory have a triumph over the MCRS theory for a perfect agreement with the simulation results in the entire density. It is easy to be understood: The contribution of Mon’s correction becomes gradually larger with inceasing density; thus, Mon’s correction is more important at higher densities and it gives a fundamental improvement of the HS variational prediction for the excess entropy of a fluid.

Fig. III.9 shows the density dependence of the packing fraction η of the HS fluid estimated by various perturbation theories. For either the WCA or the Lado-WCA value of σHS, the η increases monotonically with increasing density; however, their η values are over 0.7 at reduced system density ρ = 3.305, which is close to the density of realistic liquid Ga. The large packing fraction, which is in the solid phase of the HS system, indicate that the overestimation of the effective HS diameter makes the reference HS fluid too much

0 1 2 3

0 0.2 0.4 0.6

EMCRS MCRS WCA L_WCA

Figure III.9: Density dependence of packing fraction η of the HS fluid with HS diameter of the E-MCRS (solid line), MCRS (dashed line), WCA (dotted line) and L-WCA (dot-dashed line) values shown in Fig .III.4.

estimated by the MCRS and E-MCRS theories, the packing fraction (η = πρσ_HS³ /6) of the MCRS theory is reasonably smaller than that of the E-MCRS theory. In following, we explain that why the anomaly in ˜S predicted by the MCRS and E-MCRS theories occurs in our model. In Fig. III.9, the density curves of MCRS and E-MCRS theories also have apparently a anomaly, which shows a negative slope from ρ = 1.3 to 1.7. We maintain that the negative slope corresponds to the substantial shrinkage in the size of the effective HS with increasing density at this region as shown in Fig. III.4; thus, the anomaly appears in the packing fraction and the same physical picture cause the occurrence of anomaly in excess entropy for MCRS and E-MCRS theories. Beyond this anomaly region, the η continuously increases and ˜S continuously decrease, because the contribution of density effect more than that of particle size effect. In the E-MCRS and MCRS theories, the η of the HS reference systems at ρ = 3.305 is close to the physcially reasonable values: 0.473 and 0.5, which are in the fluid-solid coexistence region of the HS system.

III.5 Conclusions

The interatomic pair potential of liquid Ga model has a ledge-shape repulsive core and a long-range oscillatory part induced by conduction electrons. We have investigated thermodynamic properties of this model fluid to investigate the applicability of some HS perturbation theories which include the WCA, Lado-WCA, MCRS and E-MCRS approx-imations. These approximate methods have been used for the Lennard-Jose,

inversed-power fluids and liquid Na. The MD simulations of our model fluid at constant NVT ensemble from low to high densities are performed, and the highest density corresponds to realistic liquid Ga at T = 323K. There is an interesting above the density variation of the model fluid: The main peak of the radial distribution function is located in the repulsive core region of the interatomic pair potential at high density, the position of the main peak suddenly jumps to a larger distance at some intermediate density and the main peak moves outwardly to the first attractive well of the interatomic pair potential as the density is continuously decreased.

Although the WCA and Lado-WCA methods are accurate for the Lennard-Jones flu-ids, the predicted HS diameters are too large for our model at high densities. The packing fractions in the WCA and Lado-WCA theories are over the freezing point of the HS ref-erence system so that they are failure for realistic liquid Ga. Based on a variational approach, the MCRS and E-MCRS theories take into account the effect of the anomaly of packing fraction in a region of intermediate density, and the HS reference systems for the model fluid are physically reasonable at entire density. According to the deficiency of the HS perturbation theory, Mon gives a correction which is related to the configuration space of the hard-sphere fluid. Considering this correction, the E-MCRS theory indeed improves the predictions on the thermodynamics properties of the model fluid: the estimated Hel-mohltz free energy and entropy by E-MCRS method are closer to the simulation results than those by others mothods. An anomalous region, associated with the shrinkage of the effective HS diameter by increasing system density, is found in the excess entropy of our

model varied with density. For the predictions of the thermodynamic properties of liquid Ga, our results suggest that the E-MCRS perturbation theory is the most accurate one in all HS perturbation theories.

Chapter IV Conclusions

In this thesis, the model of simulation is a gallium interatomic pair potential, which consists of a legde-shape repulsive core and the long-range Friedel oscillations induced by the conduction electrons, is obtained from the first-principles GEINMP theory at the thermodynamic conditions of liquid gallium close to the melting point. The repulsive core of the pair potential varies continuously from an exponential-decay inner core, through a ramp-like intermediate region, to a LJ-like outer core. By MD simulation, we have successfully reproduced the static structure factor and dynamic structure factor which agree well with the results of experimrnts; further, the structures, dynamic properties, and thermodynamic properties of this pair potential have been investigated for studying the realistic liquid gallium. We summarize the conclusions in the following.

In chapter II, a shoulder appears on the high-q side of the first peak in the static structure factor and an anomaly appears on the linewidth of dynamic structure factor, in good agreements with the experimental data of liquid gallium. According to our results,

the structure anomaly is determined by the ledge-shape repulsive core and first three attractive well. In local structure, a modulation produced by the Friedel oscillation makes the apearance of some solid-like cluster, which are, in concept, more or less like the Peierls-distortion mechanism proposed for liquid-Arsenic before [90, 91]. It is suggested that the shoulder in the static structure factor is associated with the appearance of these solid-like clusters.

Also, the anomaly in the linewidth of the dynamic structure factor is determined by the ledge-shape repulsive core and the first two attractive wells. The dynamic anomaly, which can be well described by the revised Enskog theory with the HS diameter determined by the first-peak distance of the radial distribution function, is interpreted by cage diffusion.

The structure and dynamic anomaly are confirmed to occur at the same wavenumber, which suggests that the ion-density fluctuation at this wavevector has the same wavelength with the Friedel oscillations induced by the conduction electrons. Thus, at this wavevector, the density waves of ions and electrons are coherent and the attractions between the two systems would be more rigid so that the cage around each ion becomes harder.

In chaptr III, the radial distribution function of the model has an interesting variation with number density: At high density, the main peak of the radial distribution function is located in the inner region of the repulsive core with an exponential decay. At some critical density, the position of the main peak makes a discontinuous jump to the outer core of the LJ-like region. As the density is further decreased, the main peak moves

To predict the thermodynamic behaviors of our model, the WCA, Lado-WCA, MCRS and E-MCRS methods are used. For the prediction of HS diameter, the packing fractions, which are estimated by the WCA and the Lado-WCA methods, are over the freezing point of the HS system and this makes the WCA and Lado-WCA theories breakdown for our model at high densities. In MCRS and E-MCRS theories, due to the estimated

在文檔中液態鎵的結構、動力學及熱力學性質之研究 (頁 56-77)