Stochastic Pseudobands

Stochastic pseudobands is a method to dramatically accelerate many-body perturbation theory (MBPT) calculations such as GW and GW-BSE calculations performed by BerkeleyGW. The method reduces the scaling of the GW calculations from quartic to *quasi-quadratic, and reduces the algorithmic prefactor, with *speedups of 10-1000 $\times$ ** for typical system sizes of 10's to 100's of atoms.

Stochastic pseudobands is a general technique — it is agnostic to the mean-field code that generated the wavefunctions and to the MBPT code that will use the output compressed wavefunctions. For more details on the method and implementation, please refer to the following work, and please cite this work if using stochastic pseudobands within BerkeleyGW or anywhere else:

A. R. Altman, S. Kundu, F. H. da Jornada, Mixed Stochastic-Deterministic Approach for Many-Body Perturbation Theory Calculations, Phys. Rev. Lett. 132, 086401 (2024)

Overview

The method works by approximating high-energy contributions to the Green's function in its Lehmann representation with effective stochastic vectors. This approach is fundamentally based on the stochastic resolution of the identity. In practice the code achieves this by compressing high-energy mean-field states into effective stochastic vectors, stochastic pseudobands. Here, "high-energy" means anything beyond the energy scale you care about, but this depends on the calculation (see below). The low-energy region that you care about, centered around the Fermi level, is referred to as the protection window. Everything outside the protection window is compressed, while states within the protection window are computed exactly with the usual deterministic formalism. The protection window is necessary for computing quasiparticle energies in Sigma or quantities related to Kernel or Absorption.

In general, for the computation of the dielectric function with a plasmon-pole model, no protected states are forally required. However, quasiparticle energies and BSE matrix elements must be computed on deterministic states, so these calculations require protected states. Additionally, if the dielectric function is computed explicitly in frequency, then a protection window must also be used (see the python documentation below). It is often more convenient to use a protection window for all calculations, so that a single wavefunction file can be reused.

The stochastic pseudobands approach converges rapidly with stochastic parameters and convergence testing is almost never needed if the well-tested default parameters are used. The error in quasiparticle energies when using the default values for stochastic pseudobands is typically on the order of 50 meV as compared to a highly-converged deterministic calculation.

Applicable Systems

The use of stochastic pseudobands is applicable to any system that can be treated with MBPT, including any dimension, screening response, and size. The approach has been explicitly tested for: - 3D, 2D, 1D, and 0D systems - Semiconducting and metallic systems - Spin-polarized and fully-relativistic spinor calculations - All flavors of dielectric calculations: - Plasmon-pole models - Full-frequency calculations - Coulomb truncation schemes - Subsampling schemes (NNS) - System sizes up to 400k G-vectors - This is equivalent to several hundreds to low-thousands of atoms depending on the dimension and wavefunction cutoff - Larger systems can also be handled, if the mean-field states are obtained

How Stochastic Pseudobands are Constructed

An important aspect to using the method is to understand how the code takes mean-field states and compresses them into stochastic pseudobands. First, the states in the protection window are chosen by the user with the flag protected_cond_bands. This flag is not a convergence parameter, and is free for the user to choose. Then, the code partitions the remainder of the states based on their energy relative to the Fermi level into subspaces/slices. The subspaces $\{S\}$ are chosen such that the relation \begin{align} \dfrac{\Delta E_{S}}{\bar E_{S}} &\equiv \mathcal{F} = \text{const} \end{align} is uniformly satisfied, where $\bar E_{S}$ is the average energy of the mean-field states in each subspace $S$ (referenced to the Fermi level) and $\Delta E_{S}$ is the energy range spanned by the states in $S$ . This is Eq 4 from the reference paper above. $\mathcal{F}$ is a convergence parameter for the stochastic pseudobands method, and isset by the flag accumulation_window. Finally, from the states in each subspace $S$ , we construct $N_\xi$ random complex linear combinations to obtain $N_\xi$ stochastic pseudobands. This is to make use of the stochastic resolution of the identity. All stochastic pseudobands in a given subspace $S$ are assigned the energy $\bar E_{S}$ and are normalized to have weight $\dfrac{|S|}{N_\xi}$ . $N_\xi$ is a convergence parameter and is set by the flag num_per_subspace.

Basic Usage

There are two implementations of stochastic pseudobands in BerkeleyGW, in Fortran and Python. Both versions require HDF5 support, handle all systems from above, and require complex input files.

See the workflow diagram for how the two versions interact and can be used. Below the diagram are descriptions of the codes.

Recommended Workflow

Typical ParaBands/Pseudobands workflow

Blue text refers to pseudobands-related quantities and WFN_{v/c} refers to valence/conduction pseudobands. Advanced users may prefer to perform a full diagonalization in ParaBands and then both conduction and valence pseudobands using the Python implementation of stochastic pseudobands. Such an approach has greater flexibility but comes with several drawbacks including much larger I/O and memory usage.

Typically, valence pseudobands are not required/helpful for systems with fewer than ~30 occupied bands.

Fortran

Stochastic pseudobands is bundled with the ParaBands code. This version has the following features: - Fully parallel distributed matrix manipulations of the full WFN already loaded in memory. - Negligible I/O due to much smaller output WFN. - Compression of unoccupied states only, - WFN without valence compression is required for Sigma calculations (conduction compression is recommended for Sigma, however). - Python version can be used for valence compression for Epsilon calculations.

The Fortran implementation also has the following restrictions: - The full WFN, although computed from exact diagonalization of the Hamiltonian, is not stored. Only the compressed WFN with stochastic pseudobands is stored. Thus, it is not possible to increase the protection window after-the-fact. This is potentially problematic if you realize you need to compute more quasiparticle energies half-way through your workflow, as Sigma, Kernel, and Absorption cannot be computed on the compressed pseudobands states, and only on the deterministic protected states. - Convergence testing with respect to stochastic parameters num_per_subspace and accumulation_window is difficult for the same reason. - Because ParaBands currently only takes input from QuantumEspresso, this implementation of stochastic pseudobands is restricted to this too. - Output WFN file is not compatible with full-frequncy dielectric calculaitons.

See parabands.inp for input parameters.

Python

Stochastic pseudobands is also implemented as a standalone Python script pseudobands.py. This version has the following features: - Agnostic to the mean-field input, as long as input WFNs are in [BerkeleyGW format](wfn_h5_spec.md). If you use one of the wrappers from BerkeleyGW (*2bgw.x`) then this requirement will be satisfied. - Compression of both occupied and unoccupied states. This is the recommended approach for Epsilon calculations (see the workflow diagram below). - Simultaneously compresses WFN and WFNq if using valence pseudobands to guarantee consistency. - Optional input for subsampled WFNq (e.g. from NNS). - If WFNq is not desired, a dummy file may be provided. - Extra parameters change stochastic partitioning scheme to make output WFN compatible with full-frequency dielectric calculations. - This is an advanced use case. See the Supplemental Material of the reference paper for explanation. - Simple convergence testing and adjustment of protection window if needed.

The Python version also has the following restrictions: - Serial implementation (only uses implicit threading from numpy) - If used in conjunction with an exact diagonalizer like ParaBands, I/O can be slow and can run into memory issues due to extremely large input WFN - Dependencies: - h5py - numpy - scipy

Input parameters to `pseudobands.py`

Required Keywords

--fname_in - Description: Input WFN.h5, in HDF5 format.

--fname_in_q - Description: Input WFNq.h5, in HDF5 format.

--fname_out - Description: Output WFN.h5 with pseudobands, in HDF5 format.

--fname_out_q - Description: Output WFNq.h5 with pseudobands, in HDF5 format.

(Strongly) Recomended Keywords

--N_P_val - Description: Number of protected valence bands counting from VBM. - Type: int - Default: -1

--N_P_cond - Description: Number of protected conduction bands counting from CBM. - Type: int - Default: 100 - Equivalent to protected_cond_bands in Fortran version.

--N_S_val - Description: Number of subspaces spanning the total energy range of the valence bands. - Type: int - Default: 10

--N_S_cond - Description: Number of subspaces spanning the total energy range of the conduction bands. - Type: int - Default: 100 - Inversely related to accumulation_window in Fortran version.

--N_xi_val - Description: Number of stochastic pseudobands per valence slice. Must be at least 2. - Type: int - Default: 2

--N_xi_cond - Description: Number of stochastic pseudobands per conduction slice. Must be at least 2. - Type: int - Default: 2 - Equivalent to num_per_slice in Fortran version.

Optional Keywords

--fname_in_NNS - Description: Input WFNq_NNS.h5, in HDF5 format.

--fname_out_NNS - Description: Output WFNq_NNS.h5 with pseudobands, in HDF5 format.

--NNS - Description: Using separate NNS WFNq?

--max_freq - Description: Maximum energy (Ry) before standard slicing kicks in. Used for full-frequency dielectric calculations. This should be at least the maximum frequency for which you plan to evaluate epsilon. If uniform_width is not set then this flag is ignored. - Type: float - Default: 0.0

--uniform_width - Description: Constant width accumulation window (Ry) for slices with |energy - E_F| <= max_freq. - Type: float - Default: None

--verbosity - Description: Set verbosity level. - Type: int - Default: 2

Example Run Command:

python pseudobands.py --fname_in WFN.h5 --fname_in_q WFNq.h5 --fname_out WFN_SPB.h5 --fname_out_q WFN_SPB_q.h5 --N_P_val 10 --N_P_cond 10 --N_S_val 10 --N_S_cond 150 --N_xi_val 2 --N_xi_cond 2