Source-linked AI summary
BerkeleyGW: A Massively Parallel Computer Package for the Calculation of the Quasiparticle and Optical Properties of Materials and Nanostructures
Jack Deslippe, Georgy Samsonidze, David A. Strubbe, Manish Jain, Marvin L. Cohen, Steven G. Louie
TL;DR
Excited-state properties require methods that capture interactions missed by standard ground-state theories, but GW-BSE calculations are computationally expensive and many implementations remain small-scale. BerkeleyGW combines GW quasiparticle calculations with Bethe-Salpeter optical calculations in a massively parallel package, achieving good agreement with silicon experiment when both effects are included and scaling to very large CPU counts.
Problem
Standard ground-state theories fail to correctly capture excited-particle interactions, while GW-BSE calculations are costly and many implementations are limited to tens of atoms and about 100 CPUs.
Method
BerkeleyGW performs ground-state calculations, solves Dyson’s equation within GW for quasiparticles, and solves the Bethe-Salpeter equation for correlated electron-hole excitations.
Results
Good agreement with silicon experiment is reached only when both GW quasiparticle and Bethe-Salpeter excitonic effects are included, while epsilon and sigma computations scale nearly linearly to approximately 100,000 and 25,000 CPUs, respectively.
Takeaways & Limitations
BerkeleyGW supports quasiparticle and optical-property calculations for diverse materials, including nanostructures, using massively parallel computing beyond the scale of many earlier implementations.
Takeaways & Limitations
Practical material size is constrained by computational resources, and BSE diagonalization is more limiting than kernel construction because it typically saturates at about 1000 CPUs.
Abstract
from arXiv · showhide
BerkeleyGW is a massively parallel computational package for electron excited-state properties that is based on the many-body perturbation theory employing the ab initio GW and GW plus Bethe-Salpeter equation methodology. It can be used in conjunction with many density-functional theory codes for ground-state properties, including PARATEC, PARSEC, Quantum ESPRESSO, OCTOPUS and SIESTA. The package can be used to compute the electronic and optical properties of a wide variety of material systems from bulk semiconductors and metals to nanostructured materials and molecules. The package scales to 10,000's of CPUs and can be used to study systems containing up to 100's of atoms.
1. Program Summary
BerkeleyGW addresses excited-state electronic and optical properties using GW and Bethe-Salpeter methods. It supports parallel Linux/UNIX computing, with practical system size limited by available computational resources.
- The distributed package contains 80,000 lines and requires 200MB including test data, with BLAS, LAPACK, FFTW, and optionally ScaLAPACK and MPI.
- Ground-state theories such as standard density-functional methods fail to correctly capture interactions affecting excited-particle energies, wavefunctions, and lifetimes.
- BerkeleyGW constructs and solves Dyson’s equation within GW for quasiparticle properties, then solves the Bethe-Salpeter equation for correlated electron-hole excitations.
- Materials with up to 500 atoms per periodic cell can be studied on large HPC systems, subject to available computational resources.
- Running time is 1-1000 minutes depending strongly on system size and processor number, while memory use is 50-2000MB per CPU.
2. Introduction
The paper presents BerkeleyGW as a scalable implementation of accurate GW-BSE calculations for systems beyond the size and CPU limits of many existing implementations. It integrates with multiple ground-state codes to study quasiparticle and optical properties across diverse materials.
- GW and BSE methods have demonstrated accurate quasiparticle and optical-property predictions across bulk materials, nanostructures, and molecules.
- GW-BSE calculations are typically an order of magnitude or more costly than DFT, while many implementations handle only tens of atoms and about 100 CPUs.
- BerkeleyGW implements ab initio GW and BSE methods, alleviates small-system restrictions, and scales beyond thousands of CPUs.
- The package works with PARATEC, Quantum ESPRESSO, SIESTA, PARSEC, Octopus, and an included empirical pseudopotential code for mean-field calculations.
3. Theoretical Framework
The GW-BSE workflow starts from ground-state electronic structure, computes quasiparticle corrections, and then obtains correlated electron-hole excitations. BerkeleyGW uses this framework to calculate optical spectra, with silicon agreeing well with experiment only when both GW and excitonic effects are included.
- Theoretical Framework: GW-BSE calculations proceed through ground-state properties, GW quasiparticle energies and wavefunctions, and Bethe-Salpeter electron-hole excited states.
- Theoretical Framework: GW quasiparticle energies and wavefunctions are obtained by solving Dyson’s equation with the GW self-energy operator.
- Theoretical Framework: Bethe-Salpeter equations calculate electron-hole excitation states, using an electron-hole interaction kernel and the Tamm-Dancoff approximation.
- Theoretical Framework: The imaginary dielectric function provides an optical response based on exciton transition matrix elements and excitation energies.
- Theoretical Framework: Only the combined inclusion of GW quasiparticle and Bethe-Salpeter excitonic effects gives good agreement between BerkeleyGW’s silicon absorption spectrum and experiment.
4. Computational Layout
BerkeleyGW organizes ground-state inputs, dielectric screening, self-energy corrections, and Bethe-Salpeter calculations into interoperable executables. Its layout supports convergence analysis, frequency-dependent dielectric calculations, and extensive parallelization, while certain limits remain computationally demanding.
- Workflow: Mean-field orbitals, eigenvalues, and charge density are converted from supported DFT codes into BerkeleyGW’s plane-wave format before excited-state calculations.The package includes wrappers for several mean-field codes.
- Dielectric screening: The epsilon executable computes static or frequency-dependent RPA polarizability and inverse dielectric matrices from mean-field eigenvalues and eigenvectors.Its main outputs include eps0mat and epsmat for q → 0 and q ≠ 0.
- Quasiparticle corrections: The sigma executable constructs W and G, evaluates the self-energy Σ = iGW, and writes quasiparticle energies EQP for subsequent analysis.The screened interaction comes from the inverse dielectric matrix, while G comes from mean-field eigenvalues and eigenfunctions.
- Optical properties: The BSE kernel uses the dielectric matrix to screen the attractive electron-hole interaction, after which absorption calculations use the resulting kernels and quasiparticle data.Example calculations and post-processing utilities are included for several material systems.
- Scaling: The matrix elements needed for the polarizability scale as N^3 log N, while ideal unlimited-CPU wall time approaches N log N through parallelization.The FFTs scale as N log N and the matrix-element computation can be parallelized to tens of thousands of CPUs.
- Frequency dependence: Full-frequency dielectric calculations cost about twice a static calculation for 10-50 frequencies, then scale linearly with the number of frequencies.The frequency summation becomes the bottleneck because matrix elements are computed only once.
5. Parallelization and Performance
BerkeleyGW parallelizes its major GW-BSE stages across bands, matrix elements, reciprocal-space indices, and independent sampling points. This design reduces per-CPU memory and delivers near-linear scaling to tens of thousands of CPUs for representative nanotube calculations.
- epsilon: Matrix-element computation in epsilon scales nearly linearly up to Nv · Nc processors and can reach 100,000 CPUs for large systems.Each processor handles an approximately equal fraction of valence–conduction pairs and performs FFTs for its assigned pairs.
- epsilon: The epsilon code switches from band-based parallelization to block-cyclic distribution over G,G′ for polarizability and dielectric-matrix operations.This layout supports distributed matrix multiplication and ScaLAPACK-based inversion while controlling memory per processor.
- epsilon: Memory for χ(G,G′) decreases linearly with CPU count, although dielectric-matrix inversion saturates at hundreds of CPUs and summation at a few thousand.The CPU count for the block-cyclic distribution can therefore be tuned to balance memory reduction and computation time.
- epsilon: Nearly linear epsilon scaling reaches 3,200 CPUs per q-point and approximately 100,000 CPUs when 32 q-points are parallelized independently.The reported scaling concerns both memory and computation time for the (20,20) single-walled carbon nanotube example.
- sigma: Sigma parallelizes outer-state pairs across pools and inner states within pools, yielding N log N wall-time scaling for matrix elements and N^2 scaling for summations with unlimited resources.The matrix-element and summation stages use separate levels of band-based parallelism.
- sigma and BSE: The example sigma calculation scales nearly linearly to 1,600 CPUs per k-point and 25,000 CPUs across 16 irreducible k-points, while the kernel scales nearly linearly to 1,024 CPUs.The independent k-point and kernel parallelism provides additional scaling beyond the per-sample distributions.
6. Coulomb Interaction
BerkeleyGW supports bulk and geometry-specific truncated Coulomb interactions, with cell averaging used to regularize q → 0 behavior and improve convergence. Grid uniformity must also be checked because strongly anisotropic sampling can produce incorrect integral limits.
- Cell averaging: Cell averaging replaces q-point interaction values with reciprocal-space volume averages, optionally including the inverse dielectric function when screened W is relevant.The implementation uses Monte Carlo integration with 2,500,000 random points per cell.
- Coulomb-interaction choices: BerkeleyGW offers bulk, slab, wire, box, and spherical Coulomb interactions, with truncation geometry matched to the system.Slab truncation targets slab-like systems, wire truncation targets wire-like systems, and box or spherical schemes truncate in all three directions.
- Small-q treatment: For divergent small-q interactions, averaging W over the q = 0 cell yields finite values, including for metallic systems where W metal(q →0) = C.For metals, averaging W directly is important because ǫ−1(q) cancels the small-q dependence of v(q).
- Small-q treatment: In the (14,0) nanotube, ǫ−1(0) = 1 but decreases nearly by half at the first non-zero q-grid point, motivating dielectric-function inclusion in cell averaging.BerkeleyGW uses this W-averaging procedure by default for truncated Coulomb interactions.
- Convergence: Using an extension of the averaging procedure beyond q, GG′ = 0 can speed sigma or absorption convergence by replacing finite q sums with integrals.Users can apply cell-averaged W below an input energy cutoff.
- Grid uniformity: If the longest effective grid vector exceeds the shortest by more than a factor of 2, BerkeleyGW warns that non-uniform sampling may produce incorrect answers.Users should use a more uniform grid or check convergence against the cell-averaging cutoff; truncated directions are excluded from this uniformity check.
7. Symmetry and degeneracy
BerkeleyGW uses crystal symmetries to expand reduced-zone wavefunctions and dielectric matrices, reduce selected self-energy sums, and preserve degeneracy-aware results. Incomplete treatment of degenerate subspaces can break symmetry and yield qualitatively incorrect optical results.
- Symmetry reduction: Reduced-zone wavefunctions are expanded to the full Brillouin zone using reciprocal-space rotations and fractional translations.The full-zone wavefunctions support sums in epsilon, sigma, kernel, and absorption.
- Symmetry reduction: The shifted dielectric-matrix grid is generated by unfolding reduced points, reducing under the q-vector subgroup, and applying a small q-shift.The resulting shifted grid can contain more points than the original uniform grid.
- Dielectric matrices: Symmetry maps dielectric matrices from a reduced q-zone to the full zone using transformed q- and reciprocal-lattice vectors.This reduces the number of dielectric matrices that epsilon must calculate explicitly.
- Degeneracy: Truncating unoccupied states through degenerate subspaces can break full system symmetry, although the effect is small for large numbers of bands.The degeneracy override flag is available for testing or cases with overlapping degenerate subspaces.
- Degeneracy: For diagonal self-energy elements over a complete degenerate subspace, irreducible-zone q sums with symmetry weights can replace full-zone sums and are symmetrized by averaging.Off-diagonal elements generally still require the full q sum.
- Degeneracy: BSE active spaces containing only part of a degenerate occupied or unoccupied subspace can break symmetry and produce qualitatively incorrect absorption results.A symmetry-forbidden dark excitation may acquire nonzero oscillator strength when only one transition is included.
- Degeneracy: The degeneracy utility identifies acceptable band counts consistent with degeneracy across multiple wavefunction files.It can check shifted, unshifted, coarse, and fine grids used in epsilon and BSE calculations.
8. Computational Issues
BerkeleyGW includes build, memory-diagnostic, testing, and continuous-integration infrastructure for parallel scientific calculations. Its test suite checks runtime behavior, reproducibility, serial/parallel consistency, and compatibility across supported platforms and libraries.
- Runtime diagnostics: Each run estimates available and required memory per CPU, warning users when allocation demands exceed available memory and failure is likely.The estimate is based on the largest arrays after system parameters are read.
- Build and installation: The build uses GNU Make with cross-directory dependencies, selectable real or complex flavors, system configuration files, and a required test-suite check.Installation selects flavor.mk, links arch.mk, runs make all, and then runs make check or a parallel equivalent.
- Verification: The test suite runs mean-field, epsilon, sigma, BSE, PlotXct, and utility calculations, then compares selected outputs with reference data within specified tolerances.It detects runtime errors and generated warnings.
- Verification: The test suite verifies builds, reproducibility, serial/parallel consistency, real/complex consistency, spin consistency, and compatibility with new compilers or libraries.It can run serially or in parallel through scheduler-specific jobscripts.
- Verification: BuildBot-based continuous integration triggers tests after repository commits to help maintain code integrity during development.A PBS submission script can submit jobs, monitor status, capture output, and determine success or failure.
- Supported environments: The package reports extensive testing across Linux, AIX, and MacOS with multiple Fortran, C/C++, MPI, LAPACK/BLAS, and ScaLAPACK/BLACS implementations.MPI is optional, while LAPACK/BLAS is required and ScaLAPACK/BLACS is required by BSE when MPI is used.
9. Utilities
BerkeleyGW provides visualization, structure-manipulation, format-conversion, image-charge, and quasiparticle band-structure utilities. These tools connect electronic-structure outputs to visualization, molecular-surface, transport, and interpolation workflows.
- Visualization and structures: The Surface code generates isosurfaces from volumetric wavefunctions, charge densities, or local potentials and exports them for POV-Ray rendering.It reads Gaussian Cube or XCrySDen XSF files and uses marching-cubes or marching-tetrahedra triangulation.
- Visualization and structures: Matter manipulates periodic atomic structures, including translations, rotations, supercell generation, fragment assembly, and format conversion.Supported formats include mat, VASP, Espresso, SIESTA, XYZ, XSF, Wien, and POV-Ray.
- Band-structure interpolation: The sig2wan utility replaces mean-field eigenvalues with quasiparticle eigenvalues in Wannier90 data for band-structure interpolation.It is unsuitable when only some eigenvalues of entangled bands are replaced.
- Band-structure interpolation: The inteqp utility instead interpolates quasiparticle corrections directly to construct band structures along arbitrary directions.The corrections are described as smoother functions of k and energy than quasiparticle eigenvalues.
- File utilities: The mf utility converts binary and ASCII wavefunction, density, and exchange-correlation-potential files, with format detection based on file headers or wrapper invocation.ASCII files improve portability, while binary files are more compact.
- Image-charge model: The icm.x utility models molecular self-energy corrections near metallic surfaces as isolated-molecule corrections plus an image-charge screening term.The approximation is intended for scanning-tunneling spectroscopy and molecular-junction transport calculations.
11. Appendix
The appendix specifies BerkeleyGW’s shared file formats, required inputs, and executable-specific wavefunction grids. It also documents header contents, data organization, and utilities for inspecting and generating these files.
- Executable inputs: BerkeleyGW executables require different wavefunction files, including unshifted and shifted grids, while sigma additionally uses RHO and VXC or precomputed vxc.dat matrix elements.epsilon uses WFN and WFNq; kernel uses WFNco; absorption and plotxct use WFNco, WFNqco, and WFNqfi.
- Common file format: These input files share a common header format containing spin, reciprocal-space, symmetry, atomic, k-point, band, energy, and occupation information.The header also records FFT and k-point grids, lattice and metric tensors, symmetry operations, atomic positions, k-point weights, and band ranges.
- Utilities and interoperability: The wfninfo.x utility reports file-header information, while shared read/write routines support mean-field codes in producing BerkeleyGW-compatible output.The library routines are used by Octopus and PARATEC.
- Common file format: G-vectors are stored as integer reciprocal-lattice coordinates with full-sphere coverage, and data is organized by G-vector and spin records.G-vector components use [−n/2, n/2), where n is the FFT-grid size; records may be split into groups for file writing.
- Exchange-correlation data: vxc.dat stores exchange-correlation matrix elements by k-point, spin, diagonal, and offdiagonal entries, with real and imaginary components recorded in eV.The recommended precomputed-matrix-element scheme applies when VXC cannot represent the exchange-correlation treatment directly; selected hybrid calculations can compensate using bare-exchange fractions in sigma.inp.