Twenty five years of Finnis–Sinclair potentials

In 1984, the paper, entitled A simple empirical N-body potential for transition metals by M.W. Finnis and J.E. Sinclair appeared in Philosophical Magazine A 1. Twenty five years later, with more than 1500 citations, it is the most cited article in the Philosophical Magazine A series (1978–2002). In this Special Issue, to celebrate 25 years of Finnis–Sinclair potentials, we have invited papers that illustrate the range of simulations these potentials have enabled, as well as reviewing more recent developments in large-scale simulations in metals.

Why have ‘Finnis–Sinclair potentials’ had such an enormous impact? Prior to 1984, atomistic simulations of metals almost invariably used pair potentials to describe atomic interactions. In nearly-free-electron metals, such as the alkalis and aluminum, pair potentials representing effective interactions between screened ion cores could be derived from second-order perturbation theory. However, these potentials accounted only for ∼10% of the cohesive energy of the metal and could even be anti-cohesive, as in the case of aluminum; the remainder was in a density-dependent, structure-independent energy. The potentials themselves were also dependent on density. If there were large local changes in density, such as at a free surface or even at a vacancy in a high-electron-density metal such as aluminum, there was little justification for ignoring the density dependence of the potentials. The use of potentials dependent on the local density was considered problematic, not least due to the difficulty in rigorously defining a local density. In transition metals, where no theoretical justification could be given to pair potentials, they were nevertheless constructed by fitting them to experimental data. If all the cohesion is represented by a sum of density-independent pair potentials, then the Cauchy relation is obeyed in cubic metals: c₁₂ = c₄₄. This is not obeyed by any real metal – a fact that highlights the inadequacy of the pair-potential model. Although there were some earlier attempts to deal with the local density dependence of atomic interactions in metals, until 1984 almost all large-scale atomistic simulations in metals used pair potentials, many of which were constructed empirically.

Finnis–Sinclair potentials solved the local density dependence of atomic interactions in metals at a stroke. The binding energy of the solid is described by two terms. The first is a short-range repulsive interaction described by a pair potential acting between neighboring atoms – nothing new there. However, the attractive interaction is now described by the square root of a sum of pair interactions to neighbors, centered on each atom in the solid. The novelty is entirely in the square root. The square root captures the dependence of atomic interactions on the local density: as the number of neighbors of an atom decreases, the strength of the remaining bonds increases. This immediately predicts an inward relaxation at metallic free surfaces with a tensile surface stress, both of which are widely observed but not predicted by models in which the cohesive energy is just a sum of pair potentials. The breakthrough came as supercomputers were becoming accessible and, crucially, the computer time required to compute Finnis–Sinclair forces is similar to that for pair potentials, so there was no need to compromise system size if researchers adopted these potentials.

The origin of the square root may be found in the simplest tight-binding theory of metallic cohesion. The story begins with the concept of a local density of states introduced by Friedel in 1954 2. In 1964, Friedel 3 argued that the parabolic variation of the cohesive energy across the transition-metal 4d and 5d series may be explained by assuming it is due to the filling of the d-band, the density of states of which may be approximated by a rectangular shape. The attractive bond energy contribution to the cohesive energy then scales with the band width. Like all distributions, the width of the local density of states is characterized by the square root of its second moment. The second moment may be calculated in real space from the sum of squares of electron hopping integrals between the site of interest and its neighbors. As Olivier Hardouin-Duparc explains in this issue, this idea goes back to two independent studies of vibrations in disordered systems in 1959. Bringing these observations together, it is clear that the bonding contribution to the energy per atom varies as the square root of a sum of pair-wise functions between the atom of interest and each of its neighbors.

It appears to have been known through the 1970s to a community centered in Paris and Cambridge that the cohesive contribution to the binding energy per atom varies with the square root of the coordination number, z, rather than with z itself. In this issue, Hardouin-Duparc traces possibly the first appearance of this idea in print (in French) to François Ducastelle in 1970 4. When Volker Heine wrote down the same idea 10 years later in Equation (11.19) of his influential review 5, the idea had already become very well established at Cambridge. In August 2009, Heine wrote in a personal letter to Adrian Sutton:

‘I have looked up the papers from the 1970s in my personal reprint file and do not find it used there. Yes, the idea was known but not taken seriously enough to publish papers on cohesion and relative phase energy.’

Indeed, the limitations of the second-moment model, in which all higher moments of the local densities of states are ignored, were fully appreciated in 1971 in the powerful moments theorem of Ducastelle and Cyrot-Lackmann 6. Their theorem showed that the relative stabilities of competing crystal structures depend on higher moments of the density of states than the second.

Mike Finnis and Jim Sinclair took a pragmatic approach and fitted their potentials to experimental data. By extending the range of the potentials to beyond the first neighbors and through the empirical fitting process, they were able to circumvent the fundamental limitations of the second-moment model concerning structural stability. Although their potentials are many-body in nature, they are not directional, and the force acting between two atoms acts along the vector separating them. The modern view 7,8 is that directional bonding is significant in transition metals, and that the non-directional Finnis–Sinclair approach is better suited to noble metals. The square root of Finnis–Sinclair potentials is equivalent to the embedding function 9, glue 10 and effective medium 11 of similar potentials. The physical picture underlying those potentials is quite different from the second-moment model, but that is of no practical significance when the potentials are constructed by fitting to empirical data.

Mike Finnis has provided the following account of how the Finnis–Sinclair potentials came into being in 1983:

Since I finished my PhD with Volker Heine in 1974 I had been employed in the Theoretical Physics Division at AERE Harwell. An important mission was to carry out research in support of the nuclear power industry and it enjoyed a tradition and a reputation for developing an understanding of materials at the atomic scale, exploiting the most powerful computers for simulation. When I started we had access to an IBM360, and soon afterwards a 370 and a Cray1. John Hubbard was calculating band structures, Mike Norgett was developing and using empirical potentials for ionic materials, Ron Bullough (leader of the Group to which I was allocated) was modeling structural evolution of metals under irradiation, Roy Perrin was modeling deformation of fuel pins in reactor cores, and Marshall Stoneham was modeling everything else! Suffice it to say that nowadays that Division would have had multiscale in its title. Among the 30 or so theoreticians there, Jim Sinclair, two doors up the corridor, talked with me quite often about his multiscale modeling approach (published in 1978), in which an atomistic region was embedded in a continuum elastic region by means of a lattice Green function, and how nice it would be to have a ‘toy’ model, an interatomic potential whose predicted elastic constants could be varied in an arbitrary way by turning knobs on the potential, i.e. by varying parameters. That is not possible of course if the Cauchy constraint of pure pair potentials has to be satisfied. At the time there was no such model.

I knew about tight-binding models from my PhD days, when Volker Heine, Roger Haydock, my one-time office mate Chris Nex and my contemporary as a PhD student Mike Kelly, were talking about them and about the newly invented, or rediscovered, recursion method. They had taught me that the total bonding energy of a transition metal scaled rather with the square root of the coordination z than linearly with z. More to the point, this is a convenient way to say that the energy scales as the square root of a sum over bonds of functions that could be identified as hopping integrals squared. I had also seen Francois Ducastelle's PhD thesis in the seventies with similar ideas. So that much seemed to be common knowledge. Although the limitations of pair potentials were clear, until the CECAM Workshop in Paris in 1983 I had tended to think that the next best thing would be full, real space tight-binding rather than the ‘second-moment approximation’ that underlies the Finnis–Sinclair approach. The empirical ‘effective medium’ or ‘embedded atom’ approaches, which end up looking formally similar to our potentials, were also unattractive to me at the time, and I didn’t really understand their origin until I had carefully read Jacobsen, Norskov and Puska's physical derivation of them, published in 1987. However, Jim's needs for a toy model, and discussions at the CECAM Workshop, as well as simply being away from my desk, must have triggered me to work out the consequences of the second moment model with analytic functions that could be fitted to data. Since I was staying in a quiet, borrowed apartment in the Rue du Bac, and CECAM Workshops lasted longer in those days, I had plenty of time to derive expressions for elastic constants, the Cauchy relation and unrelaxed vacancy formation energies. When I got home I discussed these with Jim, who proposed using cubic splines to parameterize the functions, and wrote some code to generate them. I then fitted the elastic constants and lattice parameters to the simple BCC transition metals and we wrote the paper.

We called it an empirical N-body potential because the functions were completely empirical, the only physical content was the square root. It was N-body rather than three- or four-body, which would have been an alternative way to go beyond the pair-wise description.

Even more empirical, although ultimately more useful in practice, was the next step, to fit parameters also for the fcc metals Cu and Ni. This seemed like a reasonable project for a student, and there was a candidate from the Physics Department in Oxford who took it on as my first PhD student, Graeme Ackland. He was successful, and quite a lot of people since have been generating similar toy models to describe metallic bonding.

A year or two later, a discussion with Vasek Vitek and Graeme pointed me to a troublesome paper concerning self-consistency. In metals, self-consistency tends to enforce local charge neutrality, and in tight-binding one can assume that the parameters should be adjusted to do this. Simple tight-binding models had been published in which local charge neutrality was achieved by varying the on-site Hamiltonian matrix elements. Other publications ignored self-consistency altogether. I showed that, fortunately, our simple square-root model for the electronic energy retains what validity it has if we assume that, together with local charge neutrality, the density of states varies from atom to atom by being stretched or compressed along the energy axis, according to the local coordination. If these effects – local charge neutrality and the local change in the density of states – were neglected, then the electronic energy in a distorted metal could no longer be represented by the same function everywhere.

Jim Sinclair has added the following:

Although I have been out of this field for many years, I can add just a few recollections. First, I can reinforce the memory of how fertile an environment was provided by the Theoretical Physics Division at Harwell in those years. Other employees sometimes saw us through the window taking hours over our tea breaks, and complained what a lot of time-wasters we were. But the fact was that most of those long conversations were really about work. Or about some crazy idea that someone had thrown out for discussion, and which might well inspire someone to have a useful line of investigation later. If it wasn't in the tea-room, you could easily wander along the corridor and find several open doors and open minds willing to give their constructive ideas on your latest difficulties.

It was on such a mission that Mike Finnis wandered into my office to toss about some of his ideas on transition-metal potentials. Neither then nor since have I had any expertise in solid state physics or quantum mechanics, so I think I must have tried to ignore the detail in Mike's description of his motivations, and just tried to concentrate on the practical difficulties. Given this limited basis to my input, I thought then that Mike had been over-generous in adding my name to his eventual publication – and even more so now that the paper has achieved its remarkable status. Whatever my contribution to the discussion was, it must have appeared like a significant spark to Mike for him to consider me a co-inventor of what followed.

As I recall, my input arose from my experience of using valence-coordinate potentials to represent diamond or similarly structured semiconductor crystals in defect modeling calculations. I had seen valence coordinates (bond-lengths, inter-bond angles) used in harmonic potentials fitted to phonon dispersion data, and had invented a way – purely empirical – to extrapolate such harmonic potentials to cover large displacements. The idea was to modify the terms for changes in bond-lengths and inter-bond angles by substituting non-linear functions whose gradients faded to zero as the bond lengths increased. The modifying functions could be made partly respectable by fitting the potentials to data such as sublimation energy or surface energy.

I believe, then, that my sole contribution to the genesis of the Finnis–Sinclair paper was suggest to Mike that, if he wanted pair-wise interactions that varied in strength with the density of a solid, he could use for the effective density of the environment of one atom a sum of terms depending on the distances to all its neighbors within a certain range. Taking this germ of an idea through to the much more soundly based FS potential and its progeny was entirely up to Mike.

The work Mike refers to in his final paragraph above was published in 12. As Mike explains above, the functional form of the potential in the second-moment approximation remains the same in distorted atomic environments if, and only if, metallic screening is taken into account. In addition, it revealed the crucial role of metallic screening in ensuring that the square-root functional form in the second-moment approximation is independent of band filling. These are subtleties in the second-moment model that have been largely overlooked, although the importance of screening and self-consistency play a central role in the more complete tight-binding bond model 13.