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

: The ﬁrst 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 ﬂuid 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 ﬂuids.


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 energyvolume 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 couplingparameter 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.

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: 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 ∼ 2V 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 1/3 0 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 0K, 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: Here,h 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 = −dE/dV 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: ] 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].

Lattice Inversion for Potential
The lattice inversion method [20] is a convenient tool for extracting effective interparticle potentials from E 0 (V ), which is re-expressed in dimensionless functional form as E(x). More specifically, 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): 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) On separating the first shell-contribution, Equation (6) is expressed as: 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 Tables 1 and 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: 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: where {I m } are weight factors for different shells. Substituting this expansion in Equation (8) gives: 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: 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: 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 , · · · , for 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 semigroup property. This procedure is not unique and several virtual shells are introduced, as seen in Tables 4 and 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.

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: 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 emb 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 densitȳ ψ, Equation (13) is reducible to the pair-potential form [26], however, with an effective potential given by U e f f (r) = U(r) + 2F (ψ)ψ a (r) + F (ψ)ψ 2 a (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 min = −0.07143 occurs at r min = 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 min and position r min 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 min 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.

Coupling-Parameter Expansion
The CPE in thermodynamic perturbation theory [11] is based on dividing the interparticle 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 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 Here, g 0 (r) is the reference system's pair distribution function and g n (r) denotes its n th 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: 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: 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 where the Fourier transform and its inverse, for any function q(r), are defined as A 'bar' is used throughout to denote Fourier transforms. It is useful to rewrite these equations in terms of the functions [rq(r)] and [kq(k)] and the symmetric kernel [sin(kr)], 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.

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 + 2y 0 − y0 − 1, so the closure is expressed as The solution to the OZE in terms ofȳ 0 (k) is given bȳ 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: where B 0 is the derivative with respect to y 0 . Taking Fourier transform of ∆c and substituting in the second equation yields The linear operator A defining this equation involves the following operations: (i) Fourier inverse taking ∆ȳ to ∆y, (ii) multiplication by a = h 0 + g 0 B 0 , (iii) Fourier transform of the resulting product, and (iv) multiplication byb = ρ(2c 0 +ȳ 0 )/(1 − ρc 0 ). This may be symbolically written as: A ∆ȳ = ∆ȳ −bF [aF −1 [∆ȳ]], 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: where Q(r) = rq(r) andQ(k) = kq(k  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 ∆ȳ 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 ∆Ȳ is less than a prescribed value. If the norms satisfy the condition ∆Ȳ 1 < (1 − ) ∆Ȳ 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 ∆Ȳ would diverge. Accordingly, a smaller step size µ∆Ȳ with 0 < µ < 1 is to be used. Note that the norms ∆Ȳ 1 correspond to µ 1 = 1, while ∆Ȳ 0 to µ 0 = 0. Armijo's rule restricts the new step size to µ∆Ȳ with a value of µ ≤ 1/2. For choosing that µ, a sub-iteration yielding ∆Ȳ 2 with µ 2 = 1/2 is performed. A quadratic fit of ∆Ȳ 2 versus µ is then used to obtain µ 3 , which corresponds to the minimum of ∆Ȳ 2 in the interval [0.1µ 2 , 0.5µ 2 ]. Another sub-iteration providing ∆Ȳ 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 ≤ ξ ∆Ȳ 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.

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) where {C n m } 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: 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: shows thatȳ n = ρ −1s n −c n for n ≥ 1. From now on, the arguments k and λ of all Fourier functions are also omitted. Repeated differentiation of the relations = 1 + ρcs yields where the first and last terms are separately written from the summation. Simplifying the expression and substituting inȳ n gives the recursion y n = [s 2 − 1]c n +q n , n ≥ 1,q n =s Here, the definition ofs is used ands m is substituted as ρ(ȳ m +c m ) in the summation. Note thatq n only contains lower order derivativesȳ m andc m for m < n. Taking Fourier transform of Equation (26) and substituting in Equation (29) yields a linear equation forȳ n , which is expressed as 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 .

Derivatives of Bridge Function
The density-dependent bridge function for a general potential with attractive component is expressed as [28] 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: 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 .

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 7 th 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 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: Here, B 2v 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 min . 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 min , 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 ∂ 2 v 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 factors(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 thats(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 transformh(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.

'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]. 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.
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. 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.

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]  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.

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 −1 0 , 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 −1 0 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].

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.