GPAW: An open Python package for electronic-structure calculations

I. INTRODUCTION

• GPAW uniquely supports three wave-function representations—plane waves, real-space grids, and numerical atomic orbitals—and can combine them within one calculation.This multi-basis design makes GPAW versatile because the modes provide complementary accuracy, boundary-condition flexibility, scalability, and system-size coverage.Source: §I · p. 3, A. User’s perspective · p. 4
• GPAW provides extensive functionality beyond ground-state DFT, including GW, BSE, magnetic and spin–orbit calculations, TDDFT, charged defects, localized excitations, and nonlinear optical spectra.The package also supports Berry-phase properties, Wannier-function construction, electron–phonon couplings, and real-time or variational excited-state calculations.Source: §I · p. 3

A. User’s perspective

• Python and direct ASE integration make GPAW calculations easy to script, combine, inspect, and modify, including access to internal quantities and the self-consistent-field loop.Users can insert custom logic during SCF iterations for diagnostics or termination, while ASE/GPAW scripts support flexible multi-stage workflows.Source: A. User’s perspective · p. 4, B. Developer’s perspective · p. 4, Figure 1 · p. 5
• Plane-wave, finite-difference, and LCAO modes trade accuracy, boundary conditions, scalability, memory use, and system size, so the appropriate choice depends on resources and required precision.PW is suited to precise small or medium systems, FD to flexible-boundary and large parallel calculations, and LCAO to very large systems exceeding ten thousand atoms; some features remain mode-specific.Source: A. User’s perspective · p. 4

B. Developer’s perspective

• GPAW’s development ecosystem includes GPL licensing, approximately 1600 unit and integration tests, continuous integration, nightly MPI testing, and automated documentation-example runs.The documentation examples currently require approximately 5200 core-hours on a local supercomputer each weekend.Source: B. Developer’s perspective · pp. 4–5
• GPAW’s newer ground-state design improves extensibility, while production calculations are still recommended with the older code until the new implementation becomes feature-complete.GPU-accelerated PW calculations, wave-function reuse during cell optimization, and spin-spiral calculations have already been implemented on the new code.Source: B. Developer’s perspective · p. 5
• The code combines readable Python with C implementations and optimized numerical libraries for computational hot spots, preserving scripting flexibility while keeping intensive operations efficient.C and library-backed routines handle operations such as finite-difference stencils, grid interpolation, projector overlaps, FFTs, linear algebra, and MPI.Source: B. Developer’s perspective · p. 5

A. Projector augmented-wave method

• The PAW method maps smooth pseudo wave functions to all-electron wave functions through a linear transformation, retaining all-electron information while avoiding rapid nuclear-region oscillations.PAW corrections are assembled from partial waves, projectors, atomic density matrices, and explicitly evaluated exchange–correlation contributions, with key terms stored in setup files.Source: A. Projector augmented-wave method · p. 6

1. Wave-function representations

• GPAW supports plane-wave, linear-combination-of-atomic-orbitals, and finite-difference representations, with the latter representing kinetic energy on a uniform Cartesian grid.Source: §1 · p. 7

2. All-electron quantities

• PAW tools interpolate pseudo wave functions onto fine real-space grids and add corrections to recover all-electron quantities, including cusp and core-orthogonality features.Source: §2 · p. 7
• GPAW computes all-electron electrostatic potentials for TEM simulations, enabling abTEM to use ab initio scattering potentials instead of isolated-atom approximations.This supports investigation of valence-bonding effects in TEM simulations [70].Source: §2 · p. 7

3. Solving the Kohn–Sham equation

• PW and FD modes use iterative diagonalization with density mixing, whereas LCAO uses full Hamiltonian diagonalization; direct minimization is an alternative.Source: §3 · p. 7
• Each SCF step diagonalizes and iteratively solves the Hamiltonian, updates eigenvalues and occupations, then mixes and symmetrizes the density, with LCAO exceptions.For PW and FD initial wave functions come from diagonalizing a small-basis LCAO Hamiltonian built from bound-atomic pseudo partial waves.Source: §3 · p. 7
• GPAW provides direct Fourier-based and iterative multigrid Poisson solvers, with direct solving by default and iterative methods available for larger systems.Multipole corrections improve molecular and charged-system boundary conditions, while dipole layers treat slab work-function differences.Source: §3 · p. 7
• Occupation numbers can be calculated using Fermi–Dirac, Marzari–Vanderbilt [74], Methfessel–Paxton, tetrahedron, or improved-tetrahedron methods.Source: §3 · p. 7

4. Updating wavefunctions in dynamics

• ~15% speedup in realistic structure optimisations results from rigidly updating wavefunctions when atoms move slightly, despite approximate double-counting.The method targets FD and PW modes with iterative eigensolvers; LCAO overlap information could further prevent double-counting.Source: §4 · pp. 7–8

5. Direct minimization

• Direct minimization avoids Hamiltonian diagonalization at each step, is generally less costly and more robust than conventional eigensolvers and density mixing, but requires a finite band gap.GPAW’s formulation supports all orbital representations, finite and extended systems, and unitary- or nonunitary-invariant functionals.Source: §5 · p. 8

6. Convergence criteria

• Convergence is controlled through modular criteria that inspect calculation context and terminate the SCF loop only when all selected criteria report convergence.Built-in criteria include energy, forces, density, and work function, while users can implement custom Python criteria such as bandgap convergence.Source: §6 · pp. 8–9

7. PAW data sets and pseudopotentials

• GPAW supports PAW-XML datasets for multiple functionals, shared collections with Abinit, customizable dataset generation, and norm-conserving pseudopotentials.Dataset modifications can alter accuracy and computational cost, while norm-conserving pseudopotentials approximate PAW with environment-independent nonlocal terms.Source: §7 · p. 9

8. Parallelization

• Parallelization spans k-points, grids, plane waves, bands, and electron–hole pairs, using MPI with optional OpenMP and ScaLAPACK or ELPA support.Band parallelization is especially effective for real-time TDDFT, where bands propagate independently; k-point parallelization typically requires little communication.Source: §8 · p. 9
• GPAW mitigates Python’s parallel filesystem-import overhead by having one process load modules and broadcasting the imported data through MPI.This addresses the heavy load caused when all parallel processes independently read the same Python files.Source: §8 · p. 9
• FD mode scales best at high core counts, whereas PW and LCAO are limited by all-to-all communication and multi-center integrals; GPAW scales to tens or hundreds of nodes.Ground-state calculations distribute work across k-point, band, and domain communicators, with density redistribution available for efficient XC and Poisson operations.Source: §8 · p. 9

9. GPU implementation

• GPAW supports NVIDIA and AMD GPUs through CUDA and HIP, combining custom kernels, vendor libraries, and CuPy for NumPy-like GPU arrays.GPU total-energy calculations are currently supported for LDA/GGA in plane-wave mode.Source: §9 · p. 10
• A header-only porting layer maps generic GPU identifiers to CUDA or HIP equivalents, avoiding duplicated kernels while targeting both hardware platforms natively.Multiple GPUs are parallelized with MPI, assigning one GPU to each rank; GPU-aware MPI enables direct device-to-device communication.Source: §9 · p. 10

C. XC-functionals

• GPAW offers numerous exchange–correlation approximations because the exact exchange–correlation energy is unknown in practice, despite being formally exact in principle.The exchange–correlation term accounts for energy contributions missing from the kinetic, external-potential, and classical Coulomb terms.Source: C. XC-functionals · p. 10

1. Libxc and Libvdwxc

• Libxc supplies LDA, GGA, and MGGA functionals, while libvdwxc evaluates fully non-local van der Waals functionals using fast Fourier-transform convolution [31] [66] [93].GPAW also provides native implementations including TPSS, PBE, and LDA with Perdew–Wang correlation.Source: §1 · p. 10

2. GLLB-sc

• GPAW implements the solid-state GLLBsc potential, which has been shown to improve band-gap descriptions and d-electron states in noble metals [52].Source: §2 · p. 10

3. Hubbard U

• DFT+U applies a Hubbard-like energy penalty to localized orbitals, commonly improving oxide band gaps, formation energies, and magnetic states.GPAW uses a single U_eff parameter, supports occupation-matrix normalization, and permits simultaneous corrections to multiple orbital types on one atom; U_eff is often chosen semi-empirically.Source: §3 · pp. 10–11

4. Hybrids

• Hybrid functionals combine semilocal and Hartree–Fock exchange to address asymptotic-potential, charge-delocalization, and long-range charge-transfer errors.Range-separated functionals use short- and long-range mixing; LCY-PBE applies Hartree–Fock exchange at long range, whereas HSE applies it at short range [5, 116, 120, 121].Source: §4 · p. 11
• Hybrid implementation differs by representation: finite-difference hybrids are generally limited to molecules without forces, whereas plane-wave hybrids support k-points, symmetries, and forces.Source: §4 · p. 11

5. SIC

• PZ-SIC is fully self-consistent and variational, supports atomic forces in all three basis sets, and improves results for systems where common Kohn–Sham functionals fail.Applications include the Mn dimer, substitutional Al in α-quartz, and molecular excitation energies.Source: §5 · pp. 11–12

6. BEEF

• BEEF, mBEEF, and mBEEF-vdW provide GPAW error ensembles for estimating uncertainty in adsorption scaling relations, microkinetics, and materials selection.The approach uses Bayesian model distributions calibrated to benchmark data, but estimates can be unreliable when all ensemble functionals poorly describe the property.Source: §6 · pp. 12–13

IV. ION DYNAMICS

• GPAW operates as a black-box calculator through ASE, supplying energies and forces for geometry optimization, reaction paths, molecular dynamics, and constrained ionic or unit-cell dynamics.Externalizing these algorithms keeps GPAW modular while allowing ASE to provide diverse constraints and optimization methods.Source: §IV · p. 13

A. Structure relaxation

• ASE provides GPAW with quasi-Newton, MDMin, FIRE, SciPy, and preconditioned optimizers for structure relaxation, while GPmin often converges about three times faster than the best quasi-Newton methods.Preconditioning often significantly improves performance, and GPmin fits a Gaussian process to calculated energies and forces.Source: A. Structure relaxation · p. 13

B. Reaction paths and barriers

• ASE-integrated NEB and Dimer methods enable GPAW to locate reaction barriers and characterize complete atomistic pathways, including surface and molecular systems.Surrogate machine-learning models can substantially improve reaction-path optimization performance.Source: B. Reaction paths and barriers · p. 13

C. Global structure optimization

• Machine-learning-accelerated global optimization with GPAW has achieved speed-ups of several orders of magnitude over conventional methods when finding optimal structures.Applications include surfaces, clusters, and crystal structures, with Gaussian-process surrogate potentials, Bayesian optimization, and CMA-ES integration.Source: C. Global structure optimization · pp. 13–14

D. Molecular dynamics and QM/MM

• GPAW supports ASE-based ab initio NVE, NVT, and NPT dynamics and QM/MM simulations with low-overhead user-supplied electrostatic potentials and forces on MM charges.Applications include water, electrochemical interfaces, excited-state bond formation, electron transfer, and coupled solute–solvent dynamics.Source: D. Molecular dynamics and QM/MM · p. 14

V. MAGNETISM AND SPIN

• GPAW supports non-collinear spin, external magnetic fields, spin–orbit coupling, and spin-spiral calculations within the generalized Bloch theorem.These tools enable analysis of magnetic properties and complement methods for magnetic excitations described elsewhere in the paper.Source: §V · p. 14

A. Spin–orbit coupling

• Non-self-consistent spin–orbit coupling rapidly yields eigenenergies by diagonalizing the full Hamiltonian in scalar-relativistic orbitals, requiring only projector overlaps.For WS2, the valence band splits by 0.45 eV at K, while spin degeneracy remains along Γ-M because of two non-commuting mirror symmetries.Source: A. Spin–orbit coupling · p. 14, Figure 3 · p. 15
• The magnetic force theorem estimates orientation-dependent spin–orbit energy changes from rotated Kohn–Sham eigenvalues, enabling high-throughput evaluation of magnetic anisotropy.This is particularly relevant for two-dimensional magnets, where easy-axis anisotropy is decisive for magnetic order.Source: A. Spin–orbit coupling · pp. 14–15

B. Self-consistent non-collinear magnetism

• GPAW implements fully non-collinear spin self-consistency through a spin-density-matrix Kohn–Sham framework, without imposing intra-atomic collinearity in the PAW formalism.The present self-consistent non-collinear implementation is restricted to LDA because GGA generalizations face formal and numerical problems; spin–orbit coupling can be included self-consistently.Source: B. Self-consistent non-collinear… · p. 15

C. Orbital magnetization

• The PAW orbital-magnetization approximation evaluates angular-momentum contributions from all-electron partial waves, avoiding cutoff radii while excluding unbounded partial waves.Finite orbital magnetization requires magnetic order plus spin–orbit coupling or non-coplanar spin textures.Source: C. Orbital magnetization · pp. 15–16
• For bcc-Fe, fcc-Ni, fcc-Co, and hcp-Co, PAW-ACA orbital magnetizations can improve on MT-ACA and agree decently with modern-theory results.The calculations use non-self-consistent spin–orbit coupling, and the comparison includes measurements from experiment.Source: C. Orbital magnetization · p. 16, Table II · p. 16

D. Constant B-field

• A constant external magnetic field couples to electronic spins through a Zeeman term added to the Kohn–Sham Hamiltonian.This provides the field-dependent energy framework used to analyze the Cr2O3 spin-flop transition.Source: D. Constant B-field · p. 16
• At 6 T, Cr2O3 changes from Sz to Sx alignment in LDA+U (U = 2 eV), reproducing the experimental spin-flop field; without U, the critical field doubles.The transition requires fully non-collinear calculations with self-consistent spin–orbit coupling, and the zero-field energy difference gives the magnetic anisotropy.Source: D. Constant B-field · p. 16, Figure 4 · p. 16

E. Spin spirals

• GPAW evaluates planar spin-spiral energies within the chemical unit cell using Generalized Bloch’s Theorem, allowing non-collinear magnetization rotations for fixed q.The Kohn–Sham equations are solved self-consistently for the wave vector q, after which the spin-spiral energy is evaluated.Source: E. Spin spirals · pp. 16–17
• For monolayer NiBr2, the calculated ground state is an incommensurate spin spiral with Q ≃[0.1, 0.1, 0], weak q-dependent local moments, and a spiral-plane normal perpendicular to Q.Including spin–orbit coupling identifies the spiral orientation and shows that the ground state is cycloidal.Source: E. Spin spirals · p. 17, Figure 5 · p. 17

VI. RESPONSE FUNCTIONS AND EXCITATIONS

• GPAW supports response-function calculations alongside GW band structures, BSE optical excitations, magnetic susceptibilities, Raman spectra, and quadratic optical-response tensors.The transverse magnetic susceptibility can describe magnetic excitations, parameterize classical Heisenberg models, and estimate magnetic transition temperatures.Source: §VI · p. 17

A. Linear-response TDDFT

• In linear-response TDDFT, the Kohn–Sham susceptibility is obtained from orbitals, eigenvalues, and occupations, then coupled through a Dyson-like equation containing the Hartree-XC kernel [202].The kernel incorporates the Coulomb interaction and the exchange-correlation kernel of TDDFT.Source: A. Linear-response TDDFT · pp. 17–18
• The ALDA kernel is limited to the long-wavelength metallic regime, cannot describe excitons in extended systems, and produces a divergent XC-hole.A simple renormalization regulates the on-top XC-hole and improves local correlations relative to ALDA and RPA [37].Source: A. Linear-response TDDFT · p. 18

1. Implementation for periodic systems

• For periodic systems, GPAW Fourier-transforms the susceptibility so the Dyson equation becomes a numerically invertible plane-wave matrix equation diagonal in q.The implementation centers on reciprocal-space pair densities and Fourier-transformed XC kernels.Source: §1 · p. 18
• PAW corrections combined with FFT-evaluated pseudo quantities provide reciprocal-space pair densities and XC kernels, while kp-perturbation theory removes the optical-limit Coulomb divergence.The divergent Coulomb term cancels the q-dependence of the pair densities, leaving χ0vc finite.Source: §1 · p. 18

2. Spectral representation

• GPAW evaluates Kohn–Sham susceptibility frequency dependence either explicitly at selected frequencies or through a spectral function followed by a Hilbert transform.The explicit approach suits calculations requiring only a few frequencies.Source: §2 · p. 18
• The spectral-function approach uses a nonlinear frequency grid and is usually faster despite higher memory use, while linear interpolation enables tetrahedron-based k-point convergence improvement.The grid spans eigenenergy differences included in Eq. (52).Source: §2 · p. 18

B. Dielectric function

• GPAW derives macroscopic dielectric functions including local-field corrections and extracts optical absorption and electron energy-loss spectra from the dielectric matrix.The susceptibility can be evaluated within RPA or with ALDA, rALDA [37], and bootstrap kernels.Source: B. Dielectric function · p. 19
• Tetrahedron integration removes unphysical higher-frequency dielectric-function oscillations and is especially advantageous for slowly convergent metallic Ag calculations, compared with point integration.For semiconducting Si, both methods agree as ω →0; increasing point-integration k-point sampling can eventually approach tetrahedron results but at higher computational cost.Source: Figure 6 · p. 19

1. Screening in low-dimensional systems

• For low-dimensional systems, GPAW uses dimensional polarizabilities and reciprocal-space Coulomb-kernel truncation to improve unit-cell-size convergence in dielectric-property calculations.The macroscopic dielectric function depends on the unit cell and tends to unity as the non-periodic cell dimension grows, whereas the d-dimensional polarizability measures induced dipole per length or area.Source: §1 · p. 19

C. Adiabatic-connection fluctuation-dissipation theorem

• ACFDT provides an exact route to electronic correlation energy from the interacting susceptibility, avoiding the exchange–correlation error cancellation used by regular XC functionals.The correlation treatment can be combined with exact exchange.Source: C. Adiabatic-connection… · pp. 19–20
• RPA neglects fxc within ACFDT and captures nonlocal correlation across bonding types including van der Waals interactions, while rALDA renormalization improves short-range correlations and total correlation energies [37].Local adiabatic kernels such as ALDA can produce divergent pair-correlation characteristics, motivating renormalization.Source: C. Adiabatic-connection… · p. 20

D. Magnetic response

• GPAW generalizes LR-TDDFT to four-component charge–spin susceptibilities, whose Dyson equation couples charge and spin except that transverse components decouple for collinear systems without spin–orbit coupling.The four-component density is nµ(r), with µ ∈ {0, x, y, z}.Source: D. Magnetic response · p. 20

1. Transverse magnetic susceptibility

• For collinear nonrelativistic systems in ALDA, the transverse χ+− Dyson equation becomes scalar and uses spin-raising and spin-lowering density operators.GPAW evaluates transverse susceptibility in a plane-wave basis; metals require finite, carefully chosen broadening η depending on k-point sampling.Source: §1 · p. 20

2. The spectrum of transverse magnetic excitations

• The transverse spectral function governs energy dissipation and separates majority and minority excitations according to whether Sz is lowered or raised by one spin-angular-momentum unit.Its plane-wave representation permits direct comparison with inelastic neutron-scattering cross-sections.Source: §2 · pp. 20–21
• GPAW extracts magnon dispersions from spectral peaks and supports studies of nonanalytic itinerant-ferromagnet dispersions, correlation-driven magnetic transitions, and collective modes inside antiferromagnetic Stoner continua.Magnons are collective spin waves, while the Stoner continuum consists of opposite-spin electron–hole pairs.Source: §2 · p. 21
• The spectrum can be decomposed into majority and minority eigenmodes, separating individual magnon lineshapes from the many-body Stoner continuum.The full spectral weight and acoustic, optical, and Stoner contributions can be isolated through eigenvalue analysis, as illustrated for hcp-Co in Fig. 7.Source: §2 · p. 21, Figure 7 · p. 21

3. Liechtenstein MFT

• Liechtenstein MFT maps spin degrees of freedom onto a classical Heisenberg model, calculating exchange parameters from the reactive static Kohn–Sham susceptibility and exchange-correlation magnetic field.GPAW supports variable-size spherical, cylindrical, and parallelepipedic site configurations because magnetic site volumes lack an a priori definition.Source: §3 · pp. 21–22
• MFT and LR-TDDFT give excellent agreement for hcp-Co magnon dispersions, except that MFT underestimates the optical mode along K–M–Γ and misses Stoner-induced fine structure.The ALDA spectrum’s fine structure and Kohn anomalies arise from magnon–Stoner-continuum overlap, which localized-spin MFT cannot generally capture.Source: §3 · p. 22

E. GW approximation

• GPAW supports standard spin-polarized G0W0 quasiparticle calculations using first-order perturbation theory and evaluates the GW self-energy in a plane-wave basis with full real-axis frequency integration.This approach is time-consuming but numerically accurate and provides the full spectral function.Source: E. GW approximation · p. 22
• For low-dimensional systems, GPAW uses a truncated Coulomb interaction to remove periodic-image interactions and applies special q = 0 treatment to improve k-point convergence.For bulk crystals, the q → 0 head and wings of W are obtained using dense k-grid sampling around Γ and numerical or analytical integration.Source: E. GW approximation · p. 22
• GPAW and Yambo produce virtually indistinguishable real and imaginary GW self-energy matrix elements for valence and conduction states at Γ.Both calculations use full frequency integration with 0.1 eV broadening, a 12×12×12 k-point grid, 200 eV cutoff, and 200 bands.Source: E. GW approximation · p. 22, Figure 9 · p. 24

F. Bethe–Salpeter Equation (BSE)

• BSE calculations for WS2 show two spin–orbit-split excitonic peaks below the band gap and triplet-like dark excitons 70 meV below the lowest bright exciton.The valence-band splitting is 0.45 eV, while RPA places the absorption edge at the band gap because it neglects the screened interaction.Source: F. Bethe–Salpeter Equation (BSE) · p. 23, Figure 10 · p. 24
• GPAW’s BSE implementation supports finite-q calculations for plasmon dispersions from EELS and magnon dispersions from transverse magnetic susceptibility.The implementation supports spin-paired and spin-polarized systems, nonselfconsistent spin–orbit coupling, and truncated Coulomb interactions for low-dimensional systems.Source: F. Bethe–Salpeter Equation (BSE) · p. 23

G. Electron–phonon coupling

• Electron–phonon coupling in GPAW provides access to deformation potentials used for electronic and hole transport, while addressing properties including conductivity and superconductivity.The coupling measures interactions between phonon branches and electronic states, and translational-mode matrix elements vanish at q = 0 by the acoustic sum rule.Source: G. Electron–phonon coupling · pp. 23–24
• GPAW computes electron–phonon derivatives by finite differences in supercells, projects them onto LCAO orbitals, Fourier-transforms to arbitrary q, and maps the result to primitive-cell bands and phonon modes.The procedure can use any available GPAW wave-function representation.Source: G. Electron–phonon coupling · p. 24

H. Raman spectrum

• Raman intensities describe inelastic light scattering from vibrational modes, with Stokes spectra determined by phonon frequencies, Bose-Einstein occupations, light polarizations, and Raman tensors.The Raman tensor is associated with each phonon mode at q = 0.Source: H. Raman spectrum · pp. 24–25
• GPAW’s perturbative Raman implementation computes Raman tensors and can extend beyond first-order processes, whereas the more efficient KHD method is limited to first-order Raman processes.The approaches are equivalent when local-field effects are negligible [56]; KHD is implemented in ASR [22], while perturbative calculations are implemented in GPAW.Source: H. Raman spectrum · p. 25
• Bulk 2H-MoS2 exhibits a polarization-resolved Raman spectrum computed at 488 nm using the resonant term, with peaks labeled by irreducible representation.The calculation used a 700 eV plane-wave cutoff, a 2×2×2 k-point mesh, and a 3×3×2 supercell.Source: Figure 11 · p. 25, H. Raman spectrum · p. 25
• GPAW and ASR produce closely aligned Raman peak positions and dominant peaks for materials including MoS2, with minor differences attributed to implementation details and approximations.ASR uses phonopy for phonons and eigenvectors and enforces polarizability-tensor symmetry more rigorously than GPAW’s ASE phonon implementation.Source: H. Raman spectrum · p. 25

I. Quadratic optical response functions

• GPAW’s quadratic optical-response implementation is restricted to second-order dipole response without local-field effects and uses the independent-particle approximation, excluding excitonic collective behavior [57].The implementation currently requires time-reversal symmetry, limiting applications to non-magnetic systems.Source: §I · pp. 25–26
• The implemented quadratic susceptibility yields second-harmonic generation and shift-current tensors, while intrinsic permutation symmetry limits χ(2) to at most 18 independent elements before further point-group reductions.The corresponding quadratic conductivity tensor follows directly from χ(2).Source: §I · p. 26
• With sufficiently many bands, the two quadratic-susceptibility implementations produce identical results; finite Lorentzian broadening and tolerances are used to avoid numerical instabilities.The default tolerance is 10^-6 eV for energy differences and 10^-4 for Fermi-level differences.Source: §I · p. 26

• RT-TDDFT is implemented in finite-difference and LCAO modes, requires non-periodic boundaries, and can model molecules in strong fields or molecule–plasmon interactions.The LCAO implementation is newer and supports versatile analyses, while some functionality is also available in finite-difference mode.Source: §VII · pp. 26–27
• Starting from a ground state, GPAW numerically propagates the time-dependent Kohn–Sham equation, updating the density-dependent Hamiltonian after each step.Users can define external potentials, including delta-kicks and Gaussian pulses; the formulation assumes fixed ions and a time-independent overlap matrix.Source: §VII · p. 27
• Recorded time-dependent variables can be post-processed into physical and chemical observables, such as photoabsorption spectra obtained by Fourier-transforming dipole responses.Implemented observers can write dipole and magnetic moments, frequency-space Kohn–Sham density matrices, and wave functions.Source: §VII · p. 27
• Continuation from a previous propagation and wall-time limits enable short, efficiently scheduled time-propagation chunks on shared high-performance resources.Calculations can begin from the ground state or the last state of an earlier propagation.Source: §VII · p. 27
• The LCAO implementation casts the equation into a matrix problem solved with ScaLAPACK, while Hartree and exchange-correlation potentials are updated on the real-space grid.This implementation contributes to efficient modeling of large systems [52].Source: §VII · p. 27

A. Kohn–Sham decomposition

• The Kohn–Sham density matrix enables efficient observable calculations, frequency-domain analyses, and transition-contribution maps for interpreting localized surface plasmon resonances.Its frequency-space form provides information similar to Casida eigenvectors, while electron–hole decompositions support visualization of excitation contributions.Source: A. Kohn–Sham decomposition · p. 27

B. Hot-carrier analysis

• Weak-perturbation density-matrix post-processing enables spatial and energetic analyses of hot-carrier generation in nanostructures and nanoparticle interfaces.GPAW-related codes and workflows for these analyses are provided in Refs..Source: B. Hot-carrier analysis · p. 27

C. Circular dichroism for molecules

• GPAW’s FD and LCAO implementations calculate molecular electronic circular dichroism from three directional δ-kick propagations, with LCAO enabling nanoscale metal-organic nanocluster calculations [55].The rotatory strength is obtained from frequency components of the time-dependent dipole and magnetic dipole moments.Source: C. Circular dichroism for molecules · p. 28

D. Radiation reaction potential

• Radiation-reaction potentials incorporate analytically known electromagnetic Green’s tensors into TDDFT without additional computational cost for free space, one-dimensional waveguides, and dielectric spheres.Extensions to three-dimensional free space and collective Fabry–Pérot-cavity interactions are under development for polaritonic chemistry.Source: D. Radiation reaction potential · p. 28

E. Ehrenfest dynamics

• Ehrenfest TDDFT/MD accounts for explicit electronic dynamics beyond the Born–Oppenheimer approximation and has been applied to stopping, charge transfer, diffraction, bond breaking, and electrified cathodes.The implementation currently operates in FD mode, while an LCAO version is under development.Source: E. Ehrenfest dynamics · p. 28

A. Improved virtual orbitals

• Improved virtual orbitals accelerate convergence with hybrid and range-separated functionals and permit excited-state properties such as long-range charge-transfer energetics from ground-state calculations.They address spurious unoccupied orbitals and slow convergence affecting linear-response TDDFT for charge-transfer, Rydberg, and doubly excited states.Source: A. Improved virtual orbitals · pp. 28–29

B. Variational excited-state calculations

• Variational excited-state calculations optimize orbitals to non-ground-state stationary points at ground-state-like cost, enabling Hellmann–Feynman forces, geometry optimization, excited-state dynamics, and state-specific QM/MM dynamics.The approach also supports coupling to external MM potentials without additional implementation efforts.Source: B. Variational excited-state… · p. 29
• Maximum-overlap occupation helps prevent variational collapse in non-aufbau excited-state calculations, but collapse and convergence problems remain common with ground-state SCF eigensolvers and density mixing.Variationally optimized states are typically saddle points on the electronic energy surface, making them prone to lower-energy symmetry-preserving solutions.Source: B. Variational excited-state… · p. 29

1. Direct orbital optimization

• GPAW offers DO-MOM and DO-GMF across LCAO, real-space-grid, and plane-wave bases, using quasi-Newton saddle-point optimization and the limited-memory SR1 algorithm.DO-GMF is more robust but computationally costlier because it requires partial Hessian diagonalization.Source: §1 · pp. 29–30
• DO-GMF targets a specified Hessian saddle-point order, eliminating variational collapse and enabling robust excited-state following through bond-breaking configurations.In twisted ethylene, it yields an avoided crossing between the lowest doubly excited state and the ground state where other methods fail.Source: §1 · p. 30

2. Example applications of direct optimization

• Direct-optimization excited-state calculations apply to molecules in gas phase or solution and to solids, including photoinduced structural and solvation dynamics in a copper complex.QM/MM dynamics for [Cu(dmphen)2]+ linked ligand deformation with surrounding-solvent rearrangement after metal-to-ligand charge transfer.Source: §2 · p. 30
• Variational direct optimization substantially reduces charge-transfer excitation errors: 0.15–0.3 eV versus −2.0 eV for matched linear-response TDDFT and −0.85 eV with PBE0.The comparison concerns twisted N-phenylpyrrole and uses the theoretical best estimate of 5.58 eV.Source: §2 · p. 30
• For the negatively charged nitrogen-vacancy center in diamond, r2SCAN reproduces excited-state ordering with an error below 0.06 eV against high-level quantum-embedding calculations.The calculations use plane waves and supercells containing up to 511 atoms.Source: §2 · p. 30

A. Electric polarization

• GPAW evaluates electronic polarization from neighboring-k-point Bloch-orbital overlaps, enabling Berry-phase polarization calculations in the PAW formalism.Formal polarization is multivalued, while polarization changes along adiabatic paths are well defined.Source: A. Electric polarization · pp. 30–31
• For 2D ferroelectric GeS, GPAW obtains a spontaneous polarization of 490 pC/m, in excellent agreement with previous calculations.The framework also supports Born effective charges, piezoelectric response, and lattice polarizability from atomic displacements, forces, and strain.Source: A. Electric polarization · p. 32

B. Berry phases and band topology

• Berry-phase evolution across the Brillouin zone diagnoses band topology, with time-reversal-enforced degeneracies in 1T’-MoS2 indicating a quantum spin Hall state and gapless edge states.GPAW evaluates Berry phases using a smooth gauge constructed by parallel transport, and the Berry-phase spectrum is illustrated for 1T’-MoS2.Source: B. Berry phases and band topology · p. 32, Figure 15 · p. 32
• Berry-phase eigenvalues also yield formal polarization, although this route requires smooth-gauge construction unlike the direct polarization expression.The Berry-phase method has been applied to high-throughput searches for two-dimensional topological materials.Source: B. Berry phases and band topology · p. 32

C. Wannier functions

• Wannier functions are localized electronic representations obtained by unitary transformation of Bloch states that minimizes orbital spread through optimization of a spread functional.GPAW’s ASE Wannier module supports partially occupied Wannier functions for entangled bands and non-periodic systems, with a penalty term improving spread uniformity.Source: C. Wannier functions · pp. 32–33

D. Point defect calculations with hybrid functionals

• Hybrid-functional defect calculations determine formation energies and charge-transition levels using total energies, chemical potentials, Fermi levels, and electrostatic and potential-alignment corrections.The hBN example uses HSE06, an 800 eV plane-wave cutoff, a 128-atom supercell, three charge states, and Freysoldt–Neugebauer–Van de Walle corrections.Source: D. Point defect calculations with… · p. 33
• CB is lower under N-rich conditions, CN is favorable under N-poor conditions, and their charge transitions occur at 3.73 eV and 3.26 eV above or below the VBM, respectively.The results agree well with HSE06 plane-wave calculations using VASP, with minor discrepancies attributed to supercell size and non-local mixing.Source: Figure 16 · p. 33, D. Point defect calculations with… · p. 33

E. Point-group symmetry representations

• GPAW automates point-group symmetry representations for molecular and extended-system Kohn–Sham wavefunctions, enabling analysis of band degeneracies and dipole-transition selection rules.The method projects rotated wavefunctions onto a real-space grid using cubic interpolation and returns irreducible-representation weights for each eigenstate.Source: E. Point-group symmetry representations · p. 34
• Point-group analysis accepts a symmetry center, selected point group, and optional cutoff radius, allowing local symmetry purity to be assessed when distant Hamiltonian symmetry is broken.Implemented groups include C2, C2v, C3v, D2d, D3h, D5, D5h, I, Ih, Oh, Td, and Th.Source: E. Point-group symmetry representations · p. 34

F. Band-structure unfolding

• GPAW unfolds supercell band structures into the primitive-cell Brillouin zone in real-space grid, plane-wave, and LCAO modes without explicitly calculating supercell–primitive-cell wavefunction overlaps.This enables comparison despite supercells having more bands in a smaller Brillouin zone.Source: F. Band-structure unfolding · p. 34

G. The QEH model

• The QEH model couples monopole and dipole dielectric-response building blocks of individual 2D layers to obtain heterostructure screening and excitations.GPAW computes the ab initio RPA building blocks, while the coupled response is obtained by solving a Dyson equation.Source: G. The QEH model · pp. 34–35
• The static dielectric function of stacked N MoS2 and N WSe2 layers increases significantly with layer number and approaches a bulk limit.This layer dependence can inform studies of interlayer and intralayer excitons and band-edge renormalization.Source: G. The QEH model · p. 35, Figure 17 · p. 35

H. Solvent models

• GPAW models solvent effects through a position-dependent permittivity in the Poisson solver, retaining nearly vacuum-level computational cost and enabling solvation free energies and electrolyte-interfaced surfaces.The solvent is represented as a polarizable continuum.Source: H. Solvent models · p. 35

I. Charged electrochemical interfaces

• The Solvated Jellium Method enables constant-charge or constant-potential simulations by compensating excess electrons with a jellium slab embedded in implicit solvent and localizing charge on one surface side.Purely implicit and hybrid explicit-implicit solvent models are supported.Source: §I · pp. 35–36
• Grand-potential energy Ω provides force-consistent energetics for constant-potential relaxations, saddle-point searches, and molecular dynamics.A damped iterative procedure varies the electron number until the target electrode potential is reached.Source: §I · p. 36

J. Constrained DFT

• GPAW implements real-space constrained DFT in FD and LCAO modes for charge- or spin-localized states, with flexible periodicity from isolated molecules to one-, two-, and three-dimensional periodic systems.An auxiliary potential imposes predefined regional charge or spin constraints through Lagrange multipliers and weight functions.Source: J. Constrained DFT · p. 36
• The cDFT implementation provides tools for Marcus-theory transfer parameters across bulk, surface, and molecular systems and supports molecular dynamics and grand-canonical extensions.Analytical forces from the cDFT potential can be used in geometry optimization and molecular dynamics.Source: J. Constrained DFT · p. 36

K. Orbital-free DFT

• A variationally derived Kohn–Sham-like equation applies to the square root of the density, providing the basis for GPAW’s orbital-free formulation.The formulation approximates the Pauli principle and incorporates quantum effects in an averaged way.Source: K. Orbital-free DFT · p. 37
• 10 meV per atom mean absolute error is achieved against reference all-electron values while retaining linear-scaling system-size complexity in GPAW’s PAW real-space OFDFT implementation [338].OFDFT models kinetic energy directly from the density, approximately enforces the Pauli principle, and enables assessment of density functionals using all-electron values, including in large-Z and semiclassical limits.Source: K. Orbital-free DFT · p. 37

L. Zero-field splitting

• GPAW computes zero-field-splitting D tensors from reciprocal-space double integrals using pseudo densities with compensation charges.The splitting arises from magnetic dipole–dipole interactions between the two electrons of a localized triplet state in the absence of an external magnetic field.Source: L. Zero-field splitting · p. 37

M. Hyperfine coupling

• GPAW evaluates hyperfine tensors from pseudo spin densities with compensation charges, including isotropic Fermi-contact and anisotropic dipole–dipole contributions.The isotropic term depends on spin density at the nuclear center, whereas the anisotropic term describes nuclear–electronic magnetic dipole interactions.Source: M. Hyperfine coupling · p. 37

X. OUTLOOK

• GPAW’s ongoing refactoring is intended to improve developer friendliness and facilitate implementation of new functionality, alongside continued efforts to improve performance and applicability.The review emphasizes that these developments depend on contributions from the GPAW developer team.Source: §X · p. 37
• Planned many-body developments include band- and k-point parallelization for plane-wave hybrids, LCAO hybrid functionals with resolution of identity, LCAO GW, self-consistent and vertex-corrected GW, and constrained RPA.The accuracy of LCAO GW with a limited basis remains uncertain, while constrained RPA targets partially screened interactions for low-energy model Hamiltonians.Source: §X · pp. 37–38
• The outlook includes new soft and norm-conserving PAW potentials, LCAO Ehrenfest dynamics, and broader GPU compatibility beyond standard ground-state calculations.Planned GPU targets include linear response and GW, while LCAO Ehrenfest dynamics is being developed to replace comparatively slower grid-mode support.Source: §X · p. 38