How does one extract many-body interatomic potentials from ab-initio band structure calculations

An original approach to the ab-initio deriving many-body interatomic potentials in metals is discussed. It is based on calculating the Kohn-Sham total energy functional using the pseudopotential method. The local density approximation (LDA) is shown to be applied, within the pseudopotential concept, to the analysis of the valence electron kinetic energy. Utilizing the LDA and linear superposition assumption for total electron density enables us to treat the exchange-correlation and kinetic energies of valence electrons in terms of contributions to indirect many-body interactions. Equations for the pair and for the triplet potentials in metals are given. Relationship between the method developed and other approaches is analyzed.


Introduction
The paper concerns a problem of derivation from the first principles multi-ion interparticle potentials in metals.The approach suggested is based on the densityfunctional and pseudopotential versions of the Kohn-Sham theory [1,2].The focal points of the method developed are as follows: a) an accurate real-space representation of the metal total energy in terms of widely transferable, structure-independent interatomic potentials [3]; b) treatment of the valence electron kinetic energy, determined through the abinitio band structure calculations in terms of appropriate contributions to the manyion indirect interactions [4].
There are two motivating factors behind this work.First, calculating the inter-atomic potentials is one of the central problems in the modern condensed matter physics having a long history, see review of literature in [3][4][5][6].A role of angular forces, caused by the multi-ion interactions, is the subject of special interest in describing metal properties microscopically [7].
A second motivating factor in this work is the prospect of doing molecular dynamics (MD) simulation of various metal properties having derived potentials [5,7].Comparing the results of MD simulation with the corresponding experimental data we get a convincing verification of a theoretical background used at the interatomic potential evaluation.

History of the problem
One can distinguish several trends in deriving interatomic potentials in metals.2.1.The rigorous development of interatomic potentials from quantum mechanics is based on the pseudopotential method [8,5,6].Working within the framework of the perturbation theory in pseudopotential one can represent metal total energy in the form of expansion in the effective multi-ion interactions [5][6][7][8]: The first term E 0 depends on atomic volume Ω 0 only and includes the energy of the interacting electron gas as well as all one particle intraatomic contributions to the E tot .The interatomic potentials V 2 , V 3 etc.are implicitly volume dependent but explicitly structure independent and thus rigorously transferable at a given volume Ω 0 to the all bulk structures, either ordered or disordered.To put it another way the potentials themselves are independent of the absolute ion positions R i and depend only on relative separations R ij = |R i − R j | [6,7].For example, the angular-force triplet potential V 3 is the three-dimensional function: [7].The prime on each summation in (1) denotes exclusion of all the self-interaction terms where two indices are equal.The interatomic potentials V 2 , V 3 etc.present on the right hand side of equation ( 1) are completely transferable at a fixed volume because of their independence on metal structure.The structure dependence of the total energy (1) appears through the summations in (1) over all N ion positions.One should stress the important feature of the equation (1) obtained within the pseudopotential theory.Expansion (1) is rapidly convergent in the sense that the four-ion quadruplet potential V 4 is smaller than V 2 , V 3 and higher potentials (V 5 , V 6 ...) appear to be negligible [6,7].2.2.The embedded atom method (EAM), see [9] and literature cited there, is a semi-empirical approach.It is based on the assumption that the total energy of a system can be represented as the sum of two parts: one of them having its origin from the band structure is a function of the electron density and the other is a short-ranged pairwise potential describing interactions between the atomic cores.Thus, The embedding function F (ρ), the pairwise potential ϕ(R ij ) and the charge density ρ i at the site R i are modeled in the EAM.For example, the authors of [9] proposed the following forms for ϕ(r) and ρ i : where R 1 and R 5 are the nearest and fifth neighbour distances, respectively, K n , α and ρ 0 are parameters to be determined [9].The effects of the electron density Friedel oscillations are introduced in ρ i through the function f (r) [9] Here K F is the Fermi vector, ∆ and δ are free parameters.One can notice that there are a lot of fitting parameters in the EAM.Besides, the introduced forms for F (ρ i ) and ϕ(R ij ) modelling functions demand theoretical justification.2.3.The "pure phenomenological" approaches were exploited formerly [10][11][12] with representation of the total energy entirely in terms of empirical pair potentials.One should put E 0 = V 3 = V 4 = ... = 0 in equation ( 1) to compare the ab-initio and phenomenological approaches.Arbitrary forms for V 2 (R ij ) containing a number of adjustable parameters are used typically [10].Such potentials are implicitly structure dependent and hence their transferability is always strongly doubted.The authors of [11] tried to determine a pair potential more systematically by inverting the energyband calculation of E tot as a function of volume Ω 0 .To put another way round they obtained V 2 ≡ V 2 (R ij , Ω) [11].Some attempts were also made to extract a pair potential from experimental phonon spectrum [12].
In some cases an adequate theoretical description of metal and alloy properties could not be achieved within the pair interatomic interaction model.It means the many-body interactions (three-, four-particles) play essential role [13].Their being taken into account is especially important for transition metal investigations [7,14].
This brief review of literature shows clearly that the ab-initio approach is the most consistent and promising on the problem of interatomic potential evaluation.A somewhat different way of deriving the effective interatomic potentials in metals within the ab-initio approach has been advanced in [3,4].It is based on solution of the Kohn-Sham equations as well as utilizing the pseudopotential concept and the local density approximation [3,4].The approach developed has the following features: a) the Kohn-Sham total energy functional is treated in terms of the multi-ion interactions; b) the many-body potentials are evaluated without using the perturbation theory in pseudopotential.

Many-body potentials in metal
Let us examine an electron-ion system within the pseudopotential concept.It means the core electrons are excluded from the explicit consideration in the all electron problem and the true electron-ion potential is replaced by a pseudopotential [15].Thus, the pseudoions are placed in lattice sites R and their interactions with the valence electrons are described by a pseudopotential w(r − R) which is an operator in a general case.
We start from the Kohn-Sham total energy functional of such a system [2] Here T [ρ] and E xc [ρ] are functionals of the kinetic and exchange-correlation energies, respectively.The second and the third terms in ( 6) describe the energy of the electron subsystem in the external field of the pseudoions and the Hartree energy, correspondingly.The last term E i−i in ( 6) is an energy of the ion-ion direct interactions.The total energy E[ρ] (6) is the universal functional of an electron density with Ψ k (r) the valence electron pseudowave functions and k the wave vectors characterising electron states.Summation in ( 8) is performed over all occupied states.One can find the pseudowave functions Ψ k (r) solving a system of the Kohn-Sham equations with the pseudopotential [4].As a result the pseudodensity ρ(r) (8) of valence electrons that minimizes the total energy functional E[ρ] (6) is known.The pseudodensity differs from the true density of valence electrons ρ v (r) in the regions of ion cores and coincides with it outside cores.Unlike ρ v (r) the pseudodensity is a smooth function in the ion core region.The importance of this circumstance will be demonstrated below.
In [3] the known quantity of ρ(r) (8) has been presented in the form where ρ 0 = z/Ω 0 (Ω 0 = Ω/N) is a density of the uniform electron distribution, z an ion valency, Ω 0 the atomic volume and ρ i (r − R) is an electron density related to the pseudoion at the site R.
The focal point of the linear superposition assumption ( 9) is an accurate realspace treatment of the metal total energy in terms of well-defined interatomic potentials [4].The local density approximation (LDA) [2] and equation ( 9) makes it possible to represent the T [ρ] and E xc [ρ] functionals in terms of appropriate contributions to the many-body interatomic interactions and get the total energy in the form of equation (1).Let us prove this statement.
The LDA is almost universally used for the E xc [ρ] functional in total energy pseudopotential calculations [6,7,16] with ε xc (ρ(r)) the exchange-correlation energy per electron at point r which is equal to the exchange-correlation energy of the homogeneous electron gas possessing the same density ρ(r) as the electron subsystem under investigation.Consider the kinetic energy functional T [ρ].Within the pseudopotential concept the unknown Ψ k (r) functions of the valence electrons are searched on the plane-wave basis [5,6,8], that is Here a G (k) are the expansion coefficients and G, the reciprocal lattice vectors.Solution of the Kohn-Sham equations provides information on coefficients a G (k) as eigenvectors of this problem.The valence electron kinetic energy with allowance for ( 11) is equal to The G M vector is determined by the cutoff energy |k + G M | 2 /2.Let us represent the known quantity T (12) as follows that is in the form of LDA like E xc [ρ] (10).One can regard (13) as a constraint on the unknown function t (ρ(r)).
In the case of simple metals the following form could be used for the t (ρ(r)) function [17] The value C = 1 corresponds to the homogeneous electron gas LDA expression.We propose to treat C as a free coefficient fitting the right hand side of equation ( 13) to the valence electron kinetic energy value (12).The formula (14) has a physical justification.The valence electrons behave like free electrons in nontransitive metals [5,8,13,15,16] and the Fermi-surface is close to the spherical one.That is why, in simple metals the kinetic energy of valence electrons can be approximated by functionals like (14), see [18] and literature cited there.The results of [17] confirm this statement.Certainly, the form ( 14) fails in the all electron problem when both the core and valence electrons are considered, see review of literature in [18].There are convincing reasons for explaining this failure: behaviour (distribution) of the core electron density ρ c (r) does not satisfy the conditions of the Hohenberg-Kohn formalism applicability [1] (ρ c (r)-function is not smooth).
The pseudopotential method permits to avoid these difficulties.The core electrons are excluded from the explicit consideration but the true electron-ion potential is replaced by a nonlocal operator that describes an effective electron-ion interaction [5,8,15].The pseudoelectron density ρ(r) (8), see (11) is used in the total electron density calculations within the pseudopotential method [16].It satisfies the conditions of the Kohn-Hohenberg theory applicability better than the true electron density because the potentials generated by ρ(r) are much smoother than those caused by the true density [1].We should emphasize that the approximation ( 13) with ( 14) is used just for interpreting the valence electron kinetic energy in terms of appropriate contributions to the many-ion indirect potentials in (1).Such an interpretation is the most natural within the pseudopotential concept.
Let us expound briefly the main ideas of the approach advanced in [3,4].The regular metal of the valency z with one atom per primitive cell is considered.The functions t(r) and ε xc (r) are periodical ones in the regular crystal because of the electron density ρ(r) periodicity.Therefore, one should consider them in the region of the Wigner-Seitz cell only.Let us define as the electron density of a pseudoatom at the site R 1 and expand t (ρ(r)) and ε xc (ρ(r)) in ρ ps (r − R 1 ) for r ∈ Ω 0 (1).Symbol r ∈ Ω 0 (1) denotes that the running variable r takes values in the unit cell (the Wigner-Seitz cell) centred at the site R 1 .
For function f (ρ(r)) = t (ρ(r)) and f (ρ(r)) = ε (ρ(r)) such an expansion reads with r ∈ Ω 0 (1) and It is seen from ( 15) and ( 9) that the total pseudoelectron density ρ(r) in the vicinity of each atom is expressed as a sum of the density ρ ps contributed by the pseudoatom in question plus contributions from the surrounding ions that is The condition provides convergence of the series for ρ(r)t (ρ(r)) and ρ(r)ε (ρ(r)) [4].Thus, it is the main trick to represent the known total electron density ρ(r) in each unit cell Ω 0 (j) in such a manner (see equation (18)) that the condition (19) would take place in the whole Ω 0 (j) region [3].
Owing to the ( 16) and ( 18) one can treat the E xc [ρ] and the T [ρ] in terms of contributions of the exchange-correlation and the kinetic energy effects to the manybody interatomic indirect interactions [4] and represent the sum and do not depend on ion configuration.Functions V ind 2 (R 1 , R 2 ), V ind 3 (R 1 . . .R 3 ) . . .describe indirect interactions between two, three. . .ions caused by the band structure and exchange-correlation effects: where One can derive explicit expressions for the indirect many-ion interactions in metal using equations ( 13), ( 16), ( 18) and (20) and C m n the binomial coefficients.The following notations are introduced in (24) to (26): The first summand in (24) and the first term in (25) represent irreducible indirect pair and triplet interactions, respectively, in metal.The next summands in ( 24) and (25) are contributions to the ) arising from the n-particle potentials (23) when (n − 2) or (n − 3) indices of ion coordinates coincide accordingly.
They should emphasize that the present stage of the metal microscopic theory permits one to estimate only the reducible terms of the third and forth order in pseudopotential contributing to the V ind 2 (R 1 , R 2 ) [13].The problem of correct calculation of irreducible contributions to the V ind 3 (R 1 . . .R 3 ), V ind 4 (R 1 . . .R 4 ) has not been solved completely yet.Unlike the perturbation theory in pseudopotential the approach developed enables us to work out a general procedure of renormalizing the n-ion potential by the (n + 1)-, (n + 2)-, etc. reducible interionic interactions.An elegant analytical form for such a renormalization is obtained in the cases of the pair and triplet potentials, see equations ( 24) and (25).That is why the proposed method looks like a promising one.

Conclusions
A scheme of deriving the effective interatomic potentials is presented.Information on the total valence electron density ρ(r) is needed only.The ρ(r) expansion in contributions of the individual atoms is the essential point of the method suggested.The linear superposition assumption (9), see also (18), for the total electron density enables us to get the analytical expressions for the indirect interactions.The n-particle potentials (n > 2) are formed by the kinetic (band structure) and exchange-correlation effects only.Equation for the total pair potential is given in [3].Evaluation of the derived potentials as well as utilizing them in metal property investigations are highly desirable for providing a convincing verification of the approach proposed.
these terms of series for T [ρ]/N and E xc [ρ]/N which depend on the coordinates of n-ions.