Effective Pair Interactions and Structure in Liquid Noble Metals within Wills-Harrison and Bretonnet-Silbert Models

Abstract

Recently, for calculating the effective pair interactions in liquid transition metals, we have developed an approach which includes the Wills-Harrison and Bretonnet-Silbert models as limit cases. Here, we apply this approach to noble liquid metals. The dependencies of pair potentials and corresponding MD-simulated pair correlation functions in pure liquid Cu, Ag and Au on the portion of the non-diagonal (with respect to the magnet quantum number) d-d-electron couplings in the metal under consideration are studied. The model provides a good agreement with experimental and ab initio data for pair correlation functions, structure factors and velocity autocorrelation functions.

1. Introduction

At present, theoretical investigations of liquid transition metals and their alloys are being carried out intensively [1,2,3,4,5,6,7,8,9,10,11,12]. The use of ab initio simulations as well as accurate many-body potentials and force fields is one of the trends in this field [13,14]. On the other hand, the effective pair potentials providing general relations between the interatomic interactions, structure, properties and phase transitions can explain the features of the surface energy behavior, the formation of complex crystal and quasicrystal phases [15,16,17,18] and are relevant in computer simulations which require high computational efficiency [19,20]. Therefore, the development of effective pair interactions for metals and their alloys is still an actual task.

Among different approaches for transition-metals’ melts, the Wills-Harrison (WH) approximation [21] and its various modifications are well suited for studying the thermodynamics of the aforementioned substances [22,23,24].

The WH model is based on the Harrison-Froyen (HF) approximation [25] which introduces some elements of the muffin-tin orbital theory [26] into transition-metal pseudopotential theory [27,28]. Wills and Harrison used the HF approximation in conjunction with the simple-metal pseudopotential theory to describe the internal energy and effective pair interaction in transition metals. The WH model assumes two contributions into the pair interatomic interaction, due to the valence s electrons and the external d electrons, respectively. The first contribution is calculated in the framework of the nearly-free-electron (NFE) approximation [29] via the pseudopotential method. For this aim, Wills and Harrison used the simplest local model pseudopotential suggested by Ashcroft [30]. In another work [31], we proposed to use the Bretonnet-Silbert (BS) local model pseudopotential [32] instead of the Ashcroft one. Note that the BS pseudopotential taking into account the external d electrons in the metal by means of the replacement of the s-electron valence by the effective s-electron valence was developed especially for the transition metals and sometimes still is used for this purpose independently on the WH approach [33].

In works [21,25], only the diagonal matrix elements with respect to the magnet quantum number, m, between d states of neighboring atoms are taken into account. It is correct in the case of the rotational symmetry with axis along the interatomic distance. Since the d-electron contribution to the potential energy is non-pairwise, the mentioned above symmetry should be disturbed and non-diagonal d-d couplings may arise. Recently, we suggested the modification of the WH model taking into account these non-diagonal d-d-electron couplings between neighboring atoms [34].

In [35], it was observed that in the case when all possible non-diagonal couplings are taken into account, the d-electron contribution to the WH effective pair potential vanishes and latter becomes equal to the NFE (i.e., in our case, to the BS) contribution only. Thus, our modification of the WH model has two limit cases—the original WH model and the BS model.

Earlier, we studied the behavior of the pair potential, $\phi (r)$, between these limit cases in liquid Fe, Co, Ni and Fe-Co, Fe-Ni and Co-Ni melts [34,35,36]. Here, the same investigation is performed for the liquid noble metals: Cu, Ag and Au. Additionally, we analyze the pair correlation functions (or “radial distribution functions”—RDFs), $g(r)$, obtained from these potentials by means of molecular dynamics (MD) simulations.

2. Theory

2.1. Bretonnet-Silbert Model

In the Bretonnet-Silbert model, the ion–s-electron potential for the unscreened ion is the following [32]:

$${\omega }_{\mathrm{BS}}(r)=\{\begin{array}{ll}{\displaystyle \sum _{n=1}^2{B}_n\exp \left(\frac{r}{na}\right)},& r\le R_{\mathrm{C}}\hfill \\ -z_s/r,\hfill & r\ge R_{\mathrm{C}}\hfill \end{array}$$

(1)

where the parameters $R_{\mathrm{C}}$ and $a$ are the ion radius and the softness parameter, respectively; $z_s$ is the effectives-electron valence; and $B_1$ and $B_2$ are the coefficients which have to be chosen so as to make the potential and its first derivative continuous at r = $R_{\mathrm{C}}$:

$$B_n=(n{z}_s/R_{\mathrm{C}}){(-1)}^n\left(\frac{2a-n{R}_{\mathrm{C}}}{n{R}_{\mathrm{C}}}\right)\exp \left(\frac{R_{\mathrm{C}}}{na}\right)$$

(2)

The form factor of ${\omega }_{\mathrm{BS}}(r)$ is equal to

$${\omega }_{\mathrm{BS}}(q)=4\pi \rho a^3\left[\frac{B_1{J}_1(q)}{{(1+a_{}^2{q}^2)}^2}+\frac{8B_2{J}_2(q)}{{(1+4a_{}^2{q}^2)}^2}\right]-(4\pi \rho z_s/q^2)\cos (q{R}_{\mathrm{C}})$$

(3)

where $\rho$ is the mean atomic density;

$$J_n(q)=2-\exp \left(-\frac{R_C}{na}\right)\left\{\left[R_{\mathrm{C}}\frac{1+n^2{a}^2{q}^2}{na}+1-n^2{a}^2{q}^2\right]\frac{\sin (q{R}_{\mathrm{C}})}{naq}+\left[2+\frac{R_{\mathrm{C}}(1+n^2{a}^2{q}^2)}{na}\right]\cos (q{R}_{\mathrm{C}})\right\}$$

(4)

Hereafter, all quantities are given in atomic units (a.u.).

The NFE pair potential is expressed as follows:

$${\phi }_{\mathrm{NFE}}(r)=\frac{z^2}{r}+\frac{1}{8{\pi }^2\rho }{\displaystyle \underset{0}{\overset{\infty }{\int }}F(q)}\frac{\sin (qr)}{qr}{q}^2\mathrm{d}q$$

(5)

where $F(q)$ is the energy-wavenumber characteristic:

$$F(q)=-\frac{q^2{\omega }^2(q)}{\pi \rho [{(\epsilon (q)-1)}^{-1}+1-f(q)]}$$

(6)

Here $\epsilon (q)$ is the Hartree dielectric function, and $f(q)$ is the exchange-correlation function, which is calculated in the approximation of Vashishta and Singwi [37].

2.2. Wills-Harrison Model

The effective pair potential within the Wills-Harrison model consists of two parts [21]:

$${\phi }_{\mathrm{WH}}(r)={\phi }_{\mathrm{NFE}}(r)+{\phi }_d(r)$$

(7)

where:

$${\phi }_d(r)={\phi }_b(r)+{\phi }_c(r)$$

(8)

Here, ${\phi }_b(r)$ is the contribution due to the existence of the d-band width, and ${\phi }_c(r)$ is the contribution due to the shift of the d-band center:

$${\phi }_b(r)=-z_d\left(\frac{10-z_d}{10}\right){\left(\frac{12}{\nu }\right)}^{1/2}\frac{r_d{}^3}{r^5}{K}_b$$

(9)

$${\phi }_c(r)=z_d\frac{r_d{}^6}{r^8}{K}_c$$

(10)

where $z_d=z-z_s$ is the effective d-electron valence; $z$ is the total valence; $r_d$ is the d-state radius; $\nu$ is the coordination number; and $K_b$ and $K_c$ are the coefficients taking into account the d-d-electron couplings between neighboring atoms. In [21], only the couplings diagonal with respect to the magnet quantum number, m, are considered (Figure 1). In this case, the expressions for these characteristics, at the orbital quantum number, $l$, equal to 2, are the following [25]:

$$K_b={\left({\displaystyle \sum _{m=-2}^2\frac{y_m^2}{5}}\right)}^{\frac{1}{2}}$$

(11)

$$K_c=-2{\displaystyle \sum _{m=-2}^2\frac{y_m{x}_m}{5}}$$

(12)

where $y_m$ and $x_m$ are the combinatoric coefficients:

$$y_m=y_{\left|m\right|}=-\frac{{(-1)}^{\left|m\right|}180}{\pi (2+\left|m\right|)!(2-\left|m\right|)!}$$

(13)

$$x_m=x_{\left|m\right|}=-\frac{1}{8}\left(1+\frac{4m^2-1}{9}\right)y_m$$

(14)

From Equations (11)–(14), one can get:

$$K_b=\frac{28.06}{\pi }$$

(15)

$$K_c=\frac{225}{{\pi }^2}$$

(16)

2.3. Modification of the Wills-Harrison Model

If all 25 possible d-d-electron couplings are realized between two atoms (Figure 2), then the expressions for $K_b$ and $K_c$ are rewritten as follows [34]:

$$K_b=\frac{1}{5}{\left(y_0^2+4y_1^2+4y_2^2+4y_0{y}_1+4y_0{y}_2+8y_1{y}_2\right)}^{\frac{1}{2}}$$

(17)

$$K_c=-\frac{2}{25}\left(y_0{x}_0+2y_0{x}_1+2y_0{x}_2+2x_0{y}_1+2x_0{y}_2+4y_1{x}_1+4y_1{x}_2+4y_2{x}_1+4y_2{x}_2\right)$$

(18)

Equations (13), (14), (17) and (18) lead to $K_b=0$ and $K_c=0$.

In [34], we introduced the probability $p$ that all 25 d-d couplings are realized in the metal and probability $(1-p)$ that only diagonal couplings are possible. Then

$$K_b=\{(1-p)\left(\frac{y_0^2+2y_1^2+2y_2^2}{5}\right)+p(\frac{y_0^2}{25}+\frac{4y_1^2}{25}+\frac{4y_2^2}{25}+\frac{4y_0{y}_1}{25}{+\frac{4y_0{y}_2}{25}+\frac{8y_1{y}_2}{25})\}}^{\frac{1}{2}}$$

(19)

$$\begin{gathered} K_c=-2\{(1-p)\left(\frac{y_0{x}_0+2y_1{x}_1+2y_2{x}_2}{5}\right)+p(\frac{y_0{x}_0}{25}+\frac{2y_0{x}_1}{25}+\frac{2y_0{x}_2}{25}+ \\ +\frac{2x_0{y}_1}{25}+\frac{2x_0{y}_2}{25}+\frac{4y_1{x}_1}{25}+\frac{4y_1{x}_2}{25}+\frac{4y_2{x}_1}{25}+\frac{4y_2{x}_2}{25})\} \end{gathered}$$

(20)

After simple operations, Equations (19) and (20) can be respectively transformed to the following ones:

$$\begin{array}{l}{K}_b=\frac{1}{5}{\left[\left(5-4p\right)y_0^2+\left(10-6p\right)(y_2^2+y_1^2)+4p{y}_0(y_1+y_2)+8p{y}_1{y}_2\right]}^{\frac{1}{2}}\end{array}$$

(21)

$$K_c=-\frac{2}{5}[\left(1-\frac{4p}{5}\right)y_0{x}_0+\left(2-\frac{6p}{5}\right)(y_1{x}_1+y_2{x}_2)+\frac{2p}{5}(y_0(x_1+x_2)+x_0(y_1+y_2))+\frac{4p}{5}(y_1{x}_2+y_2{x}_1)]$$

(22)

At $p=0$, we have $K_b=28.06/\pi$ and $K_c=225/{\pi }^2$, and at $p=1$, $K_b=K_c=0$. In other words, the first limit case corresponds to the traditional Wills-Harrison model and the second one corresponds to the Bretonnet-Silbert model.

3. Results and Discussion

The calculations for pure liquid Cu, Ag and Au were performed near their melting temperatures using experimental values of the mean atomic volumes, $\mathsf{\Omega }=1/\rho$, [38] (

Table 1. Input experimental information.

Property	Cu	Ag	Au
$T$ (K)	1423	1273	1423
$\mathsf{\Omega }$ (a.u.)	89.34	130.38	128.31

). Other input data used for the calculations are listed in

Table 2. Input data used for calculations.

Metal	${\mathit{r}}_{\mathit{d}} (\mathbf{a}.\mathbf{u}.)$ [25]	${\mathit{R}}_{\mathbf{C}} (\mathbf{a}.\mathbf{u}.)$ [39]	$\mathit{a} (\mathbf{a}.\mathbf{u}.)$ [39]	${\mathit{z}}_{\mathit{s}}$ [39]	${\mathit{z}}_{\mathit{d}}$	$\mathit{\nu }$ [40]
Cu	1.266	1.2	0.247	1.4	9.6	12
Ag	1.682	1.1	0.213	1.4	9.6	12
Au	1.909	1.0	0.185	1.4	9.6	12

. The values of the parameter $r_d$ are taken from [25] where they were obtained by fitting the bandwidths using the atomic-sphere approximation. The BS parameters $R_{\mathrm{C}}$ and $a$ were fitted in [39] to reproduce numerically the experimental [38] isothermal compressibility and the oscillations of $g(r)$ for the metals being investigated, respectively.

To address how the obtained pair potentials describe the structure of the melts under consideration, we calculated the corresponding RDFs using MD simulations and compared the results with those obtained experimentally [38]. For MD simulations, we used LAMMPS code [41]. The systems of 6912 particles were simulated under periodic boundary conditions in Nose-Hoover NVT ensemble. The MD time step was 1 fs. Initial configurations were created as fcc lattices, which were melted at high temperatures and then equilibrated at a target temperature for 50 ps (50,000 MD steps). Then, data used to calculate $g(r)$ were collected for another 50 ps. Simulations were performed at $T$ and $\rho$ corresponding to experimental conditions presented in [38]. RDFs were calculated using visual molecular dynamics (VMD) code.

Results for $\phi (r)$ and $g(r)$ are presented in Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8.

As seen from Figure 3, Figure 5 and Figure 7, calculated pair potentials demonstrate a non-monotonous behavior with an increase of the portion of the non-diagonal d-d couplings. Indeed, at low values of $p$, the position of the first potential’s minimum, $r_{\min }$, shifts to the left-hand side, and its depth, $\phi (r_{\min })$, increases with increasing $p$ in all three metals. Then, at some value of $p$, these tendencies can change to the opposite ones depending on the metal. This behavior of the pair interactions is not intuitively expectable and it has no clear explanation yet. It is obvious that there is some d-electron characteristic whose dependence on the portion of the non-diagonal d-d couplings has an extremum over $p$ which causes the non-monotonous behavior of pair potentials discussed above.

However, as seen from Figure 4, Figure 6 and Figure 8, the observed non-monotonicity in behavior of pair potentials does not affect the structure of the melts, at least on the level of pair correlation functions. The reason is that the $g(r)$ behavior at studied densities is mostly determined by the repulsive parts of the potentials. As seen from the insets in Figure 3, Figure 5 and Figure 7, RDF correlates with the softness of the repulsive part of $\phi (r)$; the bigger the softness, the smaller the height of the first peak of $g(r)$.

The comparison of calculated and experimental [38] RDFs reveals that the experimental height of the first peak of $g(r)$ is better described at $p=0.25$ for Cu and at $p=0.75$ for Ag and Au. At the same time, for subsequent oscillations of $g(r)$, in Ag and Au the better agreement is observed in the case of the BS model ($p=1$), while for Cu the theoretical results demonstrate weak dependence on $p$. Thus, it is possible to conclude that on the whole the best agreement with experimental RDFs for Cu, Ag and Au are achieved at $p=0.25$, $p=0.75$ and $p=1$, respectively.

At these optimal values of $p$, for each metal under consideration, we calculated structure factor, $S(q)$, and found a quite satisfactory agreement with experimental one [38] except the height of its first peak (Figure 9, Figure 10 and Figure 11). The biggest difference between calculated and experimental first-peak’s heights of $S(q)$ takes place in the case of Cu for which the agreement with experiment of oscillations of $g(r)$ is the worst among the metals under study (Figure 4, Figure 6 and Figure 8).

To demonstrate the reliability of the developed potentials, we compare the results on $g(r)$ with the results of ab initio molecular dynamics (AIMD) simulations performed at the same thermodynamic conditions [42]. Moreover, we carried out MD simulations using available embedded-atom-method (EAM) potentials for Cu [43], Ag [44] and Au [45] at the same protocol as described above. One can see that all calculated $g(r)$ curves agree well with each other and with the experiment (Figure 12, Figure 13 and Figure 14).

Thus, the developed potentials provide a good description of the structural properties of liquid noble metals. Note that the calculation of dynamical properties requires much more precision in the description of interatomic forces than in the case of structural characteristics. To address this point, we used the LiquidLib code [46] to calculate velocity autocorrelation function (VAF, normalized on its initial value) for each metal using both the developed pair potentials (at optimal values of $p$) and EAMs [43,44,45] and compared the results with those obtained by AIMD [42] (see Figure 15, Figure 16 and Figure 17). We see that the pair potentials provide satisfactory agreement with AIMD data. For Ag and Au, the agreement is even better than for EAM-based simulations.

4. Conclusions

Here we utilize the modified Wills-Harrison model, which takes into account the non-diagonal d-d-electron couplings between neighboring atoms, to generate effective pair interactions in noble metals, namely Cu, Ag, and Au. The results obtained show that the pair potentials and structure of liquid noble metals change substantially with an account of the non-diagonal d-d couplings. MD simulations with the developed potentials reveal that the best accuracy of simulated correlation functions is achieved when the used model is closer to the original WH model for Cu, to the BS model for Au, and in an intermediate case for Ag. To analyze in more detail how the portion of the non-diagonal d-d couplings affects observable properties of the metals considered, it is necessary to study its thermodynamic properties that will be the subject of our future investigations.

The results obtained show that reliable effective pair interactions can be used for quantitative description of the structure and properties of liquid metals and alloys. Although more sophisticated potentials and force fields exist in the literature, we have shown that pair potentials developed in [34] can provide comparable accuracy but, at the same time, have smaller complexity and larger computational efficiency.