Many-body perturbation theory calculations using the yambo code

yambo is compiled using the standard autotools procedure: ./configure; make all will generate the main executables listed in figure 2. Since the first release the configure script has been wholly upgraded to reflect the widespread availability of high performance software libraries and to aid portability across a wider range of system architectures.

By running ./configure; make, the list of possible executables is returned

[all projects] all

[project-related suite] project

     (core, rt-project, ...)

[core] yambo

[core] ypp

[core] a2y

[core] p2y

[ph-project] yambo_ph

[ph-project] ypp_ph

[rt-project] yambo_rt

[rt-project] ypp_rt

[nl-project] yambo_nl

[nl-project] ypp_nl

[kerr-project] yambo_kerr

While yambo and ypp are the main code components, a series of projects appear in the form of yambo_PJ/ypp_PJ with PJ being the specific project identifier (ph,rt,nl,kerr). These projects correspond to pre-processor flags that, during the compilation, activate lines of code and procedures that are project-specific. In yambo several different codes coexist in the same source.

In many cases configure will manage to detect the compilation environment and external libraries automatically. For more control, a flexible list of options is available (see ./configure –help). A wide range of optional features can be activated via –enable-FEATURE[=ARG] flags, e.g.

./configure –enable-open-mp

including options controlling serial/parallel linear algebra, timing/memory profiling, type of fast Fourier transform (FFT) library, etc. External libraries can be linked to by specifying either the installation directory including the 'libs' and 'include' folders,

–with-libname-path  =  <path  >  

or the 'libs' and 'include' paths directly

–with-libname-libdir  =  <path  >  

–with-libname-includedir  =  <path  >  

or finally the libraries and the include command

–with-libname-libs  =  <libs  >  

–with-libname-incs  =  <include command  >  

This is an important improvement for allowing installation on machines with non-standard system directories.

Choice of compilers and preprocessors can be overridden via the environmental variables FC, CPP etc. Finally, the generated config/setup file can be tweaked by hand prior to compilation.

yambo can make use of several external libraries for improving performance and portability (see table 1). In addition to standard MPI (openmpi, Intel MPI, etc) and OpenMP for parallel computation, these include standard scientific computation libraries such as BLAS and LAPACK (including the Intel MKL and IBM ESSL), scalable versions of these (BLACS, ScaLAPACK; –enable-par-linalg), as well as advanced parallel numerical libraries (SLEPc, PETSc; –enable-slepc-linalg). Use of the latter in yambo is discussed in detail in section 8.6. Heavy use is made of FFTs. yambo supports many FFT implementations: Goedecker (–enable-internal-fftsg), FFTW (internal default) and 3D or standard FFT implementation of Quantum ESPRESSO (–enable-3D-fft or –enable-internal-fftqe) can be compiled while MKL and ESSL can be externally linked. Regarding internal I/O, linking to NetCDF or HDF5 format libraries is a requirement. The exchange-correlation functional library libxc is also required. Interfacing with the yambopy and AiiDA platforms is explained thoroughly in section 9. Libraries related to porting data from DFT codes are discussed in the following section.

2.2.1. External libraries.

yambo is interfaced with two widely used plane-wave first-principles codes: pwscf from the Quantum ESPRESSO (QE) distribution [39, 40] and Abinit [5, 41, 42]. The two interfaces have been introduced in [4] (sections 5.1 and 5.3). Both work with norm conserving pseudo-potentials and import Kohn–Sham (KS) eigen-energies $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\kk}{{\bf k}} \newcommand{\e}{{\rm e}} \epsilon_{n\kk}$ and eigen-functions $\newcommand{\kk}{{\bf k}} \psi_{n\kk}$ as well as information needed to compute the non-local part of the pseudo-potential $\newcommand{\xx}{{\bf x}} V^{nl}(\xx,\xx')$. Since the publication of [4] both interfaces have been largely improved and extended. All interfaces are now able to deal with both collinear and non-collinear spin systems. All interfaces take advantage of the XC library [43, 44], thus a very broad class of functionals is supported. A more detailed summary of the changes follows.

2.3.1. Interface with Quantum ESPRESSO.

2.3.2. Interface with Abinit.

ypp is the yambo postprocessing and preprocessing tool. It has several capabilities which can be used to prepare yambo simulations (preprocessing) or subsequently analyse (postprocessing) the outcome.

As one of the preprocessing options, ypp can generate random grids of $\newcommand{\kk}{{\bf k}} \kk$-points to be used as input for a DFT code to compute the corresponding KS energies. The same ypp can then generate an auxiliary database with a map linking the KS energies on the random grid to the uniform grid used to compute, for example, spectral properties. The approach is useful to speed up convergence as discussed in section 5.1.1. Another preprocessing option is the removal of a specific set of symmetries and thus the expansion of the wave-functions from the IBZ associated to the full set of symmetries to the resulting IBZ. This is needed to perform real-time simulations as described in section 7. Finally preprocessing can also be used to map DFT calculations with and without spin–orbit coupling (SOC) to compute absorption spectra with SOC corrections included in a perturbative way as described in the supplemental material of [48].

Most of the postprocessing features involve data analysis. ypp can prepare readable ascii files to plot several single-particle properties such as wave-functions, charge density, density of states (DOS), magnetization, current and band structures. In particular, it can be used to obtain the QP-DOS and to interpolate QP-energies to plot the resulting band-structure along high-symmetry paths. A mixed feature (i.e. which can be used both for preprocessing and for postprocessing) is the ability of ypp to manipulate QP-databases (ndb.QP). Indeed, this is useful both for QP plots or for using ndb.QP files as input in the BSE calculations. Finally, it can be used as a tool to analyse the excitonic wave-function. As examples of postprocessing, we discuss in detail (i) how to plot the QP band structure starting from calculations on a regular grid in section 4.4 and (ii) how to plot the excitonic wave-function in section 5.3.

yambo relies on a powerful and user friendly command line interface for generating and modifying input files as well as for launching the executables. The basic functionality is unchanged from that described in CPC2009; however, some flags have been changed since the initial release. Several new options have been added to aid usage or debugging on parallel clusters or cross-compiled architectures. For instance, yambo -M and yambo -N switch off the MPI and OpenMP functionalities, respectively, yambo -Q stops the text editor from launching, and yambo -W  <opt  >  places an internal wall clock limit on the runtime. Launching yambo -H shows the fully updated list of command options: see table 2.

Despite the computational cost, the numerical algorithm to compute the terms in equation (2) is straightforward. Since absorption is defined as the macroscopic average of the density-density response function, $\newcommand{\qq}{{\bf q}} \chi(\qq\to 0)$, the knowledge of $\newcommand{\kk}{{\bf k}} \newcommand{\qq}{{\bf q}} \rho_{nm\kk}(\qq\to 0,0)$ is also needed. To this end, the dipole matrix elements $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \newcommand{\kk}{{\bf k}} \newcommand{\rr}{{\bf r}} {\rr_{nm\kk}=\langle n\kk |{\mathbf{r}} | m \kk \rangle}$ are commonly computed [1, 49] within periodic boundary conditions (PBC) using the relation $\newcommand{\rr}{{\bf r}} {[\rr,H]=\mathbf{p}+[\rr,V_{nl}]}$. Explicitly, this gives

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rr}{{\bf r}} \newcommand{\kk}{{\bf k}} \newcommand{\ra}{\rangle} \newcommand{\la}{\langle} \newcommand{\ep}{\ege_{p}} \newcommand{\ege}{\varepsilon} \displaystyle \langle n\kk |\rr | m \kk \rangle=\frac{\langle n\kk | \,{\mathbf{p}}+[\rr,V_{nl}] \,| m \kk \rangle} {\epsilon_{n\kk}-\epsilon_{m\kk}}\, . \label{eq:dip} \nonumber \end{align} \tag{ 3 }$$

The direct evaluation of equation (3) (G-space v approach) is quite demanding due to the $\newcommand{\rr}{{\bf r}} [\rr,V_{nl}]$ term, and is evaluated from the KB form factors loaded by the interfaces, see also section 2.3. This implementation has been strongly optimized and extended to account for projectors with angular momentum l  >  2.

We have also made available alternative strategies for computing the dipoles. The shifted grids approach is based on the idea of numerically evaluating $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\kk}{{\bf k}} \newcommand{\qq}{{\bf q}} \newcommand{\e}{{\rm e}} \rho_{nm \kk}(\qq_\epsilon,0)$ for a very small $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\qq}{{\bf q}} \newcommand{\e}{{\rm e}} q_\epsilon=|\qq_\epsilon|$. Thus the wave-function at $\newcommand{\kk}{{\bf k}} \kk$ and the wave-functions at $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\kk}{{\bf k}} \newcommand{\qq}{{\bf q}} \newcommand{\e}{{\rm e}} \kk-\qq_\epsilon$ are needed. Since the $\newcommand{\qq}{{\bf q}} \qq\rightarrow 0$ limit may be direction dependent, this is done in practice by means of wave-functions computed on four different grids in $\newcommand{\kk}{{\bf k}} \kk$-space, i.e. a starting $\newcommand{\kk}{{\bf k}} \kk$-grid plus three grids with $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\kk}{{\bf k}} \newcommand{\cc}{{\bf e}} \newcommand{\e}{{\rm e}} \kk+q_\epsilon\,\cc_i$ slightly shifted along the three Cartesian directions $\newcommand{\cc}{{\bf e}} \cc_x, \cc_y, \cc_z$. Such approach is computationally more efficient, although it requires to generate a larger set of wave-functions. However, there exists a random phase associated to the wave-functions on each of the four $\newcommand{\kk}{{\bf k}} \kk$-grids, since they are obtained by independent diagonalizations of the KS Hamiltonian. Because of this, shifted grids dipoles have inconsistent phases among different directions and it is not possible to use them when the dipole matrix elements are needed (instead of their square modulus only) as for example in the evaluation of the Kerr effect (see section 5.2.1) or for non-linear optics (see section 7.2).

The G-space v approach assumes that the only non-local terms in H are the kinetic energy and the pseudopotentials. There are however cases, for example when the Hamiltonian contains non-local hybrid functionals, Hubbard U terms, or nonlocal self-energies, in which the evaluation of the commutator may become very cumbersome. To solve this issue one could in principle use the shifted grids approach. However, this approach may also become impractical because of the calculation of wave-functions on the shifted grids.

For those cases we have implemented in yambo two alternative strategies, one for extended and one for isolated systems. For extended systems the Covariant approach exploits the definition of the position operator in $\newcommand{\kk}{{\bf k}} \kk$ space: $\newcommand{\kk}{{\bf k}} \newcommand{\rr}{{\bf r}} \hat{\rr} = {\rm i} \partial_\kk$. The dipole matrix elements are then evaluated as finite differences between the $\newcommand{\kk}{{\bf k}} \kk$-points of a single regular grid. A five-point midpoint formula is used, with a truncation error $\newcommand{\kk}{{\bf k}} O(\Delta \kk^4)$. The shifted grids and the Covariant approach are very similar, however in the latter the arbitrary phase of the wave-functions at different $\newcommand{\kk}{{\bf k}} \kk$-points is correctly accounted for. To this aim, $\newcommand{\kk}{{\bf k}} i \partial_\kk$ is implemented as a covariant derivative which cancels the relative phase factor (see appendix D for details). For these reasons the Covariant approach overcomes the limitations of the shifted grids approach. The main drawbacks of the Covariant approach is that the numerical value of the dipoles needs to be converged against the size of the $\newcommand{\kk}{{\bf k}} \kk$ grid and the present implementation does not work for metals. However, in practice the convergence of dipole matrix elements is usually faster than that of the absorption spectrum.

For finite systems, finally, the dipole matrix elements can be directly evaluated in real space (R-space x approach).

We underline that in the case of a local Hamiltonian all approaches are equivalent. The desired strategy can be selected via the input variable:

DipApproach="G-space v"

#[Xd] [G-space v/R-space x/Shifted grids/Covariant]

(G-space v being the default value).

The Coulomb interaction enters in many sections of the yambo code, such as linear response, self-energy, and BSE kernel calculation. In reciprocal space, the bare Coulomb interaction for bulk systems is defined as $v({\bf q+G})=4\pi/|{\bf q+G}|^2$. For the calculation of quantities requiring integration over transferred momenta in the Brillouin zone (BZ), such as the self-energy, the integrals are evaluated by summations over regular $\newcommand{\qq}{{\bf q}} \qq$-grids. In order to remove divergencies in systems of reduced dimensionality, i.e. in the presence of a 2D or 1D sampling of $\newcommand{\kk}{{\bf k}} \kk$-points, or to speed up the convergences in 3D systems, yambo offers the possibility to evaluate Coulomb integrals by using the random integration method (RIM), which consists of evaluating these integrals by Monte Carlo sampling (as already discussed in detail in section 3.1 of [4]), dividing the full BZ in small regions around each $\newcommand{\kk}{{\bf k}} \kk$-point of the chosen uniform grid.

In order to avoid spurious interaction between replicas when dealing with low-dimensional materials such as clusters, slabs, or wires, yambo can also use Coulomb cutoff truncation techniques. These consists of truncating the Coulomb interaction beyond a certain region (depending on the chosen geometry):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{v}(\mathbf{r}) = \left\{\begin{array}{@{}lll@{}} 1/r & {{\rm if~}} \it \mathbf{r} \in D \nonumber \\ 0 & {{\rm if~}} \it \mathbf{r} \notin D. \end{array}\right. \nonumber \end{align} \tag{ 4 }$$

Different geometrical choices are available. Spherical and cylindrical shapes, suitable to treat 0D and 1D systems, respectively, have been already described in details in [50]. In addition, a box-like cutoff obtained by performing a numerical Fourier transform of the real space expression is available for 0D systems. By defining only one or two sides of the box, it is possible to treat 2D or 1D systems within the same numerical approach. It is important to stress that, as the construction of such potential requires integration over the BZ, the RIM method discussed above must be activated.

Finally, a Wigner–Seitz truncation scheme, similar to the one discussed in [51] is also available. In this scheme the Coulomb interaction is truncated at the edge of the Wigner–Seitz super-cell compatible with the $\newcommand{\kk}{{\bf k}} \kk$-point sampling. This truncated Coulomb potential turns out to be suitable for finite systems as well as for 1D and 2D systems, provided that the supercell size, determined by the adopted $\newcommand{\kk}{{\bf k}} \kk$-point sampling, is large enough to get converged results [52].

The independent particle polarizability $\chi_{{\bf{G}}{\bf{G'}}}^0 ({\bf{q}}, \omega)$, equation (1), and the correlation part of the GW self-energy $\Sigma_c(\omega)$, equation (7) in section 4.1, are evaluated through sum-over-states (SOS) expressions obtained by applying an energy cutoff to the infinite sum over virtual states. These expressions are, however, slowly convergent and, especially for large systems, require the inclusion of a large number of empty states (N_b). This condition makes GW calculations computationally demanding, both in terms of time-to-solution and memory requirements. In order to overcome this limitation, a number of approaches have been proposed to reduce [53–56] or remove [3, 8, 57] sum over states; among them, we have implemented in yambo the extrapolar correction scheme proposed by Bruneval and Gonze (BG) [53].

This scheme, here referred as X-terminator, permits to accelerate GW convergence by reducing of a sensible amount the number of virtual orbitals necessary to calculate both polarizability and self-energy. In this procedure extra terms, whose calculation implies a small computational overhead, are introduced to correct both polarizability and self-energy by approximating the effect of the states not explicitly taken into account. The method consists in replacing the energies of empty states that are above a certain threshold, and that are not explicitly treated, by a single adjustable parameter defined as extrapolar energy. When the method of terminators is applied, the independent-particle polarizability can be written as [53]:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\qq}{{\bf q}} \newcommand{\la}{\langle} \newcommand{\ep}{\ege_{p}} \newcommand{\ege}{\varepsilon} \displaystyle \chi_{{\bf{G}}{\bf{G}'}}^0 (\qq,\omega)= \chi^{0,{\rm trunc}}_{{\bf{G}}{\bf{G}'}} (\qq,\omega) + \Delta\chi_{{\bf{G}}{\bf{G}'}} (\qq,\omega, \bar{\epsilon}_{\chi_0}) \label{eq:X_term_def} \nonumber \end{align} \tag{ 5 }$$

where the first term on the rhs is truncated at the $N_b^{\prime}$ state (in general $N_b^{\prime} \ll N_b$) and the second term depends on the extrapolar energy for the polarizability $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\e}{{\rm e}} \bar{\epsilon}_{\chi_0}$. The explicit expression for $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\qq}{{\bf q}} \newcommand{\e}{{\rm e}} \Delta\chi_{{\bf{G}}{\bf{G}'}}(\qq,\omega, \bar{\epsilon}_{\chi_0})$ is provided in appendix C.

In the present implementation of yambo, the input parameter governing the use of the terminator corrections on the response function (X-terminator) is

XTermKind  =  "none" # [X] X terminator ("none","BG")

(default: "none"). When the variable is set to none (default option), the X-terminator is not applied. On the contrary when XTermKind assumes the value BG, the extrapolar corrective term is calculated. The extrapolar energy $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\e}{{\rm e}} \bar{\epsilon}_{\chi_0}$, see equation (C.3), is defined by the input variable (default: 1.5 Ha)

XTermEn  =  1.5 Ha # [X] X terminator energy

The value means $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\kk}{{\bf k}} \newcommand{\e}{{\rm e}} \bar{\epsilon}_{\chi_0}=\epsilon_{N_{b'} \kk}+1.5\ {{\rm Ha}}$, with $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\kk}{{\bf k}} \newcommand{\e}{{\rm e}} \epsilon_{N_{b'} \kk}$ the highest energy state included in the calculation.

For demonstration purposes, in figure 3 we report the calculated QP corrections for the valence band maximum (VBM) and the conduction band minimum (CBM) of a bulk Si described in a supercell (36 Si atoms, 72 occupied states). Results are obtained by increasing the number of bands explicitly included in the calculation of the response function $\chi$ and by imposing a very high number of bands in the self-energy, that is therefore converged. Empty circles connected with solid lines denote the results obtained for the VBM and CBM states without applying any correction. Improvements induced by the use of the X-terminator are depicted by solid circles connected with solid lines that have been obtained imposing XTermEn  =  1.5 Ha. We can observe that the X-terminator leads to a relevant reduction in the number of bands necessary to converge the polarizability and thus the GW corrections.

Zoom In Zoom Out Reset image size

Within the GW approximation, the matrix elements of the correlation self-energy over the KS basis are expressed as:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\qq}{{\bf q}} \newcommand{\kk}{{\bf k}} \newcommand{\ra}{\rangle} \newcommand{\la}{\langle} \newcommand{\ep}{\ege_{p}} \newcommand{\ege}{\varepsilon} \displaystyle \langle n\kk |\Sigma_c (\omega)| n^\prime \kk \rangle = \sum_{m\qq} \int \frac{{\rm d}\omega'}{2\pi {\rm i}} I^{nn'\kk}_{m\qq} (\omega^\prime) \left[ \frac{f_{m \kk-\qq}\theta(\omega^\prime)}{\omega-\omega'-\epsilon_{m \kk-\qq} -{\rm i}\eta}+ \frac{(1-f_{m \kk-\qq})\theta(-\omega^\prime)}{\omega-\omega'-\epsilon_{m \kk-\qq} +{\rm i}\eta} \right]. \label{eq:sigma_c} \nonumber \end{align} \tag{ 7 }$$

I is linked to the self-energy spectral function. From a computational point of view the definition of I is really critical as, in equation (7), it is connected to the self-energy via a complex Hilbert transformation. In yambo I is defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\GG}{{\bf G}} \newcommand{\qq}{{\bf q}} \newcommand{\kk}{{\bf k}} \newcommand{\la}{\langle} \displaystyle I^{nn'\kk}_{m\qq} (\omega^\prime)= -\frac{1}{N_k \Omega}\sum_{\GG \GG'} W^{\delta}_{\GG \GG'}(\qq,\omega') \times \rho_{nm\kk}( \qq, \GG) \rho^*_{n'm\kk}(\qq, \GG'). \label{eq:I_c} \nonumber \end{align} \tag{ 8 }$$

In equation (8), $W^{\delta}$ is the delta-like part of the screened interaction. This is defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\GG}{{\bf G}} \newcommand{\qq}{{\bf q}} \newcommand{\la}{\langle} \displaystyle W^{\delta}_{\GG \GG'}(\qq,\omega)=\left[\frac{1}{2}\Im\left(W_{\GG \GG'}(\qq,\omega)+W_{\GG' \GG}(\qq,\omega)\right)-\frac{{\rm i}}{2}\Re\left(W_{\GG \GG'}(\qq,\omega)-W_{\GG' \GG}(\qq,\omega)\right)\right]. \label{eq:W_delta} \nonumber \end{align} \tag{ 9 }$$

In equation (9) $\newcommand{\qq}{{\bf q}} \newcommand{\GG}{{\bf G}} W_{\GG \GG'}(\qq,\omega)$ is the screened Coulomb potential defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\GG}{{\bf G}} \newcommand{\qq}{{\bf q}} \newcommand{\la}{\langle} \newcommand{\ep}{\ege_{p}} \newcommand{\ege}{\varepsilon} \displaystyle W_{\GG \GG'}(\qq,\omega)=\epsilon^{-1}_{\GG \GG'}(\qq,\omega)\frac{4\pi}{|\qq+\GG||\qq+\GG'|}. \label{eq:W_def} \nonumber \end{align} \tag{ 10 }$$

Note that, in the case of systems with both spatial and time reversal symmetry, $\newcommand{\qq}{{\bf q}} \newcommand{\GG}{{\bf G}} W_{\GG \GG'}(\qq,\omega)=W_{\GG' \GG}(\qq,\omega)$ and $\newcommand{\qq}{{\bf q}} \newcommand{\GG}{{\bf G}} W_{\GG \GG'}(\qq,\omega)$ reduces to the imaginary part of W.

In order to take into account the frequency dependence of the self-energy, two different strategies are implemented in yambo. As already described in [4], it is possible to adopt the plasmon-pole approximation (PPA) in order to model the dynamic screening matrix. This approximation essentially assumes that all the spectral weight of the dielectric function is concentrated at a plasmon excitation. Among different models present in the literature yambo implements the Godby–Needs construction [65] where the parameters of the model are chosen in such a way that $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\e}{{\rm e}} \epsilon^{-1}_{{\bf G G^\prime}}({\bf q},\omega)$ is reproduced at two different frequencies: the static limit $\omega=0$ and another imaginary frequency $\omega={\rm i}\omega_p$ given in the input file by PPAPntXp (default: 1 Ha). Quasi-particle energy levels calculated within this approximation have been shown to agree to a large extent with numerical integration methods for materials with different characteristics including semiconductors and metal-oxides [59, 66]. Moreover, it has the great advantage to avoid the computation of the inverse of the dielectric matrix for many frequency points and to make the frequency integral of equation (7) expressible in an analytic form. Nevertheless the assumption made for the PPA breaks down in certain situations as when dealing with metals [67–69] or interfaces [70] and the frequency integral needs to be solved numerically. In yambo the integral is solved on the real-axis which implies the knowledge of the full frequency dependence of $\newcommand{\qq}{{\bf q}} \newcommand{\GG}{{\bf G}} W_{\GG \GG'}(\qq,\omega)$. In practice, first the inverse dielectric function $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\qq}{{\bf q}} \newcommand{\GG}{{\bf G}} \newcommand{\e}{{\rm e}} \epsilon^{-1}_{\GG \GG'}(\qq,\omega)$ is evaluated for a number of frequencies set by the variable ETStpsXd, and uniformly distributed in the energy range given by the maximum electron–hole pairs included in the response function defined in equation (1). Next, the summation over $\newcommand{\GG}{{\bf G}} \GG$ and $\newcommand{\GG}{{\bf G}} \GG^{\prime}$ is performed computing $\newcommand{\kk}{{\bf k}} \newcommand{\qq}{{\bf q}} I^{nn'\kk}_{m\qq} (\omega^\prime)$ defined in equation (8), and finally the correlation part of the self-energy is computed via a Hilbert transform defined in equation (7).

In this scheme, the evaluation of equation (10) is the most time consuming step due the computation of the inverse dielectric matrix for a large number of frequencies (order of 100) in order to have converged results. Nevertheless as the calculations for each frequency are independent from each other, parallelization over frequencies provides a linear speedup.

Quasi-particle energies calculated by using the real-axis method have been demonstrated to provide the same level of accuracy of other beyond plasmon-pole techniques such as the contour deformation scheme [71].

The ability of yambo to calculate the real-axis GW self-energy allows direct access to the quasi-particle electron-mediated lifetimes. Indeed if we define $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} \newcommand{\kk}{{\bf k}} \newcommand{\e}{{\rm e}} \Gamma^{\rm e{{\rm \mbox{--}}}e}_{n\kk}(\omega)\equiv \Im\left(\langle n\kk |\Sigma_c (\omega)| n \kk \rangle\right)$, from equation (7) it is easy to see that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\qq}{{\bf q}} \newcommand{\kk}{{\bf k}} \newcommand{\la}{\langle} \newcommand{\ep}{\ege_{p}} \newcommand{\ege}{\varepsilon} \displaystyle \Gamma^{\rm e{{\rm \mbox{--}}}e}_{n\kk}(\omega) &= \frac{1}{2}\sum_{m\qq} I^{nn\kk}_{m\qq} (\omega-\epsilon_{m \kk-\qq}) \nonumber \\ &\quad \times \Big[ \theta(\omega-\epsilon_{m \kk-\qq})\,f_{m \kk-\qq} -\theta(\epsilon_{m \kk-\qq}-\omega)\left(1-f_{m \kk-\qq}\right) \Big]. \label{eq:ee_life} \nonumber \end{align} \tag{ 11 }$$

yambo can evaluate the quasi-particle lifetimes $\newcommand{\kk}{{\bf k}} \tau_{n\kk}(\omega)$, proportional to the inverse of $\newcommand{\kk}{{\bf k}} \Gamma^{\rm e{{\rm \mbox{--}}}e}_{n\kk}(\omega)$, either in the on-the-mass-shell approximation (OMS) or in the full GW approximation. The difference between the two is the inclusion of the renormalization factors, $\newcommand{\kk}{{\bf k}} Z_{n\kk}$. For details on the theory see, for example, [68] and references therein.

In the OMS we have that $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\kk}{{\bf k}} \newcommand{\e}{{\rm e}} \left.\Gamma^{\rm e{{\rm \mbox{--}}}e}_{n\kk}\right|_{\rm OMS}=\Gamma^{\rm e{{\rm \mbox{--}}}e}_{n\kk}(\epsilon_{n \kk})$. The e–e lifetimes of bulk copper are shown in figure 4 using several flavours of GW approximations [68].

Zoom In Zoom Out Reset image size

An important numerical property of the electron-mediated lifetimes calculation is that they depend only on the $\newcommand{\kk}{{\bf k}} \kk$-grid. Indeed, as evident from equation (11) the band summations are limited by the two theta functions that confine the scattering events in reduced regions of the BZ. This is the mechanism that, in simple metals, leads to the well known quadratic scaling of $\newcommand{\kk}{{\bf k}} \left.\Gamma^{\rm e{{\rm \mbox{--}}}e}_{n\kk}\right|_{\rm OMS}$ near the Fermi level, as a function of distance of $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\kk}{{\bf k}} \newcommand{\e}{{\rm e}} \epsilon_{n \kk}$ from the Fermi level itself.

More physical insight in the electronic lifetimes will be given in section 6.2 where the phonon-mediated case will be described.

In section 3.3 we have discussed the X-terminator procedure. A similar scheme can be adopted to study the correlation part of the GW self-energy, as from equation (16) of [53]. Also in this case the approximation implies the introduction of an extra term that takes into account contributions arising from states not explicitly included in the calculation. The input parameter governing the use of the terminator corrections on the self-energy (G-terminator) is

GTermKind  =  "none" # GW terminator ("none","BG")

(default: "none"). When the variable is set to none, the G-terminator is not applied. On the contrary when it assumes the value BG, the extrapolar corrective term is calculated. The extrapolar energy for the self-energy is defined by the tunable input variable

GTermEn  =  1.5 Ha # [X] X terminator energy

(default: 1.5 Ha).

Also in this case, the value is referenced to the highest band included in the calculation. In figure 5 we reconsider the system discussed in the example of section 3.3, figure 3. In this case, however, we study the convergence of the self-energy by exploiting the G-terminator procedure. Empty circles connected with solid lines show the usual GW convergence for the VBM and CBM states (no corrections applied). Calculations have been performed by imposing a high number of bands in the polarisability (that is therefore converged) and by increasing the number of bands included in the self-energy. We set GTermEn  =  1.5 Ha, that represents the best choice for this system. Improvements provided by the use of the G-terminator procedure are represented by solid circles connected with solid lines; it is evident that the application of this scheme accelerates the convergence by leading to a significant reduction in the number of states necessary to converge the GW self-energy and therefore the calculated QP correction.

Zoom In Zoom Out Reset image size

In order to elucidate the role played by the extrapolar parameter, we report in figure 6 a convergence study of the VBM GW correction for a ${{{\rm TiO}}}_{2}$ nanowire (NW). The black line is obtained without applying any correction. Coloured lines are instead obtained by applying both X- and G-terminators, moving the extrapolar energy from 1.0 to 3.0 Ha. Results are reported as a function of the number of states explicitly included in the calculation of both polarisability and self-energy. As pointed out in [53], the extrapolar energy for the self-energy can be safely taken equal to the extrapolar energy introduced in equation (C.3) for the polarizability; for this reason we impose XTermEn  =  GTermEn. Consistently with the study of figure 6, the convergence of the VBM without terminators is very slow and requires the inclusion of a large number of bands to be achieved; this condition makes the calculation cumbersome also on modern HPC-machines. When the terminator technique is adopted to correct both polarisability and self-energy, the convergence becomes much faster; especially for some values of the extrapolar energy (about 1.5 Ha), we observe a significant reduction in the number of bands necessary to converge the calculation, with a strong reduction of both the time-to-solution and the allocated memory. Noticeably, the correction is almost independent on the selected extrapolar energy (terminators are convergence accelerators and the extrapolar correction vanished in the limit of infinite bands included); this parameter therefore influences the number of bands necessary to converge the calculation (and thus the computational cost of the simulation) but not the final result.

Zoom In Zoom Out Reset image size

In DFT the eigenvalues of the Hamiltonian at every $\mathbf{k}$-point can be obtained by the knowledge of the ground-state charge density, allowing one to perform non-self-consistent calculations on an arbitrary set of $\mathbf{k}$-points. Instead at the HF or GW level, to obtain QP corrections for a given $\mathbf{k}$-point it is necessary to know the KS wave-functions and eigen-energies on all $(\mathbf{k+q})$-points, having chosen a regular grid of $\mathbf{q}$-points as convergence parameter. In practice yambo computes QP corrections on a regular grid. As a consequence the evaluation of band structures along high-symmetry lines can be computationally very demanding.

A simple strategy which is implemented in ypp is to interpolate the QP corrections from such regular grid to the desired high symmetry lines. The approach implemented is based on a smooth Fourier interpolation [72], which is particularly efficient for 3D grids. The interpolation scheme can also take, as additional input, the KS energies computed along the high symmetry lines to better deal with bands crossing and regions with non analytic behaviour, such as cusp-like features.

A more involved strategy is instead based on the Wannier interpolation scheme as implemented in the wannier90 [73] and WanT [74] codes, where electronic properties computed on a coarse reciprocal-space mesh can be used to interpolate onto much finer meshes at low cost [75]¹⁹. In the context of GW calculations, the Wannier interpolation scheme can be used to interpolate the QP energies and other band structure properties [74] (e.g. effective masses) from QP corrections computed only on selected $\mathbf{k}$-points. Wannier interpolation of GW band structures requires two sets of inputs: on one side quantities computed at the DFT level such as KS eigenvalues, overlaps between different KS states, and orbital and spin projections of KS states, that are imported from Quantum ESPRESSO, and on the other side the QP corrections computed by yambo. In fact, wannier90 works with uniform coarse meshes on the whole BZ, while yambo uses symmetries to compute quantities on the IBZ. In addition, converging the GW self-energy typically requires denser meshes with respect to what is needed for the charge density or Wannier interpolation. To address this issue, ypp allows one to unfold the QP corrections from the IBZ to the whole BZ, as required by wannier90 for interpolation purposes. Finally the wannier90 code yields a GW-corrected Wannier Hamiltonian and interpolates the GW band structure. A similar procedure is implemented in WanT.

For example, in monolayer WS₂ a grid of $48\times48\times1$ (or denser) is required to converge the GW self-energy. In this case, the band structure can be obtained either by explicitly computing the QP corrections on all $\mathbf{k}$-points of the $48\times48\times1$ grid, or it can be Wannier-interpolated from the QP corrections computed onto coarse subgrids, such as a $6\times6\times1$ corresponding to seven symmetry-nonequivalent $\mathbf{k}$-points only in the IBZ (see figure 7). The second approach requires substantial less CPU time.

Zoom In Zoom Out Reset image size

The computational cost of the BSE grows as a power of the number of electron–hole pairs. As this number can be as large as 10⁵–10⁶, it is crucial to devise numerically efficient algorithms for the calculation of the $V$ and W matrix elements and the solution of the BSE. The massive parallelization and memory distribution which contributes in making these calculations possible for very large systems are discussed in section 8. Here we discuss the use of the double grid for the sampling of the BZ [85], where the BSE is solved (section 5.1.1)—which aims at reducing the number of degrees of freedom involved—and the use of Lanczos-based algorithms together with the interface to the SLEPC library (Scalable Library for Eigenvalue Problem Computations) [86] (section 5.1.2)—which aims at avoiding the full diagonalization of $H^{\rm 2p}$.

5.1.1. Double-grid and the inversion solver.

5.1.2. Spectra and exciton wavefunctions via Krylov subspace methods.

Solving the BSE implies the solution of an eigenvalues problem for the two-particle Hamiltonian that in the e–h basis can result in a matrix as large as $10^6\times 10^6$. The standard dense matrix diagonalization algorithm is available in yambo through the interface with the LAPACK and the ScaLAPACK libraries [88] (section 8). Alternatively, when only the spectrum is required, yambo provides the Haydock–Lanczos solver [89]. The latter, originally developed for the Hermitian case only (see section 3.2 of CPC2009)—by neglecting the coupling between e–h at positive and negative energies within the Tamm–Dancoff approximation—has been extended to treat the full non-Hermitian two-particle Hamiltonian [90, 91]. Cases in which considering the full non-Hermitian two-particle Hamiltonian turns out to be important have been discussed in [90, 92, 93].

More recently, yambo has been interfaced with the SLEPC library [86] which uses objects and methods from the PETSC library [94] to implement Krylov subspace algorithms to iteratively solve eigenvalue problems. These are used in yambo to obtain selected eigenpairs of the excitonic Hamiltonian. This allows the user to select a fixed number of excitonic states to be explicitly calculated thus avoiding the full dense diagonalization and saving a great amount of computational time and memory. Two options are available for the SLEPC solver. The first, which is the default, uses the PETSC matrix-vector multiplication scheme; it is faster but duplicates the BSE matrix in memory when using MPI. The second, which is activated by the logical BSSSlepcShell in the input file, uses the internal yambo subroutines (the same also used for the Haydock solver); it is slower but distributes the BSE matrix among the MPI tasks. To select the part of the spectra of interest, the library allows one to use different extraction methods controlled by the variable BSSSlepcExtraction. The standard method, ritz, obtains the lowest lying eigenpairs, while the harmonic method obtains the eigenpairs closest to a defined energy. The SLEPC solver makes it possible to obtain and plot exciton wave functions (ypp -e w) in large systems where the full diagonalization might be computationally too demanding. For example, the spectrum and the wave function of the lowest-lying exciton in monolayer hBN are shown in figure 8. The BSE eigenmodes were extracted only for the two lowest-lying excitonic states, and a $\sqrt{3}\times\sqrt{3}\times 1$ supercell was used in the calculation (the SLEPC spectrum is shown in blue). The full Haydock solution is displayed with a red line for comparison.

Zoom In Zoom Out Reset image size

5.2.1. Spin–orbit coupling and Kerr.

5.2.2. Fractional occupations, gauges and more.

Other extensions have been made available. The implementation has been modified so that the excitonic matrix is now Hermitian (or pseudo-Hermitian if coupling is included) also in the presence of fractional occupations in the ground state. This is done in practice by introducing a slightly modified four-point response function $\tilde{L}$ which is divided by the square root of the occupations as discussed in equations (14)–(16) of [98]. The resulting excitonic Hamiltonian has the form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\langle} \newcommand{\ep}{\ege_{p}} \newcommand{\ege}{\varepsilon} \displaystyle \tilde{H}_{ll'}=\Delta\epsilon_l \delta_{ll'} - \sqrt{\Delta f_l} \left( v_{ll'}-W_{ll'} \right) \sqrt{\Delta f_{l'}} \label{eq:BSE_symm} \nonumber \end{align} \tag{ 12 }$$

with $\newcommand{\kk}{{\bf k}} l=\{nk\kk\}$ a super-index in the transition space, with the square root of the occupation factors appearing on the left and on the right of the BSE kernel $v-W$. This has been used to compute absorption of systems out of equilibrium, but it is also important to describe metallic systems like graphene or carbon nanotubes where excitonic effects can be non-negligible due to the reduced dimensionality of the system.

Further, it is now possible to compute the dielectric tensor starting from the different response functions, as described in [99]. Indeed, starting from the excitonic propagator L, it is possible to construct the density–density response function $\chi_{\rho,\rho}$, or the dipole–dipole response function $\newcommand{\dd}{{\bf d}} \chi_{\dd,\dd}$ at $\newcommand{\qq}{{\bf q}} \qq=0$ (length gauge), and the current–current response function $\newcommand{\jj}{{\bf j}} \chi_{\jj,\jj}$ (velocity gauge). This can be controlled by setting Gauge="length" or Gauge="velocity" in the input file (the length gauge is the default). In case the velocity gauge is chosen the conductivity sum rule is imposed unless the flag NoCondSumRule is activated in the input file. At zero momentum, changing response function is equivalent to change gauge. At finite $\newcommand{\qq}{{\bf q}} \qq$ instead the use of $\newcommand{\jj}{{\bf j}} \chi_{\jj,\jj}$ allows for the calculation of both the longitudinal and the transverse components of the dielectric function. The finite-q BSE has been implemented and it is currently under testing before its final release.

Another extension is connected to the output of a BSE run, which also generates a file with the joint density-of-states, at the IP level, and the excitonic density of states, at the BSE level. These can be used for example to visualize dark or very small intensity peaks as shown in figure 9.

Zoom In Zoom Out Reset image size

Once a BSE calculation is performed using an algorithm which explicitly computes the excitonic eigenvectors $\newcommand{\la}{\langle} \newcommand{\kk}{{\bf k}} A^{\lambda}_{cv\kk}$, several properties of the excitons can be analyzed as shown in section 5.2.2 (see figure 8). First of all, the excitonic eigenvalues $\newcommand{\la}{\langle} E_{\lambda}$ can be sorted and plotted. The so-called amplitudes and weights can also be calculated to inspect which are the main contributions in terms of single quasi-particles to a given excitonic state. The weights are defined as the squared modulus of the excitonic wavefunction $\newcommand{\la}{\langle} |A^{\lambda}_{\rm eh}|^2$ (by default only electron–hole pairs that contribute to the exciton more than 5% are considered; the threshold can be tuned by modifying the input file 'MinWeight'). The amplitudes are defined as $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\la}{\langle} \newcommand{\kk}{{\bf k}} \newcommand{\e}{{\rm e}} \sum_{cv\kk} | A^{\lambda}_{cv\kk} |^2 \delta (\epsilon_{c\kk} -\epsilon_{v\kk} -\hbar \omega)$.

Moreover the excitonic wavefunction written in real-space $\newcommand{\la}{\langle} \newcommand{\kk}{{\bf k}} \newcommand{\rr}{{\bf r}} {\Psi_\lambda(\rr_e,\rr_h)=\sum_{cv\kk} A^{\lambda}_{cv\kk} \psi^*_{v\kk}(\rr_h)\psi_{c\kk}(\rr_e)}$ can be computed. $\newcommand{\la}{\langle} \newcommand{\rr}{{\bf r}} {\Psi_\lambda(\rr_e,\rr_h)}$ is a two-body quantity or joint-correlation function. Fixing the position of the hole $\newcommand{\rr}{{\bf r}} \rr_{h}=\bar{\rr}_h$, $\newcommand{\la}{\langle} \newcommand{\rr}{{\bf r}} {|\Psi_\lambda(\bar{\rr}_h,\rr)|^2}$ provides the conditional probability of finding the electron somewhere in space. This quantity is clearly nonperiodic and its spatial decay can change from material to material, marking the distinction between Frenkel and Wannier excitons. As an alternative it is also possible to plot $\newcommand{\la}{\langle} \newcommand{\rr}{{\bf r}} {|\Psi_\lambda(\rr,\rr)|^2}$ which is instead Bloch-like.

In figure 10 we focus on two interlayer excitonic states of bilayer hexagonal boron nitride ($\newcommand{\la}{\langle} \lambda=3$ and $\newcommand{\la}{\langle} \lambda=8$). We first plot $\newcommand{\la}{\langle} \newcommand{\rr}{{\bf r}} {|\Psi_\lambda(\bar{\rr}_h,\rr)|^2}$ (top frames), then we proceed to extract more information by analyzing the phase of $\newcommand{\la}{\langle} \newcommand{\rr}{{\bf r}} \Psi_\lambda(\bar{\rr}_h,\rr)$ (lower frames). By comparing the phase for two positions of the hole related by inversion symmetry ($\newcommand{\rr}{{\bf r}} \bar{\rr}_h$ and $\newcommand{\rr}{{\bf r}} \mathcal{I}(\bar{\rr}_h)$), we see that the first exciton (figure 10(a)) is even with respect to inversion symmetry, while the second one is odd (figure 10(a)). Symmetry analysis of the wave function permits us to conclude that $\newcommand{\rr}{{\bf r}} \Psi_3(\rr_h,\rr_e)$ and $\newcommand{\rr}{{\bf r}} \Psi_8(\rr_h,\rr_e)$ transform as the A_1g and A_2u representations of point group D_3h of the lattice, respectively.

Zoom In Zoom Out Reset image size

The HAC approach is based on the static Rayleigh–Schrödinger perturbation theory, allowing one to calculate the temperature-dependent correction of the bare electronic energies, owing to the phonon field perturbation. In the QP approximation, the bare energy is instead renormalized because of the virtual scatterings described by the real part of the self-energy and it also acquires a finite lifetime due to the imaginary part of the self-energy. The eigenvalues $E_{n\bf k}(T)$ are then complex and depend on the temperature. The more the QP approximation is valid the more the renormalization factors $Z_{n\bf k}$ are close to 1, analogously to the GW method.

If the QP approximation holds, the spectral function $A_{n\bf k}(\omega,T)=\Im \left[G_{n\bf k}(\omega,T)\right]$ is a single-peak Lorentzian function centered at $\Re[E_{n\bf k}]$ with width $\newcommand{\kk}{{\bf k}} \Gamma_{n\bf k}=\Im{[E_{n\kk}]}$. In case of strong EP interaction it has been proven that the spectral function spans a wide energy range [114, 115] and the QP approximation is no longer valid.

Figure 11 shows the spectral functions (SFs) of trans-polyethylene and trans-polyacetylene, calculated at 0K. In figure 11(a) multiple structures appear in the SFs. SFs are then spread over a large energy range. In figure 11(b) a two-dimensional plot of the SFs reveals a completely different picture with respect to the original electronic band structures. Since SFs are featured by a multiplicity of structures, each carries a fraction of the electronic charge $\newcommand{\kk}{{\bf k}} Z_{n\kk}$ depriving the dominant peak of its weight. A crucial aspect is that some SFs overlap, like in the case of trans-polyacetylene, and in the end it is impossible to associate a well defined energy to the electron and to state which band it belongs to.

Zoom In Zoom Out Reset image size

By following the same strategy used in the electronic case, the phonon-mediated contribution to the electronic lifetimes can be easily calculated from equation (13). Indeed if we define $\newcommand{\kk}{{\bf k}} \newcommand{\e}{{\rm e}} \Gamma^{\rm e{{\rm \mbox{--}}}p}_{n\kk}(\omega,T)\equiv \Im \Sigma^{\rm FAN}_{n\kk}(\omega,T)$ it is easy to see that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\qq}{{\bf q}} \newcommand{\kk}{{\bf k}} \newcommand{\la}{\langle} \newcommand{\ep}{\ege_{p}} \newcommand{\ege}{\varepsilon} \displaystyle &\Gamma^{\rm e{{\rm \mbox{--}}}p}_{n\kk}(\omega,T) = \frac{\pi}{N_{\qq}} \sum_{n'\qq\lambda} |g^{\qq \lambda}_{n n'\kk}|^2 \nonumber \\ &\quad \Big[ \delta(\omega-\epsilon_{n' \kk-\qq}-\omega_{\qq\lambda}) (N_{\qq\lambda}(T)+1-f_{n' \kk-\qq})\nonumber \\ &\quad +\delta(\omega-\epsilon_{n' \kk-\qq}+\omega_{\qq\lambda}) (N_{\qq\lambda}(T)+f_{n' \kk-\qq}) \Big]. \label{eq:ep_life} \nonumber \end{align} \tag{ 15 }$$

In perfect analogy with the electronic case, within the OMS, we have that $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\kk}{{\bf k}} \newcommand{\e}{{\rm e}} \left.\Gamma^{\rm e{{\rm \mbox{--}}}p}_{n\kk}(T)\right|_{\rm OMS}=\Gamma^{\rm e{{\rm \mbox{--}}}p}_{n\kk}(\epsilon_{n \kk},T)$. Like in the electronic case the most important numerical property of the lifetimes calculation is that they depend only on the $\newcommand{\qq}{{\bf q}} \qq$-grid.

It is very instructive to compare $\newcommand{\kk}{{\bf k}} \left.\Gamma^{\rm e{{\rm \mbox{--}}}p}_{n\kk}(T)\right|_{\rm OMS}$ and the $\newcommand{\kk}{{\bf k}} \left.\Gamma^{\rm e{{\rm \mbox{--}}}e}_{n\kk}(T)\right|_{\rm OMS}$ for a paradigmatic material like bulk silicon. This is done in figure 12 in the zero temperature limit. The very different nature of the two lifetimes appear clearly. By simple energy conservation arguments, the electronic linewidth are zero by definition in the two energy regions $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\e}{{\rm e}} \epsilon_{\rm VBM}-E_{\rm g}$ and $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\e}{{\rm e}} \epsilon_{\rm CBM}+E_{\rm g}$, with $E_{\rm g}$ the electronic gap (in silicon $E_{\rm g}\approx1.1$ eV) and $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\e}{{\rm e}} \epsilon_{\rm CBM}$ $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\e}{{\rm e}} (\epsilon_{\rm VBM})$ the conduction band minimum (valence band maximum). In these energy regions the e–p contribution is stronger and the corresponding linewidths are larger then the e–e ones. The quadratic energy dependence of the e–e linewidths inverts this trend at higher energies.

Zoom In Zoom Out Reset image size

While the e–e contributions grow quadratic in the energy dependence, the e–p ones follow the electronic density of states profile. This property is confirmed by figure 13(b) (see also [117]) and remains accurate when the temperature increases. Figure 13(a) shows instead the EP correction in single-layer MoS₂ of the valence and conduction band states for several temperatures, together with the widths and the density of states (DOS). In general, the EP correction tends to close the bandgap. This is visible in figure 13 (panels (a) and (b)), the conduction state energy decreases with temperature, while the valence one increases. Only in a few cases do we find an opening of the bandgap when temperature increases [118].

Zoom In Zoom Out Reset image size

Once the temperature-dependent corrections to electron and hole states have been calculated, they constitute the key ingredients of the finite temperature excitonic eigenvalue equation. Since the electron and hole eigenvalues are complex numbers the resulting excitonic eigenvalues have a real part (the exciton binding energy) and an imaginary part (the exciton lifetime). The dielectric function then depends explicitly on T, $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\la}{\langle} \newcommand{\e}{{\rm e}} \epsilon_2(\omega,T)= -(8\pi/V) \sum_\lambda |S_\lambda(T)|^2\Im{[\omega-E_\lambda(T)]^{-1} }$, where $\newcommand{\la}{\langle} S_{\lambda}(T)$ are the excitonic optical strengths and $\newcommand{\la}{\langle} E_\lambda(T)$ are the complex excitonic energies. As shown in figure 14 for a single layer MoS₂, the main effect of the temperature on the optical spectra is the renormalization of the energy transitions along with a broadening of the spectrum lineshape related to the finite lifetime of the underlying excitonic states which increases with T [119, 120]. This picture is also valid when $T\rightarrow 0$ because of the zero-point vibrations. A remarkable effect of the exciton-phonon coupling has been observed in hexagonal BN. It has been proven that the optical brightness turns out to be strongly temperature-dependent such as to induce bright to dark (and viceversa) transitions [121].

Zoom In Zoom Out Reset image size

A technically relevant issue is the slowing down of energy correction convergence at some high-symmetry points. Some EP matrix elements might be zero by symmetry and are not representative of the discretization of an integral. yambo deals with this issue by computing the energy shift corrections on a random $\bf q$-wavevector grid of transferred momenta. The numerical evaluation of the EP self-energy on a dense $\bf q$-grid is a formidable task (see equation (5) of [110]). The reason is that such dense grids of transferred momenta are inevitably connected with the use of equally dense grids of $\bf k$-points. The solution implemented in yambo is a double grid approach: matrix elements are calculated for a fixed $\bf k$-point grid while energies are integrated using a larger grid of random-points in the whole BZ.

To speed-up the convergence with the number of random points, the BZ is divided in small spherical regions $\newcommand{\qq}{{\bf q}} R_{\qq}$ centered around each $\newcommand{\qq}{{\bf q}} \qq$ point of the regular grid and the integral is calculated using a numerical Monte-Carlo integration technique. Furthermore, the divergence at $\newcommand{\qq}{{\bf q}} \qq \rightarrow 0$ of the $\newcommand{\la}{\langle} \newcommand{\kk}{{\bf k}} \newcommand{\qq}{{\bf q}} |g^{\qq \lambda}_{n'n\kk}|^2$ matrix elements is explicitly taken into account for the 3D case for which the $\newcommand{\qq}{{\bf q}} \qq$ integration compensates the $\newcommand{\qq}{{\bf q}} \qq^{-2}$ divergence. In the case of 2D materials the divergence of the EP matrix elements would not be lifted by the surface element $2\pi q$. In principle, an analytic functional form for the EP matrix elements can be envisaged as reported by [123].

The EOM for the density matrix projected in the space of the single particle wave-functions, $\newcommand{\rhom}{\underline{\rho}} \rhom$, is derived from non-equilibrium (NEQ) many-body perturbation theory and reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Um}{\underline{U}} \newcommand{\hm}{\underline{h}} \newcommand{\rhom}{\underline{\rho}} \newcommand{\kk}{{\bf k}} \newcommand{\la}{\langle} \displaystyle \label{eq:td_rho} {\rm i}\hbar \partial_t\, \rhom_\kk(t) = \left[ \hm^{\rm rt}_\kk[\rhom]+\Um^{\rm ext}_\kk(t), \rhom_\kk(t) \right] -{\rm i} \Gamma_\kk \rhom_\kk . \nonumber \end{align} \tag{ 16 }$$

Here we underline quantities which are vectors in the transition space (and we will underline twice matrices in transition space). $h^{\rm rt}$ contains the equilibrium eigenvalues $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\kk}{{\bf k}} \newcommand{\e}{{\rm e}} \epsilon_{n\kk}$ plus the variation of the self-energy $\newcommand{\rhom}{\underline{\rho}} \newcommand{\Sm}{\underline{\Sigma}} \Delta\Sm^{\rm Hxc}[\rhom]$; $\newcommand{\ege}{\varepsilon} \newcommand{\ep}{\ege_{p}} \newcommand{\kk}{{\bf k}} \newcommand{\e}{{\rm e}} \epsilon_{n\kk}$ can be the KS energies or the QP corrected energies. QP corrections can be loaded either from a previous calculation or by adding a scissor operator from input. For $\newcommand{\rhom}{\underline{\rho}} \newcommand{\Sm}{\underline{\Sigma}} \Delta\Sm^{\rm Hxc}[\rhom]$ different levels of approximation can be chosen, setting the HXC_kind input variable. $\newcommand{\efield}{\boldsymbol{\mathcal{E}}} \newcommand{\Um}{\underline{U}} \newcommand{\rmat}{\underline{r}} \newcommand{\e}{{\rm e}} \Um^{\rm ext}=-e\efield\cdot {\bf \rmat}$ represents the external potential written in the length gauge; shape, polarization, intensity (and eventually frequency) of the field $\newcommand{\efield}{\boldsymbol{\mathcal{E}}} \newcommand{\e}{{\rm e}} \efield$ can be selected in input. $\newcommand{\rmat}{\underline{r}} {\bf \rmat}$ is the position operator. The coupling to the external field is exact up to first order. From the knowledge of the density matrix, the first order polarization $\newcommand{\PP}{{\bf P}} \PP(t)=-e\sum_{i\neq j{\bf k}} {\bf r}_{ij{\bf k}} \rho_{ij{\bf k}}$ is computed at each time step. The spectrum of the system can then be obtained by the Fourier transform of the polarization, which can be done as a post-processing step. Absorption is thus obtained via the dipole–dipole response function (equivalent to the length gauge in linear response). A delta-like external field is convenient to obtain the spectrum for all frequencies.

The implementation of the external field in the velocity gauge (equivalent to the velocity gauge in linear response) is currently under testing before its final release. Equation (16) represents a set of equations, one for each $\newcommand{\kk}{{\bf k}} \kk$-point in the BZ, coupled via the functional dependency of $\newcommand{\Sm}{\underline{\Sigma}} \Delta\Sm^{\rm Hxc}$ on the whole $\newcommand{\rhom}{\underline{\rho}} \rhom$. Different options of the self-energy are available, by setting the HXC_kind variable to: IP, Hartree, DFT, Fock, Hartree  +  Fock, SEX, or Hartree  +  SEX. For HXC_kind  =  IP one has $\newcommand{\Sm}{\underline{\Sigma}} \Delta\Sm^{\rm Hxc}=0$. For local HXC_kind like Hartree and DFT, $\newcommand{\Sm}{\underline{\Sigma}} \Delta\Sm^{\rm Hxc}$ can be computed on the fly from the real-space density $\newcommand{\rr}{{\bf r}} n(\rr,t)$ and the approach is in practice equivalent to TD-DFT to linear order in the field. For non-local HXC_kind like Hartree  +  Fock (HF) and Hartree  +  SEX (HSEX), the self-energy is written in the form $\newcommand{\rhom}{\underline{\rho}} \newcommand{\Sm}{\underline{\Sigma}} \newcommand{\Kmm}{\underline{\underline{K}}} \Delta\Sm^{\rm Hxc}[\rhom]=\Kmm^{\rm Hxc}\cdot\rhom$ with $\newcommand{\Kmm}{\underline{\underline{K}}} \Kmm^{\rm Hxc}$ computed before starting the real-time propagation. The calculation of $\newcommand{\Kmm}{\underline{\underline{K}}} \Kmm^{\rm Hxc}$ can be either done in a preliminary run, with the matrix-elements stored on disk and then reloaded, or on-the-fly before starting the real-time propagation. In case the HSEX approximation is used the resulting spectrum is equivalent to a BSE calculation in the limit of small perturbations, as shown both analytically and numerically in [126]. Thus, to linear order, TD-SEX is is able to properly capture excitons, which can be hardly described within TD-DFT. The comparison between the two approaches is reported in figure 15.

Zoom In Zoom Out Reset image size

When local self-energies are computed directly from the real space density, the numerical cost is mainly due to the projection of the potential on the KS-basis set at each iteration. This step is avoided in real-space TD-DFT where, however, the wave-functions on the real-space grid need to be propagated. The relative computation cost of the two strategies depends critically on the size of the real-space grid versus the number of KS function used. Instead when non local self-energy are used, the computational cost is mostly due to the preliminary calculation of the kernel $\newcommand{\Kmm}{\underline{\underline{K}}} \Kmm^{\rm Hxc}$. This step has, roughly, the same computational cost of a standard BSE run and requires to store $\newcommand{\Kmm}{\underline{\underline{K}}} \Kmm^{\rm Hxc}$ in memory (disk or RAM). The subsequent real-time propagation is instead very fast. In some cases it is convenient to use this scheme also for local self-energies. However, regardless of the self-energy used, only variations of the self-energy which are linear in the density matrix are described when using $\newcommand{\Kmm}{\underline{\underline{K}}} \Kmm^{\rm Hxc}$.

To run simulations and compute the spectra as described in the present section the yambo_rt and ypp_rt executables need to be used.

7.1.1. Double-grid in real time.

Alternative to the time-evolution of the density matrix, it is possible to perform the time-evolution of the Schrödinger equation for the periodic part of the Bloch states projected in the eigenstates of the equilibrium Hamiltonian: $\newcommand{\ra}{\rangle} |v_{m{\bf k}} \rangle$. Here we briefly present the actual implementation in yambo and how it can be used to obtain non-linear optics response, for more details see [129]. The EOM for the valence band states reads:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\efield}{\boldsymbol{\mathcal{E}}} \newcommand{\kk}{{\bf k}} \newcommand{\ra}{\rangle} \newcommand{\la}{\langle} \displaystyle {\rm i}\hbar \partial_t\, | v_{m\kk} \rangle = \left( h^{{\rm rt}}_{\kk}[\rho] +{\rm i} \efield \cdot \tilde \partial_\kk\right) |v_{m\kk} \rangle \label{eq:td_wfs} \; \nonumber \end{align} \tag{ 17 }$$

where the effective Hamiltonian $\newcommand{\kk}{{\bf k}} h^{{\rm rt}}_\kk$ has been introduced in section 7.1 and $\rho(t)$ is constructed starting from $\newcommand{\ra}{\rangle} \newcommand{\kk}{{\bf k}} | v_{m\kk} \rangle$. The second term in equation (17), $\newcommand{\efield}{\boldsymbol{\mathcal{E}}} \newcommand{\e}{{\rm e}} \efield \cdot \tilde \partial_{\bf k}$, describes the coupling with the external field $\newcommand{\efield}{\boldsymbol{\mathcal{E}}} \newcommand{\e}{{\rm e}} \efield$ in the dipole approximation. As we imposed Born–von Kármán periodic boundary conditions, the coupling takes the form of a ${\bf k}$-derivative operator $\tilde \partial_{\bf k}$. The tilde indicates that the operator is 'gauge covariant' and guarantees that the solutions of equation (17) are invariant under unitary rotations among occupied states at ${\bf k}$ (see [130] for more details).

Propagating the single particle wave-functions presents advantages and disadvantages with respect to the density matrix. The major advantage is that the coupling of electrons with the external field, within the length gauge, is now written in terms of Berry's phase, which is exact to all orders also in extended systems [131]. Moreover, from the evolution of $\newcommand{\ra}{\rangle} | v_{m{\bf k}} \rangle$ in equation (17) also the time-dependent polarization [132] $\newcommand{\PP}{{\bf P}} \PP_\parallel$ along the lattice vector $\mathbf a$ can be computed in terms of the Berry phase:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\PP}{{\bf P}} \newcommand{\la}{\langle} \displaystyle \PP_\parallel = -\frac{ef |\mathbf a| }{2 \pi \Omega_c} {{\rm Im~log}} \prod_{{\bf k}}^{N_{{\bf k}}-1}\ {{\rm det}} S\left({\bf k} , {\bf k} + \mathbf q\right), \label{berryP} \nonumber \end{align} \tag{ 18 }$$

where $\newcommand{\la}{\langle} \newcommand{\ra}{\rangle} S({\bf k} , {\bf k} + \mathbf q) = \langle v_{n{\bf k}}|v_{m{\bf k} + {\bf q}}\rangle$ is the overlap matrix between the valence states, $\Omega_c$ is the unit cell volume, f  is the spin degeneracy, $N_{{\bf k}}$ is the number of ${\bf k}$ points along the polarization direction, and $\mathbf q = 2\pi/(N_{{\bf k}} {\mathbf a})$. The resulting polarization can be expanded in a power series of the field $\newcommand{\efield}{\boldsymbol{\mathcal{E}}} \newcommand{\e}{{\rm e}} \efield_j$ as:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\efield}{\boldsymbol{\mathcal{E}}} \displaystyle {\bf P}_i = \chi^{(1)}_{ij} \efield_j + \chi^{(2)}_{ijk} \efield_j \efield_k + \chi^{(3)}_{ijkl} \efield_j \efield_k \efield_l + O(\efield^4) \; , \nonumber \end{align} \tag{ 19 }$$

where the coefficients $\chi^{(i)}$ are functions of the frequency of the perturbing fields and of the outgoing polarization. From the Fourier analysis of the ${\bf P}_i$ it is possible to extract all the non-linear coefficients (see [129] for more details). As in section 7.1, the level of approximation of the so-calculated susceptibilities depends on the effective Hamiltonian that appears in the right hand side of equation (17). Different choices are possible, namely, the independent particle approximation (IPA), the time-dependent Hartree, the real-time Bethe–Salpeter equation (RT-BSE) framework, or TD-DFT. This approach has been successfully applied to study second-, third- harmonic generation and two-photon absorption in bulk materials and nanostructures [129, 133, 134]. As before, in the limit of small perturbation equation (17) reproduces the optical absorption calculated with the standard GW  +  BSE approach [135].

Since the exact polarization is available, the approach based on equation (17) not only reduces to TD-DFT when local functionals are considered, it also includes TD density polarization functional theory (DPFT) as a special case. Thus specific approximations for both the microscopic and the macroscopic part of $\Delta\Sigma^{\rm Hxc}$ are available, which, within TD-DPFT, are expressed as functionals of $\rho$ for the microscopic part ($v^{\rm Hxc}[\rho]$) and of $\newcommand{\PP}{{\bf P}} \PP$ for the macroscopic part ($\newcommand{\efield}{\boldsymbol{\mathcal{E}}} \newcommand{\e}{{\rm e}} \efield^{\rm Hxc}$) as discussed in [136, 137].

A comparison between TD-DPFT in the real-time framework and the solution of the Bethe–Salpeter equation for different zinc-blende compounds has been recently published by Riefer et al [138]. While for linear response the different functionals give a satisfactory result [137], for the second harmonic generation the situation is less clear. This is probably due to the fact all exchange-correlation kernels implemented in yambo and tested in the previous papers were derived in the linear response regime.

In non-linear optics simulations the system is excited with a laser at given frequency $\omega$ and dephasing term $\newcommand{\la}{\langle} \lambda_{\rm deph}$ is added to the Hamiltonian. After a time $\newcommand{\la}{\langle} T \gg \lambda_{\rm deph}$, sufficient to damp out the eigenmodes of the system, the signal is analyzed to extract the non-linear response functions, see figure 16 and [129]. To run simulations and compute the spectra as described in the present section the yambo_nl and ypp_nl executables need to be used.

Zoom In Zoom Out Reset image size

The multilevel MPI structure of yambo is reflected in the input file where, for each computational kernel (runlevel here) there are two related input variables: the first one, runlevel_ROLEs, sets on which parameters the user wants to distribute the MPI workload, while the second runlevel_CPU defines how many MPI tasks will be associated to such parameters. As an example

X_finite_q_ROLEs="q.k.c.v.g"

X_finite_q_CPU="2.3.5.2.1"

is a possible input for running on $2\times 3\times 5\times 2 \times 1= 60$ MPI tasks, with the $\newcommand{\qq}{{\bf q}} \qq$-points distributed on two tasks, the $\newcommand{\kk}{{\bf k}} \kk$-points on three, the conduction and valence bands on five and two respectively. One more level of parallelism (g) is present, acting and distributing the response matrix over space degrees of freedom (plane waves). The order of the parameters in the runlevel_ROLEs variables is irrelevant. On top of that, more input variables are available to handle parallel linear algebra (e.g. via scalapack and blacs libs) and to select the number of OpenMP threads (on a runlevel basis if needed). Such hierarchical organization makes it possible to have MPI communication only within the subgroups, thus avoiding, whenever possible, to deal with all2all communications.

If the user does not wish to deal with the complexity of such multi-level parallelization a default layout is provided. However, the fine-tuning of the MPI/OpenMP related variables can (further) reduce the load imbalance, improve the memory distribution or decrease the total time-to-solution. For this reason, in the sections 8.3–8.6 specific suggestions for best parallel exploitation of each runlevel are provided.

yambo stores binary data using the netcdf library. Depending on the configuration flags, data can be stored in classic netcdf format (file size limit of 2 GB, activated with –enable-netcdf-classic), 64-bit netcdf format (no file size limit, default) or HDF5 format (requires at least netcdf v4 linked to HDF5, activated with –enable-netcdf-hdf5). Since version 4.4, in case the HDF5 format is specified, parallel I/O can also be activated (–enable-hdf5-par-io) to store the response function in $\newcommand{\GG}{{\bf G}} \GG$ space or the kernel of the BSE. For the $\newcommand{\GG}{{\bf G}} \GG$-space response function, parallel I/O avoids extra communication among the MPI tasks and also reduces the amount of memory allocated per task. For the BSE case, parallel I/O makes it possible to load the kernel computed from a previous calculation using a different parallel scheme and/or a different number of MPI tasks. Indeed the calculation of the kernel matrix elements is very time consuming but has a very efficient memory and load distribution. In contrast, the solution of the BSE eigenproblem is less time consuming but also less efficiently distributed. It is thus suggested to first compute the kernel matrix on a large number of cores and then to solve the BSE on fewer tasks as a follow-up step.

According to equation (1), the computation of the response function in G-space can be distributed over five different levels: $\newcommand{\qq}{{\bf q}} \qq$-points, $\newcommand{\kk}{{\bf k}} \kk$-points, conduction and valence bands ($c,v$), and $\newcommand{\GG}{{\bf G}} \GG$-vectors. The distribution over the $\newcommand{\qq}{{\bf q}} \qq$-points would be the most natural choice, since the response functions at different $\newcommand{\qq}{{\bf q}} \qq$ are completely independent. However, it may lead to significant memory duplication (multiple sets of wavefunctions are managed at the same time) and load imbalance since the number of possible transitions varies from point to point. It is usually not recommended unless a large number of $\newcommand{\qq}{{\bf q}} \qq$-points has to be considered. Instead, it is usually more effective to parallelize over $\newcommand{\kk}{{\bf k}} \kk$, and bands (c,$v$) indexes. This requires slightly more MPI communication (due to a MPI reduction at the end of the calculation) but is very efficient in terms of speedup (almost linear) and in terms of memory distribution (especially for $c,v$, since usually wavefunctions are the leading memory contribution). While transitions are evenly distributed (balanced workload), yambo sorts and groups transitions which are (almost) degenerate in energy (see appendix B) thus, in practice, small imbalances can still be present.

Further, the distribution over the $\newcommand{\GG}{{\bf G}} \GG$-vectors is the most internal one, and requires more communication among MPI tasks (unless parallel I/O is activated). It can be useful for systems with a very large number of $\newcommand{\GG}{{\bf G}} \GG$-vectors (such as low dimensional systems or surfaces) to distribute the response function and ease memory usage. Finally, the computation of $\chi^0$ can also benefit from OpenMP parallelism. The distribution over threads has been implemented at the same level of the MPI parallelism (i.e. over transitions), resulting in a very good scaling while reducing the memory usage per core. Note however that some memory duplication (a $\newcommand{\GG}{{\bf G}} M(\GG,\GG')$ workspace matrix per thread) has to be paid to make the implementation more efficient. The OpenMP parallelism of $\chi^0$ (including dipoles) is governed by the input variables:

X_Threads  =  8

DIP_Threads  =  8

both defaults being set to 0, i.e. controlled as usual by the OMP_NUM_THREADS environment variable. Once the independent-particle response function $\newcommand{\qq}{{\bf q}} \newcommand{\GG}{{\bf G}} {\chi^0_{\GG\GG'}(\qq,\omega)}$ has been computed, a Dyson equation is solved for each frequency to construct the RPA response function. This can be done either by distributing over different frequencies or by using parallel linear algebra (see section 8.6).

Following equation (7), the HF and GW correlation self-energies can be parallelized with MPI over three different layers: $\newcommand{\qq}{{\bf q}} \qq$-points (q); bands in the Green's function (b) [see m in equation (15)]; and quasiparticle corrections $\newcommand{\kk}{{\bf k}} \Sigma_{n n' \kk}$ to be computed (qp). OpenMP parallelism here acts at the lowest level, dealing with sums over $\newcommand{\GG}{{\bf G}} \GG$ and $\newcommand{\GG}{{\bf G}} \GG'$, i.e. spatial degrees of freedom. The following variables can be modified to fine-tune the self-energy parallelization (here shown for 60 MPI tasks and eight OpenMP threads):

SE_ROLEs="q.qp.b"

SE_CPU="1.4.15"

SE_Threads  =  8

Since the sum over $\newcommand{\qq}{{\bf q}} \qq$-points in equation (7) is over the whole BZ, the $\newcommand{\qq}{{\bf q}} \qq$-parallelism for the self-energy may be even more unbalanced than that for the response function (here every $\newcommand{\qq}{{\bf q}} \qq$-point needs to be expanded to account by symmetry for its whole star) and is recommended only when a large number of $\newcommand{\qq}{{\bf q}} \qq$-points is available. Instead, the parallelism over bands b tends to distribute evenly both memory and computation, at the price of a mild MPI communication, thereby resulting a natural choice (when enough bands are included in the calculation). qp-parallelism distributes the computation but tends to replicate memory (wavefunctions are not further distributed). In general, the OpenMP parallelism is extremely efficient for the GW self-energy without having to pay for any extra memory workspace.

In the solution of the BSE most of the CPU time is spent in building up the excitonic matrix, or more precisely, its kernel. The input flags which control the parallel distribution of the workload needed to build the kernel are eh.k.t. To distribute the workload, first all possible transitions $\newcommand{\kk}{{\bf k}} cv\kk$, i.e. from valence band $v$ to conduction band c at the k-point $\newcommand{\kk}{{\bf k}} \kk$, are split into transition groups (TGs). Then for each pair of TGs a block of the BSE matrix is created $B_{ij}=\{T_i\rightarrow T_j\}$. Defined N_t the total number of TG, then the BSE matrix will be divided into $N_b=N_t^2$ blocks. In the Hermitian case (as in the Tamm–Dancoff approximation), only $N_b=N_t(N_t+1)/2$ blocks will be computed. The parallelization flags for the BSE define both N_t and N_b, and how the resulting blocks are distributed among the MPI tasks. Indeed $N_t=n_{\rm eh}n_k^{\rm ibz}$ where $n_{\rm eh}$ is the number of MPI tasks assigned to the eh field and $n_k^{\rm ibz}$ is the number of $\newcommand{\kk}{{\bf k}} \kk$-points in the IBZ. This means that even setting eh  =  1 a minimum number of $\newcommand{\kk}{{\bf k}} \kk$-based TGs (k-TGs) is always created, which is eventually split into subgroups when $n_{\rm eh}>1$. It is important to note that k-TGs are defined using the k-sampling in the the IBZ, while the BSE matrix is defined in the whole BZ, resulting in groups of non-uniform size. However, the symmetry operations relating matrix elements within a given k-TGs are taken into account by yambo. As a consequence, in systems where $n_k^{\rm ibz}\neq n_k^{\rm bz}$ the use of $n_{\rm eh}>1$ is discouraged, as the splitting of k-TGs over different MPI-tasks implies that symmetry-related matrix-elements can be assigned to different MPI-tasks and need to be recomputed. Once N_t and hence N_b are defined, transitions and blocks are distributed among the MPI tasks as explained in the following example.

Suppose we have a system with 18 $\newcommand{\kk}{{\bf k}} \kk$-points in the IBZ, and we adopt the parallelization strategy 2.3.3 for eh.k.t in the case the BSE is Hermitian. Then $N_t=2\times 18=36$ and N_b  =  666. Thus, in our example we are using in total $2\times 3\times 3$ MPI-tasks. The eh.k fields define $2\times 3=6$ MPI-groups which split the 36 transition-groups. Thus, each MPI-group has to deal with six transition-groups. For each transition group T_n, there are N_t blocks $B^n_{ij}=T_i\rightarrow T_j$ for the Hermitian case, where the T_n appears as initial ($T_i=T_n$) or final ($T_j=T_n$) state. Most of the blocks belongs to two transition-groups (only the blocks B_ii belong to one transition-group). This means that each MPI-group builds half of the B_ij ($6 \times 35/2$) plus all B_ii (6) blocks. These 111 blocks are divided according to the t field and thus each MPI-task will be assigned to 37 blocks.

Dense linear algebra is extensively used in yambo. Among the most time-consuming tasks we have identified the inclusion of local field effects [1] in the RPA response function

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\qq}{{\bf q}} \newcommand{\la}{\langle} \displaystyle \chi_{{\bf{G}}{\bf{G}'}}^{RPA}(\qq,\omega)&= \chi^{0}_{{\bf{G}}{\bf{G}'}} (\qq,\omega) \nonumber \\&+ \chi^{0}_{{\bf{G}}{\overline{\bf{G}}}}(\qq,\omega) \frac{4\pi}{|\qq+\overline{\bf{G}}|^2} \chi_{\overline{\bf{G}}{\bf{G}'}}^{RPA}(\qq,\omega). \label{eq:X_rpa} \nonumber \end{align} \tag{ 20 }$$

The solution of equation (20) can be cast in the form of a matrix inversion. Indeed:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\GG}{{\bf G}} \newcommand{\qq}{{\bf q}} \newcommand{\la}{\langle} \displaystyle \nonumber \chi_{{\bf{G}}{\bf{G}'}}^{RPA}(\qq,\omega)&= \left[ \delta_{\GG,\overline{\GG}} - \chi^{0}_{{\bf{G}}{\overline{\bf{G}}}}(\qq,\omega) \frac{4\pi}{|\qq+\overline{\bf{G}}|^2}\right]^{-1}\nonumber \\ &\qquad \qquad \qquad \chi_{\overline{\bf{G}}{\bf{G}'}}^{0}(\qq,\omega). \label{eq:X_rpa_inv} \nonumber \end{align} \tag{ 21 }$$

Equation (21), and the solution of the BSE (diagonalization), which can be considered prototype kernels.

Once a finite basis set is adopted, the operators involved are represented as (dense) $N\times N$ matrices, with N easily reaching few-to-tens of thousands or more, making standard linear algebra tasks (such as matrix multiplication, inversion, diagonalization) quite intense. We have therefore implemented dense parallel linear algebra by exploiting the ScaLAPACK library [88] within the MPI parallel structure of yambo. Concerning the RPA response, this means that on top of the MPI parallelism over $\mathbf{q}$-vectors, multiple instances of parallel linear algebra are run at the same time (one per $\mathbf{q}$ vector) to compute $\chi^{{\rm RPA}}$.

The behavior of the yambo parallel linear algebra is governed by the variables:

runlevel_nCPU_LinAlg_INV  =  64

runlevel_nCPU_LinAlg_DIAGO  =  64

where (runlevel could be, for example, the RPA response function or the BSE). Given the relevance, the calculation of the IP response function $\newcommand{\qq}{{\bf q}} \newcommand{\GG}{{\bf G}} {\chi^0_{\GG\GG'}(\qq,\omega)}$ has also been block-distributed over $\newcommand{\GG}{{\bf G}} \GG,\GG'$ vectors (g-parallelism in section 8.3), both in terms of computation and memory-usage.

When using the SLEPC diagonalization method to obtain the BSE spectra, the memory distribution of the eigensolver (not to be confused with the memory distribution of the BSE matrix discussed in section 5.1.2) is handled by the SLEPC library itself. For more details, the reader is referred to the SLEPC specific literature [86].

Yambopy is a community project to develop Python classes and scripts to express, automate, and analyze calculations with the yambo code. A typical yambo workflow involves a few steps: generating the KS states with a DFT code, preparing the yambo databases and then running the yambo code. These workflows become more complicated when performing convergence tests or when repeating the same calculations for multiple materials. With yambopy the user can express these yambo workflows on a python script that can then be shared and reproduced among different users. We currently provide python classes to read and write the input files as well as the output of the yambo code. A lightweight python interface for the Quantum ESPRESSO Suite is also provided in the qepy package distributed along with yambopy. For more comprehensive python interfaces for the DFT codes see ASE [142], Abipy [143] for the Abinit code, or AiiDA [141].

The YamboIn class provided by yambopy is used to read and write the base yambo input file generated by yambo and modify it in a programmatic way. With this tool it is possible, for example, to create a set of input files by changing a single variable inside a for loop. We also provide classes to read the yambo NetCDF databases in Python (for a complete list see the on-line documentation [144]). These classes provide methods to manipulate and represent the data using the matplotlib [145] library giving great flexibility for the interpretation and analysis of the results. Running these workflows on a HPC context requires to write job submission scripts for different job schedulers, this can be done using the schedulerpy package also accompanying yambopy.

For quick access to some features from the command line, we provide the yambopy shell command. The script is automatically installed with yambopy in such a way that some functionalities of yambopy are directly callable from the command line. This script has features to plot the convergence tests for GW and BSE calculations, excitonic wave functions, and dielectric functions among others.

To ensure software quality and usability we provide yambopy as an open-source code along with documentation and automatic testing. A detailed documentation of the classes, features, and a tutorial are available on the yambopy website [144]. We keep a public git repository hosted on Github where the users can download the latest version of the code as well as contribute with patches, features, and workflows. Sharing the workflows among users allow us to avoid repeated technical work and greatly simplifies the use of the yambo code. Continuous integration tests are done using the Travis-CI platform [146] leading the code to be tested at each commit and thus enforcing its reliability. yambopy is a project under development, and will be described on its own in a future work.

The AiiDA platform is a materials' informatics infrastructure which implements the so-called ADES model (Automation, Data, Environment and Sharing) for computational science [141]. The AiiDA plugins and workflows for yambo are publicly available on GitHub [147], while online documentation and tutorials are available on Read the Docs [148].

9.2.1. The yambo-AiiDA plugin.

9.2.2. AiiDA workflows: automated GW.

A new important tool introduced to improve and stabilize the development of the yambo code is the test-suite. The yambo test-suite is stored in a dedicated repository (yambo-tests) on GitHub and contains a series of tests which can be run in an automated manner. The repository is freely accessible after registering as a 'yambo user' on the GitHub. While the test-suite is mainly aimed at developers, users can also benefit from accessing its input and reference files and automatically checking if their compiled version of yambo works properly.

The test-suite is governed by using a Perl script, driver.pl. This script uses internal Perl modules to perform several tasks: it automatically compiles the yambo code (a precompiled version can also be used), it runs the code and checks the output against reference files stored in the test-suite repository.

The code can be run in serial, parallel with OpenMP threads and checking parallel I/O and/or parallel linear algebra. At least two different groups of tests are available: smaller (and faster) tests which are run on a daily basis and longer tests which are used for a deeper testing of the code before a release.

The same driver can also be used to run yambo benchmarks. Benchmarks tests are a particular group of materials that, describing complex nano-structures (a 1D polymer or carbon-based ribbon) or a water cell, require a large number of reciprocal space vectors and/or $\newcommand{\kk}{{\bf k}} \kk$-points. As a consequence these systems are suitable to be executed using a large number of cores on parallel machines. In this case the test-suite can collect the results and loop on different parallel configurations testing their performances. More importantly the test-suite organizes the results in machine dependent folders that can be, eventually, post-analyzed.

The results of the night runs of the test-suite are publicly available on the web-page [150] and can be inspected without having access to the machines that run the tests. This is very useful in order for any development to reproduce a specific error to be fixes.