SPHInX

h1. History of SPHInX

_by Christoph Freysoldt_

h2. Initial version 2000-2003

In 2000, Sixten Boeck started his diploma Thesis in the group of Jörg Neugebauer at the "Fritz-Haber-Institut (FHI)":http://www.fhi-berlin.mpg.de on an implementation of density-functional perturbation theory in the FHI plane-wave code fhi98md. Coming from a computer science background, he quickly realized that FORTRAN90, the dominant computer language used for high-performance computing (HPC) in the physics community at that time, would not allow him to use post-1970 programming paradigms such as modularization, encapsulation, dynamic memory allocation, etc., efficiently.

For his "PhD,":https://d-nb.info/99806761X/34 Sixten suggested to write a modern DFT code in a modern language (C++). However, it was generally believed that C++ programs are much slower than FORTRAN programs due to the abstraction overhead. Sixten quickly demonstrated that this was a misconception: using conventional C++ paradigms such as a high degree of abstraction for all objects, including vectors for number crunching, indeed was much slower than the conventional FORTRAN paradigm of properly ordered loops. However, C++ as a programming language was very well capable of reaching and sometimes even exceeding FORTRAN performance, if the programming style of core routines was adapted to the specific challenges of high-performance computing. In contrast to FORTRAN, however, it was very easy to encapsulate these core functionalities into classes with very intuitive interfaces that allowed a very concise coding style. The same strategy was also applied to integrating numerical libraries that performed specific tasks (Fast Fourier transforms, FFTs, and linear algebra) in an hardware-specific, severely optimized manner.

The surprising outcome was soon that his programs would often run faster than usual FORTRAN code, because the simple interfaces would seduce the developer to use optimized routines systematically, while FORTRAN seduced the developer to prefer straight-forward explicit implementations over complicated function calls that could not exploit the actual hardware. The programming efforts were then rapidly joined by others (notably Alexey Dick and Lars Ismer, later Abdallah Qteish, Matthias Wahn), and resulted in a very fast plane-wave DFT code in 2003, named sfhingx 1.0.

h2. SPHInX 1: Growing capabilities. 2004-2010

One of the key innovations in sfhingx 1.0 was the idea to combine a mathematical vector with the metadata required for the physical context (e.g. the plane-wave basis set for a vector of plane-wave coefficients) into a single object. This allowed to write non-algebraic transformations such as FFTs in an elegant, yet transparent way within the code. When the analogy to the use of Dirac's bra-ket notation became clear, Sixten decided to redesign the entire program around this concept, removing along the way also other design mistakes of the original version. This version became operational in 2005, and was named S/PHI/nX. /PHI/ stands for the Greek letter φ , the first letter of physics. Sphinx was also the Egyptian goddess of wisdom in ancient times.

However, the appointment of J. Neugebauer at "Max-Planck Institute für Eisenforschung GmbH":https://www.mpie.de in 2005 and the subsequent relocation of the group to Düsseldorf slowed down the progress. Sixten Boeck became main responsible for building up the new Computer Center at MPIE, and in addition to his ongoing PhD thesis, found little time to work on the functionality. Moreover, the change in the group's focus from semiconducting to transition metals required to go beyond norm-conserving pseudopotentials. The projector-augmented wave (PAW) formalism was the obvious choice, but required a significant change to the implementation.

It was only in 2009, that the PAW implementation was finally addressed by Christoph Freysoldt, who had joined the SPHInX developer team a few years earlier. Thanks to the modular design of the SPHInX library, he put together a first working version of PAW within about 30 working days, while initial estimates had been 6-12 months. The first implementation closely followed the original formalism by Peter Bloechl, and also used his PAW potential formats. However, to ensure that results could compared more directly to other codes, readers for other PAW potential formats (atompaw, abinit, vasp) were added later on. During the time that the PAW implementation was gradually stabilized and improved, a new generation of developers joined in, notably Björn Lange and Gernot Pfanner. 

A little earlier, other people had started on expanding the code's capabilities in other directions: Hazem Abu-Farsakh introduced a self-consistent tight-binding model, Oliver Marquardt created a plane-wave based multi-band *k·p* code, Blazej Grabowski and later Klaus Reuter added MPI parallelization, which soon developed into an easy-to-use interface for MPI, that once again follows the SPHInX coding paradigm: to provide access to the most relevant advanced high-performance routines with simple usage rules.

h2. SPHInX 2: Maturing. 2010-2020

In 2010, the main developer S. Boeck decided to leave academia for founding his own IT company, and he passed the scepter on to Christoph Freysoldt, who is in charge of supervising SPHInX development until today. Yet, the work that had went into creating a practical interface to a variety of standard computer programming tasks was attractive for bringing the new company to speed. Therefore, they decided to split the code base into a basic computer-science part,  "SxAccelerate,":https://gitlab.com/sphinxlib/sxaccelerate and a more physics-related part, and continue collaboration on the former part.

In the following years, a large number of smaller projects addressed adding new functionality to the code (such as spin constraints, geometry optimization, Hubbard-U corrections for molecular orbitals) and improving the performance on modern CPU-based multi-core computer architecture. Today, our core routines run 2-3x faster than the corresponding "naive" algorithms based on external libraries by combining close-to-theoretical-peak-performance kernel routines with dispatcher algorithms that exploit the specific characteristics of our data and computational tasks.

h2. SPHInX 3: 2021-

While the Dirac vector concept proved to be a central innovation in SPHInX's code design, it's actual implementation had gradually grown out of Sixten's original vector library from 2000, and it was severly limited by the compiler capabilities of those times. As a result, the vector implementation showed a few pecularities that created a barrier for new developers. For instance, the use of reference counting as a key technique to avoid data copies could lead to unexpected duplicate references to the same data in different variables. Also, the underlying scalar type was always wrapped in a datafree type mapper class. Last, purely mathematical vectors and basis&metadata-aware Dirac vectors were treated as incompatible classes, while they shared 99% of their implementation.  First ideas how to remove such flaws appeared in 2012, but it was only in late 2019 that Christoph decided to take this crucial step towards the future. 

The new vector class separates the control over memory containing the data from its possible access patterns, by implementing references to such data as the central Dirac vector concept: the vector reference. Each reference has its independent set of metadata and access patterns (such as column and row strides for non-contiguous subsets of the data). This is supplemented by the "standard" vector, a contiguous specialization of the vector reference with "copy-on-assignment" semantics that resolves the frequent confusion between (efficient) referencing and actual copying when passing data around. Thanks to the capabilities of modern C++ (since the 2011 standard), many unecessary copies can be avoided by the compiler. Likewise the distinction between metadata-free math vectors and meta-data containing Dirac vectors could be dropped, because we had realized in the mean time that the overhead for unconditionally copying a few dozen bytes of (possibly unused) metadata was negligible in all use cases, and definitely cheaper than making any distinction. Christoph finalized the implementation of the vector class in june 2020, but it took him until the end of the year to adapt the physics part to the subtle differences between the two implementations, notably when the distinction between copy and reference did matter. SPHInX 3.0 was declared functional in april 2021, after rehauling the k dot p part and incorporating the latest developments by Olli Marquardt.