text:
# Screened Coulomb interactions in metallic alloys: II Screening beyond the single-site and atomic sphere approximations

## I Introduction

One of the main problems of modern alloy theory is to establish a quantitatively accurate description of the configuration-dependent part of the free energy, i.e., the difference in the total energies of alloys with different atomic arrangements on the underlying lattice, in terms of effective cluster interactions which may subsequently be used in statistical thermodynamics simulations . Even without lattice relaxation effects (which are not considered here, although they play an important role in the phase equilibria of many alloy systems) a solution to the problem is still a challenge especially in the case of inhomogeneous systems such as surfaces in the presence of long-range and multi-site interactions which cannot be neglected.

The challenge originates from the fact that quantitatively accurate and reliable (within the accuracy of the approximation for the exchange-correlation part of the total energy of the electronic subsystem in density functional theory (DFT) ) effective cluster interactions can be obtained only by the Connolly-Williams (CW) or structure inversion method  on the basis of the total energies of a set of specifically chosen ordered structures calculated by the so-called full potential (FP) methods, which have no restrictions on the form of the one-electron potential and density. However, the structure inversion methods become practically unusable in the case of an inhomogeneous system, not only because of the large number of basic structures which must be calculated to extract position- (layer-, for instance) dependent interactions, but mainly because of the large size of those basic structures (supercells) which are needed to factorize a specific interaction inside homogeneous parts of the systems, e.g., inside the layers parallel to a surface.

In this situation there appears to be only one alternative to the structure inversion methods: The so-called generalized perturbation method (GPM) proposed by Ducastelle et al.  on the basis of the coherent potential approximation (CPA)  and formulated within tight-binding (TB) theory. Later the GPM was generalized in a straightforward manner  for use in ab initio calculations based either on the Korringa-Kohn-Rostoker (KKR) method or the linearized muffin-tin orbitals (LMTO) method in the atomic sphere approximation (ASA). The main idea behind the GPM is to calculate perturbatively the total energy difference between the alloy in the initial state, which is completely random, and in a final state in which only one specifically chosen atomic distribution correlation function or short-range order parameter is different from that in the random state. This makes the GPM very efficient and convenient to use as it directly yields the needed effective cluster interactions.

However, it is known, although rarely mentioned in the literature, that the interactions obtained by the GPM yield a quantitatively poor description of the ordering in real alloys (see, for instance, Ref. ) in those cases where there is a substantial size mismatch between the alloy components. This failure may only partly be attributed to lattice relaxation effects. Rather, it originates not from the GPM method itself but is a consequence of inappropriate approximations in the basic methods underlying the GPM calculations. This is so, because, as has been demonstrated by Bieber et al.  in parametrized tight-binding calculations and by Singh et al.  in ab initio KKR-CPA calculations, the GPM interactions may provide (under certain conditions) a consistent description of the ordering or configurational energy. That is, the ordering energy obtained from the GPM interactions, calculated in the framework of a particular technique, agrees reasonably well with the ordering energy obtained directly from the total energy calculations by the same technique.

The ab initio techniques underlying GPM calculations are usually the KKR-ASA, KKR-ASA-CPA, and LMTO-CPA methods  which are based on a number of approximations, such as the CPA, the single-site (SS) approximation for the electrostatic part of the DFT problem, and the spherical approximation for the form of the potential which, depending on the geometry, is called either the muffin-tin (MT) or atomic sphere approximation (ASA). The question is which of these approximations is the most severe in the cases where the alloy components have a substantial size difference. To answer the question, we note that a size difference leads to so-called "charge transfer effects" or, to be more precise, to a non-zero net charge for each alloy component inside their atomic spheres chosen to be of equal size.

Although there are systems, where the CPA may lead to substantial errors, it is clear from a general point of view that the CPA cannot be responsible for the errors in the case of pronounced charge transfer effects, because the error of the CPA is mainly related to specific features of the electronic structure of the individual alloy components such as the difference in the position and overlap of the energy bands . Moreover, there is a number of different calculations which show that in such systems as, for instance, CuPd and CuAu, where the alloy components have similar electronic structures but different atomic size, the CPA works fairly well  as a method for obtaining the average electronic structure of random alloys.

As far as the underlying KKR-CPA or LMTO-CPA methods are concerned, much larger errors may in fact come from the use of the single-site approximation in the self-consistent DFT part of the calculations as this yields no information about the distribution of the charge outside the individual atomic spheres of the alloy components. In fact, the effective medium outside the individual atomic spheres of the alloy components is electroneutral, and therefore, if the net charge of an atomic sphere is non-zero, Poisson’s equation cannot be solved properly. A number of different models have been proposed to include the missing screening charge in the solution to Poisson’s equation  and, most recently, a general formalism of screened Coulomb interactions (SCI) based on the knowledge of the spatial distribution of the screening charge around an impurity has been developed in Ref.  together with a formalism for the SCI contribution to the GPM interactions . Although the SCI may now be included in SS-DFT-CPA calculations, this does not solve all problems connected with the description of the energetics of alloys.

It is not surprising that the main source of inaccuracy in the KKR-CPA and LMTO-CPA methods is the spherical approximation, MT or ASA, for the form of the electron density and potential (in the following we will consider only the ASA, since the difference between the ASA and MT is unimportant for the later discussion and results). For instance, in the extreme case, where one of the alloy components is a vacancy, the error due to the use of the ASA is about 100% (or several eV in absolute values) for the vacancy formation energy . As has been shown by Korzhavyi et al. , this kind of error originates from the oversimplified description of the non-spherical electrostatic contribution to the one-electron potential and energy mainly from the charge density on the atoms next to the vacancy.

This is similar to the case of surfaces where the quite large ASA error may be substantially reduced by the inclusion of the multipole moments of the electron charges inside the atomic spheres . The so-called ASA+M approach significantly improves vacancy and defect formation energies , surface energies  and alloy energetics . Recently, Finnis et al.  have included the multipole moments in their self-consistent tight-binding model which allowed them to obtain quite accurate description of the energetics of zirconia.

In this paper we show that the use of the ASA+M approach leads to a representation of the configurational part of the total energy, which is very close to the full-potential results. Since the polarization of the atoms in an alloy is almost entirely determined by their closest local environment, it is obvious that the effect of polarization cannot be described properly in the single-site approximation. Nevertheless, the SS-DFT-CPA methods may be still used for the electronic structure and total energy calculations of random alloys if the definition of the SCI is modified. It is the main purpose of the present paper to demonstrate how this may be done in the cases of ordinary bulk homogeneous random alloys and inhomogeneous systems such as surfaces.

The paper is organized as follows. In section II we introduce the ASA+M approximation and outline some details of our calculations. In section III we compare the ordering energies of NiPt alloys, calculated by the KKR method in different approximations and by the Vienna ab initio simulation package (VASP) . In section IV we define the on-site screening Madelung potential, which should be added to the one-electron potential in the SS-DFT calculations. The SCI and the Madelung energy of a random alloy are defined in section V. In section VI we calculate the intersite SCI in NiPt in different approximations and by different methods. In section VII the screened generalized perturbation method interactions are calculated and compared with the Connolly-Williams interactions. In section VII we show how the formalism for the SCI should be modified in the case of inhomogeneous systems where there are several non-equivalent sublattices, like partially ordered alloys or surfaces.

## II Beyond the ASA

In a previous publication , in the following referred to as I, we presented a consistent and variational, within DFT, approach to the electrostatic screening effects in random alloys, and within the ASA we found that these screening effects were almost independent of alloy composition, lattice spacing, and crystal structure. However, in those cases where the alloy components have a substantial size difference one cannot obtain a quantitatively correct description of the configurational part of the total energy of metallic alloys within the ASA. One must therefore go beyond the spherical approximation at least for the density. In the following we will show how this may be done. All other detail of our approach may be found in I.

### A Multipole correction to the atomic sphere approximation: ASA+M

The idea behind the multipole correction is simply to include those contributions to the electrostatic multipole-moment expansion of the intercell or Madelung part of the one-electron potential and energy, which are neglected in the ASA. If the multipole moments of the electron charge, $`q_𝐑^L`$, inside an atomic sphere centered at R are defined as

$$q_𝐑^L=\frac{\sqrt{4\pi }}{2l+1}_{S_𝐑}\left(\frac{r}{S_𝐑}\right)^ln_𝐑(𝐫_𝐑)Y_L(\widehat{r}_𝐑)𝑑𝐫_𝐑Z_𝐑\delta _{0l},$$
(1)
where $`L`$ is short-hand for the $`(l,m)`$ quantum numbers, $`S_𝐑`$ the radius of the atomic sphere, $`n_𝐑`$ the non-spherical charge density, and $`Y_L`$ a real harmonic, the Madelung contribution to the one electron potential is given by

$$v_{M𝐑}^L=\frac{1}{S}\underset{𝐑^{},L^{}}{}M_{𝐑,𝐑^{}}^{L,L^{}}q_𝐑^{}^L^{},$$
(2)
while the Madelung energy which now includes the multipole-multipole electrostatic interactions between different lattice sites may be written

$$E_M=\frac{1}{2S}\underset{𝐑,L}{}q_𝐑^L\underset{𝐑^{},L^{}}{}M_{𝐑,𝐑^{}}^{L,L^{}}q_𝐑^{}^L^{}.$$
(3)
In these expressions, $`M_{𝐑,𝐑^{}}^{L,L^{}}`$ is the multipole Madelung matrix which is equivalent to the conventional LMTO structure constants and the number of multipoles included in the $`L,L^{}`$ summations is determined by the angular momentum cutoff $`l_{max}`$ in the basis set used in the Green’s functions calculations. Owing to the properties of the Gaunt coefficients non-zero multipole moments of the charge density may be generated for $`l`$-values up to $`l_{max}^M=2l_{max}`$.

We note that, since in the ASA+M the one-electron potential is still kept spherically symmetric inside each atomic sphere, the only term which contributes to the one-electron potential is the L=(0,0) or monopole term. This simple restriction on the form of the-one electron potential violates the variational connection between the Madelung potential and energy and, in turn, between the one-electron potential and the total energy, i.e.,

$$v_{M𝐑}v_{M𝐑}^{00}\frac{\delta E_M}{\delta n_𝐑}.$$
(4)
However, since this is just a consequence of the model, but not of theory in general, it does not create any problems. On the other hand, the reinstatement of the variational connection between the one-electron potential and the total energy by keeping only the monopole-multipole term in (3) may lead to a substantial error in the total energy, as the multipole-multipole interactions will not be accounted for.

### B Details of calculations

The Green’s function technique has been used in both the KKR-ASA and the locally self-consistent Green’s function (LSGF) calculations in the scalar relativistic and atomic sphere approximations. This part of the techniques is described in Refs. . The basis functions have been expanded up to $`l_{max}=3`$ ($`spdf`$-basis) inside the atomic spheres, while the multipole moments have been calculated up to $`l_{max}^M=6`$. We have also performed a number of calculations in the $`spd`$-basis, i.e., $`l_{max}=2`$ and $`l_{max}^M=4`$. The integration of the Green’s function over energy was performed in the complex plane over 16 energy points on a semicircular contour using a Gaussian technique. We have used the generalized gradient approximation (GGA) of Perdew and co-workers . For each structure the integration over the Brillouin zone (BZ) has been done by using equally spaced k-points in the irreducible part of the appropriate BZ and the number of k-point has been chosen to be equivalent to 500 - 1000 uniformly distributed k-points in the irreducible part of the BZ of the fcc structure. Core states were recalculated at each DFT iteration.

For benchmark calculations we applied the Vienna ab initio simulation package (VASP) described in detail in Refs. . These calculations were performed in a plane-wave basis, utilizing fully non-local Vanderbilt-type ultrasoft pseudopotentials (US-PP)  which allow the use of a moderate cutoff in the construction of the plane-wave basis for the transition metals. In the actual calculations the energy cutoff was set to 302 eV, exchange and correlation were treated in the framework of the GGA , and the integration over the Brillouin zone was performed on a Monkhorst-Pack k-mesh . Test calculations showed that, depending on structure, the required convergence was reached for 35 to 275 k-points in the irreducible wedge of the BZ.

## III Ordering energies in Ni-Pt

We start by demonstrating the accuracy of the various approximations which are usually used in KKR(LMTO)-ASA-like calculations. For this purpose we have calculated a set of ordered fcc NiPt alloys by the KKR-ASA method and by the US-PP for a fixed lattice constant without any local or anisotropic relaxations. The lattice constant has been chosen to be $`a`$ 3.791 Å, which corresponds to an atomic Wigner-Seitz radius of 2.8 a.u.. The ordered structures include: L12 (Cu3Au-type), DO22 (TiAl3-type), Z3, $`\beta `$ , $`\gamma `$ (Pt2Mo-type), L10 (CuAu-type), CH or "40" (NbP-type), Z2 , L11 (CuPt-type), and the so-called SQS-16 .

To simplify the comparison we present in Table I the values of the calculated "mixing" energies of the above mentioned ordered structures

$$E_{mix}^{Ni_mPt_n}=E_{tot}^{Ni_mPt_n}\frac{mE_{tot}^{Ni}+nE_{tot}^{Pt}}{m+n},$$
(5)
where $`E_{tot}^{Ni}`$ and $`E_{tot}^{Pt}`$ are the total energies of the pure components calculated at the same lattice constant. All total energies are per atom. The mixing energies of the random alloys have been obtained on the basis of all the energies included in Table I, except SQS-16, plus the energies of Ni7Pt and Pt7Ni (CuPt7-type) ordered alloys (not presented in the Table) by the Connolly-Williams method in which the total energy expansion included pair interactions at the first four, 7th and 10th coordination shells (these are the largest pair interactions in this as well as in many other fcc transition metal alloys), the first four triangle interactions and the two tetrahedron interactions corresponding to the tetrahedra of nearest neighbors, and the straight line along  direction (the last being quite substantial in many systems).

The SQS-16 is a so-called special quasirandom structure  which consists of 8 atoms of one type and 8 atoms of another type distributed in the unit cell in such a way that the first 7 pair, the nearest neighbor triangle, and the tetrahedron atomic distribution correlation functions are the same as in the random alloy. Hence, the fact that the values of $`E_{mix}^{SQS16}`$ and $`E_{mix}^{rand}`$ are nearly eaqual indicates that: 1) the SQS-16 provides a good model for the random NiPt alloy, and 2) the convergence of the CW method is reasonably good. Part of the convergence of the CW method is provided by the use of total energies on a fixed lattice, whereby the volume dependent contribution to the total energy is not expanded in terms of cluster interactions, which is an ill-defined procedure in metallic systems and usually leads to very bad convergence of the CW method (see, for instance, ).

In Table I we present KKR-ASA results in the $`spdf`$\- as well as the $`spd`$-basis and in both cases we show results in the pure ASA, i.e., without multipole moment contributions, and in the ASA+M. We also include the results of neutral sphere calculations, in which the atomic sphere of Pt is chosen to be larger than that of Ni in order to provide zero net charges of the atomic spheres. Although the comparison of the mixing energies should be done with some caution, because the ground state properties of the alloys are different in different approximations and because all the calculations have been performed at the same fixed lattice constant, it is clear that the ASA+M approach in the $`spdf`$-basis leads to values of $`E_{mix}`$ which are in considerably better agreement with the US-PP results than any of the other approaches.

Using the results of the Connolly-Williams method for the total energy of random alloys one can calculate the ordering energies, defined as the difference between the total energies of an ordered and a random alloy at the same composition. In Table II we compare the ordering energies of different structures calculated in different approximations for Ni3Pt, NiPt, and Pt3Ni. Such a comparison makes sense since the ordering energies are much less volume-dependent than the total energies themselves. Again, it is seen that, relative to the US-PP results, the ASA+M approach in the $`spdf`$-basis gives not only the best, but also a quite accurate description of the ordering effects in NiPt. It is also seen that the ordering energies in the ASA+M but without multipole-multipole contribution to the Madelung energy appear to be half way beteeen ASA and ASA+M results. As we will see below, this is in fact the limit of accuracy which can be reached in consistent single-site mean-field calculations.

It is also obvious from the table that the KKR-ASA does not yield reasonable values for the ordering energies in the case of transition metal alloys unless $`f`$-states are included in the basis. These $`f`$-states are needed to supply a better augmentation of the basis functions at the atomic sphere and a better interstitial charge density . The neglect of $`f`$-states can only be partly compensated by the use of the so-called combine-correction term in the LMTO method .

Another important conclusion, which can be drawn from the results in Table II, is the fact that the use of neutral spheres leads to a substantial underestimate of the ordering effects in KKR-ASA(+M) calculations. In other words, although the neutral-sphere approach formally solves the problem of the electrostatic interaction in an alloy, the electrostatic contribution to the one-electron potential and energy being zero by definition, it introduces errors which are unacceptable in a quantitative description of the configurational energetics.

The reason for this failure is the following. If we compare the values of $`E_{mix}`$ from Table I obtained with neutral spheres with those obtained with equal spheres, we find that the neutral sphere approach leads to substantial lowering of the total energy of the ordered alloys. However, the amount of lowering is structure dependent: In the case of "phase separated"-like structures, where there is a certain clustering of atoms of one type, like Z3, Z2, the lowering is much greater, than it is in some of the more "normal" structures. Through the CW procedure the exaggerated lowering of the energies of the Z3 and Z2 structures leads to an exaggerated lowering of the total energy of the random alloy, making, in the end, all the ordering energies much smaller in absolute value, than they should be. In the case of the "phase-separated"-like structures this kind of error originates in the poor filling of the crystal space provided by the non-equal atomic spheres. The division of a crystal into atomic spheres is in itself quite a crude approximation, the error of which is not generally known. Now, if such division is done differently for different ordered structures by means of non-equal atomic spheres, it is obvious that the resulting errors will be structure dependent and render a comparison of total energies meaningless.

Based on the above considerations it appears that the use of equal atomic spheres for the alloy components provides the only consistent and correct way of dividing space for configurational or ordering energy calculations, if the underlying lattice is a simple Bravais lattice, such as fcc and bcc. In that case the packing of space by the atomic spheres of the alloy components is homogeneous and independent of the alloy configuration which makes at least part of the ASA error systematically compensated when the energy difference of different structures is calculated. Of course, the situation is different in alloys where the underlying lattice has a more complicated structure. However, in the configurational energy calculations one should follow the recipe of choosing equal sphere radii of the alloy components on the sublattice, where the alloying (ordering) is taking place (see, for instance,).

## IV The on-site screening Madelung potential in the ASA+M

Although the ASA+M approach seems a natural generalization of the ASA, that changes little in the formalism, it has a large effect on the way the SCI must be treated in random alloys. The most dramatic consequence of the ASA+M is the fact that the correct SCI can no more be obtained in single-site electronic structure or Green’s function calculations: The largest polarization effects, which give multipole-moment contributions to the SCI, actually comes from the nearest neighbors of the atom the screening of which is considered.

This is easy to see in the case of a single impurity in an otherwise perfect crystal which is the dilute limit of a random alloy. The point-group symmetry of the impurity site is exactly the same as that of the underlying lattice, while none of its neighbors has even inversion symmetry and, therefore, every atomic sphere around the impurity has a non-zero dipole moment. Such a dipole moment plus the higher multipole moments induced by the presence of the impurity can only be found in Green’s function calculations which include these neighboring sites in the perturbation part of the Dyson equation. Thus, the multipole-moment contribution to the SCI can be obtained only beyond the single-site approximation in the Green’s function calculations.

Nonetheless, we will determine a simple parameterization of the SCI for the SS-DFT-CPA method, which will allow us to obtain accurate results in SS-KKR(LMTO)-CPA calculations. Such a parameterization of the SCI in the ASA+M is purely a fitting procedure in contrast to the SCI determined by SS-ASA-Green’s function calculations in I basically for the purpose of showing the existence of a consistent theory of the SCI in random alloys in the framework of the single-site-CPA theory.

Let us first consider the on-site Madelung potential in a random alloy due to the screened Coulomb interactions which should be used in SS-DFT-CPA calculations for the random alloy. In the case of a binary random alloy, A1-cBc, this potential can be defined as (see I)

$$v_i=e^2\alpha _{scr}(0)\frac{q_i}{S},$$
(6)
where $`q_i`$ is the net charge of the $`i`$th alloy component, $`S`$ the Wigner-Seitz radius, and $`\alpha _{scr}^i(0)`$ the on-site screening constant which can be obtained from the screening charge by performing "impurity" calculations (see I), for instance, by exchanging the type of atom at a particular site of the supercell modeling the random alloy. Then $`\alpha _{scr}(0)`$ is

$$\alpha _{scr}(0)=\frac{2S}{e^2}\underset{L,𝐑}{}S_{0,L}(𝐑)Q_L(𝐑),$$
(7)
where $`S_{L,L^{}}(𝐑^{})`$ are the canonical structure constants, the first few terms of which are explicitly defined, for instance, in Ref. . Further, $`Q_L(𝐑)`$ are the normalized multipole moments of the screening charge in the atomic sphere centered at R

$$Q_L(𝐑)=\frac{\mathrm{\Delta }q_L(𝐑)}{\mathrm{\Delta }q(𝐑=0)}\frac{\mathrm{\Delta }q_L(𝐑)}{\mathrm{\Delta }q_{L=00}(𝐑=0)},$$
(8)
where $`\mathrm{\Delta }q_L(𝐑)`$ is the difference between the L-moments of the charge in the atomic sphere at $`𝐑)`$ after and before the impurity has been introduced at the site $`𝐑=`$0.

There is one important point. In the pure ASA $`\alpha _{scr}(0)`$ can be determined in an alternative although formally equivalent manner by

$$\alpha _{scr}(0)=(\alpha _{rand})=\frac{S}{e^2}\frac{<v_i>}{<q_i>},$$
(9)
where $`<V_i>`$ and $`<q_i>`$ are the average values of the Madelung potential and net charges of the $`i`$th alloy component in the self-consistent supercell calculations of the random alloy. Being practically exact in the ASA this scheme is only approximately valid in the ASA+M, where the average Madelung potential of the supercell, which is equal to the Madelung potential of the underlying lattice, is not equal to zero in general due to the presence of the non-zero multipole moments at least for $`l=`$4. In the case of inhomogeneous systems like surfaces, where the multipole contribution to the average Madelung potential is quite large, Eq. (9) cannot be used at all. Therefore, (9) should be modified by substracting the corresponding average values of the net charges and Madelung potentials of the alloy components on the (sub)lattice, $`\overline{q}=(1c)<q_A>+c<q_B>`$ and $`\overline{v}=(1c)<v_A>+c<v_B>`$

$$\alpha _{scr}(0)=(\alpha _{rand})=\frac{S}{e^2}\frac{<v_i>\overline{v}}{<q_i>\overline{q}}.$$
(10)
In Fig. 1 we compare the monopole moment of the normalized screening charge, $`Q_{L=00}(R)`$, in a Ni50Pt50 random alloy obtained by changing the type of the atom from Pt to Ni at a some particular site in a 384-atom Ni50Pt50 supercell, the first seven short-range order (SRO) parameters of which are equal to zero. The calculations have been carried out in the single-site (SS)-LSGF-ASA as well as the embedded-cluster (EC)-LSGF-ASA+M methods (see I). The local interaction zone (LIZ) in the EC-LSGF-ASA+M calculations included 43 atoms, i.e., the central atom and its first three coordination shells (LIZ=4), while in the SS-LSGF calculations the LIZ included only one atom (LIZ=1).

It may be seen from Fig. 1 that the screening is more efficient in the ASA+M than in the ASA, although the difference between the two cases is very small. Nevertheless, the effect of the multipole moments on the on-site screening constant $`\alpha _{scr}(0)`$ is quite pronounced, it increases to about 0.74 from the 0.61 in the ASA calculations. The largest multipole contribution comes from the dipoles in the first coordination shell, which contribute almost 0.1 to the $`\alpha _{scr}(0)`$, while the quadrupole and octuple moments contribute 0.036 and 0.016, respectively.

The effective charge transfer, $`\mathrm{\Delta }q=<q_{Ni}><q_{Pt}>`$, increases from 0.505 in the ASA to 0.583 in the ASA+M. However, it is still reproduced correctly in SS-DFT-CPA calculations, provided the correct, i.e., corresponding to the ASA+M, value of the on-site screening constant, $`\alpha _{scr}(0)`$, is used in (6). This is so, because the on-site screening constant is a parameter which determines a constant shift of the one-electron potential both in the ASA and ASA+M. This constant may be chosen to contain the monopole-multipole interactions of the charge at a given site with its screening cloud, thereby yielding the correct effective charge transfer. Unfortunately, the same is not the case for the Madelung energy and the SCI in general.

## V The SCI and Madelung energy of a random alloy

The screened Coulomb interactions, $`V_{scr}(R)`$, are the energies of the electrostatic interaction between the electron charge density inside an atomic sphere centered at some site (it is convenient to choose this site as the origin) and the perturbed electron density and its screening charge at another, in general, different site. In the ASA+M, the SCI can be represented in the form of the multipole expansion,

$$V_{scr}(R)=\underset{L}{}V_{scr}^L(R),$$
(11)
where $`V_{scr}^L(𝐑)`$ is the $`L`$-component contribution to the SCI at the distance $`R`$, which for a binary A1-cBc alloy on a Bravais lattice may be expressed as in I

$`V_{scr}^L(𝐑)`$ $`=`$ $`{\displaystyle \frac{e^2}{2}}\mathrm{\Delta }q^2{\displaystyle \frac{Q_L(𝐑=0)\alpha _{scr}^L(𝐑)}{S}}`$ (12)
$`=`$ $`{\displaystyle \frac{e^2}{4S}}\mathrm{\Delta }q^2Q_L(𝐑=0){\displaystyle \underset{L^{},𝐑^{}0}{}}S_{L=00,L^{}}(𝐑^{})Q_L^{}(𝐑^{}𝐑).`$ (13)
Here $`\mathrm{\Delta }q`$ is the difference of the net charges in the atomic sphere after and before the perturbation, i.e., the exchange of the type of atom at site $`𝐑`$, $`\alpha _{scr}^L(𝐑)`$ is a generalized screening constant, $`S_{L,L^{}}(𝐑^{})`$ are the canonical structure constants , $`Q_L(𝐑^{}𝐑)`$ are the normalized multipole moments of the screening charge, and $`\mathrm{\Delta }q_{L=00}=<q_A><q_B>`$ is the effective charge transfer in the alloy.

Within the single-site mean-field considerations presented in I, all the multipole moments on the alloy sites are uncorrelated, the average value of $`Q_L(𝐑=0)`$ being either very small or equal to zero, unless $`L=00`$, and thus the only non-zero SCI is $`V_{scr}^{L=00}`$, which, for instance, in the case of the fcc underlying lattices can be written in the form (see I)

$$V_{scr}(R)=\mathrm{\Delta }q^2\frac{e^2}{2}\frac{\alpha _{scr}(R)}{S},$$
(15)
where

$$\alpha _{scr}(R)=\frac{1}{2}\underset{L^{},𝐑^{}0}{}S_{L=00,L^{}}(𝐑^{})Q_L^{}(𝐑^{}𝐑).$$
(16)
The on-site term $`V_{scr}(R=0)`$ is the energy of the electrostatic interaction between the net charge of an alloy component and its screening density or, as has been shown in I, it is the screening Madelung energy of the random alloy. It is easy to see that this energy is DFT-consistent with the corresponding screening Madelung shift of the one-electron potential (6). However, in contrast to the screening Madelung potential which correctly reproduces the effective charge transfer in random alloys, $`V_{scr}(R=0)`$ underestimates the corresponding Madelung energy in the ASA+M because of the missing contribution from the multipole-multipole interactions in the single-site mean-field approximation. In the case of, for instance, a Ni50Pt50 random alloy this contribution is about $``$3 mRy/atom for the Wigner Seitz radius $`S=`$2.8 a.u..

This means that if one wants a quantitatively accurate value of the total energy of a random alloy in the SS-DFT-CPA calculations consistent with the supercell ASA+M calculations one needs to modify the definition of the Madelung energy of the random alloy by introducing a fitting parameter. The simplest way to do so is to define the Madelung energy of the random alloy as

$$E_{Mad}^{rand}=\underset{i}{}c_iE_i^{scr},$$
(17)
where $`c_i`$ is the concentration of the $`i`$th alloy component, and

$$E_i^{scr}=\frac{e^2}{2}\beta q_i^2\frac{\alpha _{scr}(R=0)}{S},$$
(18)
which means that in the case of a binary random alloy

$$E_{Mad}^{rand}=c(1c)\beta V_{scr}(R=0).$$
(19)
Here, $`\beta `$ is the renormalization coefficient which is approximately equal to 1.16 for most fcc and hcp transition metal random alloys. Thus, $`E_i^{scr}`$ and $`V_{scr}^i`$ are no more DFT-consistent (see (4)).

As discussed above, this violation of general theory is a consequence of the ASA+M, which on the other hand brings the ordering energies of the much more efficient SS-DFT-CPA approach into good quantitative agreement with the corresponding full-potential results. Although the difference between the multipole-multipole and monopole-multipole results in Table II might not look so dramatic, the omission of the multipole-multipole interactions in the ASA+M Madelung energy has much more serious consequences in the case of, for instance, surface energy anisotropy calculations, which cannot be reproduced even qualitatively without this term.

Finally, we show that the CPA itself introduces relatively small errors in the electronic structure calculations of a random alloy. In Fig. 2 we compare the total density of states (DOS) of a random Ni50Pt50 alloy and the local, Pt and Ni, contributions obtained by two different methods, SS-KKR-CPA and EC-LSGF. The SS-KKR-CPA calculations have been performed with $`\alpha _{scr}=0.74`$ for the screening Madelung potentials of the alloy components (6), and in the EC-LSGF calculations LIZ=4 has been used to calculated the SQS-16 supercell of Ni50Pt50. Since, all the SRO parameter of the SQS-16 are equal to zero up to the 7th coordination shell, beyond which the net charge at each site is practically completely screened, the EC-LSGF results may be considered as a benchmark for the SS-KKR-CPA method. The agreement between the EC-LSGF and SS-KKR-CPA results is seen to be good clearly indicating that the CPA works fairly well for this system.

## VI Intersite screening constants in NiPt

There are several ways of obtaining the intersite screening constants, $`\alpha _{scr}(R)`$, which determine the corresponding intersite SCI. First, they may be obtained directly from the normalized moments of the screening charge, $`Q_L(𝐑)`$, by menas of (16). This requires two self-consistent supercell calculations: One performed for some initial atomic configuration in the supercell and the other for the same supercell with an "impurity" at the site where the type of atom is changed (see paper I). In this manner, however, only the monopole-multipole part of the SCI can be found.

Secondly, one may take advantage of the special properties of the LSGF method, namely of the fact that the electronic structure in the LSGF method is obtained in the so-called combined cluster-effective medium approach . That is, the local environment effects are taken into account only inside the LIZ during the electronic structure calculations, while the rest of the crystal is seen by each atom as a random alloy, described by the CPA effective medium. This means that all the correlated atomic configurations attributed to the non-zero SRO parameters beyond the LIZ do not contribute to the electronic structure of the supercell. At the same time, since the Madelung problem is solved exactly, they are accounted for in the electrostatic part of the total energy.

Considering the change of the atomic configuration on the lattice as a small perturbation, one may argue on the basis of Andersen’s force theorem  that the difference in the total energies of the alloy with some non-zero SRO parameters and a completely random alloy in the LSGF calculations may be given by

$`E_{ord}^{LIZ=N}`$ $``$ $`E_{tot}^{SRO}E_{tot}^{rand}`$ (20)
$`=`$ $`{\displaystyle \frac{1}{2}}c(1c)\left[{\displaystyle \underset{i=1}{\overset{N1}{}}}\alpha _iz_iV_{oneel}(R_i)+{\displaystyle \underset{i=1}{\overset{N_{scr}}{}}}\alpha _iz_iV_{scr}(R_i)\right],`$ (21)
where $`\alpha _i`$ are the Warren-Cawley SRO parameters (see, for instance, Ref. ) at the $`i`$th coordination shell, $`N`$ the size of the local interaction zone, which is determined as the number of coordination shells around the impurity site plus 1, and $`N_{scr}`$ the coordination shell beyond which the SCI vanish. In the case of an fcc alloy it may safely be assumed that, $`N_{scr}=7`$. It is obvious that the first term in (20) may be identified with the usual GPM interactions.

It follows from (20), that if $`N1<N_{scr}`$ and $`\alpha _i`$=0 for all $`i<N_{scr}`$ except for one coordination shell, $`j`$, which is beyond the LIZ, $`j>N1`$, then the intersite SCI for this specific coordination shell can be determined in two LSGF calculations: One performed for a supercell which corresponds to the random alloy ($`\alpha _i=0`$ for all $`i<N_{scr}`$) and the other for a supercell which satisfies the above described conditions with $`\alpha _j0`$. In this case

$$\alpha _{scr}(R_j)=\frac{4SE_{ord}}{e^2c(1c)\mathrm{\Delta }q^2z_j\alpha _j}.$$
(22)
In the ASA all the SCI can be determined by using (22) in the corresponding SS-LSGF-CPA calculations, while $`\alpha _{scr}(R_2)`$ is the first SCI which can be determined in the ASA+M in this way, since there is no multipole contribution to the SCI in the single-site approximation (N(LIZ)=1). It is also clear that if the multipole-multipole Madelung energy is included in the corresponding LSGF calculations, then, in principle, $`\alpha _{scr}(R)`$ determined from (22) will also contain the multipole-multipole contribution.

However, in this approach the higher order atomic distribution correlation functions (multisite SRO parameters) of the supercell should be also optimized. This is so, since all the multisite interactions for the figures inscribed in the LIZ with an vertex located at the central atom of the LIZ, also contribute to the ordering energy. Moreover, in the case of the ASA+M, in principle, there is a non-zero contribution from the SCI to the multisite interactions themselves. In our calculations, we have not optimized the multisite SRO parameters of the supercells and, therefore, considering the quite small values of the SCI beyond the first coordination shell we have not used (22) in the ASA+M calculations.

In Table III we compare $`\alpha _{scr}(R)`$ for the first four coordination shells calculated either by direct summation of the normalized multipole moments (16) or by means of (22) in the SS-LSGF total energies. The LSGF calculations have been performed for 384-atom supercells of an equiatomic Ni50Pt50 alloy with the corresponding sets of SRO parameters. Although we have not optimized the higher order SRO parameters, they turned out to be small: $``$ 0.02 or less, at least for the triangle and tetrahedra of the nearest neighbors on the fcc lattice. Thus they should not affect the results, at least in the ASA. As seen from Table III the agreement between the intersite screening constants, $`\alpha _{scr}(R)`$, determined in two different calculations in the ASA, is quite reasonable.

As shown in I, $`\alpha _{scr}(R)`$ is practically a universal function for metallic alloys on simple Bravais lattices in the ASA and the single-site approximation. In the ASA+M approach, this is not the case any more. However, our calculations for the dilute limit of NiPt alloys and for different lattice constants show that $`\alpha _{scr}(R)`$ changes very little in these cases and, in fact, less than the difference between the values of $`\alpha _{scr}(R)`$ calculated in two different ways. For instance, in Table III we show the values of $`\alpha _{scr}(R)`$ obtained for a Ni impurity in pure Pt at $`S=`$ 3 a.u., which are very close to those of the random equiatomic alloy at $`S=`$ 2.8 a.u.

## VII Screened GPM interactions in NiPt

As demonstrated in I the intersite screened Coulomb interactions, $`V_{scr}(R)`$, given by (15) must be added to the corresponding one-electron term given by the GPM in order to satisfy the force theorem . In the case of a binary random alloy AcB1-c  the screened generalized perturbation method (SGPM) interactions are defined as in Ref. :

$`V^{SGPM}(R)`$ $`=`$ $`V^{GPM}(R)+V_{scr}(R)`$ (23)
$`=`$ $`V^{GPM}(R)+{\displaystyle \frac{e^2}{2}}\mathrm{\Delta }q^2{\displaystyle \frac{\alpha _{scr}(R)}{S}},`$ (24)
where $`V^{GPM}(R)`$ is the usual GPM interactions obtained as the change of the one-electron energies due to specifically induced alloy configurations on the alloy’s underlying lattice. In fact, $`V^{GPM}(R)`$ should be renormalized due to the intersite SCI. However, this problem, as well as the complete SGPM formalism will be considered elsewhere .

In Table IV we compare the first four most important effective pair interactions (the rest of the pair interactions and the multisite interactions are less than 0.1 mRy) in the NiPt fcc alloys at $`S=`$2.8 a.u. obtained by three different techniques, i) the SGPM for equiatomic alloy composition, ii) the Connolly-Williams method on the basis of the total energies of the ordered alloys described above, and iii) direct calculation from the EC-LSGF total energies of Ni50Pt50 alloys similar to the case of the intersite screening constant calculation described in the previous section. That is, in iii) the effective interactions have been obtained from the EC-LSGF total energy calculations for a completely random alloy, $`\alpha _i=`$0 for all $`i<N_{scr}`$, and for an alloy with one non-zero $`\alpha _j`$. However, now $`j<N1`$, where $`N`$ is the size of the LIZ in the LSGF calculations, and therefore the local environment effects attributed to the non-zero SRO parameter, are included in the electronic structure calculations. In this case the effective pair interactions can be determined as

$$V_j=\frac{2E_{ord}^{LIZ=N}}{c(1c)z_j\alpha _j}.$$
(25)
It is important, that (25) is not based on any additional approximations and therefore constitutes a direct way of determining the effective interactions with an accuracy which, in principle, is restricted only by the approximations used in the LSGF calculations, that is, mainly by the ASA+M, since the CPA yields very small relative errors in the EC-LSGF calculations with LIZ$`>2`$.

Although the SGPM interactions as well as the interactions determined from (25) are concentration dependent, while the Connolly-Williams interactions are concentration independent, they can be compared since the concentration dependent interactions obtained for an equiatomic alloy composition are equal to those of the concentration independent interactions, at the same fixed volume . Therefore, if the basis in the Connolly-Williams method includes all the important interactions for a given system, the Connolly-Williams interactions obtained from the KKR-ASA+M calculations of the ordered alloys at a fixed lattice constant should be equal to those obtained from (25) for the equiatomic alloy composition at the same lattice constant.

It is clear from Table IV, that in general the agreement between the Connolly-Williams and the direct calculations is quite good. Further, the Connolly-Williams interactions obtained by the KKR-ASA+M and by the US-PP (CW-KKR-ASA+M and CW-US-PP) agree well with each other, except for $`V_3`$ which is a little larger in the CW-US-PP calculations than in the CW-KKR calculations. This, then, confirms our point that the ASA+M approach allows us to obtain quantitatively accurate configurational energies of metallic alloys.

It is interesting to note that although the neutral-sphere approach yields quite large errors for the ordering energies, it seems to work remarkably well for the GPM interactions, except for $`V_1`$, which is more than twice as small as it should be. It is probably a coincidence that they come out very close to the CW-KKR-ASA+M interactions, since the SCI are quite small beyond the first coordination shell. Unfortunately, the CW interactions obtained in the neutral sphere calculations do not seem to be convergent, the three- and four-site interactions being of the same order of magnitude as the interaction at the first coordination shell and therefore, they are not given in the table.

The SGPM interactions, obtained in the ASA+M and in the ASA (in these two cases we have used the intersite SCI calculated by (15)) are not very different, except for the interactions at the first coordination shell. This is most probably due to the missing multipole-multipole contribution to the intersite SCI in the ASA+M. One may also see from Table Table IV, that the SGPM-KKR-ASA+M interactions are in fact quite close to those of the CW-KKR-ASA+M obtained without multipole-multipole electrostatic interactions (0-L). These interactions have recently been used in Monte Carlo simulations of the ordering in NiPt and reproduced quite well the order-disorder transition temperature for an equiatomic alloy composition and the values of the SRO parameters in a random alloy at T=1200 K .

Finally, in Table V we show the ordering energies of the four equiatomic ordered alloys, obtained from the first 20 SGPM interactions, although the contribution from the interactions beyond the fourth coordination shell is only a few percent of the total ordering energy. Comparing these energies with those from the direct total energy calculations, presented in Table II, we find reasonable agreement for the ASA+M and the ASA results and very good agreement in neutral-sphere approach.

## VIII The SCI at alloy surfaces

The generalization of the SCI formalism to inhomogeneous systems, such as partially ordered alloys or surfaces, is straightforward. In the latter case the SCI also become inhomogeneous and, therefore, definition (15) should be rewritten as

$$V_{scr}^{\lambda \lambda ^{}}(R)=\frac{e^2}{2}\mathrm{\Delta }q_\lambda \mathrm{\Delta }q_\lambda ^{}\frac{\alpha _{scr}^{\lambda \lambda ^{}}(R)}{S},$$
(26)
where $`\mathrm{\Delta }q_\lambda `$ is the effective charge transfer on the $`\lambda `$ sublattice, and $`\alpha _{scr}^{\lambda \lambda ^{}}(R)`$ the screening constant, which is defined as

$$\alpha _{scr}^{\lambda \lambda ^{}}(R)=\frac{1}{2}\underset{L^{},𝐑^{}0}{}S_{L=00,L^{}}(𝐑^{})Q_L^{}^\lambda (𝐑^{}𝐑).$$
(27)
Here, the vector $`𝐑`$ connects the site on sublattice $`\lambda `$, where the perturbation of the charge density is induced to the site on sublattice $`\lambda ^{}`$, at which the SCI is determined (see paper I for details), and $`Q_L^\lambda (𝐑)`$ are the multipole moments of the screening charge in the atomic sphere centered at $`𝐑`$ normalized by $`\mathrm{\Delta }q_\lambda (R=0)`$.

It is clear that $`\alpha _{scr}^{\lambda \lambda ^{}}(R)`$ depends on the direction of $`𝐑`$, i.e., $`\alpha _{scr}^{\lambda \lambda ^{}}(R)\alpha _{scr}^{\lambda ^{}\lambda }(R)`$, and therefore the SGPM interactions, which in the inhomogeneous systems should be invariant under a sublattice index interchange, are defined as

$`V_{\lambda \lambda ^{}}^{SGPM}(R)`$ $`=`$ $`V_{\lambda \lambda ^{}}^{GPM}(R)+{\displaystyle \frac{1}{2}}(V_{scr}^{\lambda \lambda ^{}}(R)+V_{scr}^{\lambda ^{}\lambda }(R))=`$ (28)
$`=`$ $`V_{\lambda \lambda ^{}}^{GPM}(R)+{\displaystyle \frac{e^2}{4}}\mathrm{\Delta }q_\lambda \mathrm{\Delta }q_\lambda ^{}{\displaystyle \frac{\alpha _{scr}^{\lambda \lambda ^{}}(R)+\alpha _{scr}^{\lambda ^{}\lambda }(R)}{S}},`$ (29)
where $`V_{\lambda \lambda ^{}}^{GPM}(R)`$ is the usual GPM interaction.

As an example of the inhomogeneous system we have chosen three low-index fcc surfaces: (111), (100) and (110). The calculations has been performed by the EC-LSGF method with an inhomogeneous effective medium, which was fcc(111), fcc(100), and fcc(110) slabs consisting of 3 vacuum and 6 atomic layers, 4 vacuum and 6 atomic layers, and 6 vacuum and 10 atomic layers, respectively, parallel to the surface. The actual supercells for the impurity calculations were built on the basis of the effective medium supercells by $`N_x\times N_y`$ translations in the plane parallel to the surface, which were $`6\times 6`$ for the (111) and (100) surfaces, and $`6\times 4`$ for the (110) surface. In the case of the (110) surface $`N_x`$=6 was the period in the closed-packed  direction, while $`N_y`$=4 was the period in the  direction. Since the screening is insensitive to the alloy composition and the lattice parameter, the screening density was obtained for a Ni impurity in pure Pt at $`S=`$3 a.u..

In Table VI we present the on-site screening constant, $`\alpha _{scr}^\lambda (0)\alpha _{scr}^{\lambda \lambda }(R=0)`$, in the first three layers of these surfaces (the impurity is in the $`\lambda `$ layer). It is clear that the surface makes the screening more efficient, since the larger value of the on-site screening constant means a closer position of the screening charge to the impurity, and, in fact, $`\alpha _{scr}^{\lambda =1}(0)`$ increases from the most close-packed (111) to the most open (110) surfaces. One can also see that the surface influences the screening mainly in the first layer for the fcc(111) and (100) surfaces, and very little in the second layer of the fcc(110) surface.

The intersite screening constants , $`\alpha _{scr}^{\lambda \lambda ^{}}(R)`$, are presented in Table VII. Again, one can see that the surface has a quite substantial effect on the screening constants $`\alpha _{scr}^{1\lambda ^{}}(R)`$, which is due to the perturbed electron density in the first layer, $`\lambda =`$1. However, the values for $`\alpha _{scr}^{\lambda \lambda ^{}}(R)`$ are already very close to the corresponding bulk values (see Table III) for $`\lambda =`$2.

## IX Conclusion

The polarization of the electron density of the alloy components due to their size mismatch makes a substantial contribution to the electrostatic energy of the alloy. This contribution, which is missing in the pure ASA, may be accounted for in the ASA+M approach through the multipole-moment interactions in the Madelung part of the electrostatic problem, and, as we have shown, it plays a crucial role in obtaining the correct ordering energetics.

We have also demonstrated that the neutral sphere approach based on the use of different atomic spheres for the alloy components on the corresponding underlying lattice (or subllatice, in general) leads to unacceptable, quantitative errors, and therefore the only consistent way of obtaining correct configurational energetics is to use spheres of equal radii for the alloy components on the sublattice where the alloying is taking place.

Since the multipole moments due to polarization effects originate from the specific local atomic configuration around each site, they may in principle be accounted for only by methods which go beyond the single-site approximation in the electronic structure (Green’s function) calculations. However, the use of a simple parametrized form for the on-site Madelung potential and energy in the SS-DFT-CPA calculations still allows one to obtain a reasonably accurate description of the electronic structure (if the CPA works for a given system) and total energy, although, obviously, such a parameterization is possible only on the basis of the calculations by more accurate methods.

The monopole-multipole intersite SCI have been obtained for NiPt fcc bulk and surface alloys. The SGPM interactions, which are the usual GPM interactions plus the SCI, reproduce the corresponding monopole-multipole KKR-ASA+M results, which give a semiquantitatively correct description of the ordering in NiPt.

## X Acknowledgements

Valuable discussions with Dr. I. A. Abrikosov, are greatly acknowledged. Center for Atomic-scale Materials Physics is sponsored by the Danish National Research Foundation. S.I.S. acknowledges support from the Swedish National Research Council (NFR) and Materials Consortium "ATOMICS". The work of P.K. is financed by the SKB AB, The Swedish Nuclear Fuel and Waste Management Company.