GPAW: An open Python package for electronic structure calculations Open Access

The electronic-structure (ES) problem, i.e., the solution to the time-independent Schrödinger equation for a collection of electrons and atomic nuclei, forms the starting point for the quantum-mechanical treatment of matter. Indeed, all chemical and physical properties of any substance (solid, molecule, surface, etc.) can, in principle, be obtained from the energies and wave functions that constitute the solution. A pioneering step towards solving the many-body ES problem was the formulation and formal proof of density functional theory (DFT) by Hohenberg and Kohn in 1964¹ and a practical scheme for its solution by Kohn and Sham in 1965.² At present, most codes solving the ES problem from first principles are based on DFT. Such codes are extremely powerful and allow one to determine the atomic structure of solids and molecules containing hundreds of atoms with a relative error below 1%.^3–5 Once the atomic structure of the compound has been solved, its properties (electronic, magnetic, optical, topological, etc.) can, in principle, be determined. The evaluation of properties often involves theories beyond the formal DFT framework to account for effects such as temperature and lattice vibrations,^6,7 many-body interactions in excited states,^8,9 or time dependence.^10,11 As such, first-principles atomistic calculations often involve two successive phases: the solution of the ground-state ES problem (including ion dynamics) and the subsequent evaluation of physical properties. This review is structured accordingly, as Secs. III and IV deal with the first phase while Secs. V–IX are devoted to the second.

In recent years, the scientific significance of ES codes has shifted from a useful tool to describe and understand matter at the atomic scale to an independent driver of the discovery and development of new materials.^12–16 This change in scope has been fueled by the exponential increase in computer power accompanied by improved numerical algorithms^17,18 as well as the use of workflow management software for high-throughput computations^19–22 and the adoption of machine-learning techniques to leverage the rapidly growing data generated by ES codes.^23–25 In parallel with these capacity-extending developments, continuous progress in the fundamental description of exchange–correlation effects has advanced the predictive power of ES calculations to a level where they rival experiments in terms of accuracy for many important properties.^26–31 

The grid-based projector-augmented wave (GPAW) code was originally intended as a Python-based multigrid solver of the basic DFT equations within the projector-augmented wave (PAW) formalism.³² The name GPAW, accordingly, was an abbreviation for “grid-based projector-augmented wave.” At present, other choices than regular grids for representations of the wave functions exist in GPAW, but the name has stuck. During the years 2005–2010, GPAW evolved to a full-blown DFT package³³ supporting most of the functionality expected from a modern ES code, in addition to a few more specialized features, including real-time propagation of wave functions³⁴ and an alternative basis of numerical atomic orbitals [referred to as the linear combination of atomic orbitals (LCAOs) basis]³⁵ to supplement the real-space grid. In 2011, a plane-wave (PW) basis set was also implemented. At present, the possibility to use three different types of basis sets and even combine them within a single run remains a unique feature of GPAW, rendering the code very versatile.

The implementation of the PW basis set laid the groundwork for GPAW’s linear-response module, which at present supports the calculation of linear response functions,³⁶ total energies from the adiabatic connection fluctuation–dissipation theorem,^29,37 the GW self-energy method for quasiparticle band structures,³⁸ the Bethe–Salpeter Equation (BSE) for optical excitations,³⁹ and more. The code also supports a wide range of features related to the calculation of magnetism and spin–orbit effects. Examples include spin-spiral calculations using the generalized Bloch theorem,⁴⁰ external magnetic fields, orbital magnetization, magnetic anisotropy,⁴¹ adiabatic magnon dispersions from the magnetic force theorem,⁴² and dynamic magnetic response from time-dependent density functional theory (TDDFT).⁴³ For solids, the k-space Berry phases can be computed directly from the Bloch orbitals and may be used to obtain spontaneous polarization,⁴⁴ Born effective charges, piezoelectric response tensors,⁴⁵ and various indices characterizing the band topology.⁴⁶ 

In addition, GPAW can compute the localization matrices forming the basis for the construction of Wannier functions with, e.g., the Atomic Simulation Environment (ASE)⁴⁷ or Wannier90.⁴⁸ Electrostatic corrections to the formation energies of charged point defects in insulators are implemented, as are calculations of the hyperfine coupling and zero-field splitting for localized electron spins. GPAW also offers the possibility to perform time-independent, variational calculations of localized electronic excitations in, e.g., molecules or crystal point defects, using direct orbital optimization strategies implemented for all three types of basis sets.^49–51 This provides an efficient and robust alternative to traditional “ΔSCF” approaches. GPAW can also be used to describe ultrafast electron dynamics within time-dependent density functional theory (TDDFT), with wave functions represented either on a real-space grid³⁴ or in the LCAO basis.⁵² The latter can provide a significant speed-up due to the relatively small size of the basis.^53–55 The LCAO representation also forms the basis for the calculation of electron–phonon couplings as well as non-linear optical spectra such as Raman scattering⁵⁶ (which can alternatively be obtained in the PW mode as a finite difference of the dielectric tensor), second-harmonics generation,⁵⁷ and shift currents⁵⁸ using higher-order perturbation theory.

There are dozens of electronic-structure codes available for interested users. The codes differ in their license (in particular, whether open or proprietary), the underlying programming language (e.g., Fortran, C, and Python), their treatment of core electrons (all-electron vs pseudopotentials), the employed representations of the wave functions (plane waves, atom-centered orbitals, and 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. Section II B focuses on its possibilities for more advanced users, who may 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 the common user also 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-sized systems where high precision is required. Convergence is easily and systematically controlled by tuning the cutoff 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 graphics processing unit (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 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 (LCAOs) mode. The size of the basis set can be varied through the 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 the 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 to a project. For larger systems, it is, for example, possible to divide 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 Sec. VII) and electron–phonon coupling (see Sec. VI G). Conversely, some new features only work for the PW mode, which was added after the real-space modes. Examples are RPA total energies (see Sec. VI C) and the 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 webpage has a table showing up-to-date information about which features work with which modes.⁵⁹ 

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

An advantage of having 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 of having 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 implementing 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 be performed with the old code and that work on new features be performed 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 Sec. III B 9), reuse of the wave functions after unit-cell changes during cell optimization, and spin-spiral calculations (see Sec. V E).

GPAW uses pytest⁶¹ for its test suite, which currently consists of ∼1600 unit and integration tests (see Table I). 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⁶² require resources [time and number of central processing units (CPUs)] beyond what would make sense to run as part of the pytest test suite. For that, we use MyQueue²¹ to submit those scripts as jobs to a local supercomputer every weekend. At the moment, this amounts to ∼5200 core-hours of calculations.

As can be seen from Table I, 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 into 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 calculating overlaps between projector functions and wave functions. In addition, we have Python interfaces to the numerical libraries FFTW,⁶³ ScaLAPACK,⁶⁴ ELPA,¹⁷ BLAS, Libxc,^65,66 libvdwxc,⁶⁷ and MPI. Finally, GPAW makes heavy use of the NumPy⁶⁸ and Scipy⁶⁹ 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 performed 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 of LCAO mode are low memory usage and high efficiency for large systems; for small unit cells with many k-points, the PW mode is the 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 $⟨Φμ|ṽ|Φν⟩$ 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.⁷⁰ 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 smooth convergence.³² At the moment, we only provide double-ζ 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⁷¹ 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 (|ψ_n⟩) to pseudo-wave functions $(|ψ̃n〉)$⁠.

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

The beauty of the PAW method is that you never need to transform the pseudo-wave functions to all-electron wave functions, but you can do it if you want to. GPAW has tools for interpolating the pseudo-wave functions to a fine real-space grid and adding the PAW corrections. A fine real-space grid is needed to properly represent the cusp and all the oscillations necessary for the wave function to be orthogonal to all the frozen core states.

GPAW also has tools for calculating the all-electron electrostatic potential. This is useful for transmission electron microscopy (TEM) simulations.⁷² Most TEM simulations have relied on the so-called independent atom model (IAM), where the specimen potential is described as a superposition of isolated atomic potentials. While this is often sufficient, there is increasing interest in understanding the influence of valence bonding.⁷³ This can be investigated by a TEM simulation code such as abTEM,⁷⁴ which can directly use ab initio scattering potentials from GPAW.

The default method for solving the Kohn–Sham equation for PW and FD modes is to do iterative diagonalization combined with density mixing; for LCAO mode, we do a full diagonalization of the Hamiltonian. Alternatively, one can do direct minimization as described in Sec. III B 5.

For PW and FD modes, we need an initial guess of the wave functions. For this, we calculate the effective potential from a superposition of atomic densities and diagonalize an LCAO Hamiltonian in a small basis set consisting of all the pseudo-partial waves corresponding to bound atomic states.

Each step in the self-consistent field (SCF) loop consists of the following operations: (1) diagonalization of the Hamiltonian in the subspace of the current wave functions (skipped for LCAO); (2) one or more steps through the iterative eigensolver (except for LCAO, where a full diagonalization is performed); (3) update of eigenvalues and occupation numbers; and (4) density mixing and symmetrization. See previous work^32,33 and Ref. 75 for details.

GPAW has two kinds of Poisson equation solvers: direct solvers based on Fourier transforms or Fourier-sine transforms, and iterative multi-grid solvers (Jacobi or Gauss–Seidel). The default is to use a direct solver, whereas iterative solvers may be chosen for larger systems where they can be more efficient.

For 0-, 1-, and 2-dimensional systems, the default boundary conditions are to have the potential go to zero on the cell boundaries. This becomes a problem for systems involving large dipole moments. The potential due to the dipole is long-ranged and, therefore, the converged potential requires large vacuum sizes. For molecules (0D systems), the boundary conditions can be improved by adding multipole moment corrections to the density so that the corresponding multipoles of the density vanish. The potential of these corrections is added to the obtained potential. The same trick is used to handle charged systems. For slabs (2D systems), a dipole layer can be added to account for differences in the work functions on the two sides of the slab.

Methods for calculating occupation numbers are the Fermi–Dirac, Marzari–Vanderbilt,⁷⁶ and Methfessel–Paxton distributions, as well as the tetrahedron method and the improved tetrahedron method.⁷⁷ 

Simulations commonly move the atoms without changing other parameters. If an atom moves only slightly, we would expect most of the charge in its immediate vicinity to move along with it. We use this to compute an improved guess for the wave functions in the next self-consistency loop with FD or PW mode, where the eigensolver is iterative.

The method could be further improved by using the LCAO basis set and the overlap matrix to prevent double-counting.

The GPAW formulation of DM is applicable to all representations of the orbitals available in GPAW, as well as Kohn–Sham (unitary invariant) and orbital density-dependent (nonunitary invariant) energy functionals, and can be used for both finite and extended systems. In LCAO calculations, the reference orbitals are expressed as a linear combination of the atomic basis functions Φ, Ψ₀ = ΦC₀, where the matrix of coefficients C₀ is fixed. Therefore, only a minimization with respect to the elements of matrix A is required.⁸² For plane wave and real-space grid representations, a minimization in the space tangent to the reference orbitals is sufficient if the functional is unitary invariant. Otherwise, if the functional is nonunitary invariant, such as when self-interaction correction is used (see Sec. III C 5), an inner loop minimization in the occupied–occupied block of the matrix A is performed to make the energy stationary with respect to the unitary transformation of the occupied orbitals.⁵⁰ 

The DM method avoids diagonalization of the Hamiltonian at each step, and as a result, it usually involves a smaller computational effort. The DM method has also been shown to be more robust than conventional eigensolvers and density mixing in calculations of molecules and extended systems.⁸² However, the current implementation does not support a finite-temperature distribution of occupation numbers and thus can only be used for systems with a finite bandgap.

The modular architecture of GPAW allows the user to have precise control over how the SCF loop decides that the electronic structure has converged to sufficient precision. GPAW contains simple keywords for common convergence criteria, such as “energy,” “forces,” (electron) “density,” and “work function,” which are sufficient for the most common use cases.

Internally, all convergence criteria are Python objects that are instances of convergence classes. For each step through the SCF loop, each convergence object is called. When any convergence object is called, it is passed a context object that contains the current state of the calculation, such as the wavefunctions and the Hamiltonian. The criterion can thus pull relevant data from the calculation to decide if it is convergent. Because the convergence criterion itself is an object, it can store information such as previous values of the energy for comparison to the new value. When all convergence criteria report that they are converged, the calculation as a whole is considered to be converged and terminates.

This modular nature gives the user full control over how each convergence criterion operates. For example, the user can easily ask the energy criterion to check the differences in the last four values of the energy rather than the last three. If a convergence criterion itself is expensive to calculate, it can make sense to not check it until the rest of the convergence criteria are met. This can be accomplished by activating an internal “calculate_last” flag within the convergence criterion.

Users can easily add their own custom convergence criteria to the SCF loop. If a user would like to use a criterion not included by default with GPAW, it is straightforward to write a new criterion as a Python class and pass this to the convergence dictionary of GPAW. For example, if one wanted to be sure the bandgap of a semiconductor was converged, the criterion could check the bandgap at each iteration, compare it to stored values from previous iterations, and report that the calculation is converged when the peak-to-peak variance among the last n iterations is below a given threshold.

GPAW can read PAW datasets from (possibly compressed) XML files following the PAW-XML specification.⁸³ Dataset files for most of the periodic table can be downloaded from the GPAW webpage or installed with the gpaw install-data command-line tool. The datasets are available for the local-density approximation (LDA), Perdew–Burke–Ernzerhof (PBE),³ revPBE,³⁴⁶ RPBE,³⁴⁵ and GLLBSC⁹⁷ functionals. The electronic structure code Abinit⁸⁴ also reads the PAW-XML format, allowing GPAW and Abinit to share PAW dataset collections such as the Jollet–Torrent–Holzwarth collection.⁸⁵ 

Specialized datasets can be generated with the gpaw dataset command-line tool. This allows one to tweak the properties of a dataset. Some examples could be: (1) add more/less semi-core states; (2) increase/decrease the augmentation sphere radius to make the pseudo-wave functions more/less smooth; (3) add/remove projector functions and corresponding pseudo and all-electron partial waves; or (4) base the PAW dataset on a different XC-functional. These changes will affect the accuracy and cost of the calculations.

GPAW is also able to use norm-conserving pseudopotentials (NCPP) such as Hartwigsen-Goedecker-Hutter (HGH)⁸⁶ and pseudopotentials in the Unified Pseudopotential Format (UPF), such as SG15.⁸⁷ Non-local NCPPs can be considered an approximation to PAW: in the PAW description, the non-local part of the Hamiltonian [the term containing $ΔHija$ in Eq. (6)] will adapt to the environment, whereas for a NCPP, $ΔHija$ will be diagonal and have a fixed value taken from a reference atom. Because of norm-conservation, NCPPs will have $ΔSija=0$⁠.

GPAW can parallelize over various degrees of freedom, depending on the type of calculation, and it implements multiple algorithms for achieving good parallel performance and scalability. In calculations involving k-points, parallelization over them is typically the most efficient, as little communication is needed during the summation of wave functions to calculate the density and any derived quantities. As the number of k-points is often limited, especially in large systems, parallelization is also possible over real-space grids in the FD and LCAO modes, as well as over plane waves in the PW mode. All modes also support parallelization over electronic bands, which is particularly efficient for real-time TDDFT where the time-propagation of each band can be carried out independently. Additional parallelization possibilities exist depending on calculations, such as over electron–hole pairs in linear-respose TDDFT calculations.

Parallelization is performed mainly with MPI. In the FD and LCAO modes, it is possible to use OpenMP within shared-memory nodes, which can improve performance when the number of CPU cores per node is large. Dense linear algebra, such as matrix diagonalization and Cholesky decomposition, can be carried out with the parallel ScaLAPACK or ELPA libraries. This applies to both the direct diagonalization in the LCAO mode as well as to subspace diagonalizations in the iterative Davidson and Residual Minimization Method with Direct Inversion in the Iterative Subspace (RMM-DIIS) methods in the FD and PW modes.

For ground-state calculations, GPAW will divide the cores into three MPI communicators: k-points, bands, and domain. When parallelizing over k-points and/or bands, all the cores of the k-point and/or band communicators will have a copy of the density (possibly distributed over the domain communicator). GPAW has the option to redistribute the density from the domain-communicator to all cores so that operations such as evaluating the XC-energy and solving the Poisson equation can be performed more efficiently.

For each k-point, all the cores in the band and domain communicators will cooperate on calculating the matrix elements of the Hamiltonian and the overlap operators. Dense-matrix linear-algebra operations on those matrices can be performed on a single core (most efficient for small systems) or with ScaLAPACK, where the matrices are distributed over either all or some of the cores from the pool of cores in the band and domain communicators.

One drawback of Python is that in large parallel calculations, its import mechanism may incur a heavy load on the filesystem because all the parallel processes are trying to read the same Python files. GPAW tries to alleviate this by using a special “broadcast imports” mechanism: during initial module imports, only a single process loads the modules; afterward, MPI is used to broadcast the data to all the processes.

Parallel scalability depends strongly on the calculation mode and the system. The FD mode offers the best scalability for high core counts, as only nearest-neighbor communication is needed across domains over a domain decomposition. In PW mode, the limiting factor is all-to-all communication in parallelization over plane waves. In LCAO mode, communications arise from multi-center integrals of basis functions across domains. At best, GPAW scales to tens or hundreds of nodes in supercomputers.

The GPU implementation of GPAW works both on NVIDIA and Advanced Micro Devices (AMD) GPUs, targeting either CUDA or Heterogeneous-Compute Interface for Portability (HIP) backends, respectively. GPAW uses a combination of manually written GPU kernels, external GPU libraries (such as cuBLAS/hipBLAS), and CuPy.⁸⁸ CuPy offers an easy-to-use Python interface for GPUs centered around a NumPy-like GPU array and makes many hardware details completely transparent to the end-user.

In the manually written GPU kernels, both GPU backends (CUDA and HIP) are targeted using a header-only porting approach⁸⁹ in which generic GPU identifiers are translated to vendor-specific identifiers at compile time. For example, to allocate GPU memory, the identifier gpuMalloc is used in the code, which is then translated either to cudaMalloc or hipMalloc depending on which GPU backend is targeted. This allows us to avoid unnecessary code duplication and still target multiple hardware platforms natively.

An earlier GPU implementation of GPAW^90,91 served as the starting point for the recent work on a new GPU code based on the rewritten ground-state code. The objects that store quantities like $ψ̃n(r)$⁠, $p̃ia(r)$⁠, $⟨p̃ia|ψ̃n⟩$⁠, $Dii′a$⁠, $ñ(r)$ and $ṽ(r)$ use Numpy arrays for the CPU code and CuPy arrays when running on a GPU. At the moment, GPUs can be used for total-energy calculations with LDA/generalized gradient approximation (GGA) in the PW mode.

Parallelization to multiple GPUs is performed using MPI. Each MPI rank is assigned a single GPU, and communication between the GPUs is handled by MPI. Support for GPU-aware MPI makes it possible to do direct GPU-to-GPU communication without unnecessary memory copies between the GPU device and the host CPU.

The libxc library⁶⁶ provides implementations of several (semi-)local variants of the XC functional, given by the LDA, GGA, and meta-GGA (MGGA)⁹¹ families. These are available in GPAW by a combination of their names from libxc, e.g., “GGA_X_PBE+GGA_C_PBE” for PBE.³ At the moment, using MGGA functionals from libxc in GPAW requires libxc to be compiled with --disable-fhc. Most MGGA functionals in libxc do not depend on the laplacian of the density—those that do are not supported.

Additionally, GPAW provides its own implementation of several (semi-)local functionals called by their short names, e.g., TPSS,⁹³ PBE,³ and LDA, the latter with the correlation of Perdew and Wang.⁹⁴ Several hybrids (see below) are implemented in GPAW with the support of the libxc library for their local parts.

For fully non-local van der Waals functionals, like the van der Waals (vdW)-DF functional,³¹ GPAW uses the efficient fast Fourier-transform convolution algorithm by Román-Pérez and Soler⁹⁵ as implemented in the libvdwxc library.⁶⁷ 

All LDA, GGA, and common MGGA functionals, such as Strongly Constrained and Appropriately Normed semilocal density functional (SCAN)²⁸ and TPSS, from libxc can be used in GPAW.

GPAW has an implementation of the solid-state modification of the Gritsenko–van Leeuwen–van Lenthe–Baerends exchange–correlation model potential (GLLB-sc).^96,97 This potential has been shown to improve the description of the bandgap¹⁶ and the d-electron states in noble metals.^52,98–100

DFT+U calculations using a Hubbard-like correction can be performed in GPAW to improve the treatment of Coulombic interactions of localized electrons. This correction is most commonly applied to the valence orbitals of transition metals to assist with obtaining experimental bandgaps of oxides,¹⁰¹ which may otherwise be underestimated.^102,103 In addition, formation energies and magnetic states are often improved due to a more correct description of the electronic structure. It may also be applied to main-group elements such as N, O, P, or S, but this is less commonly performed.¹⁰⁴ 

GPAW supports the normalization of the occupation matrix, accounting for the truncated portion of the wave function outside the augmentation sphere. To maintain consistency with other codes that do not support normalization, this normalization can be disabled, but large disagreements are expected when applied to p orbitals. GPAW is one of the few codes that currently supports multiple simultaneous corrections to orbitals on the same atom; this is useful when two types of orbitals, such as p and d orbitals, are nearly degenerate but both are partially occupied.

There is no U^eff that is strictly correct, but methods such as RPA¹⁰⁷ or linear response^108,109 can allow for a value to be calculated from first principles. More commonly, U^eff is chosen semi-empirically to fit experimental properties such as formation energy,¹¹⁰ bandgap¹¹¹ or, more recently, machine learning predictions.¹¹² 

Hybrid functionals, especially range-separated functionals (see below), correct problems present in calculations utilizing (semi-)local functionals such as the wrong asymptotic behavior of the effective potential leading to an improper description of Rydberg excitations,¹¹³ the improper calculation of the total energy against fractional charges leading to the charge delocalization error,¹¹⁴ and the wrong description of (long-range) charge transfer excitations due to the locality of the exchange hole.^115–117 

Details on the FD-mode implementation of long-range corrected RSF can be found in Refs. 124 and 125. In general, the FD-mode implementation of hybrids is limited to molecules, and forces have not been implemented. The PW-mode implementation of hybrids handles k-points, exploits symmetries, and can calculate forces.

PZ-SIC has been shown to give accurate results in cases where commonly used KS functionals fail. This includes, for example, the Mn dimer, where the PBE functional gives qualitatively incorrect results while the corrected functional gives close agreement with high-level quantum chemistry results as well as experimental measurements.¹³³ Another example is the defect state of a substitutional Al atom in α-quartz.¹³⁴ In addition, PZ-SIC has been shown to improve the values of the excitation energy of molecules obtained in variational calculations of excited states^50,135 (see Sec. VIII B).

A great strength of Kohn–Sham DFT and its extensions is that a reasonably high accuracy of physical, material, and chemical properties can be obtained at a relatively moderate computational cost. DFT is thus often used to simulate materials, reactions, and properties where, de facto, there is no “better” alternative. Even though more accurate electronic structure methods in principle might exist for a given application, the poor scaling of systematically more accurate methods often makes these computationally infeasible at the given system size that is being studied using DFT. Therefore, one is often in the situation where the accuracy of a given DFT calculation of some materials or chemical properties cannot be verified against, e.g., a more accurate solution of the Schrödinger equation, even on the biggest available supercomputers.

On the other hand, the wealth of available XC functionals naturally allows one to look at how sensitive a DFT result is to the choice of functional, and often accuracy is, therefore, being judged primarily by applying a small set of different XC functions, especially if no accurate benchmark theoretical simulation or experimental measurement is available. A challenge, however, is that the different available functionals are often known to be particularly good at simulating certain properties and poor at others. It is thus not at all clear how much one should trust a given functional for a given simulated property. The Bayesian Error Estimation (BEE) class of functionals¹³⁶ attempts to develop a practical framework for establishing an error estimate from a selected set of functional “ingredients.”

Assume that the XC model M(a) is a function of a set of parameters, a, that can be varied freely. If a benchmark dataset, D, of highly accurate properties established from experiments or higher-accuracy electronic structure simulations is available, we can attempt to identify the ensemble of models given by some distribution function, P(M(a)), such that the most likely model in the ensemble, M(a₀), makes accurate predictions for the benchmark dataset, while the spread of the ensemble reproduces the spread between the predictions of the most likely model and the benchmark data.

Figure 2 shows an ensemble of exchange enhancement factors F_x(s) from the BEEF-vdW functional, where s = |∇n|/(2k_Fn) is the reduced density gradient. The approach has given rise to several functionals, BEEF-vdW,¹³⁷ meta-BEEF (mBEEF),¹³⁸ and mBEEF-vdW,¹³⁹ which all include error estimation and are readily available in GPAW. Through the ASE interface, one can, for example, utilize the error ensembles to establish error estimates on Python-implemented models that are using DFT simulations for their parametrization. An example of this is the application of error estimates to adsorption scaling relations, microkinetics, and material selection.¹⁴⁰ 

One risk to the approach of establishing error estimates from a small selected “ensemble” of XC-functionals is clearly that if the simulated property in question is poorly described by all functionals in the ensemble, then the error estimate might also become poor. This could, for example, be the case if one tried to establish an error estimate for bandgaps in oxides or van der Waals bonding of adsorbates on surfaces based on an ensemble of GGA XC-functionals, since no GGA functional may be accurate for simulating either property.

GPAW can be employed as a “black-box” calculator, supplying energies and forces to other programs such as, e.g., ASE, which then perform optimization of ground state geometries and reaction paths or carry out molecular dynamics. In fact, this is a key design principle behind GPAW. Methodological developments and general implementations that are not directly dependent on fast access to detailed electronic structure information should preferably be implemented externally by GPAW. This led to the maximal simplicity of the GPAW code itself while also allowing for the external code to be utilized and tested with other electronic structure codes and atomistic potentials. A key to the high efficiency of GPAW simulations involving ionic displacements is the versatile implementation of restraints in ASE. Here, many types of restraints are readily accessible, from the simple removal of degrees of freedom to more exotic restraints allowing for rigid molecule dynamics,¹⁴¹ harmonic restoring forces,¹⁴² spacegroup preservation, and combined ionic-unit cell dynamics.¹⁴³ Many algorithms are available for various structure optimization and molecular dynamics tasks.

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^144–147 and limited-memory BFGS,¹⁴⁸ or Newtonian dynamics-based algorithms such as MDMin and FIRE.¹⁴⁹ The fact that the optimizers have been implemented externally in GPAW provides benefits in terms of simple procedures for restarting long simulations and monitoring their convergence. Some optimizers from SciPy⁶⁹ 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,¹⁵⁰ 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 the potential within this model. 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.¹⁵¹ 

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¹⁵² 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¹⁵³ and a range of later improvements are available.^154–159 Special care has to be taken in selecting the optimizer for a NEB algorithm, as this choice can have a 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.^160–162 GPAW has been used to drive NEB calculations on both surface systems^163–166 as well as in molecules.^166,167

GPAW is integrated with various tools for global optimization of structures, compositions, as well as material morphologies. Some quite generally applicable global optimization tools available are basin hopping,¹⁶⁸ minima hopping,¹⁶⁹ and genetic algorithms.¹⁷⁰ 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,¹⁷¹ and crystal structures generally.¹⁷² 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 the optimal structures of the systems under investigation.^173–175 The method has been augmented by introducing extra (hyper)dimensions by interpolating between chemical elements, which speeds up the efficiency of the global search.¹⁷⁶ 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.¹⁷⁷ 

Ab initio molecular dynamics (MD) can be performed and analyzed through the ASE interface to GPAW. This includes 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¹⁷⁸ and the water/Au(111) electrochemical interface.¹⁷⁹ 

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.¹⁸⁰ GPAW has been central in a range of 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^181–184 and electron transfer, as well as coupled solute–solvent conformational dynamics in photosensitizer systems.^167,185 Work on polarizable embedding QM/MM within GPAW is ongoing.^186,187

Many important technological applications utilize magnetic order or the manipulation of spin in materials. GPAW has a wide range of functionalities that facilitate the analysis of magnetic properties. This includes calculations with non-collinear spin, the inclusion of external magnetic fields, spin–orbit coupling, and spin spiral calculations within the generalized Bloch theorem. The implementation of these features is described below, while additional methods for calculating magnetic excitations are described in Sec. VI D.

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.

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 wave is included, although the all-electron atomic expansion is only formally exact inside the PAW spheres. Additionally, this means that the unbounded, i.e., non-normalizable, 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 its real counterpart. Practically, this means that a finite orbital magnetization requires magnetic order and either the 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. The orbital magnetization has been calculated using Eq. (35) for the simple ferromagnets bcc-Fe, fcc-Ni, fcc-Co, and hcp-Co, where the spin–orbit interaction is included non-self-consistently. The PAW-ACA results are displayed in Table II along with MT-ACA results, modern theory results, and measurements from experiments, demonstrating mainly that the PAW-ACA can be an improvement over the MT-ACA and secondly that there is a decent agreement between the PAW-ACA and the modern theory for these systems.

Using the GPAW functionality to compute the spin-spiral energy as a function of q, one can then search for the ground-state wave vector Q that minimizes the energy. In Fig. 5, we show that the monolayer NiBr₂ has an incommensurate spin-spiral ground state and that the local magnetic moment on the Ni atom depends only weakly on the wave vector q.

Furthermore, the orientation of the ground state spin spiral can be obtained by including spin–orbit coupling non-self-consistently in the projected spin–orbit approximation²⁰² and search for the orientation that minimizes the energy. The ordering vector Q and the normal vector n of the spiral plane thus constitute a complete specification of the magnetic ground state within the class of single-q states. The normal vector is of particular interest since spin spirals may lead to spontaneous breaking of the crystal symmetry, and the normal vector largely determines the direction of magnetically induced spontaneous polarization.⁴⁰ In Fig. 5, we show that n is perpendicular to the wavevector Q in the monolayer NiBr₂ ground state, corresponding to a cycloidal spin spiral.

Linear response functions are the bread and butter of condensed matter physics. Their applications include the description of dielectric screening, optical and electron energy loss spectra, many-body excitations, and ground state correlation energies. In this section, we describe the methods available in GPAW for calculating electronic response functions as well as GW quasiparticle band structures and optical excitations from the Bethe–Salpeter Equation (BSE). In addition to the electronic response function, GPAW can also calculate the transverse magnetic susceptibility, which holds information about the magnetic excitations, e.g., magnons, and can be used to derive parameters for classical Heisenberg spin models and to estimate magnetic transition temperatures. Finally, we present methods to calculate Raman spectra of solids and quadratic optical response tensors for describing second harmonics generation and the Pockels electro-optical effect.