GPAW: An open Python package for electronic-structure calculations

There are dozens of electronic-structure codes available for the interested user. The codes differ in their license (in particular, whether open or proprietary), the underlying programming language (e.g. Fortran, C, Python), their treatment of core electrons (all-electron versus pseudopotentials), the employed representations of the wave functions (plane waves, atom-centered orbitals, real-space grids), and the beyond-DFT features they support. Why should one choose GPAW?

In this section, we describe some of the features that make GPAW interesting from the point of view of a common user who wants to perform electronic-structure calculations. The next section focuses on its possibilities for more advanced users, who perhaps want to modify the code or implement completely new functionalities.

A first point to note is that GPAW is written almost exclusively in Python and is directly integrated with the Atomic Simulation Environment. This integration with ASE makes the setup, control, and analysis of calculations easy and flexible. The programming language is of course a key issue for developers, but also the common user benefits from Python and the ASE/GPAW integration. A typical stand-alone program only offers a fixed (though of course possibly large) set of tasks that it can perform, while Python scripting allows for a more flexible use of the code. This could for example mean combining several different GPAW calculations in new ways. Another advantage is that “inner parts” of the code like the density or the Kohn–Sham eigenvalues are directly accessible in a structured format within Python for further analysis. It is even possible to “open up the main loop” of GPAW and have access to, inspect, and also modify key quantities during program execution (see Fig. 1).

As already mentioned in the introduction, GPAW distinguishes itself from other available ES codes by supporting three different ways of representing the wave functions. The most commonly used basis set is plane waves (PW), which is appropriate for small or medium-size systems, where high precision is required. Convergence is easily and systematically controlled by tuning the cut-off energy. A large number of advanced features and “beyond-DFT” methods are available in the PW mode. These include the calculation of hybrid functionals, RPA total energies, linear-response TDDFT, and many-body perturbation theory techniques like GW and the Bethe–Salpeter equations. The new GPU implementation also uses the PW mode.

The wave functions can alternatively be represented on real-space grids, which was the original approach in GPAW. The implementation of this so-called finite-difference (FD) mode relies on multi-grid solutions of the Poisson and Kohn–Sham equations. The FD mode allows for more flexible boundary conditions than the PW mode, which is restricted to periodic supercells. The boundary conditions may for example be taken to reflect the charge distribution in the unit cell. Calculations in the FD mode can be systematically converged through lowering of the grid spacing, but the approach to full convergence is slower than in the PW mode. The FD mode is particularly well suited for large systems because the wave-function representation allows for large-scale parallelization through real-space decomposition. Furthermore, it is possible to perform time-propagation TDDFT including Ehrenfest dynamics in this mode.

The third representation of the wave functions is a basis of numerical atom-centered orbitals in the linear combination of atomic orbitals (LCAO) mode. The size of the basis set can be varied through inclusion of more angular momentum channels, additional orbitals within a channel, or polarization functions. GPAW comes with a standard set of orbitals, but a basis-set generator is included with the code so that users may construct different basis sets depending on their needs and requirements. The LCAO mode is generally less accurate than the PW and FD modes, but it allows for the treatment of considerably larger systems – more than ten thousand atoms. It is also possible to study electron dynamics through a fast implementation of time-propagation DFT, and Ehrenfest dynamics is under development.

As explained, the different modes have different virtues and limitations, and it can therefore be an advantage to apply several modes in a project. For larger systems, it is for example possible to divide a structure optimization into two steps. First, an optimization is performed with the fast LCAO basis leading to an approximately correct structure. This is then followed by an optimization in either the PW or FD mode, which now requires much fewer steps because of the good initial configuration. Due to the ASE/Python interface this combined calculation can easily be performed within a single script.

Since GPAW was originally created with the FD mode only, and the LCAO mode was added next, some features have been implemented for only those modes. Examples are real-time TDDFT (see section VII) and electron–phonon coupling (see section VI.7). Conversely, some new features only work for the PW mode, which was added after the real-space modes. Examples are RPA total energies (see section VI.3) and calculation of the stress tensor. To summarize, given the limitations just mentioned, users should most of the time use PW or LCAO mode and the choice will depend on the accuracy needed and the resources available.

The GPAW source code is written in the Python and C languages and is hosted on GitLab [59] licensed under the GNU General Public License v3.0. This ensures the transparency of all features and allows developers to fully customise their experience and contribute new features to the community.

An advantage of having a Python code is that the Python script you write to carry out your calculations will have access to everything inside a GPAW calculation. An example showing the power and flexibility this affords is the possibility to have user-code inserted inside the self-consistent field (SCF) loop as demonstrated in Fig. 1.

At the time of this writing (July 2023), GPAW has two versions of the ground-state DFT code in the main branch of the code. There is the older version that has grown organically since the birth of GPAW: it has many features, but also a lot of technical debt that makes it harder to maintain and less ideal to build new features on top of. The newer ground-state code addresses these issues by having a better overall design.

The new design greatly improves the ease of implementation of new features. The goal is to have the new code be feature-complete so that it can pass the complete test suite and then delete the old code once that is achieved. At the moment, we recommend that all production calculations are done with the old code and that work on new features is done on top of the new code even though certain features are not yet production-ready. Three new features, not present in the old code base, have already been implemented based on the new code: GPU implementation of PW mode calculations (see section III.2.9), reuse of the wave functions after unit-cell changes during cell optimization, and spin-spiral calculations (see section V.5).

GPAW uses pytest [60] for its test suite that currently consists of approximately 1600 unit and integration tests (see Table 1). A subset of those tests runs as part of GitLab’s continuous integration (CI) thereby checking the correctness of every code change. Unfortunately, the full test suite is too time-consuming to run as part of CI, so we run that nightly both in serial as well as in parallel using MPI.

Many of the code examples in GPAW’s documentation, exercises and tutorials [61] require resources (time and number of CPUs) beyond what would make sense to run as part of the pytest test suite. For that, we use MyQueue [21] to submit those scripts as jobs to a local supercomputer every weekend. At the moment this amounts to approximately 5200 core-hours of calculations.

As can be seen from Table 1, the majority of the code is written in Python, which is an interpreted language that is easy to read, write and debug.

Interpreter-executed code will not run as efficiently as code that is compiled to native machine code. It is therefore important to make sure that the places in the code where most of the time is spent (hot spots) are in native machine code and not in the interpreter. GPAW achieves this by implementing the hot spots in C-code with Python wrappers that can be called from the Python code. Examples of such computationally intensive tasks are applying a finite-difference stencil to a uniform grid, interpolating from one uniform grid to another, or calculation of overlaps between projector functions and wave functions. In addition, we have Python interfaces to the numerical libraries FFTW [62], ScaLAPACK [63], ELPA [17], BLAS, Libxc [64, 65], libvdwxc [66], and MPI. Finally, GPAW makes heavy use of the NumPy [67] and Scipy [68] Python packages. NumPy provides us with the numpy.ndarray data type which is an $N$-dimensional array that we use for storing wave functions, electron densities, potentials, matrices like the overlap matrix or the LCAO wave function coefficients, and much more. The use of NumPy arrays allows us to use the many sub-modules of SciPy to manipulate data. This also gives us an efficient memory layout allowing us to simply pass a pointer to the memory whenever we need to call the C-code from the Python code. With this strategy, we can get away with having most of the code written in a relatively slow interpreted language and still have most of the time spent in highly efficient C-code or optimized numerical libraries.

The advantage of the original FD-mode, where there are no Fourier transforms of the wave functions, is that the algorithms should parallelize well for large systems. In practice, it has turned out that FD-mode has a number of disadvantages: 1) Due to integrals over the unit cell being done as sums over grid-points, there will be a small periodic energy variation as you translate atoms and the period of the variation will be equal to the grid-spacing used (the so-called egg-box error); 2) The system sizes that are typically most interesting for applications of DFT are too small for the parallel scalability to be the decisive advantage; 3) The memory used to store the wave functions on uniform grids in real-space is significant. In contrast, the PW mode has practically no egg-box error, is very efficient for the most typical system sizes and often uses a factor of 10 less memory compared to an FD mode calculation of similar accuracy. The main advantages for LCAO mode are low memory usage and high efficiency for large systems; for small unit cells with many $\mathbf{k}$-points the PW mode is most efficient. One disadvantage of LCAO-mode is egg-box errors: LCAO-mode uses the same uniform grids as used in FD-mode for integration of matrix elements like $\langle\Phi_{\mu}|\tilde{v}|\Phi_{\nu}\rangle$ and therefore has similar egg-box energy variation. A second disadvantage of LCAO-mode is that as for any localized basis-set, reaching the complete basis-set limit is more involved compared to PW and FD modes. This can have severe consequences even for ground state calculations of difficult systems such as Cr₂ for example[wurdemann_density_2015]. In PW or FD modes the complete basis set limit is easy to reach by simply increasing the number of plane-waves or grid-points, respectively, which leads to a smooth convergence[32]. At the moment we only provide double-$\zeta$ polarized (DZP) basis sets and going beyond DZP is left for users to do themselves.

The diverging Coulomb potential causes rapid oscillations in electronic wave functions near the nuclei, and special care is required to be able to work with smooth wave functions. The projector augmented-wave (PAW) method by Blöchl [69] is a widely adopted generalization of pseudopotential methods, utilizing their strength of smooth pseudo wave functions while retaining a mapping from all-electron wave functions ($|\psi_{n}\rangle$) to pseudo wave functions ($|\tilde{\psi}_{n}\rangle$).

The crux is to define a linear transformation $\hat{T}$ from pseudo to all-electron space,

where

Here $\tilde{p}^{a}_{i}(\mathbf{r})$, $\phi^{a}_{i}(\mathbf{r})$ and $\tilde{\phi}^{a}_{i}(\mathbf{r})$ are called projectors, partial waves and pseudo-partial waves, respectively. The pseudo-partial waves and projectors are chosen to be biorthogonal, $\int{\rm d}{\mathbf{r}}\,\tilde{\phi}^{a}_{i}(\mathbf{r})\tilde{p}^{a}_{j}(\mathbf{r})=\delta_{ij}$, allowing for an approximate closure relation

which is utilized heavily to obtain an efficient, but all-electron, picture. In addition to biorthogonality, the pseudo and all-electron partial waves are chosen to be equal outside the PAW augmentation sphere cutoff radius $\tilde{\phi}_{i}(\mathbf{r})=\phi_{i}(\mathbf{r})$, $r>r_{c}.$

The basic recipe for converting an operator $\hat{O}$ is $\langle\psi_{n}|\hat{O}\psi_{n}\rangle=\langle\tilde{\psi}_{n}|\underbrace{\hat{T}^{\dagger}\hat{O}\hat{T}}_{\hat{\tilde{O}}}\tilde{\psi}_{n}\rangle$. For example, the all-electron Kohn–Sham equations

where $\hat{H}$ is the single-particle all-electron Kohn–Sham Hamiltonian operator, can be transformed to their PAW counterparts:

We have used Eq. 1, and also multiplied with $\hat{T}^{\dagger}$ from left to make its dual space the pseudo one, i.e. $\langle\tilde{\psi}|\in\mathcal{H}^{*}$ can act from left. This results in a PAW Hamiltonian and in PAW overlap operators as follows:

Terms such as $\Delta H^{a}_{ii^{\prime}}$ and $\Delta S^{a}_{ii^{\prime}}$ represent so-called PAW-corrections. In each part of the description, which handles a particular kind of operators, such as kinetic energy, spin-operators, or the electrostatic potential, respective PAW-corrections must be calculated. The most crucial ones are precalculated, such as overlap, kinetic energy, Coulomb, and stored in the ‘setup’-file, which also stores the partial waves and projectors. As an example, the overlap PAW corrections are precalculated to setup as follows:

We further define the atomic density matrices as

The atomic density matrix contains all information required to construct PAW corrections to any local all-electron expectation value:

The all-electron atomic density can be constructed as

and the corresponding equation for pseudo-densities holds with $n\rightarrow\tilde{n}$ and $\phi\rightarrow\tilde{\phi}$.

Since the exchange–correlation (xc) potential is non-linear, the PAW corrections must be evaluated explicitly. The xc PAW corrections are performed by constructing the atomic all-electron and pseudo-electron densities as given by Eq. 11:

This integral is numerically evaluated in a Cartesian product of a Lebedev angular grid and a non-uniform radial grid (denser mesh closer to the nucleus) for each atom.

Exchange–correlation (XC) functionals provide a mapping between the interacting and the non-interacting system of electrons. In Kohn–Sham DFT, the density is built from a set of occupied non-interacting single-particle orbitals $\psi_{i}$: $n=\sum_{i}f_{i}\left|\psi_{i}\right|^{2}$, with $f_{i}$ denoting the occupation numbers, leading to the same density as the interacting system. The total energy in DFT is expressed as a sum of density functionals for the different contributions:

where $T_{\mathrm{S}}[n]$ denotes the kinetic energy of the non-interacting system, $V_{\mathrm{ext}}[n]$ the energy of the density in the external potential, $U_{\mathrm{H}}[n]$ the classical Coulomb energy of the density with itself, and $E_{\mathrm{xc}}[n]$ the so called exchange–correlation energy, which collects all energy contributions missing in the prior terms and therefore provides a mapping between the interacting and the non-interacting system of electrons. While the first three terms can be calculated exactly, even the form of $E_{\text{xc}}$ is unknown and although proven to exist and be exact in principle, it has to be approximated in practice. A huge number of approaches belonging to several families exists [1, 2, 90]. Several of these approximations are available in GPAW, and an overview is given in the following.

Local structure optimization in GPAW is typically achieved through the use of an optimizer from ASE. Here a larger range of standard optimizers are available such as quasi-Newton algorithms including BFGS [142, 143, 144, 145] and limited-memory BFGS [146], or Newtonian dynamics-based algorithms such as MDMin and FIRE [147]. That the optimizers have been implemented externally to GPAW provides benefits in terms of simple procedures for restarting long simulations and monitoring their convergence. Some optimizers from SciPy [68] are also available through the open-source ASE package, which provides a simple way to interface any optimizer in the SciPy format to GPAW. Preconditioning is implemented in an accessible way [148], which often leads to significant performance improvements.

Of the classical optimization methods, the quasi–Newton algorithms are often highly competitive. Here one builds up information on the Hessian or inverse Hessian matrix from calculated forces, ultimately leading to an accurate harmonic model of the potential energy surface in the vicinity of the local minimum. Such algorithms can, however, have problems both dealing with anharmonicity in the potential energy surface and with any noise in the electronic structure simulations. It often makes sense instead to fit a Gaussian process to the calculated energies and forces, and minimize within this model potential. This is implemented as the so-called GPmin method in ASE and often converges on the order of three times faster than the best quasi–Newton optimizers. [149]

Reliable calculations of energy barriers are of key importance for determining the rate of atomistic processes. In many quantum-chemical codes utilizing accurate atom-centered basis functions, this is achieved using analytical second derivatives in the direct search for first-order saddle points. This approach is less useful in the plane-wave based or grid-based modes of GPAW. The Dimer method [150] is implemented in ASE and can be used with GPAW, but often one would like to have an overview of the entire reaction path of an atomic-scale process to verify that the correct energy barrier for a process has been determined, and to obtain a full overview of the atomistic mechanism. For this purpose, the Nudged Elastic Band (NEB) method is typically employed through the ASE interface to GPAW. Both the original method [151] and a range of later improvements are available [152, 153, 154, 155, 156, 157]. Special care has to be taken in selecting the optimizer for a NEB algorithm, as this choice can have drastic influence on the convergence rate. For optimization of the reaction path, drastic performance improvements can be obtained by carrying out the optimization in a surrogate machine-learning model fitted to the potential energy surface [158, 159, 160]. GPAW has been used to drive NEB calculations on both surface systems [161, 162, 163, 164] as well as in molecules [164, 165].

GPAW is integrated with various tools for global optimization of structures, compositions, as well as materials morphologies. Some quite generally applicable global optimization tools available are Basin Hopping [166], Minima Hopping [167], and Genetic Algorithms [168]. Some of the most powerful global optimization problems addressed using GPAW rely on machine-learning accelerated global-optimization strategies. These strategies have for example been applied to surfaces, clusters [169], and crystal structures in general [170]. In other machine-learning accelerated global-optimization routines, GPAW was used to generate initial databases of surface and bulk systems, and for later model validation in a strategy that uses Gaussian processes to generate surrogate potential-energy surfaces. These were then explored with Bayesian-optimization techniques, achieving speed-ups over more conventional methods of several orders of magnitude in finding optimal structures of the systems under investigation [171, 172, 173]. The method has been augmented through introducing extra (hyper)dimensions by interpolating between chemical elements, which speeds up the efficiency of the global search [174]. GPAW has been integrated with a Covariance Matrix Adaptation Evolution Strategy (CMA-ES) framework, providing energies and forces, generating training data, and evaluating CMA-ES candidate structures [175].

Ab initio molecular dynamics (MD) can be performed and analyzed through the ASE interface to GPAW. This include standard simulations such as constant-NVE, constant-NVT, and constant-NPT simulations. There is access to e.g. Langevin-, Andersen-, Nose-Hoover-, and Berendsen-dynamics. Due to the externalization of the dynamics, the development of novel algorithms and analysis tools becomes facile. Some examples of the use of GPAW span ab initio MD studies that have explored the liquid structure of water [176] and the water/Au(111) electrochemical interface [177].

Furthermore, GPAW is capable of working with external electrostatic potential terms from user-supplied point charge values and positions, enabling Quantum Mechanical (QM) / Molecular Mechanical (MM) simulations. Since GPAW is from the outset designed around highly efficient grid operations, the computational overhead of evaluating this potential as well as the resulting forces on the MM point charges from the QM density is kept low [178]. GPAW has been central in a range studies of ion dynamics in solution, both in and out of equilibrium. By using the molecular-dynamics functionality of ASE, researchers have performed QM/MM MD simulations of excited-state bond formation in photocatalyst model systems [179, 180, 181, 182], and electron transfer as well as coupled solute-solvent conformational dynamics in photosensitizer systems [183, 165]. Work on polarizable embedding QM/MM within GPAW is ongoing [184, 185].

The spin–orbit coupling is typically completely dominated by regions close to the nuclei where the electrostatic potential becomes strong. Within the PAW augmentation sphere of atom $a$, the spin-$\sigma$ component of orbital $n$ is given by

Assuming a spherically symmetric potential, the spin–orbit Hamiltonian for atom $a$ is then written as

where $V_{r}^{a}$ is the radial electrostatic potential of atom $a$. We evaluate $V_{r}^{a}$ as the spherical part of the XC- and Hartree-potential from the local expansion of the density given by Eq. 11. Since the partial waves $\phi^{a}_{i}$ are eigenstates of the scalar-relativistic Hamiltonian, these are independent of spin and it is straightforward to evaluate the action of $\hat{\mathbf{L}}$ on them.

The eigenenergies can be accurately calculated in a non-selfconsistent treatment of the spin–orbit coupling and may be obtained by diagonalizing the full Hamiltonian in a basis of scalar-relativistic orbitals [186]:

Here $\varepsilon^{0}_{m}$ and $\psi_{m}^{0}$ represent the scalar-relativistic eigenenergies and eigenstates, respectively. This constitutes a fast post-processing step for any scalar-relativistic calculation and only requires the projector overlaps $\langle\tilde{p}^{a}_{j}|\tilde{\psi}_{n\sigma^{\prime}}^{0}\rangle$. It should be noted that this approach in principle requires convergence with respect to the number of scalar-relativistic states included in the basis, but the eigenvalues typically converge rather rapidly with respect to the basis. In Fig. 3, we show the band structure of a WS₂ monolayer obtained from PBE with non-self-consistent spin–orbit coupling. The W atoms introduce strong spin–orbit coupling in this material and the valence band is split by 0.45 eV at the K point. The spin degeneracy is retained along the $\Gamma$-M line, which is left invariant by two non-commuting mirror symmetries.

For magnetic materials the non-selfconsistent treatment of spin–orbit coupling is convenient for evaluating the magnetic anisotropy. The magnetic force theorem [187] implies that rotating the magnetic moments away from the ground state configuration yields a contribution to the energy, which is well approximated by the change in Kohn–Sham eigenvalues. The change in energy for a given orientation of the magnetization density can thus be obtained as

where $\varepsilon_{n}^{\theta,\varphi}$ are the eigenvalues obtained from diagonalizing Eq. (29) with the spins rotated to a direction defined by the angles $(\theta,\varphi)$ and $f_{n}^{\theta,\varphi}$ are the associated occupation numbers. This corresponds to rotating the xc-magnetic field (defined below), which will lead to different eigenvalues when spin–orbit coupling is included. In two-dimensional magnets, an easy-axis anisotropy is decisive for magnetic order [188, 189, 41] and Eq. (30) is easily applied to high throughput computations of magnetic properties [190, 191].

The Kohn–Sham framework for treating non-collinear magnetism was developed by Barth and Hedin [192] and involves the spin density matrix

as the basic variable. The electronic density and magnetization are then given by $n=\mathrm{Tr}[\rho]$ and $\mathbf{m}=\mathrm{Tr}[\boldsymbol{\sigma}\rho]$ respectively. The XC-potential acquires four components, which are the functional derivatives of the XC-energy with respect to the density matrix and may be represented as a $2\times 2$ matrix acting on spinor Kohn–Sham states. These can be expressed in terms of the density and magnetization, which lead to the XC-part of the Kohn–Sham Hamiltonian

where the scalar potential and XC-magnetic field are given by

In GPAW, the self-consistent treatment of non-collinear spin is implemented within LDA where $\mathbf{B}_{\mathrm{xc}}$ is approximated as

Here $m(\mathbf{r})=|\mathbf{m}(\mathbf{r})|$ is the magnitude and $\mathbf{\hat{m}}(\mathbf{r})$ is the direction of the magnetization. The generalization of Eq. 34 to GGAs are plagued by formal as well as numerical problems [193, 194] and the self-consistent solution of the non-collinear Kohn–Sham equations is presently restricted to LDA. It should be emphasized that the PAW formalism allows for a fully non-collinear treatment that does not rely on intra-atomic collinearity, which is often imposed in other electronic structure packages.

Spin–orbit coupling may be included by adding Eq. 28 to the Kohn–Sham Hamiltonian and this constitutes a fully self-consistent framework for spin–orbit calculations within LDA.

Current ab initio methods for determining the orbital magnetization of a material involve either the modern theory, i.e. a Berry-phase formula, or the atom-centered approximation (ACA), where contributions to the expectation value of the angular momentum operator are restricted within atom-centered muffin-tin (MT) spheres with specified cutoff radii. The PAW formulation of the wave functions allows for an approximation similar to the MT-ACA where only the PAW expansion of the wave functions is assumed to contribute significantly to the expectation value of the angular momentum. The orbital magnetic moments can then be calculated through

where $D^{a}_{\sigma ii^{\prime}}$ are elements of the atomic density matrix as defined in Eq. (9), $\phi^{a}_{i}(\mathbf{r})$ are the bound all-electron partial waves of atom $a$, and $\hat{\mathbf{L}}$ is the angular-momentum operator.

We note that this expression does not entail cutoff radii for the atomic contributions to the orbital magnetic moments; instead, the entire all-electron partial waves are included although the all-electron atomic expansion is only formally exact inside the PAW spheres. Additionally, this means that the unbounded, i.e. non-normalisable, all-electron partial waves must be excluded in Eq. (35)

A prerequisite for nonzero orbital magnetization is broken time-reversal symmetry as represented by a complex Hamiltonian that is not unitarily equivalent to a real counterpart. Practically this means that a finite orbital magnetization requires magnetic order and either inclusion of spin–orbit coupling or non-coplanar spin textures. In GPAW, the spin–orbit interaction can be included either self-consistently in a non-collinear calculation or as a non-self-consistent post-processing step following a collinear calculation.