Next Article in Journal
Plausible Physical Mechanisms for Unusual Volatile/Non- Volatile Resistive Switching in HfO2-Based Stacks
Previous Article in Journal
Acknowledgment to Reviewers of Condensed Matter in 2020
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Scaling of Phase Diagram and Critical Point Parameters in Liquid-Vapour Phase Transition of Metallic Fluids

Shiv Enclave, 304, 31-B-Wing, Tilak Nagar, Mumbai 400089, India
Retired from Bhabha Atomic Research Centre, Mumbai 400085, India.
Submission received: 14 December 2020 / Revised: 23 January 2021 / Accepted: 25 January 2021 / Published: 28 January 2021

Abstract

:
The first objective of this paper is to investigate the scaling behavior of liquid-vapor phase transition in FCC and BCCmetals starting from the zero-temperature four-parameter formula for cohesive energy. The effective potentials between the atoms in the solid are determined while using lattice inversion techniques as a function of scaling variables in the four-parameter formula. These potentials are split into repulsive and attractive parts, as per the Weeks–Chandler–Anderson prescription, and used in the coupling-parameter expansion for solving the Ornstein–Zernike equation that was supplemented with an accurate closure. Thermodynamic quantities obtained via the correlation functions are used in order to obtain critical point parameters and liquid-vapor phase diagrams. Their dependence on the scaling variables in the cohesive energy formula are also determined. An equally important second objective of the paper is to revisit coupling parameter expansion for solving the Ornstein–Zernike equation. The Newton–Armijo non-linear solver and Krylov-space based linear solvers are employed in this regard. These methods generate a robust algorithm that can be used to span the entire fluid region, except very low temperatures. The accuracy of the method is established by comparing the phase diagrams with those that were obtained via computer simulation. The avoidance of the ’no-solution-region’ of the Ornstein-Zernike equation in coupling-parameter expansion is also discussed. Details of the method and complete algorithm provided here would make this technique more accessible to researchers investigating the thermodynamic properties of one component fluids.

1. Introduction

The scaling and universal features in phase transition theory were first brought out with the van der Waals equation of state [1]. Scaling the thermodynamic variables provided a universal equation of state and critical behavior that is characterized via universal critical exponents. The mean field approach implied in the van der Waals equation neglects long range fluctuations that are close to the critical point and so provide only the classical critical behavior as opposed to the renormalization group theory of critical phenomena [2]. The Ornstein–Zernike equation (OZE) defines the relationship between the short range and long range correlation functions in fluids. When supplemented with appropriate closure relations, which relate the correlation functions to inter-particle interaction, OZE provides a framework to apply the mean field theory to fluids with arbitrary inter molecular potentials [3]. The thermodynamic properties with an accuracy comparable to those obtained via simulation are presently obtained with the OZE.
The zero-temperature cohesive energy for a variety of metals follows a universal curve with suitable scaled variables [4]. It is then natural to explore whether the universal energy-volume curve can be extended to the expanded volume states to describe the liquid-vapor phase transition in metallic fluids [5]. This investigation becomes easy when used with the corrected rigid spheres (CRIS) model, which is a thermodynamic perturbation theory that is based on the cohesive energy of nearest-neighbor pairs in the fluid [6]. In comparison to other perturbation theories [7], which are based on inter-particle interaction potentials, the CRIS model provides an approximate theory of fluids based on cohesive energy curves. However, the accurate prediction of the liquid–vapor phase diagram within the CRIS model needs a tuning of the reference-hard-sphere pair distribution function [8].
The first objective in this paper is to investigate the scaling aspects of the liquid–vapor phase transition in metallic fluids, starting from the scaled cohesive energy (SCE) versus volume curves. An improved SCE formula [9], involving four-parameters, is used for this purpose. First of all, effective pair interaction potentials are derived from the SCE formula by employing lattice inversion techniques, and split into repulsive and attractive components while using the Weeks–Chandler–Anderson (WCA) prescription [10]. These components are used in an accurate thermodynamic perturbation theory, called coupling-parameter expansion (CPE) [11], for solving the Ornstein-Zernike equation (OZE) and an appropriate closure relation. The coupling parameter ( 0 λ 1 ) tunes the strength of the attractive component of the potential, and all correlation functions are expressed as Taylor’s series in λ around λ = 0 . It is important that the method offers series expansion to arbitrary orders in Taylor’s series. Therefore, all of the thermodynamic quantities are easily computed in the entire phase plane once the correlation functions are determined to sufficient accuracy. The critical point parameters and phase diagrams of metallic fluids are then obtained in terms of the scaling variables in the SCE formula. The results for critical point parameters are similar to those obtained earlier [5] while using the CRIS model; however, there are important differences. For example, the FCC and BCC lattices are found to generate different type of critical point parameters. Further, the scaling of phase diagrams was not considered earlier.
The second objective of the paper is to revisit, in some detail, the methods to be used for implementing the CPE technique. It needs the correlation functions that are given by the non-linear OZE for a reference system, and solutions to a hierarchy of linear integral equations for derivatives of the correlation functions. It is most appropriate to use the Newton–Armijo non-linear solver and Krylov space-based linear solvers [12] for this purpose. In addition to facilitating computation of correlation functions to high order perturbation theory, the CPE technique avoids the ’no-solution region’ of the OZE in the liquid-vapor transition region. Indeed, the occurrence of such regions, for the hypernetted chain (HNC) closure [13] and other closure relations [14], is a bottleneck in applying OZE with the full potential.
The following sections discuss different aspects of the paper, such as the SCE versus volume formula, corrections to be applied for small volumes, lattice inversion to derive potentials, and algorithmic details of the CPE technique. The simulation results on phase diagrams, taken from literature, for specific potentials are compared with the results that were obtained while using present method to establish its accuracy. Critical point parameters and phase diagrams are obtained for potentials that correspond to FCC and BCC lattices. The Appendix A provides details of the complete algorithm for easy implementation of CPE.

2. Scaled Cohesive Energy Formula

The Fermi-pressure in degenerate electron systems, which is a pure quantum effect, manifests as the zero-temperature isotherm in metals. This significantly contributes to total pressure in compressed solids, and it becomes the dominant component at strong compression. Density functional theories (DFT) are routinely used now [15] in order to generate energy versus specific volume (or volume per atom) tables, which are then easily incorporated into semi-empirical formula. Such a formulation for zero-temperature energy and pressure, involving four-parameters [9], is expressed in the SCE formula:
E c ( V ) = E 0 ( 1 + a + δ a 3 ) e a ,
P c ( V ) = 3 B 0 ( 1 ( V / V 0 ) 1 / 3 ) ( V / V 0 ) 2 / 3 ( 1 3 δ a + δ a 2 ) e a , η = ( 9 B 0 V 0 / E 0 ) 1 / 2 , a = η ( ( V / V 0 ) 1 / 3 1 ) , δ = ( 1 / 2 ) ( B 0 1 ) / η 1 / 3 .
Here, the variable V is atomic volume or volume per atom. The four parameters in this model are the atomic volume V 0 , the bulk modulus B 0 , its pressure derivative B 0 , and the cohesive energy per atom E 0 at equilibrium conditions. This is a refinement over Rose equation [4], and it provides accurate zero-temperature energy and pressure in compressed states up to V 0 / 2 , (which corresponds to about 150 GPa) , and expanded states up to about 2 V 0 for about forty metals [9]. If energy E c and pressure P c are scaled with E 0 and B 0 , respectively, there only remains two dimensionless parameters η and δ in these expressions. Note that the parameter a is related to the dimensionless length variable ( V / V 0 ) 1 / 3 , which is simply the (scaled) side length of the atomic volume. Unless specified explicitly, length and energy will be scaled with V 0 1 / 3 and E 0 , respectively, throughout. Negative a corresponds to compressed states while it is positive for the expanded states. Furthermore, the equilibrium atomic volume V 0 , which corresponds to 0 K , should be determined, such that the sum of zero-temperature pressure and thermal pressure of ions and electrons is just one bar at 300 K.

Correction for Strong Compression

The expressions given above are inadequate in the region of strong compression. This is evident from Equation (1), which approaches a finite value as V 0 . In this limit, energy and pressure must approach those of an electron gas around the nucleus, as given by the quantum statistical model (QSM) [16]. This model provides accurate electronic properties above ∼250 GPa of pressure, as it accounts for exchange and correlation effects in addition to incorporating corrections for electron density gradients [17]. Electron pressure in a compressed atom within the QSM model is analytically expressed as:
P q ( V ) = 2 2 m e [ 2 5 ( 3 π 2 ) 2 / 3 N s 5 / 3 2 a B 13 36 ( 3 / π ) 1 / 3 N s 4 / 3 ] , N s ( V ) = Z n V E x p [ α ( V V 0 ) 1 / 3 β ( V V 0 ) 2 / 3 ] , α = 0.1935 Z n [ 0.495 0.039 ( log 10 Z n ) ] R w a B , β = [ 0.068 + 0.078 ( log 10 Z n ) 0.086 ( log 10 Z n ) 2 ] ( R w a B ) 2 .
Here, is reduced Planck’s constant, m e electron mass, a B Bohr radius, Z n atomic number, and R w the equilibrium Wigner–Seitz radius. The quantity N s ( V ) is electron density at the atomic cell surface, and the negative contribution in pressure is due to the exchange effect. This contribution extends the range of validity of P q ( V ) to larger V . The internal energy per atom, E q ( V ) , is obtained by integrating the thermodynamic relation P = d E / d V from a suitable initial volume, say V 0 .
It is easy to use an interpolation procedure [18] to smoothly go from SCE model to QSM. Choose a volume V m , such that the SCE model is accurate for V V m , so that the corrected zero temperature pressure and energy are, respectively, P 0 ( V ) = P c ( V ) and E 0 ( V ) = E c ( V ) for V V m . The interpolated expressions for smaller volumes are expressed as:
E 0 ( V ) = [ E q ( V ) E q ( V m ) ] B i ( V ) + E c ( V m ) , V V m
P 0 ( V ) = P q ( V ) B i ( V ) + [ E q ( V ) E q ( V m ) ] B i ( V ) , V V m .
where B i ( V ) = [ 1 + b 1 V + b 2 V 4 / 3 + b 3 V 5 / 3 ] is an interpolating function. The parameters b 1 , b 2 , and b 3 are chosen, such that P 0 ( V ) and its first two derivatives are continuous at V m [18]. Note that E 0 ( V ) is continuous at V m by definition. This procedure gives a smooth transition from SCE model to the QSM, and it is expected to be better than modifying the definition of E c ( V ) arbitrarily [19].

3. Lattice Inversion for Potential

The lattice inversion method [20] is a convenient tool for extracting effective inter-particle potentials from E 0 ( V ) , which is re-expressed in dimensionless functional form as E ( x ) . More specifically, E ( x ) = E 0 ( V 0 x 3 ) / E 0 , where x = ( V / V 0 ) 1 / 3 is the scaled side length of the atomic volume. The lattice is imagined as an assembly of successive shells around a central atom, and so E ( x ) is expressed as a lattice sum over inter-particle potential U ( r ) , where r is the (scaled) nearest-neighbor distance (NND):
E ( x ) = 1 2 n = 1 γ n U ( b n z 0 x ) .
Here n is the shell index, γ n number of atoms on the shell, and b n the normalized shell-radius ( b 1 = 1 ) in units of NND. The factor 1/2 in this equation arises as U ( r ) is defined for a pair of atoms, while E ( x ) corresponds to a single atom. For numerical applications, the infinite sum is truncated at some finite shell (∼20), so that U is negligible there after. The NND and x are related as r = z 0 x , where z 0 is a constant that is dependent on the lattice type. For the FCC lattice (four atoms in unit cell) z 0 = 2 1 / 6 , while for BCC lattice (two atoms in unit cell) z 0 = 3 / 4 3 . On separating the first shell-contribution, Equation (6) is expressed as:
U ( r ) = 2 γ 1 E ( r / z 0 ) 1 γ 1 m = 2 γ m U ( b m r ) .
The subtracted term is the contribution to U ( r ) from second and higher shells. This equation is solved via iteration [21] to express U ( r ) explicitly in terms of E ( x ) [20]. However, for computational purposes, it is reasonable to assume that U ( r ) = 0 for r r for a sufficiently large value of r . Subsequently, Equation (7) is solved [21] by marching to lower values of r starting from r , as the second term is already computed. It is advantageous to define a geometrically progressing mesh over the required domain, r 1 r r , and interpolate U ( r ) for in-between values of r for computing the second term. Alternatively, starting with the approximation provided by the first term, few iterations on the second term readily provide a converged profile of U ( r ) .
Tables for the lattice dependent constants have been enumerated [22] for 50 shells. In fact, tabulations are provided for elementary sub-lattices, where atoms are placed only at cell corners (SC), faces centers (FC), and edge centers (EC). All of the cubic lattices (FCC, BCC, ECC, etc.) are decomposed into these sub-lattice types. The lattice dependent constants γ n and b n for 1 n 20 for FCC and BCC lattices are reproduced in Table 1 and Table 2, respectively, for convenience.

Chen-Möbius Formula

A remarkable explicit formula for the potential in terms of E ( x ) is obtained [23] by generalizing the lattice sum expression to the form:
E ( x ) = 1 2 n = 1 Γ n U ( B n z 0 x ) .
Here, the expanded sequence of shell radii { B n } monotonically increses ( B 1 = b 1 = 1 ) , and it contains the original sequence { b n } . Most importantly, the expanded sequence has the additional property that it forms a semi-group. That is, for any two integers n and m, there exists an integer p, such that B p = B n B m . The elements of { B n } that are not contained in the original set { b n } are radii of virtual shells, and the corresponding occupation numbers in the expanded set { Γ n } are zero. Thus, for the FCC lattice, the expanded radii are given by: B n = n . Hence, 14 is a virtual shell and Γ 14 = 0 as B 14 = 14 is not in the original set { b n } (see Table 1). For the BCC lattice, the set { B n } is to be generated appropriately from { b n } , so as to have the semi-group property. Now, motivated by Equation (7), a general expansion for U ( r ) is attempted in the form:
U ( r ) = 2 m = 1 I m E ( B m r / z 0 ) ,
where { I m } are weight factors for different shells. Substituting this expansion in Equation (8) gives:
E ( x ) = n = 1 Γ n m = 1 I m E ( B m B n x ) .
The double sum on the right-hand-side can be rewritten by grouping terms while using the semi-group prperty. Substituting B m B n = B p , and grouping terms yields:
E ( x ) = p = 1 E ( B p x ) m = 1 p I m Γ m [ p ] δ [ B p B m B m [ p ] ] .
Here, δ [ x y ] equals 1 for x = y and 0 for x y . For a given p, the indices m and m [ p ] are those satisfying the semi-group property. The second sum terminates at p, due to the monotonic property of the elements in { B p } . Now, Equation (11) reduces to an identity if { I p } are chosen, so that the second sum is 1 for p = 1 and 0 for p 2 . A recursion formula for I p , for p 2 , then follows on separating out the last term in this sum:
I 1 = 1 Γ 1 ; I p = 1 Γ 1 m = 1 p 1 I m Γ m [ p ] δ [ B p B m B m [ p ] ] , p 2 .
Note that m [ p ] = p for p = 1 . Thus, the weight factors { I p } in Equation (9) are readily computed once { Γ n } and { B n } are specified. This choice is already discussed for FCC lattice, and Table 3 provides the generated inversion-constants for 20 shells. Note that shell 14 is virtual and Γ 14 = 0 .
Generating a new set of shell radii { B n } and, consequently, { Γ n } is more involved in the case of BCC lattice. A possible approach [23] is to add virtual shells of radii { b 2 b 2 j , b 3 b 2 j , b 4 b 2 j , , f o r j = 1 , 2 , } to the existing set { b n } . The occupation numbers Γ n of these new shells are zero. Subsequently, { B n } is generated by reordering the enlarged sets to obtain monotonically increasing shell radii which also satisfy the semi-group property. This procedure is not unique and several virtual shells are introduced, as seen in Table 4 and Table 5. These tables, which cover a total of 30 shells, have 19 virtual shells. Nevertheless, the method provides an inversion formula, as its parameters for sufficient number of shells are easily generated on a computer.

4. Inter-Particle Potentials for FCC and BCC Metals

The cohesive energy formula for metals is more general than in Equation (6), as it is necessary to account for embedding energy [24] due to the presence of free electrons at lattice sites. The generalized form is given by:
E e m b ( x ) = F ( ψ ) + 1 2 n = 1 γ n U ( b n z 0 x ) .
Here, F ( ψ ) is the embedding energy that results due to the electrons, and ψ is the electron density at the central atom. This electron density is expressible as a sum of the density ( ψ a ( r j ) ) contributions from all other atoms, which is, ψ = j 0 ψ a ( r j ) . By employing suitable empirical forms for F ( ψ ) , its contribution from E e m b is readily subtracted out before determining the inter-particle potentials [25]. However, such potentials would be devoid of the binding effects modeled as embedding energy. On the other hand, approximating F ( ψ ) with a Taylor’s expansion (truncated to second order) about an average electron density ψ ¯ , Equation (13) is reducible to the pair-potential form [26], however, with an effective potential given by U e f f ( r ) = U ( r ) + 2 F ( ψ ¯ ) ψ a ( r ) + F ( ψ ¯ ) ψ a 2 ( r ) . Furthermore, it is also found that the radial distribution functions of liquid metals obtained with the effective potentials agree quite well with the simulation results using full potentials [26]. Acordingly, it will be assumed here after, although not indicated explicitly, that all of the potentials derived from cohesive energy are effective potentials.
As an illustration of the model, the effective potential U ( r ) versus nearest-neighbor distance r for Cu (curve-1) is shown in Figure 1A. Table 6 provides the parameters used in the model. The potential minimum U m i n = 0.07143 occurs at r m i n = 1.359 in reduced units. The repulsive (curve-2) and attractive (curves-3) components of the potential, as per WCA specification (see below), are also shown. The cohesive energy E ( x ) versus side length ( x ) of atomic volume is shown in Figure 1B. Keeping all of the parameters except η fixed, variations of the potential minimum U m i n and position r m i n versus η are shown in Figure 2A. These results, which correspond to FCC lattices in general, show that potential parameters are somewhat insensitive to η beyond η 6 . Therefore, critical point parameters, which are to be discussed below, also would show this weak dependence for larger values of η . Figure 2B shows similar results, obtained using data for Fe (see Table 6), which correspond to BCC lattices. The trends of variation of the parameters are similar, although the range of variation of r m i n is much less.
Having discussed the details of obtaining the inter-particle potentials, it is now appropriate to consider a general method for computing properties of metallic fluids. Critical point parameters and phase diagrams in liquid-vapor phase transition are obtained thereafter. The scaling of cohesive energy and potential would imply such a scaling of critical point parameters with respect to η and some (approximate) universal form for phase diagrams.

5. Coupling-Parameter Expansion

The CPE in thermodynamic perturbation theory [11] is based on dividing the inter-particle potential U ( r ) into a reference (repulsive) and perturbation (attractive) components. The most appropriate division is based on the WCA prescription [10], which is defined as
U R ( r ) = U ( r ) U m i n , r r m i n , U R ( r ) = 0 , r > r m i n , U A ( r ) = U m i n , r r m i n , U A ( r ) = U ( r ) , r > r m i n .
The reference part is purely repulsive and it determines the structure of dense liquids, while the attractive perturbing component is responsible for cohesion and phase changes [10]. A coupling-parameter λ is introduced into the potential and it is expressed as U ( r , λ ) = U R ( r ) + λ U A ( r ) . Its magnitude 0 λ 1 determines the strength of the perturbation; with U ( r , 0 ) being the reference part and U ( r , 1 ) the full potential. All of the correlation functions, which determine the structure of the fluid, like the pair distribution function g ( r , λ ) , are now to be treated as functions of λ . In the thermodynamic perturbation theory (TPT), these functions are expanded as Taylor’s series in λ regarding their reference part as
g ( r , λ ) = g 0 ( r ) + λ g 1 ( r ) + + 1 n ! λ n g n ( r ) + .
Here, g 0 ( r ) is the reference system’s pair distribution function and g n ( r ) denotes its n t h derivatives at λ = 0 . Similar Taylor’s series representations hold for other correlation functions as well: the direct correlation function c ( r , λ ) , the total correlation function h ( r , λ ) = g ( r , λ ) 1 , and the indirect correlation function y ( r , λ ) = h ( r , λ ) c ( r , λ ) . For all of them, the subscript 0 will denote the function that corresponds to the reference system. The effect of the perturbation is accurately determined if the derivatives, like g n ( r ) , are obtained up to sufficient orders. Traditional perturbation theory is mainly limited to the first two orders [7]; however, CPE accounts for quite higher orders by taking recourse to the OZE, which defines the relation between the short ranged direct correlation function c ( r , λ ) and the long ranged total correlation function h ( r , λ ) . The latter is the appropriate function to be used in the OZE as it tends to zero for large r, just like the pair-potential. For one-component systems with spherically symmetric potential, as in the present case, the OZE is written as:
h ( r , λ ) = c ( r , λ ) + ρ c ( | r r | , λ ) h ( r , λ ) d r .
where ρ is the number density of atoms. This relation has to be supplemented with a closure relation involving h ( r , λ ) , c ( r , λ ) and the pair-potential U ( r , λ ) . An exact closure is unfeasible; however, several approximate forms are available [3]. All of the closures are generally expressed in terms of a ’bridge function’ B ( r , λ ) in the form:
g ( r , λ ) = E x p [ β U ( r , λ ) + y ( r , λ ) + B ( r , λ ) ] ,
where β = ( k B T ) 1 , k B Boltzmann’s constant and T absolute temperature. Note that k B T is also scaled with energy unit E 0 . Approximate forms of B ( r , λ ) generate different closures; for instance, the one called HNC that was mentioned earlier corresponds to the choice B ( r , λ ) = 0 .
Being an integral equation with displacement kernel, the OZE is easily solved in the Fourier space, leading to the solution
h ¯ ( k , λ ) = c ¯ ( k , λ ) [ 1 ρ c ¯ ( k , λ ) ] 1 .
where the Fourier transform and its inverse, for any function q ( r ) , are defined as
q ¯ ( k ) = 4 π 0 [ sin ( k r ) / ( k r ) ] q ( r ) r 2 d r , q ( r ) = 1 2 π 2 0 [ sin ( k r ) / ( k r ) ] q ¯ ( k ) k 2 d k .
A ’bar’ is used throughout to denote Fourier transforms. It is useful to rewrite these equations in terms of the functions [ r q ( r ) ] and [ k q ¯ ( k ) ] and the symmetric kernel [ sin ( k r ) ] , as discussed below.
The two algebraic equations Equations (17) and (18), although defined in two different spaces, define the non-linear systems that determine the correlation functions and entire thermodynamic properties. When the inter-particle potential contains an attractive component, solutions to the OZE with most of the closures develop singularities in the liquid-vapor transition region [14]. Hence, the determination of the critical points and phase diagrams is a difficult task, and very careful strategies need to be implemented [27]. The CPE technique completely avoids this problem of dealing with singularities, as all correlation functions are expanded in Taylor’s series. The issues that are related to convergence of the series and ’no-solution-region’ are discussed later. The Newton–Armijo non-linear solver considered next is applicable to general potentials; however, it is discussed below for the case of reference potential, as required in the CPE method.

5.1. Newton-Armijo Solver for Reference System

For repulsive inter-particle potentials, as in the case of reference system, an accurate form of bridge function [28] is B 0 = 1 + 2 y 0 y 0 1 , so the closure is expressed as
c 0 ( r ) = E x p [ β U R ( r ) + 1 + 2 y 0 ( r ) 1 ] y 0 ( r ) 1 .
The solution to the OZE in terms of y ¯ 0 ( k ) is given by
y ¯ 0 ( k ) = ρ c ¯ 0 2 ( k ) [ 1 ρ c ¯ 0 ( k ) ] 1 .
These nonlinear equations are readily solved while using Newton’s method [29]. While, in this reference, these are treated as two equations and two unknowns, it is possible to take Equation (20) as providing c 0 if y 0 is given. Accordingly, only Equation (21) is to be solved by Newtons’ method. Starting with initial guess solution y 0 (and, hence, c 0 , c ¯ 0 and y ¯ 0 ), an improved solution y 0 + Δ y is obtained while using the first order correction terms:
Δ c = [ h 0 + g 0 B 0 ] Δ y a n d Δ y ¯ = [ ρ 2 c ¯ 0 + y ¯ 0 1 ρ c ¯ 0 ] Δ c ¯ + [ ρ c ¯ 0 2 1 ρ c ¯ 0 y ¯ 0 ] .
where B 0 is the derivative with respect to y 0 . Taking Fourier transform of Δ c and substituting in the second equation yields
Δ y ¯ = ρ 2 c ¯ 0 + y ¯ 0 1 ρ c ¯ 0 [ ( h 0 + g 0 B 0 ) Δ y ] ¯ + [ ρ c ¯ 0 2 1 ρ c ¯ 0 y ¯ 0 ] .
The linear operator A defining this equation involves the following operations: (i) Fourier inverse taking Δ y ¯ to Δ y , (ii) multiplication by a = h 0 + g 0 B 0 , (iii) Fourier transform of the resulting product, and (iv) multiplication by b ¯ = ρ ( 2 c ¯ 0 + y ¯ 0 ) / ( 1 ρ c ¯ 0 ) . This may be symbolically written as: A Δ y ¯ = Δ y ¯ b ¯ F [ a F 1 [ Δ y ¯ ] ] , where F denotes the Fourier transform. Thus, the action of A essentially involves two Fourier transform operations, which is re-defined with a symmetric kernel as:
Q ¯ ( k ) = 4 π 0 [ sin ( k r ) ] Q ( r ) ] d r , Q ( r ) = 1 2 π 2 0 [ sin ( k r ) ] Q ¯ ( k ) d k .
where Q ( r ) = r q ( r ) and Q ¯ ( k ) = k q ¯ ( k ) . Subsequently, the operator A is re-defined as A [ Δ Y ¯ ] = [ Δ Y ¯ ] b ¯ S [ a S [ Δ Y ¯ ] ] , where S represents integration with the symmetric kernel [ sin ( k r ) ] . The adjoint of A , as needed in Krylov space-based linear equation solvers discussed next, is then easily computed as A [ Δ Y ¯ ] = [ Δ Y ¯ ] S [ a S [ b ¯ Δ Y ¯ ] ] .
Numerical implementation of this algorithm uses a uniformly spaced discrete coordinates { r i = i δ r } , 1 i N m with N m = 2 N u where N u is, typically, 10. This choice allows the use of FFT (Fast Fourier-transform) techniques. With a typical δ r 0.025 , the discrete mesh sufficiently extends to account for the asymptotic variation of the correlation functions. In Fourier space, δ k is chosen, such that δ k δ r = π / N m , so that orthogonality of trigonometric functions is maintained in the discrete space. With these choices, the linear equation for Δ y ¯ is reduced to a matrix equation of order N m , which is most efficiently solved while using Krylov space-based methods [12].
The process of improving the solution, called Newton’s iterations, is continued until the Euclidean norm Δ Y ¯ is less than a prescribed value. If the norms satisfy the condition Δ Y ¯ 1 < ( 1 ϵ ) Δ Y ¯ 0 (subscripts ’1’ and ’0’ denote the current and previous iteration values, and ϵ 10 4 is a small number), next Newton’s iteration follows. Otherwise, it is likely that iteration with full step Δ Y ¯ would diverge. Accordingly, a smaller step size μ Δ Y ¯ with 0 < μ < 1 is to be used. Note that the norms Δ Y ¯ 1 correspond to μ 1 = 1 , while Δ Y ¯ 0 to μ 0 = 0 . Armijo’s rule restricts the new step size to μ Δ Y ¯ with a value of μ 1 / 2 . For choosing that μ , a sub-iteration yielding Δ Y ¯ 2 with μ 2 = 1 / 2 is performed. A quadratic fit of Δ Y ¯ 2 versus μ is then used to obtain μ 3 , which corresponds to the minimum of Δ Y ¯ 2 in the interval [ 0.1 μ 2 , 0.5 μ 2 ] . Another sub-iteration providing Δ Y ¯ 3 is performed and, if it satisfies the norm-reduction criterion, next Newton’s iteration follows. Otherwise, more sub-iterations are performed, each time with the replacements μ 2 μ 1 and μ 3 μ 2 , and, similarly, the respective norms. If the sub-iterations exceed a prescribed limit ( 10), it would be necessary to restart Newton’s iterations with a new guess solution.
The linear equations to be solved, at each Newton’s iteration, is of the standard form A x = b . The Conjugate Gradient (CG) method is the simplest of all iterative methods attempting a solution in the Krylov space, although it is applicable to symmetric matrices [30]. For non-symmetric matrices, as in the present case, the CG method is applicable to the normal equation A A x = A b as A A is symmetric. The algorithm starts with a guess solution x 0 = b and its error vector e 0 = b A x 0 , and auxiliary vectors p 0 = A e 0 . Afterwards, the iteration process, called CGNR (Conjugate Gradient to Normal equations to minimize Residual) (see Appendix A), generates new solutions { x n } , such that the residual error norm e n = b A x n is minimized at each iteration n. The process converges to the exact solution in utmost N m (dimension of the matrix) iterations if A is non-singular [30]. In practice, iterations are terminated when the error norm e n satisfies a prescribed criteria. There is little point in getting a high degree of convergence in the early stages of Newton’s iterations. In what are called inexact methods, the criteria used are e n ξ Δ Y ¯ with a typical value ξ 0.9 . A more common method, called GMRES (Generalized Minimum Residual), although it does not need A , requires much larger storage space [30].
The correlation functions of the reference system are determined with the methods discussed (also see Appendix A) in the entire fluid-phase-plane, except at very low temperatures, where perturbation theories are inapplicable.

5.2. Derivatives of Correlation Functions

The next step in CPE is to compute the derivatives, like g n ( r ) in Equation (15). Writing the general closure relation in Equation (17) in the form g = exp [ f ] , where f = β U + y + B , repeated differentiation with respect to λ yields
g n = g f n + m = 1 n 1 C m n 1 f n m g m , f n = β U A δ n 1 + y n + B n . n 1 .
where { C m n } denote the Binomial coefficients and δ n 1 is the Kroenecker delta. This is valid for arbitrary λ and, in particular, for λ = 0 . The arguments r and λ of all functions are omitted for simplifying the equations. Note that U A is the attractive component in the potential according to the WCA prescription. The summation in Equation (25) only contains derivatives lower than n as the term containing y n is separated out. Substituting for f n and using the relation c n = g n y n for n 1 yields the recursion relation:
c n = [ ( g 1 ) + g Λ ] y n + w n , n 1 w n = g β U A δ n 1 + g B n * + m = 1 n 1 C m n 1 [ β U A δ n m , 1 + ( 1 + Λ ) y n m + B n m * ] [ y m + c m ] , n 1 .
The derivative B n is written as B n = Λ y n + B n * , so that B n * only contains derivatives of y of an order lower than n. This is needed, as y n is yet to be determined with another equation. The definition of structure factor:
s ¯ ( k , λ ) = [ 1 ρ c ¯ ( k , λ ) ] 1 = 1 + ρ h ¯ ( k , λ ) ,
shows that y ¯ n = ρ 1 s ¯ n c ¯ n for n 1 . From now on, the arguments k and λ of all Fourier functions are also omitted. Repeated differentiation of the relation s ¯ = 1 + ρ c ¯ s ¯ yields
s ¯ n = ρ s ¯ n c ¯ + ρ s ¯ c ¯ n + ρ m = 1 n 1 C m n c ¯ n m s ¯ m , n 1 .
where the first and last terms are separately written from the summation. Simplifying the expression and substituting in y ¯ n gives the recursion
y ¯ n = [ s ¯ 2 1 ] c ¯ n + q ¯ n , n 1 , q ¯ n = s ¯ m = 1 n 1 C m n c ¯ n m ρ [ y ¯ m + c ¯ m ] n 1 .
Here, the definition of s ¯ is used and s ¯ m is substituted as ρ ( y ¯ m + c ¯ m ) in the summation. Note that q ¯ n only contains lower order derivatives y ¯ m and c ¯ m for m < n . Taking Fourier transform of Equation (26) and substituting in Equation (29) yields a linear equation for y ¯ n , which is expressed as
y ¯ n = [ s ¯ 2 1 ] [ ( h + g Λ ) y n ] ¯ + [ s ¯ 2 1 ] w ¯ n + q ¯ n , h = g 1 , n 1 .
This expression is valid for arbitrary value of λ . The linear operator in this equation (for all n) is identical to that shown in Equation (23) if the derivatives are evaluated at λ = 0 , which is, for the reference system. Thus, the Krylov space-based methods mentioned above are also applicable here; with the only change being the source term [ s ¯ 0 2 1 ] w ¯ n + q ¯ n .

5.3. Derivatives of Bridge Function

The density-dependent bridge function for a general potential with attractive component is expressed as [28]
B = 1 + 2 ζ ζ 1 , ζ = y λ ρ β U A ,
where all are functions of r and λ . Defining O = B + ζ + 1 = 1 + 2 ζ yields B n = O n ζ n for n 1 where ζ n = y n δ n 1 ρ β U A . Now, repeated differentiation of the relation O O 1 = ζ 1 readily yields the recursion formula:
O n + 1 = O 1 [ y n + 1 m = 0 n 1 C m n O n m O m + 1 ] , n 1
which provides all derivatives B n in terms of y n . The parameter Λ defined earlier is now identified as Λ = O 1 1 , and the summation term provides B n + 1 * .

5.4. Thermodynamic Functions

After computing the Taylor’s series for all functions to an arbitrary high order, setting λ = 1 provides correlation functions that correspond to the full potential. The results that are discussed below are obtained with 7th order approximation, as a further increase in order does not have any significant effect. Thermodynamic properties, like pressure ( P ) and inverse compressibility ( χ 1 ) , are expressed as
β P / ρ = 1 ( 2 / 3 ) ( π β ρ ) 0 g ( r ) U ( r ) r 3 d r ,
χ 1 = β ( P / ρ ) T = 1 4 π ρ 0 c ( r ) r 2 d r ,
Similar expressions hold for internal energy and chemical potential [27]. The argument λ = 1 is omitted from these expressions. For computational purposes, Equation (33) is re-expressed in terms of the ’cavity function’ Y = exp ( β U ) g ( r ) , as:
β P / ρ 1 = ρ B 2 v + ( 2 / 3 ) π ρ 0 [ exp ( β U ) ] ( Y 1 ) r 3 d r .
Here, B 2 v is the second virial coefficient for U, and it is separated out explicitly. An excellent first order perturbation theory [31] for the equation of state follows from this equation with the approximations: U U R and Y Y 0 , thereby terminating the integral at r m i n . The second approximation is the basis of first order theory. The first approximation is motivated from the observation that derivative [ exp ( β U ) ] is only predominant near r m i n , where Y 0 1 is small.
The OZE does not provide exact thermodynamic consistency [3] in the sense that pressure that is computed from different routes does not agree with each other. The inverse compressibility χ 1 , being expressed in terms of the short ranged correlation function c ( r ) , is the more appropriate response function to compute pressure via integration with respect to ρ . The method followed here employs a fine mesh of T and ρ over the phase plane where χ 1 is computed. Subsequently, two-dimensional (2-D) interpolation is used for intermediate values, and P is obtained via integration over ρ from zero. Solutions of the critical conditions v P = 0 and v 2 P = 0 provide the critical point parameters and Maxwell’s construction the phase diagram. Using thermodynamic relations to obtain other quantities, like energy and entropy, equation of state tables [32] can be generated in the expanded volume regions.
There are some important issues underlying the above procedure. It is known that there are multiple solutions to OZE and closure relations in the two-phase region when attractive interactions are present. This is illustrated in Baxter’s analytical solution [33] of the sticky hard sphere model. These include isotherms with negative compressibility, complex, and even diverging pair distribution function. Numerical methods also show similar features for general potentials and closure relations (see following section). Solutions that are given by van der Waals equation, first order perturbation theory, and the CPE show isotherms with negative compressibility in the two-phase region. The nonphysical pressure P ( ρ , T ) that is generated by these approaches is replaced with a flat isotherm while using thermodynamic conditions of equal pressure and chemical potential, or the equivalent Maxwell’s construction. The use of virial pressure (Equation (33)) or inverse compressibility (Equation (34)) produces similar results via Maxwell’s construction, although some differences arise due to the inherent thermodynamic inconsistency mentioned earlier. Another approach is to compute virial pressure and chemical potential, only in regions where compressibility is positive, and then derive the phase diagram while using thermodynamic conditions [27]. The structure factor s ¯ ( k ) within the two-phase region is to be obtained as a weighted average of contributions of the gas and liquid phases that correspond to any point on the isotherm. Only renormalization group approaches, such as in the the hierarchical reference theory [34], generate flat isotherms in the two-phase region without using Maxwell’s construction.
Yet another relevant issue regarding the utility of Taylor’s series expansion of g ( r , λ ) (and other correlation functions) needs clarification. Negative compressibility, at any phase point, implies that s ¯ ( 0 ) 1 is negative and, by continuity requirements, s ¯ ( k 0 ) 1 would vanish at some k 0 , as it should approach unity for large k. Therefore, s ¯ ( k ) and, consequently, h ¯ ( k ) would diverge at k 0 . Thus, negative compressibility at any phase point implies a singularity in the Fourier transform h ¯ ( k ) (The author thanks one of the anonymous reviewers of Condensed Matter (MDPI) journal for this argument). Taylor’s series expansions that are assumed in CPE do not converge to this singular solution, in fact, numerical results indicate that they converge to a different solution (see following section). This is plausible in a scenario, wherein multiple solutions exist in the two-phase region. Their utility in computing phase diagram (via thermodynamic or Maxwell’s construction) needs to be evaluated by comparing the results with those from independent methods.

5.5. ’No-Solution-Region ’

Multiple solutions and singularities of the OZE in the liquid-vapor transition region have been investigated in detail for the Lennard-Jones (LJ) potential: U ( r ) = 4 ϵ [ ( σ / r ) 12 ( σ / r ) 6 ] , with the HNC closure [13] and other closures [14]. The pseudo-arc-length continuation algorithm is needed for the direct application of Newton’s method to the full potential for tracking these singularities [13]. The no-solution-line is defined as the curve in the phase plane, along which the Jacobian is singular and so Newton’s method is inapplicable. Figure 3A displays (curve-1) this curve for the LJ-HNC model [13].
In fact, there are more characteristic features, like the ’fold-bifurcation’ of solutions in the liquid-vapor region of the nonlinear OZE. However, the CPE method does not pick up these singularities, as it is based on series expansion about the repulsive component of the potential. The situation is similar that of van der Waals equation, although free of singularities, which predicts a region of negative compressibility. Figure 3A also shows the spinodal line ( defined by χ 1 = 0 ) for the LJ-HNC model, which is obtained using 7th order CPE of this paper (curve-2).
The LJ system is also investigated while using the density-dependent bridge function [28] for comparison with simulation data [35]. This is important because the bridge function is the sole ingredient that determines the predictive power of OZE. The results for phase diagram (co-existence curve) that were obtained with CPE shown (curve-1) in Figure 3B compare well with simulation data (symbols). There is only slight disagreement along the liquid branch near the transition point. The critical point parameters obtained: ρ c σ 3 = 0.299 [ 0.312 ] , k B T c / ϵ = 1.326 [ 1.316 ] and P c σ 3 / ϵ = 0.102 [ 0.127 ] also compare well with simulation data given in square brackets. The spinodal line (curve-2), which is different from that obtained with HNC, is also shown.
As discussed in the previous section, it is necessary to check wheter the CPE generates any meaningful correlation functions, particularly in the two-phase region. Figure 4A shows the direct correlation function c ( r ) (curve-c), the reference function c 0 ( r ) (curve-0), and six successive derivative terms c n ( r ) / n ! (curve-1 to curve-6) for the LJ system using the density-dependent bridge function [28]. These correspond to the phase point ρ σ 3 = 0.25 and k B T / ϵ = 1.1 , which is in the spinodal region. It is evident that the expansion approaches a well defined function c ( r ) . The reference function c 0 ( r ) is mostly negative, but the positive part in c ( r ) for r > 1 , arising out of attractive interactions, gives rise to negative compressibility. Similarly, Figure 4B shows the pair distribution function g ( r ) (curve-g), the reference function g 0 ( r ) (curve-0), and six derivative terms g n ( r ) / n ! (curve-1 to curve-6), also corresponding to the same phase point. The reference function g 0 ( r ) with a small peak develops to g ( r ) with two well defined peaks due to attractive interactions. It will be more accurate to use the series expansions in Equation (17) (closure relation) for computing g ( r ) . These ’correlation functions’ are nonphysical, as they generate negative compressibility and van der Waals’s loop in pressure, but they may be used in order to compute phase diagram via Maxwell’s construction. Afterwards, it will be necessary to take the weighted averages (as explained earlier for structure factor) of contributions of gas and liquid phases to obtain physically acceptable correlation functions.

6. Liquid-Vapor Phase Transition in Metals

The CPE method and inter-particle potentials that are derived via lattice inversion method for metallic fluids are next evaluated as these have softer repulsive components. The simulation results using Morse potential [36] with parameters corresponding to Au and Cu are used for comparing phase diagrams. The data that are given in Table 6 and the cohesive energy formula in Equation (4) provide the parameter-free potentials for use in CPE.
The results shown (curve-2) in Figure 5A and simulation data [36] (symbols) for Au agree quite well. Critical point parameters of CPE method, viz., ρ c = 5.729 [ 5.925 ] g/cm 3 , T c = 7643 [ 7566 ] K and P c = 0.598 [ 0.525 ] GPa and simulations (given in square brackets) also agree well. The spinodal curve (defined by χ 1 = 0 ) is also shown (curve-1). A similar comparison of phase diagram for Cu is shown (curve-2) in Figure 5B with simulation data [36] (symbols). The critical point parameters ρ c = 2.309 [ 2.631 ] g/cm 3 , T c = 8800 [ 8650 ] K and P c = 0.905 [ 0.954 ] GPa also compare reasonably with the simulation results that are given in square brackets.

Scaling of Critical Parameters and Phase Diagram

Having established the accuracy of methods, the scaling behavior of critical point parameters and phase diagram is now investigated. As the cohesive energy formula in Equation (4) only contains two parameters, ( η and δ ), it is obvious that the effective potentials and, hence, all of the features that are related to co-existence also only depend on these variables. Of these two, δ is generally a small number (see Table 6). Computations utilizing CPE are performed to determine the variation of critical point parameters: ρ c , T c , P c and Z = P c / ( k B T c ρ c ) , with the scaling variable η , and results for FCC lattices are shown in Figure 6A. Note that the graph displays dimensionless quantities that are expressed in the units: V 0 1 , E 0 / k B and E 0 / V 0 , respectively. Similar results for BCC lattices are provided in Figure 6B. The parameter δ is kept constant at 0.0808 for FCC and 0.0148 for BCC lattices. Varying this parameter within its range does not produce significant changes in these curves. For FCC lattices, the variation in ρ c and T c are more pronounced for smaller values of η . These results are somewhat similar to those that were obtained earlier [5]; however, their variations and dependence on the lattice type are new.
Phase diagrams for three FCC metals, Al (curve-1), Cu (curve-2), and Au (curve-3) are displayed in reduced variables T / T c and ρ / ρ c in Figure 7A. Similar results are given in Figure 7B for two BCC metals Fe(curve-1) and W(curve-2). Panels A (curve-4) and B (curve-3) also contain van der Waals phase diagram for comparison. These results indicate that there are similarities in the co-existence region for these metallic fluids, although there is no perfect universality, as in the case of van der Waals model. This is due to the fact that the interaction potentials, scaled with V 0 1 and E 0 , do not assume a universal form. Similar observations are also found in the simulation data of phase diagrams while using Morse potential [36].

7. Summary

The main aim in this paper is to discuss the scaling aspects of liquid–vapor phase transition in metallic fluids. There are two requirements to carry out this investigation: (i) realistic (effective) inter-particle potentials and (ii) an accurate statistical mechanical model. The four-parameter formula for cohesive energy expressed in units E 0 and L 0 consists of just two dimensionless variables ( η and δ ). This is quite accurate in the compressed as well as expanded volume states, although it needs to be corrected for extreme compression while using the QSM for electrons. Accordingly, effective potentials derived from the cohesive energy formula will mainly depend on these variables η and δ . Two complementary ways to impact lattice inversion are considered: (i) direct iterative method and (ii) Chen-Möbius inversion formula. New tables for inversion parameters covering up to 20 atomic shells are generated for FCC and BCC lattices, which are easily extended to other lattice types. Variations of the potential depth and first neighbor distance versus η are then evaluated for the two lattices.
The statistical mechanical model using CPE is revisited and explicit recursion formulas for the derivatives of correlations and bridge function are derived. In fact, very high order derivatives are easily computed with the present formulation. The Newton–Armijo nonlinear solver together with Krylov space-based linear solvers are implemented within the CPE. The same linear solvers are applicable to solve the linear equations for the derivatives of correlation functions, as the matrices involved are the same. With these techniques, it is possible to cover the entire fluid region, except in the low temperature limit, where perturbation theory is inapplicable. The entire procedure is evaluated with simulation data for LJ system as well as metallic fluids. The ’no-solution-region’ concurring when OZE is applied to the full potential is also checked with CPE results for the LJ system. The variations of critical point parameters versus η are then evaluated for FCC and BCC lattices. Further, it is also shown that the phase diagrams of these systems show (approximate) universal behavior when expressed in terms of reduced variables. It is hoped that the CPE elaborated in this paper would be increasingly used to derive thermodynamic properties of one component fluids. To that end, the Appendix A (Algorithms A1–A5) provides short outlines of the main algorithmic components.

Funding

This research received no external funding.

Acknowledgments

The author thanks the reviewers of Condensed Matter Journal for critical reviews and suggestions to improve the presentation of the paper.

Conflicts of Interest

The author declare no conflict of interest.

Appendix A. Algorithm for CPE

The main components of the algorithm to implement CPE is provided below. Newton is the non-liner solver for the reference potential and CGNR is the linear solver. Other linear solvers, like GMRES, can also be used. Matrix computes matrix-vector products in the linear solver. Recursive equations for derivatives of correlation functions are in Derivative. The heart of the algorithm is fast Fourier transforms provided in Fourier. Note that h 0 * and Ψ are the global variables used in Matrix.
T y p i c a l   p a r a m e t e r s N u = 9 N m = 2 N u τ n w = 0.0001 τ i n = 0.0001 m a x N = 15 m a x I = 15 m a x D = 7 d x = . 025 C i = cos [ π ( i 1 ) / N m ] i : 1 , N m + 1 d k = π / ( N m d x )
Algorithm A1: Newton [ y 0 , c 0 ]
in: y 0
out: c 0 , y 0
i = 1 , y 0 = 0
While [ i m a x N ] Do      
g 0 = exp [ β U R + ( 1 + 2 y 0 ) 1 ]       
Λ = 1 / 1 + 2 y 0 1       
h 0 * = g 0 1 + g 0 Λ       
c 0 = g 0 1 y 0       
y ¯ 0 = F o u r i e r [ 1 , y 0 ]       
c ¯ 0 = F o u r i e r [ 1 , c 0 ]       
Ψ = ρ ( y ¯ 0 + 2 c ¯ 0 ) [ 1 ρ c ¯ 0 ] 1       
b ¯ = ρ c ¯ 0 2 [ 1 ρ c ¯ 0 ] 1 y ¯ 0       
C a l l C G N R [ b ¯ , Y ]       
y ¯ 0 y ¯ 0 + Y       
y 0 y 0 + F o u r i e r [ 2 , Y ]       
i i + 1       
if k · Y < τ n w Exit
Algorithm A2: CGNR [ b , x ]
in: b
out: y
x 0 = b
e 0 = b A x 0
p 0 = A
e 0 w 0 = p 0
n = 1
While [ n m a x I ] Do      
   q n = A p n 1       
   a = ( q n · e n 1 ) / ( q n · q n )       
   x n = x n 1 + a p n 1       
   e n = e n 1 a q n       
   w n = A e n       
   α = ( w n · w n ) / ( w n 1 · w n 1 )       
   p n = w n + α p n 1       
   n n + 1       
  if e n < τ i n Exit
Algorithm A3: Matrix [ p , x , y ]
in: x
out: y
( p = 1 : y = A x , p = 2 : y = A x )
If (p .eq. 1)       
   z = h 0 * F o u r i e r [ 2 , x ]       
   y = x Ψ F o u r i e r [ 1 , z ]
If (p .eq. 2)       
   z = Ψ x w = h 0 * F o u r i e r [ 2 , z ] y = x + F o u r i e r [ 1 , w ]
Algorithm A4: Derivative [ m a x D ]
in: m a x D
out: y n , c n
Λ = 1 / 1 + 2 y 0 1
c ¯ 0 = F o u r i e r [ 1 , c 0 ]
S ¯ 0 = ( 1 ρ c ¯ 0 ) 1
Ψ = S ¯ 0 2 1
h 0 * = g 0 1 + g 0 Λ
While [ 1 n m a x D ] Do      
   B n * = Λ β ρ U A δ n 1 ( 1 + Λ ) m = 0 n 2 C m n 1 O m + 1 O n 1 m       
   w n = g β U A δ n 1 + g 0 B n *             
   + m = 1 n 1 C m n 1 [ β U A δ n m , 1 + y n m + B n m ] [ y m + c m ]       
   q ¯ n = ρ S ¯ 0 m = 1 n 1 C m n c ¯ n m [ y ¯ m + c ¯ m ]       
   b ¯ = R F o u r i e r [ 1 , w n ] + q ¯ n       
   C a l l C G N R [ b ¯ , y ¯ n ]       
   y n = F o u r i e r [ 2 , y ¯ n ]       
   c n = w n + h 0 y n       
   c ¯ n = F o u r i e r [ 1 , c n ]       
   O n = ( 1 + Λ ) y n + B n *       
   B n B n * + Λ y n       
   n n + 1
Algorithm A5: Fourier [ p , x , y ]
in: p
out: y
c 1 = 4 π d x
c 2 = d k / ( 2 π 2 )
Do [ i : 1 , N m 1 ] :       
   p = 1 : x i c 1 r i x i       
   p = 2 : x i c 2 k i x i
l = 1
Do [ m : 1 , N u 1 ] :       
   m 1 = N m / l       
   m 2 = m 1 / 2       
   m 5 = 0       
   m 6 = m 2       
  Do [ l 1 : 1 , l ] :            
    i i = 1             
    s u m = 0             
   Do [ i : 1 , m 1 1 , 2 ] :                  
     s u m = s u m + i i × x m 5 + i                   
     i i = i i                   
     y m 6 = s u m
   Do [ i : 1 , m 2 1 ] :
     i 2 = m 5 + 2 i
     y m 5 + i = x i 2 1 + x i 2 + 1
     y m 6 + i = x i 2
   Do [ i : 1 , m 2 1 ] :
     x m 5 + i = y m 5 + i
     x m 6 + i = y m 6 + i             
    x m 6 = y m 6             
    m 5 = m 5 + m 1             
    m 6 = m 6 + m 1       
   l = l
l = 2 N u 2
Do [ m : 1 , N u 1 ] :      
   m 1 = N m / l       
   m 2 = m 1 / 2       
   m 5 = 0       
   m 6 = m 2       
  Do [ l 1 : 1 , l ] :
    m 7 = m 6 + m 2             
    i c = l             
   Do [ i : 1 , m 2 1 ] :                  
     [ t 1 = x m 5 + i / ( 2 C i c + 1 )                   
     t 2 = x m 6 + i                   
     y m 5 + i = t 1 + t 2                   
     y m 7 i = t 1 t 2                   
     i c = i c + l
   Do [ i : 1 , m 2 1 ] :           x m 5 + i = y m 5 + i
       x m 7 i = y m 7 i             
      m 5 = m 5 + m 1             
      m 6 = m 6 + m 1       
   l = l / 2
Do [ i : 1 , N m 1 ] :       
   p = 1 : y i y i / k i       
   p = 2 : y i y i / r i

References

  1. Landau, L.D.; Lifshitz, E.M. Statistical Physics; Pergamon Press: London, UK, 1959. [Google Scholar]
  2. Menon, S.V.G. Renormalization Group Theory of Critical Phenomena; Wiley Eastern Ltd.: Bengaluru, India, 1995; Available online: https://www.researchgate.net/publication/260165199 (accessed on 23 January 2021).
  3. Hansen, J.P.; McDonald, I.R. Theory of Simple Liquids, 4th ed.; Academic Press: Cambridge, MA, USA, 2013. [Google Scholar]
  4. Rose, J.H.; Smith, J.R.; Guinea, F.; Ferrante, F. Universal features of the equation of state of metals. Phys. Rev. B 1984, 29, 2963. [Google Scholar] [CrossRef]
  5. Rosenfeld, Y. Predicting the liquid-vapor critical point from the crystal anharmonicity. Contrib. Plasma Phys. 2001, 41, 183. [Google Scholar] [CrossRef]
  6. Kerley, G.I. Perturbation theory and the thermodynamic properties of fluids. II. J. Chem. Phys. 1980, 73, 478. [Google Scholar] [CrossRef]
  7. Solana, J.R. Perturbation Theories for Thermodynamic Properties of Fluids and Solids; CRC Press: New York, NY, USA, 2013. [Google Scholar]
  8. Cowen, B.J.; Carpenter, J.H. Improved reference system for corrected rigid spheres equation of state model. J. Appl. Phys. 2020, 128, 055901. [Google Scholar] [CrossRef]
  9. Li, J.H.; Liang, S.H.; Guo, H.B.; Liu, B.X. Four-parameter equation of state and determination of the thermal and mechanical properties of metals. J. Alloys Compd. 2007, 431, 23. [Google Scholar] [CrossRef]
  10. Weeks, J.D.; Chandler, D.; Andersen, H.C. Role of Repulsive Forces in Determining the Equilibrium Structure of Simple Liquids. J. Chem. Phys. 1971, 54, 5237. [Google Scholar] [CrossRef]
  11. Ramana, A.S.V.; Menon, S.V.G. Coupling-parameter expansion in thermodynamic perturbation theory. Phys. Rev. E 2013, 87, 022101. [Google Scholar] [CrossRef]
  12. Kelley, C.T. Solving Nonlinear Equations with Newton’s Method, Fundamentals of Algorithms; SIAM: Philadelphia, PA, USA, 2003. [Google Scholar]
  13. Peplow, A.T.; Beardmore, R.E.; Bresme, F. Algorithms for the computation of solutions of Ornstein-Zernike Equation. Phys. Rev. E 2006, 74, 046705. [Google Scholar] [CrossRef]
  14. Tikhonov, D.A.; Sarkisov, G.N. Singularities of solution of the Ornstein-Zernike Equation within the gas-liquid transition region. Russ. J. Phys. Chem. 2000, 74, 470. [Google Scholar]
  15. Gillan, M.J.; Alfe, D.; Brodholt, J.; Vocadlo, L.; Price, G.D. First-principles modeling of Earth and planetary materials at high pressure and temperatures. Rep. Prog. Phys. 2006, 69, 2365. [Google Scholar] [CrossRef]
  16. Kalitkin, N.N.; Kuz’mina, L.V. Curves of cold compression at high pressures. Sov. Phys. Solid State 1972, 13, 1938. [Google Scholar]
  17. More, R.M. Quantum-statistical model for high-density matter. Phys. Rev. A 1979, 19, 1234. [Google Scholar] [CrossRef]
  18. Kerley, G.I. User’s Manual for PANDA II- A Computer Code for Calculating Equation of State; Sandia Report, SAND88-229.UC-405; Sandia: Albuquerque, NM, USA, 1991. [Google Scholar]
  19. Heltemes, T.A.; Moses, G.A. BADGER v 1.0: A Fortran equation of state library. Comput. Phys. Commun. 2012, 183, 2629. [Google Scholar] [CrossRef]
  20. Carlsson, A.E.; Gelatt, C.D., Jr.; Ehrenreich, H. An ab initio pair potential applied to metals. Philos. Mag. A 1980, 41, 241. [Google Scholar] [CrossRef]
  21. Bazant, M.Z.; Kaxiras, E. Modeling of covalent bonding in solids by inversion of cohesive energy curves. Phys. Rev. Lett. 1996, 77, 4370. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. Wiley, J.D.; Seman, J.A. The enumeration of neighbors on cubic and hexagonal-based lattices. Bell Syst. Tech. J. 1970, 49, 355. [Google Scholar] [CrossRef]
  23. Nan-Xian, C.; Zhao-Dou, C.; Yu-Chuan, W. Multidimensional inverse lattice problem and a uniformly sampled arithmetic Fourier transform. Phys. Rev. E 1997, 55, R5. [Google Scholar] [CrossRef] [Green Version]
  24. Daw, M.S.; Baskes, M.I. Embedded-atom method: Derivation and application to impurities, surfaces and other defects in metals. Phys. Rev. B 1984, 29, 6443. [Google Scholar] [CrossRef] [Green Version]
  25. Xie, Q.; Zhang, W.; Chen, N. Analytical long-range embedded-atom potentials. arXiv 1997, arXiv:9610183v2. [Google Scholar]
  26. Foiles, S.M. Application of the embedded-atom method to liquid transition metals. Phys. Rev. B 1985, 32, 3409. [Google Scholar] [CrossRef]
  27. Charpentier, I.; Jakse, N. Phase diagram of complex fluids using an efficient integral equation method. J. Chem. Phys. 2005, 123, 204910. [Google Scholar] [CrossRef] [PubMed]
  28. Sarkisov, G.N. Structure of simple fluid in the vicinity of the critical point: Approximate integral equation theory of liquids. J. Chem. Phys. 2003, 119, 373. [Google Scholar] [CrossRef]
  29. Zerah, G. An efficient Newton’s method for numerical solution of fluid integral equations. J. Comput. Phys. 1985, 61, 280. [Google Scholar] [CrossRef]
  30. Kelley, C.T. Iterative Methods for Solving Linear and Nonlinear Equations; SIAM: Philadelphia, PA, USA, 1995. [Google Scholar]
  31. Song, Y.; Mason, E.A. Statistical-mechanical theory of a new analytical equation of state. J. Chem. Phys. 1989, 91, 7840. [Google Scholar] [CrossRef]
  32. Menon, S.V.G.; Nayak, B. An equation of state for metals at high temperature and pressure in compressed and expanded volume regions. Condens. Matter 2019, 4, 71. [Google Scholar] [CrossRef] [Green Version]
  33. Baxter, R.J. Percus-Yevick equation of hard spheres with surface adhesion. J. Chem. Phys. 1968, 49, 2770. [Google Scholar] [CrossRef]
  34. Parola, A.; Pini, D.; Reatto, L. The smooth cut-off hierarchical reference theory of fluids. Mol. Phys. 2009, 107, 503. [Google Scholar] [CrossRef] [Green Version]
  35. Johnson, J.K.; Zollweg, J.A.; Gubbins, K.E. The Lennard-Jones equation of state revisited. Mol. Phys. 1993, 78, 591. [Google Scholar] [CrossRef]
  36. Singh, J.K.; Adhikari, J.; Kwak, S.K. Vapor-liquid phase co-existence curves for Morse fluids. Fluid Phase Equilib. 2006, 248, 1. [Google Scholar] [CrossRef]
Figure 1. (A) Effective potential U ( r ) for Cu versus scaled nearest-neighbor distance (curve-1, with repulsive (curve-2) and attractive components (curve-3) as per WCA prescription. Potential minimum is U m i n = 0.07143 at r m i n = 1.359 in reduced units. (B) Cohesive energy E ( x ) for Cu versus scaled side length of atomic volume x. Energy minimum ( 1 ) occurs at (1) in reduced units. Length and energy units are L 0 = V 0 1 / 3 and E 0 (see text), respectively.
Figure 1. (A) Effective potential U ( r ) for Cu versus scaled nearest-neighbor distance (curve-1, with repulsive (curve-2) and attractive components (curve-3) as per WCA prescription. Potential minimum is U m i n = 0.07143 at r m i n = 1.359 in reduced units. (B) Cohesive energy E ( x ) for Cu versus scaled side length of atomic volume x. Energy minimum ( 1 ) occurs at (1) in reduced units. Length and energy units are L 0 = V 0 1 / 3 and E 0 (see text), respectively.
Condensedmatter 06 00006 g001
Figure 2. (A) Potential minimum U m i n and its position r m i n versus η for FCC lattices. All other parameters correspond to those of Cu (see text). (B) Similar results for BCC lattices. All of the parameters, except η , correspond to those of Fe (see text).
Figure 2. (A) Potential minimum U m i n and its position r m i n versus η for FCC lattices. All other parameters correspond to those of Cu (see text). (B) Similar results for BCC lattices. All of the parameters, except η , correspond to those of Fe (see text).
Condensedmatter 06 00006 g002
Figure 3. (A) The ’no-solution-line’ (curve-1) for Lennard–Jones potential and hypernetted chain (HNC) obtained with pseudo-arc-length continuation algorithm [13]. Newton’s method applied to the full potential diverges on this curve. The spinodal line (curve-2) (defined by χ 1 = 0 ) obtained using the CPE of this paper lies inside the no-solution region. (B) Comparison of co-existence line (curve-1), obtained with CPE and density-dependent bridge function [28], with simulation results (symbols) [35]. The spinodal line (curve-2), which is different from that obtained with HNC, is also shown.
Figure 3. (A) The ’no-solution-line’ (curve-1) for Lennard–Jones potential and hypernetted chain (HNC) obtained with pseudo-arc-length continuation algorithm [13]. Newton’s method applied to the full potential diverges on this curve. The spinodal line (curve-2) (defined by χ 1 = 0 ) obtained using the CPE of this paper lies inside the no-solution region. (B) Comparison of co-existence line (curve-1), obtained with CPE and density-dependent bridge function [28], with simulation results (symbols) [35]. The spinodal line (curve-2), which is different from that obtained with HNC, is also shown.
Condensedmatter 06 00006 g003
Figure 4. Convergence of coupling-parameter expansion (CPE) for the direct correlation function c ( r ) and pair distribution function g ( r ) for the LJ system (obtained with density-dependent bridge function [28]) at the phase point ρ σ 3 = 0.25 and k B T / ϵ = 1.1 , which is in the spinodal region. (A) Direct correlation function c ( r ) (curve-c), reference function c 0 ( r ) (curve-0) and successive derivative terms c n ( r ) / n ! (curve-1 to curve-6) are shown. (B) Pair distribution function g ( r ) (curve-g), reference function g 0 ( r ) (curve-0) and derivative terms g n ( r ) / n ! (curve-1 to curve-6) are shown.
Figure 4. Convergence of coupling-parameter expansion (CPE) for the direct correlation function c ( r ) and pair distribution function g ( r ) for the LJ system (obtained with density-dependent bridge function [28]) at the phase point ρ σ 3 = 0.25 and k B T / ϵ = 1.1 , which is in the spinodal region. (A) Direct correlation function c ( r ) (curve-c), reference function c 0 ( r ) (curve-0) and successive derivative terms c n ( r ) / n ! (curve-1 to curve-6) are shown. (B) Pair distribution function g ( r ) (curve-g), reference function g 0 ( r ) (curve-0) and derivative terms g n ( r ) / n ! (curve-1 to curve-6) are shown.
Condensedmatter 06 00006 g004
Figure 5. Comparison of phase diagram obtained with 7th order CPE and density-dependent bridge function [28] against simulation data for soft repulsive potentials occurring in metallic fluids: (A) The CPE phase diagram for Au (curve-2) and simulation results (symbols) using Morse potential [36] (see text). The spinodal curve (defined by χ 1 = 0 ) that shows the boundary of meta-stable states (curve-1) is also shown. (B) The CPE phase diagram for Cu (curve-2) is compared against simulation results (symbols) using Morse potential [36] (see text). The spinodal curve (curve-1) is also shown.
Figure 5. Comparison of phase diagram obtained with 7th order CPE and density-dependent bridge function [28] against simulation data for soft repulsive potentials occurring in metallic fluids: (A) The CPE phase diagram for Au (curve-2) and simulation results (symbols) using Morse potential [36] (see text). The spinodal curve (defined by χ 1 = 0 ) that shows the boundary of meta-stable states (curve-1) is also shown. (B) The CPE phase diagram for Cu (curve-2) is compared against simulation results (symbols) using Morse potential [36] (see text). The spinodal curve (curve-1) is also shown.
Condensedmatter 06 00006 g005
Figure 6. Variation of critical point parameters: ρ c , T c , P c and Z = P c / ( k B T c ρ c ) , with the scaling variable η (see text). Note that these quantities are expressed in the units: V 0 1 , E 0 / k B and E 0 / V 0 , respectively. (A) The results for FCC lattices. (B) Results for BCC lattices.
Figure 6. Variation of critical point parameters: ρ c , T c , P c and Z = P c / ( k B T c ρ c ) , with the scaling variable η (see text). Note that these quantities are expressed in the units: V 0 1 , E 0 / k B and E 0 / V 0 , respectively. (A) The results for FCC lattices. (B) Results for BCC lattices.
Condensedmatter 06 00006 g006
Figure 7. Phase diagram of metals in reduced variables T / T c and ρ / ρ c . (A) FCC metals Al (curve-1), Cu (curve-2) and Au (curve-3). van der Waals phase diagram is also shown (curve-4). (B) Similar results for BCC metals Fe (curve-1) and W (curve-2). van der Waals phase diagram is shown (curve-3) for comparison.
Figure 7. Phase diagram of metals in reduced variables T / T c and ρ / ρ c . (A) FCC metals Al (curve-1), Cu (curve-2) and Au (curve-3). van der Waals phase diagram is also shown (curve-4). (B) Similar results for BCC metals Fe (curve-1) and W (curve-2). van der Waals phase diagram is shown (curve-3) for comparison.
Condensedmatter 06 00006 g007
Table 1. Shell radii and occupation numbers for FCC lattice. (see Equation (6)).
Table 1. Shell radii and occupation numbers for FCC lattice. (see Equation (6)).
n γ n b n n γ n b n n γ n b n n γ n b n
112168 6 1124 11 1648 17
26 2 748 7 1224 12 1730 18
324 3 86 8 1372 13 1872 19
412293631448 15 1924 20
524 5 1024 10 151242048 21
Note that b 14 = 15 . There is no shell of radius 14 .
Table 2. Shell radii and occupation numbers for BCC lattice. (see Equation (6))
Table 2. Shell radii and occupation numbers for BCC lattice. (see Equation (6))
n γ n b n n γ n b n n γ n b n n γ n b n
18166 4 / 3 1112 4 2 / 3 1624 2 11 / 3
26 2 / 3 724 19 / 3 1248 35 / 3 1784
312 2 2 / 3 824 2 5 / 3 1330 2 3 1848 17
424 11 / 3 924 2 2 1424 2 10 / 3 1924 2 13 / 3
582103231524 43 / 3 2048 2 14 / 3
Note that there is no formula for radii in BCC lattice.
Table 3. Constants of Chen–Möbius inversion formula for FCC lattice.
Table 3. Constants of Chen–Möbius inversion formula for FCC lattice.
n Γ n B n I n n Γ n B n I n n Γ n B n I n n Γ n B n I n
112 1 1/1268 6 1/91124 11 1 / 6 1624 16 1 / 64
26 2 1 / 24 748 7 1 / 3 1224 12 7/721724 17 1 / 3
324 3 1 / 6 86 8 1/321372 13 1 / 2 1824 18 17 / 72
412 4 1 / 16 936 9 1/12140 14 1/31924 19 1 / 2
524 5 1 / 6 1024 10 01548 15 1/32024 20 5/24
Note that B 14 = 14 and Γ 14 = 0 are introduced.
Table 4. Constants of Chen–Möbius inversion formula for BCC lattice.
Table 4. Constants of Chen–Möbius inversion formula for BCC lattice.
n Γ n B n I n n Γ n B n I n n Γ n B n I n
1810.1250060 1.77778 0.03955110 2.17732 0.31641
26 1.15470 0.0975 70 1.88562 0.28125120 2.21108 0.5625
30 1.33333 0.07031824 1.91485 0.375 136 2.30940 0.09375
40 1.53960 0.05273 98 2.0 0.125 140 2.37037 0.02225
512 1.63299 0.1875 100 2.05280 0.02966 150 2.51416 0.31641
Note that B 3 = 1.33333 and Γ 3 = 0 , and many others are introduced.
Table 5. Constants of Chen–Möbius inversion formula for BCC lattice—continued.
Table 5. Constants of Chen–Möbius inversion formula for BCC lattice—continued.
n Γ n B n I n n Γ n B n I n n Γ n B n I n
1624 2.51661 0.375 2124 2.82843 0.375 2632 3.0 0.5
170 2.55314 0.63281 220 2.9031 0.29663 270 3.0792 0.58008
1824 2.58199 0.375 230 2.90593 0.5625280 3.16049 0.01251
190 2.66667 0.21094240 2.94811 0.632812912 3.26599 0.75
200 2.73707 0.01669 250 2.98142 0.5625300 3.35221 0.26697
Note that B 17 = 2.55314 and Γ 17 = 0 , and many others are new.
Table 6. Material Parameters for four-parameter model .
Table 6. Material Parameters for four-parameter model .
Material V 0 B 0 B 0 E 0 η δ L 0 At. No.At. Wt.
A 3 GPa eV A Z n M n
Copper11.38134.85.193.4895.06190.08082.27712963.50
Aluminum16.3579.34.373.3894.63610.02962.53821326.98
Gold16.95180.75.433.8126.7182 0.0034 2.568679196.97
Iron11.81163.04.504.2815.02700.01482.27752655.85
Tungsten15.82325.04.368.7915.732 0.0402 2.277574183.85
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Menon, S.V.G. Scaling of Phase Diagram and Critical Point Parameters in Liquid-Vapour Phase Transition of Metallic Fluids. Condens. Matter 2021, 6, 6. https://0-doi-org.brum.beds.ac.uk/10.3390/condmat6010006

AMA Style

Menon SVG. Scaling of Phase Diagram and Critical Point Parameters in Liquid-Vapour Phase Transition of Metallic Fluids. Condensed Matter. 2021; 6(1):6. https://0-doi-org.brum.beds.ac.uk/10.3390/condmat6010006

Chicago/Turabian Style

Menon, S.V.G. 2021. "Scaling of Phase Diagram and Critical Point Parameters in Liquid-Vapour Phase Transition of Metallic Fluids" Condensed Matter 6, no. 1: 6. https://0-doi-org.brum.beds.ac.uk/10.3390/condmat6010006

Article Metrics

Back to TopTop