# 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

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
```