History

My knowledge of the development of the MT-based MSW programs is fairly sketchy, deriving mostly from what I've heard by word of mouth, or snippets of information scattered through pieces of program documentation. It appears, however, that some original code from the work of Smith and Johnson passed through the hands of various programmers in several places (particularly IBM San Jose, and Stanford University in California) developing into a package consisting of two main programs, ENERGY and SCF, together with some routines to prepare an initial cluster muffin-tin potential from tabulations of free-atom electron densities, and programs to construct linear combinations of spherical harmonics at the atomic sites, which transform according to the relevant irreducible representation (irrep) of the cluster's point symmetry group.

 

The ENERGY program takes an initial potential and finds the energy levels so that they can be populated with the available electrons. All this information is passed to the rather more formidable SCF program to iterate the potential to self-consistency, at each cycle generating a new cluster electron density and projecting the resulting potential onto the muffin-tin form.

 

The application of the MT MSW method to XANES begins with the development of a version of the codes able to calculate positive-energy scattering states in the cluster potential. This was the famous CONTINUUM program of Rino Natoli, which appeared about 1980. It has been kept in development since that time and has become, in one form or another, part of several popular packages (most notably in recent years MXAN which is beginning to build up a considerable user base).

 

The actual history of the MT MSW codes up to the early 1980s is much more extensive than I have indicated here. The bound-state codes alone have been used for many calculations that have reported in the literature, occasionally implementing extensions to the basic theory (including, for example, partitioning within the cluster or molecule).

 

My own small part began in September 1983, when I joined the EXAFS group at the Hamburg Outstation of the EMBL (European Molecular Biology Laboratory) to begin my PhD work under the supervision of Robert Pettifer, as an external student of the University of Warwick. While Robert and our colleague Christoph Hermes devoted their energies to the experimental side, primarily beamline development and measurements for particular systems, frequently in collaboration with our many visitors (and with occasional help of the opposing-thumb variety from yours truly), my own task was the improvement of ab initio XAFS calculations. Although tangential to our subject, it is worth recalling here that Robert and Christoph produced some remarkable work (including absolute energy calibration of XAFS spectra, luminescence EXAFS, double crystal monochromator etc).

 

Those involved in ab initio XAFS calculations at the time will recall that the really hot topic was the inclusion of multiple-scattering effects in the calculation of EXAFS, for which, typically, only a single-scattering treatment was usually made. Evidence had been accumulating that for certain types of systems some non-collinear scattering paths could contribute significantly to the EXAFS signal. Initially, however, we were interested in calculations for the EXAFS region using standard methods with a view to experimenting later with multiple-scattering effects. Typical of the systems that we studied at the EMBL were transition-metal centres in biologically important systems like enzymes. We therefore looked at various model compounds for which good XAFS spectra were available, and which were analogous to these kinds of systems, at least from the point of view of physics. The most useful of these was chromium hexacarbonyl, Cr(CO)6, which we returned to several times. (Curiously, the chromium K-edge turned out to be outside the wavelength range of our beamline (EMBL's D2 in HASYLAB) in its then configuration, and a spectrum was kindly provided by Jürgen Röhler and his group using the Kiel Univ. beamline.) For our initial investigations Robert had provided the phase-shift calculation program MTAHARA (which included as an option the then new Dirac-Hara exchange potential) and the fitting program EXAFSFT he had developed earlier with A D Cox.

 

Looking more deeply into the calculated phase shifts and the resulting scattering amplitudes, f(θ,k), as a function of scattering angle and photoelectron momentum, we noticed some problems in the form of spurious oscillations. Fairly quickly these problems were traced to the muffin-tin approximation, more particularly the potential discontinuities which arise at the atomic-sphere surfaces. Given that XAFS is an oscillatory phenomenon, spurious oscillations of non-physical origin can hardly be welcome. That such effects were intrinsic to the method of calculating the scattering behaviour of atoms in our clusters was obviously a serious problem.

 

At this point a visitor to our lab, Lucia Incoccia, told us of some theoretical work being done on the non-muffin-tin multiple-scattering problem by some of her colleagues at the Frascati laboratories of the INFN (Istituo Nazionale di Fisica Nucleare) in Italy, and was kind enough to put us in contact with them. This was, of course, Rino Natoli and his then post-doc Maurizio Benfatto. As a consequence, I was lucky enough to be able to spend a month with them in May 1984 to learn about the MSW method; both the old MT version, and the new non-MT theory which was later to appear as the article by Natoli, Benfatto and Doniach (1986). Out of this initial visit developed a fruitful collaboration which has lasted at various levels of activity over the years.

 

Initially it was decided that I should devote my attention to an implementation of the non-MT MSW theory, with advice from Rino and Maurizio at a distance, and from Robert locally. The first task was to become familiar with the non-MT theory and also details of the existing codes, particularly the main MSW programs ENERGY, SCF and CONTINUUM. One useful byproduct of this familiarity was a handful of modifications to the codes, for example, in the construction of the MSW secular matrix (leading, in the case of SCF for one particular system, to a speed increase of a factor of 3), as well as improvements aimed at greater numerical precision in special function routines. A great many sample calculations were also performed with the MT codes for model systems, notably Cr(CO)6 again, and ZnSe. Some of this material was presented at the conference Progress in X-ray Studies by Synchrotron Radiation (1985, Strasbourg), as a poster highlighting the effects of the muffin-tin approximation on multiple-scattering based calculations of XAFS generally. (Although XAFS was only one of a variety of subjects represented at the conference, there were several particularly interesting presentations, the most impressive to my mind being that of J E Müller who demonstrated spectacular agreement with experiment for (I think) the K-edge of palladium metal based on a band-structure type code (see Müller and Wilkins (1984)). Talking with him later we learned that he had written the code from scratch, in assembly language, for a Cray - which should provoke a frisson of horror in the programmers amongst us!)

 

Beginning work in earnest on the non-MT codes the first task was to develop the routines for modelling parts of a given cluster potential in multipoles. It was decided fairly early that the codes should be written to cope with an arbitrary molecular geometry rather than implement some quick test codes for a particular system. It was also decided to follow the old MT codes in using point-group symmetry to reduce the computational burden. In developing this approach a considerable amount of work was done to implement a variety of routines for numerical integration over spherical surfaces, manipulating spherical harmonic expansions and point-group irreps. A lot of effort was put into making spherical harmonic and related routines robust for high l values. It was also necessary to get a feel for the requirements of the potential modelling parameters (l truncations, atomic sphere radii etc) for a desired accuracy of the model.

 

Some of this material was presented to the EXAFS and Near Edge Structure IV conference (1986, Fontevraud, France) (see Foulis, Pettifer and Natoli (1986)). At the same time we also presented some work which was a direct continuation of our original purpose, i.e. to include multiple-scattering paths in the standard analysis of EXAFS for a biologically relevant system (Pettifer, Foulis and Hermes (1986)). Using phase shifts derived in a way which avoided MT problems, and using Robert's implementation of a fast EXAFS algorithm, we calculated EXAFS for a model compound containing Zn coordinated by rigid imidazole rings, reproducing (for the first time I believe) the famous "camel back" feature which is sometimes seen in the spectra of analogous biological systems (such as carbonic anhydrase). (In regretful hindsight this is a good illustration of the rule of thumb which says never to publish important work first in a conference proceedings - since we never got round to publishing it more fully anywhere else!) (As an aside: this was a very useful and enjoyable conference in several ways; multiple-scattering effects and ways of including them in EXAFS calculations were much in evidence in both private and public discussions. The conference centre itself was the ancient Abbey of Fontevraud, whose significance for us Brits I only discovered by chance: wandering around the old cloister I came upon a chapel being renovated, in the centre of which were three sarcophagi ... those of Henry II of England, his wife Eleanor of Aquitaine, and, of course, their son Richard the Lionheart!)

 

With the initial phase of the non-MT MSW coding work complete I moved down to Frascati in Spring 1987 to work directly with Rino and Maurizio on the development of the main multiple-scattering codes. During a painfully intense period of activity up to Summer 1988 several major problems were overcome to complete the program development and testing, before finishing off with some calculations for our main model system, Cr(CO)6. Surprisingly, the development of the coupled radial Schrödinger equation integration codes (based on the matrix Numerov method), for the atomic T-matrices, turned out to be almost straightforward, if very complicated. The basic algorithm was outlined in the Natoli, Benfatto and Doniach (1986) article and (from numerical evidence at least) is very accurate and robust.

 

More problematic was the treatment of the interstitial region. Our initial idea was to "discretize" the interstitial volume and transform its Lippmann-Schwinger equation into a linear system which could then be solved (to give the appropriate term in the non-MT MSW secular equation) by standard methods (and, more importantly, with existing software). A number of technical difficulties with this approach, together with time pressures, led us to shelve it and use instead the first Born approximation which required "only" the evaluation of integrals of smooth functions over the interstitial volume. A volume integration method was devised based on manipulating spherical volume integrals. The code section for this part of the calculation is quite complicated nevertheless and heavily CPU and storage intensive, and proves to be the main bottleneck for the whole computational implementation of the non-MT MSW method.

 

Two main programs were developed based on the above approach: ENESHX and CNTSHX for bound and continuum states respectively (corresponding roughly to the earlier mentioned ENERGY and CONTINUUM with added "Spherical Harmonic eXpansions"). It was decided at that stage to defer the development of a full-scale self-consistent field (SCF) program since it would be too time-consuming. Thus all calculations for real systems would be based on potentials derived from overlapped free-atom (OFA) electron densities (as was then most often the case for XAFS calculations).

 

Before using the FPX codes for realistic systems it was necessary to find some way of testing unambiguously that they were correctly coded, since the non-MT algebra is quite complicated especially when symmetry is used. Thus we tried to find some simple, but non-trivial, tests. For bound states the hydrogen molecular ion, H2+, provides a particularly stringent test because results from the codes can be compared with the well-known, exact, analytic solutions. In fact, one may model various aspects of the H2+ system and the results are quite convincing. For continuum states and photoionization we could not use this system since positive energy analytic solutions appeared not to be available (although this may have changed since then). It is possible, however, to use an atomic, hydrogen-like system to make a meaningful test. (For minor technical reasons we used the Li2+ system.) Analytic photoionization cross sections are well known for this system. To test the codes one must first "break" the potential by imposing a MT form using empty spheres and a spurious interstitial region. The potential is "repaired" by allowing the potential multipoles in each of these empty regions to be taken into account. Looking at both the bound and continuum-state tests more closely one finds that overall the non-MT method as implemented in the FPX codes fixes the bulk of the errors due to the MT approximation. Indeed the coupled radial Schrödinger equations and resulting atomic T-matrices are essentially exact, implying that the remaining error is due to the Born approximation for the interstitial region calculation. One notes, however, that this latter part provides most of the correction typically for electronic states which are quite important (bonding orbitals for example).

 

It should be said that one or two of my colleagues disagree with this last point (at least as a general tendancy for real systems), being of the opinion that only taking into account the higher multipoles of the potential on the atomic spheres would be enough to give most of the non-MT correction (in the context of XANES at any rate). I have to say that my experience to date leads me to disagree strongly. Certainly for the analytic tests there can be no argument. My experience of calculations for real systems so far (for example, Cr(CO)6 which results have been reproduced in the Powerpoint file mentioned below) has always shown that the interstitial region contribution is the largest. On the other hand I note that many colleagues use MT-based MSW codes to calculate XANES and allow the atomic spheres to overlap (increasing radii by as much as 20% to 30% from just touching). This sometimes appears to improve agreement with experiment, probably because the volume of the interstitial region is less and the discontinuities at the atomic sphere surfaces are smaller. Strictly speaking the overlap of atomic spheres is not allowed in the MSW method as certain expansions may not converge. However there have been some attempts in the literature to justify this ad hoc measure on physical grounds. Although such a measure should be anathema in ab initio calculations (not least because it introduces adjustable parameters into the situation), some plausible justification may be found in the observation that a large part of the potential variation (or anisotropy associated with a particular atom) that is lost in the MT approximation, tends to occur just outside the atomic sphere surfaces when these are just touching, and therefore that increasing radii will pick up some of the effect of this variation.

 

With the codes written and tested it was time to turn to some calculations for a realistic system (Cr(CO)6) and thus complete the first phase of the non-MT MSW work. When this was all done I returned to the UK to complete my PhD (manually, and excruciatingly, typing up my thesis (see Foulis (1988)).

 

I was fortunate to be able to return to Italy for a short post-doc in 1989 to do some further work on the codes with Rino and Maurizio. During this time we published the tests of the codes (Foulis, Pettifer, Natoli and Benfatto (1990a)). The results for Cr(CO)6 were presented to the 2nd European Conference on Progress in Synchrotron Radiation Research (1989, Rome) (see Foulis, Pettifer, Natoli and Benfatto (1990b)). As far as we know these calculations were the first non-MT MSW XAFS calculations. (Of related interest at this conference was the presentation by John Rehr of his method for curved-wave muffin-tin based multiple-scattering EXAFS calculations using a separable representation of the free propagator, which became the basis for the popular FEFF code.) Also around this time Rino and colleagues (Natoli, Benfatto, Brouder, Ruiz Lopez and Foulis (1990)) presented a large theoretical work on a multi-channel extension of the non-MT MSW method. Given that the non-MT MSW method is already complicated this seemed at the time a rather formidable challenge to implement, but I have been interested recently to see that some progress has been made in this direction (see P Krüger and C R Natoli (2004)).

 

When one considers in detail the Cr(CO)6 results several issues come to light. In the absence of SCF electron densities we used the OFA approximation, using several different final-state potentials, both neutral and ionized, with variants of the so-called "Z+1" appromixation to model the 1s core hole. Although a very crude way of moving some electric charge around, it served to reveal that the full non-MT calculations were indeed sensitive to details of the potential that the MT approximation washed out almost entirely. In fact the MT results (as well as those including higher multipoles of the atomic potentials) are almost identical except in the region very close to the absorption edge. More importantly the full-potential (non-MT) results were still some way from good agreement with experiment. Plainly there were other things lacking from our treatment, and the main suspects were the lack of SCF densities (including using a final-state density suitably relaxed around a core hole), a better treatment of exchange, and the inclusion of some treatment of inelastic loss.

 

In 1990 I spent a useful and pleasant year at LURE (Laboratoire pour l'Utilisation du Rayonnement Electromagnétique) near Paris, working with Pierre Lagarde and Anne-Marie Flank, using the codes to tackle some problems of interest to XAFS experimentalists there. One interesting project involved calculations for small SiOx clusters in rare gas matrices (Foulis and Flank (1991)). While the results were interesting and suggestive, it became clearer as time went by that it was necessary to deal with some of the outstanding remaining problems of the potential model, before being able to tackle effectively some real problems. Highest on the list of targets was the use of SCF densities.

 

An opportunity to treat these problems came up in 1991 with a 3-year post-doc back at the Warwick University Physics Department working with Robert Pettifer again, with a general brief to improve aspects of the potential model and theoretical treatment of the MSW problem.

 

As it turned out the first job to be tackled was that of finding a better treatment of exchange. Until this point we had used the well-known Xα form for the exchange potential, this being the standard approach in EXAFS calculations at the time, as well as MSW XANES calculations. Not being experts in this area we were content to experiment with various approaches that had been developed by others for XAFS as well as related fields such as LEED and electron-molecule scattering. These included such things as the energy-dependent Dirac-Hara potential, and the real part of the Hedin-Lundqvist potential (the complex form of which is now a popular choice for the EXAFS region), as well as simple approaches such as merely using Xα with a reduced α (something whose effect had been noted before by other workers). We also decided to base our survey on results for a well-characterised atomic system (the K-edge of argon in this case) using our cluster codes; relying on the logic that if the results are not good for an atom, they can hardly be expected to be better for a molecule. From our many hundreds of calculations it appeared that the Hedin-Lundqvist potential was the best overall, while, surprisingly, the simple reduced-α potential (with an a priori value of 0.55) turned out to be reasonably good at least in the near-edge region.

 

At this time we also looked at the basic non-MT MSW method to see if there was any way of improving our implementation, particularly regarding the treatment of the interstitial region. To this end I spent a short period in September 1992 with Rino Natoli in Frascati exploring some possibilities. Although nothing of immediate use was developed, some interesting ideas were discussed, some of which bore fruit in later years. In particular, Rino suggested a hybrid method in which the atomic spheres were reduced in size so that they contained only the roughly spherical atomic cores (whose scattering may be then simply treated); and Schrödinger's equation for the remaining smoothly varying interstitial-region potential is solved on a cubic grid using standard methods (Lanczos ?), solutions then being matched across the atomic-sphere boundaries. I believe that this was the approach later used by Yves Joly and colleagues (although I hope they will correct me if I have got the story wrong ...) which I have mentioned in the Alternatives section. We also looked at the possibility of using the cluster version of the cellular method proposed by Brown (1988), an approach which became, after some reworking by Rino, part of a wider collaboration mentioned below. Around this time some work on multicentre spherical-harmonic expansions was published (Foulis (1993)) which may (one day) be the basis of an alternative MSW method.

 

The exchange potential work was not published until long afterwards (Foulis, Pettifer and Jennings (2002)), and then only in outline. At the time, however, we presented some of the work to the X-ray and Inner-Shell Processes Conference, X'93 (1993, Debrecen, Hungary), where we were fortunate to able to discuss its significance with real experts in the field of atomic photoionization processes, such as T Åberg who made some encouraging comments. We also had the pleasure of meeting Bo Wästberg (Chalmers UT, Göteborg) who showed considerable interest in our non-MT work, and described some of his own efforts in the area. He was kind enough to send me afterwards a letter with a preprint of his remarkable work (later to be published as Wästberg (1994)). (I take the opportunity here belatedly to apologise to Bo for not having found the time to reply to his letter: I can only plead extenuating circumstances - with two young children and job hunting providing major distractions at that period.) Once I found the time to sit down and work through the details of his method I was rather impressed. Using an interesting mathematical technique he had developed an accurate approach to the matrix solution of the interstitial region Lippmann-Schwinger equation (see above), and implemented a non-MT MSW code for bound states of diatomic systems (although the method is generally applicable). It appeared to me quite a practical way of approaching the interstitial T-matrix problem, and would have been one of the main improvements to the FPX codes had the opportunity ever arisen.

 

At this point we returned to the important problem of using more realistic electron densities to calculate the molecular potential. Since no non-MT MSW SCF program was available we had to obtain electron densities from another source. Via Martyn Guest of Daresbury Laboratory, we entered into collaboration with Paul Sherwood who worked on the GAMESS-UK code. He helped us run the code and interface it with specially written routines in FPX to generate more realistic potentials. It was important for the comparison with experiment that we choose a system with as few extraneous complications as possible. This pointed to a gaseous molecular system (even though obtaining good XAFS spectra for such systems is far from trivial), and we chose diatomic chlorine (and later hydrogen chloride) since some spectra for these systems had then not long been published (S Bodeur et al. (1990) Z. Phys. D 17, 291) which were particularly useful since absolute cross-sections were available. The results for chlorine (see Foulis, Pettifer and Sherwood (1995a); and, for HCl, see Foulis, Pettifer and Sherwood (1995b)), comparing simultaneously the experiment with MT and non-MT, as well as non-SCF and SCF versions of the theory, showed very clearly (and for the first time in the context of MSW XANES calculations) that both realistic SCF-based potentials and a full non-MT MSW treatment are necessary for good agreement with experiment; and that if both are included in the calculations, one may approach such agreement.

 

We presented our chlorine and hydrogen chloride results to the XAFS VIII conference (1994, Berlin), where we received some favourable and interesting (and, in one case, unprintable but gratifying!) commentary. The results seem to have been quietly and well digested by workers in the field and appear to have informed later work on the problem of ab initio XANES calculations. To this extent (at least) I am convinced that the FPX codes have made a useful contribution to work on this problem. It had become then clear also what was the likely future shape of an upgraded version of the codes, although the task would be non-trivial and there was no immediate prospect of it being realised.

 

Strictly speaking however the development of FPX within the context of XAFS finished at this point. I became involved for a few months in 1994 in a European network ("Multiple Scattering", led by Valérie Briois of LURE), involving many of my former colleagues, one of whose sub-topics was the development of a cellular-method version of the molecular non-MT MSW method (as mentioned above). This was primarily the effort of Tom Harrigan working with Robert Pettifer and Rino Natoli, and some promising results were obtained (and later presented at the X-96 meeting (1996, Hamburg)). In October 1994 I started working in plasma physics, primarily with Steve Rose, then head of theory at the Central Laser Facility at RAL (Rutherford Appleton Laboratory), where we explored the use of the codes to investigate the effect of transient electronic structure on the opacity of dense plasmas. It was interesting to use the codes in such an unexpected way, and some suggestive results were obtained (Foulis, Rose and Beynon (1996, 1997), Rose, Foulis and Gauthier (1998)). In a way this work "pushed the envelope" of the codes and the non-MT MSW method, but highlighted also various shortcomings which would be worth addressing in some future version of the codes.