Many-body perturbation theory calculations using the yambo code

Motivation

• GW-BSE calculations provide quantitatively accurate charged and neutral excitations, but their complexity and relatively poor scaling hinder application to increasingly challenging materials.The project responds to this barrier through physical, algorithmic, computational, and software-development advances.Source: §I · p. 2

Algorithms

• Algorithmic developments improve technical precision and numerical efficiency through convergence acceleration and alternatives to standard GW-BSE approximations.Examples include approaches addressing screening and exciton-coupling approximations, alongside dipole alternatives, empty-state terminators, full-frequency GW, and improved BSE solvers.Source: §I · pp. 2–3

Software ecosystem

• Modern software practices support yambo through modularity, standard libraries, documentation, testing, benchmarking, automation, yambopy analysis scripts, and AiiDA workflow management.The project is distributed through GitHub and includes tools for convergence, benchmarking, reproducibility, and automated workflows.Source: §I · pp. 2, 4

Project overview

• yambo is an open-source many-body perturbation theory code for electronic and optical excitations within a high-performance computational environment.The project has expanded substantially since its first public release and has been applied across several classes of novel materials.Source: §I · pp. 2–3, Figure 1 · p. 3

New physical capabilities

• yambo adds electron-phonon and exciton-phonon interactions, real-time propagation for linear and nonlinear optics, and spin-orbit coupling with Kerr effects in noncollinear BSE.These capabilities target temperature-dependent electronic and optical properties, nonlinear optical properties, and spin-orbit-related phenomena.Source: §I · pp. 2–3

High-performance computing

• The code adopts a mixed MPI-OpenMP parallel strategy and external parallel libraries to exploit massively parallel architectures across calculations with different performance and memory behaviors.The parallel approach is user-tunable and accommodates differing scalability and memory-distribution characteristics of linear response, GW, and BSE calculations.Source: §I · pp. 2–3

Software structure

• DFT outputs are converted into NetCDF database files by a2y and p2y before yambo performs linear-response, GW, and BSE calculations.The standard executable and project-specific executables reflect the general order in which the software is run.Source: §II · p. 4
• Project-specific executables include real-time, nonlinear-optics, phonon, and Kerr variants activated through project-specific compilation flags.The project identifiers include ph, rt, nl, and kerr, with multiple related codes coexisting in the source.Source: §II · p. 4, A. Installation & projects · p. 4

Configuration

• The autotools build system supports portability through configurable executables, optional features, compiler overrides, and direct control over external library paths and include commands.These options support installations on machines with non-standard system directories.Source: A. Installation & projects · p. 4, B. Configuration · pp. 4–5

Executable organization

• The package contains core executables yambo, ypp, a2y, and p2y alongside project suites selected during compilation.The configuration procedure lists available core and project-related suites.Source: A. Installation & projects · p. 4

External libraries

• yambo uses MPI and OpenMP together with BLAS, LAPACK, ScaLAPACK, BLACS, SLEPc, and PETSc for parallel computation and numerical performance.The available libraries include serial and scalable linear-algebra options, while FFT, NetCDF or HDF5, and libxc support computation and I/O.Source: B. Configuration · p. 5

Configuration

• Required libraries can be automatically downloaded, configured, compiled, and reused across subsequent yambo compilations.The installation procedure can place downloaded libraries in a persistent user-selected directory.Source: §1 · p. 5

Interface with Quantum ESPRESSO

• Quantum ESPRESSO interfaces support broader pseudopotential, exchange-correlation, hybrid-functional, spin, and input-file formats.p2y supports UPF version 2, multi-projector pseudopotentials, LDA/GGA and hybrid functionals, collinear and noncollinear spin, qexml/qexsd, and HDF5 files.Source: C. Interfaces with DFT codes · p. 6, §1 · p. 6

Interface with Abinit

• The new Abinit a2y interface reads NetCDF WFK wave-function files directly instead of relying on the discontinued KSS-file workflow.This strategy is described as more numerically efficient, lower-I/O, compatible with recent Abinit developments, and portable to future releases.Source: §2 · pp. 6–7

Data post- (and pre-) processing

• The ypp tool supports preprocessing and postprocessing, including random k-point grids, symmetry manipulation, quasiparticle data handling, band interpolation, and excitonic-wave-function analysis.It can prepare DFT inputs and analyze wave functions, charge densities, DOS, magnetization, currents, band structures, QP-DOS, and excitonic states.Source: D. Data post- (and pre-) processing · p. 9

A. Dipole matrix elements

• Dipole matrix elements can be computed using G-space v, shifted grids, or Covariant approaches, with the latter using finite differences on a regular k-grid.The G-space v method is optimized for nonlocal pseudopotential commutators, while shifted grids and Covariant methods address cases where that evaluation is cumbersome.Source: A. Dipole matrix elements · pp. 9–10
• For local Hamiltonians, the dipole approaches are equivalent and the desired strategy is selected through DipApproach, whose default is G-space v.The input variable offers G-space v, R-space x, shifted grids, and Covariant options.Source: A. Dipole matrix elements · p. 10

Coulomb interaction

• Coulomb-cutoff schemes reduce spurious interactions between periodic replicas in clusters, slabs, wires, and other reduced-dimensional systems.Available geometries include spherical, cylindrical, box-like, and Wigner-Seitz truncations, with convergence depending on supercell size where specified.Source: B. Coulomb interaction · p. 10

C. Sum-over-states terminators in IP Linear Response

• GW sum-over-states expressions converge slowly and can require many empty states, increasing time-to-solution and memory demands.The response function and correlation self-energy are both affected by this computational burden.Source: C. Sum-over-states terminators in IP… · p. 11
• The X-terminator approximates omitted virtual-state effects through corrective terms for the polarizability and self-energy, reducing the bands needed for GW convergence.For bulk Si, solid-circle results with XTermEn=1.5 Ha show reduced band requirements relative to uncorrected convergence.Source: C. Sum-over-states terminators in IP… · p. 11

• The GW self-energy is decomposed into bare exchange and dynamic correlation terms to obtain quasiparticle corrections to Kohn–Sham energies.Source: §IV · p. 11
• yambo supports both plasmon-pole and real-axis strategies for treating the self-energy’s frequency dependence.The Godby–Needs plasmon-pole model matches the inverse dielectric function at static and imaginary frequencies, while real-axis integration retains its full frequency dependence.Source: A. Full frequency GW · p. 12
• Frequency-point calculations can be parallelized independently, providing linear speedup for the computationally intensive inverse-dielectric-matrix evaluation.The calculation may require roughly 100 frequency points for convergence.Source: A. Full frequency GW · p. 12
• Real-axis GW quasiparticle energies achieve the same accuracy as contour-deformation calculations, according to the cited demonstration.Source: A. Full frequency GW · p. 12

Electron-mediated lifetimes

• Real-axis GW calculations provide direct access to quasiparticle electron-mediated lifetimes through the imaginary part of the correlation self-energy.Source: B. Electron-mediated lifetimes · pp. 12–13
• yambo evaluates lifetimes in either the on-the-mass-shell or full GW approximation, with their difference determined by inclusion of renormalization factors.Source: B. Electron-mediated lifetimes · p. 13
• Electron-mediated lifetimes depend only on the k-grid because theta functions restrict scattering events to reduced Brillouin-zone regions.In simple metals, this restriction produces the well-known quadratic scaling of Γe−e near the Fermi level.Source: B. Electron-mediated lifetimes · p. 13

Terminators

• The G-terminator accelerates GW self-energy convergence by significantly reducing the number of states needed to converge quasiparticle corrections.Source: C. Reducing the number of empty states… · pp. 13–14
• Terminator corrections reduce the states necessary to converge the calculation and computational cost, but not the final result.Source: C. Reducing the number of empty states… · p. 14

• GW quasiparticle corrections require wavefunctions and energies on neighboring k+q points, making dense-grid evaluation a convergence challenge.Source: D. Interpolation of the QP band… · p. 14
• Fourier interpolation transfers quasiparticle corrections from a regular grid to high-symmetry lines and can incorporate Kohn-Sham energies near band crossings.Source: D. Interpolation of the QP band… · p. 14
• Wannier interpolation transfers quasiparticle energies and related properties from selected coarse-grid points to much finer meshes at low cost.Source: D. Interpolation of the QP band… · p. 14
• 48x48x1 k-point sampling is required to converge the GW self-energy in monolayer WS2.Source: D. Interpolation of the QP band… · p. 14
• The 6x6x1 Wannier approach uses 7 symmetry-nonequivalent points and requires substantially less CPU time than explicit 48x48x1 quasiparticle calculations.Source: D. Interpolation of the QP band… · p. 15

BSE numerical strategies

• The BSE implementation combines dense electron-hole bases with shifted or symmetry-preserving grids and iterative solvers for large excitonic problems.Source: §V · p. 15, A. Numerical efficiency · p. 15, §1 · p. 15

Spin and magneto-optical effects

• Non-collinear BSE calculations include spin-orbit coupling and expand the excitonic matrix to 4N × 4N because Sz is not a good quantum number.Source: §1 · pp. 16–17
• The BSE matrix yields off-diagonal dielectric-tensor elements for describing the magneto-optical Kerr effect.Source: §1 · p. 17

Fractional occupations

• Occupation-weighted BSE kernels produce Hermitian or pseudo-Hermitian excitonic Hamiltonians for fractional-occupation, nonequilibrium, and metallic systems.Source: §2 · p. 17

Response functions and gauges

• yambo constructs dielectric, dipole, and current response functions and selects length or velocity gauge through the Gauge input.Source: §2 · p. 17

Exciton analysis

• Excitonic and independent-particle density of states help visualize dark or very weak-intensity peaks.Source: §2 · p. 17
• Real-space exciton wavefunctions provide conditional electron probabilities, whose spatial decay distinguishes Frenkel from Wannier excitons.Source: C. Analysis of excitonic wavefunctions · p. 17

Exciton symmetry

• In bilayer hBN, inversion-related phase comparisons classify two interlayer excitons as even and odd under inversion symmetry.Source: C. Analysis of excitonic wavefunctions · pp. 17–18, Figure 10 · p. 18

Scope and implementation

• Electron-phonon interactions affect superconducting critical temperatures, band gaps, carrier mobility, optical spectra, metallic anomalies, and carrier relaxation rates.Source: §VI · p. 18
• yambo ph obtains electron-phonon coupling quantities from Quantum ESPRESSO density-functional perturbation theory and stores them through ypp ph.Source: §VI · p. 18

Electronic structure and spectral functions

• The quasiparticle treatment represents phonon-induced energy renormalization through the real self-energy and finite lifetimes through its imaginary part.Source: A. Temperature-dependent electronic… · p. 18
• Strong electron-phonon interaction can spread spectral functions over a wide energy range, invalidating the single-peak quasiparticle approximation.Source: A. Temperature-dependent electronic… · p. 19

Lifetimes

• In bulk silicon, electron-electron linewidths vanish within energy regions set by the band gap, whereas electron-phonon scattering is stronger there.Source: B. Phonon–mediated electronic lifetimes · p. 19, Figure 12 · p. 19

Temperature dependence

• In single-layer MoS2, electron-phonon corrections generally close the bandgap as temperature increases, with conduction states shifting downward and valence states upward.Source: B. Phonon–mediated electronic lifetimes · p. 20, Figure 13 · p. 20

Finite-temperature BSE

• Finite-temperature BSE uses complex electron and hole energies, producing excitonic binding energies from real parts and lifetimes from imaginary parts.Source: C. Finite Temperature Bethe-Salpeter… · p. 20
• In single-layer MoS2, increasing temperature renormalizes transition energies and broadens optical spectral lines through longer-lived excitonic effects.Source: C. Finite Temperature Bethe-Salpeter… · p. 20, Figure 14 · p. 20

Electron-phonon integration

• Electron-phonon energy corrections converge slowly at some high-symmetry points because symmetry-zero matrix elements may not represent the integration discretization.Source: D. Double grid in the electron-phonon… · p. 20
• yambo addresses electron-phonon integration using random transferred-momentum grids, spherical regions, and Monte Carlo integration.Source: D. Double grid in the electron-phonon… · pp. 20–21

Real-time propagation

• Real-time propagation solves a time-dependent equation of motion under an external laser pulse within the space of equilibrium Kohn-Sham wavefunctions.Source: §VII · p. 21, A. Time–dependent Screened Exchange · p. 21
• The Fourier transform of the time-dependent polarization provides the optical spectrum and absorption through the dipole-dipole response function.Source: A. Time–dependent Screened Exchange · p. 21, Figure 15 · p. 22
• The screened-exchange real-time scheme describes only self-energy variations linear in the density matrix and requires storing the screened-exchange kernel.Source: A. Time–dependent Screened Exchange · pp. 21–22

Double-grid sampling

• The double-grid approach computes dipoles and KHxc on a coarse grid while assigning energies and occupations on a fine k-point grid.Coarse-grid matrix elements are extrapolated to the fine grid using a nearest-neighbour technique.Source: §1 · p. 22

Real-time propagation

• Real-time propagation evolves the electronic system under an external laser pulse and is an unusually distinctive MBPT capability of yambo.The implementation provides two schemes: density-matrix evolution and propagation of periodic Bloch states.Source: B. Nonlinear optics · p. 22

Bloch-state formulation

• Bloch-state propagation couples the external field through a gauge-covariant k-derivative, preserving invariance under unitary rotations among occupied states.In the length gauge, this approach also enables time-dependent polarization through Berry’s phase.Source: B. Nonlinear optics · p. 22

Nonlinear optics

• Nonlinear susceptibilities are extracted by Fourier analysis after laser excitation and dephasing have damped the system’s eigenmodes.The response can be evaluated within IPA, time-dependent Hartree, RT-BSE, or TD-DFT approximations.Source: B. Nonlinear optics · p. 23

Parallel execution

• The parallel implementation combines multilevel MPI with coarse-grained OpenMP to exploit many-core nodes and limited memory per core.Measured performance is demonstrated in a run using up to 1000 Intel KNL nodes.Source: §VIII · p. 24

General structure

• MPI workload is configured per computational kernel through ROLE variables that select distributed parameters and CPU variables that assign task counts.A further g level distributes response matrices over plane-wave degrees of freedom.Source: A. General structure · p. 24
• The multilevel parallel structure distributes work across q-points, k-points, bands, and G-vectors, with additional OpenMP and linear-algebra parallelism.The response-function computation can use five distribution levels, including conduction and valence bands.Source: C. Linear Response · p. 25, D. Self-Energy: HF-Exchange and GW · p. 25

Response functions

• Response-function parallelism trades memory duplication and load balance against communication, making q-point, band, and G-vector choices system-dependent.G-vector distribution can ease memory usage for systems with very large reciprocal-space bases.Source: C. Linear Response · p. 25

Self-energy

• Self-energy parallelism spans q-points, Green-function bands, and quasiparticle corrections, with band distribution generally balancing memory and computation effectively.q-point parallelism may become unbalanced because symmetry expands each q-point to its full star.Source: D. Self-Energy: HF-Exchange and GW · p. 25

Bethe-Salpeter equation

• BSE workload is divided into transition groups and matrix blocks, while symmetry-related transitions and non-uniform group sizes constrain some k-point parallelization choices.When few irreducible k-points are available, splitting k-transition groups across MPI tasks is discouraged because symmetry-related elements may need recomputation.Source: E. Bethe Salpeter Equation · pp. 25–26

F. Linear Algebra

• Dense matrix operations in response and BSE calculations are parallelized with ScaLAPACK within yambo’s MPI structure.The matrices can reach dimensions of a few to tens of thousands or more, making multiplication, inversion, and diagonalization intensive.Source: F. Linear Algebra · p. 26
• For RPA response calculations, parallel linear algebra runs concurrently across q-vectors while response data can also be block-distributed over G and G′.This distributes both computation and memory usage for χ0.Source: F. Linear Algebra · p. 26

Workflow support

• yambopy and AiiDA address the large intermediate data volumes and repetitive multi-code workflows produced by GW and BSE calculations.The platforms organize, simplify, and accelerate data transfer, convergence studies, and reproducible analysis.Source: §IX · pp. 26–27

Yambopy

• Yambopy provides Python classes and scripts for expressing, automating, sharing, and analyzing yambo workflows.It can generate parameterized inputs, read NetCDF databases, and visualize results through matplotlib.Source: A. Yambopy · p. 27
• Yambopy is distributed with documentation, automatic testing, and a public repository, but remains a project under development.Continuous integration tests run at each commit.Source: A. Yambopy · p. 27

B. yambo within the AiiDA platform

• AiiDA stores code-agnostic inputs, converts them into yambo files, submits and monitors jobs, retrieves outputs, and parses relevant data into a database.The workflow can target remote workstations or HPC clusters.Source: §1 · p. 27
• The yambo-AiiDA plugin supports quasiparticle corrections and optical properties, with workflows starting from Quantum ESPRESSO, p2y, or previous yambo runs.Some parsing functionality is provided by yambopy.Source: §1 · p. 27
• The workflow hierarchy ranges from YamboRestart to YamboFullConvergence, with increasing complexity and interdependence.YamboRestart handles GW or BSE calculations, while YamboWf manages the complete DFT-to-yambo path.Source: §2 · pp. 27–28
• YamboConvergence automates quasiparticle-gap convergence through a one-dimensional parameter search and comparisons among recent calculations.By default, convergence requires changes across four successive calculations to remain below the tolerance.Source: §2 · p. 28
• The yambo test-suite automatically compiles or runs code and compares outputs with repository reference files across serial, OpenMP, parallel-I/O, and parallel-linear-algebra modes.Daily short tests and longer pre-release tests provide different testing depths.Source: C. Test-suite and benchmark-suite · p. 28
• The project combines many-body and optical capabilities with native parallel execution, automation interfaces, and tools spanning preparation, calculation, and post-processing.The paper reports that yambo can now study complex materials and interfaces, including unit cells with 10 to 20 atoms and systems reaching hundreds of atoms on HPC facilities.Source: §X · pp. 28–29

Glossary

• The glossary expands abbreviations including BSE, DFT, DFPT, GW, HPC, IPA, and MBPT used throughout the paper.It also defines terms such as conduction-band minimum, Hartree–Fock, Kohn–Sham, and generalized-gradient approximation.Source: Appendix A · p. 30

Response-function evaluation

• The response function in G space is evaluated by separating an internal frequency-independent transition sum from an external frequency-dependent sum over transition groups.This splitting is introduced to compute the response function efficiently.Source: Appendix B · p. 30
• Degenerate poles contribute through the internal sum, while the external sum retains one member from each degenerate group.The grouping is defined over transitions indexed by bands and k-points.Source: Appendix B · p. 30
• The CGrdSpXd input controls the degeneracy threshold, with its default value 100 corresponding to 10^-5 Hartree.Setting CGrdSpXd to zero makes every group size one, so the external sum includes all transitions.Source: Appendix B · p. 30

Response corrections

• yambo implements sum-over-states terminator expressions that provide a correction to the independent-particle response function χ.The appendix identifies these expressions as introduced in an earlier reference and implemented in yambo.Source: Appendix C · p. 30
• In the terminator expression, the extrapolar energy belongs to the polarizability, N′b counts included conduction-band states, fs is the spin occupation factor, and n and m are band indices.These parameters specify the energy extrapolation, conduction-band truncation, spin weighting, and band labels used in the expression.Source: Appendix C · p. 30

Optical matrix elements

• For extended systems, Berry-phase coupling replaces the dipole operator with the k-space derivative r = i ∂/∂k; finite k-point sampling uses a finite-difference representation.In linear response, this representation yields a formula for dipole matrix elements.Source: Appendix D · p. 31