SIMPLE code: optical properties with optimal basis functions

TL;DR

SIMPLE is a computational code that efficiently calculates optical properties of materials using an optimal basis method, reducing computational costs while maintaining accuracy for metals and insulators.

Contribution

The paper introduces SIMPLE, a new code implementing the optimal basis method for optical property calculations, with systematic improvability and ease of use.

Findings

Accurate optical properties obtained with reduced computational cost.

Effective for both metallic and insulating systems.

Validated on bulk silver and silicon.

Abstract

We present SIMPLE, a code developed to calculate optical properties of metallic and insulating extended systems using the optimal basis method originally proposed by E. L. Shirley in 1996. Two different approaches for the evaluation of the complex dielectric function are implemented: (i) the independent-particle approximation considering both interband and intraband contributions for metals and (ii) the Bethe-Salpeter equation for insulators. Since, notably, the optimal basis set is systematically improvable, accurate results can be obtained at a strongly reduced computational cost. Moreover, the simplicity of the method allows for a straightforward use of the code; improvement of the optimal basis and thus the overall precision of the simulations is simply controlled by one (for metals) or two (for insulators) input thresholds. The code is extensively tested, in terms of verification and performance, on bulk silver for metals and on bulk silicon for insulators.

Tables

Table 1. Table 1: Input keywords for simple.x . The input file must contain the string &inputsimple before inserting any input flag.

Input flag	Description
prefix	Same as in pw.x
outdir	Same as in pw.x
calc_mode	Calculation mode (0 for BSE, 1 for IP)
s_bands	Threshold for the construction of the OB [bohr3]
num_nbndv	Occupied bands
num_val	Valence bands (starting from HOMO)
num_cond	Conduction bands (starting from LUMO)
nkpoints(1)	k-grid ($n_1$)
nkpoints(2)	k-grid ($n_2$)
nkpoints(3)	k-grid ($n_3$)
s_product	Threshold for the construction of the OPB [bohr3] (BSE)
w_type	Screened Coulomb interaction from GWL calculation (0) or
	from model microscopic dielectric function (1) (BSE)
nonlocal_commutator	Flag to calculate the non-local commutator (IP)
epsm	Macroscopic dielectric constant ${ϵ}_m$ from model microscopic
	dielectric function (if w_type=1) (BSE)
lambdam	Screening length $\overline{\lambda }$ from model microscopic
	dielectric function (if w_type=1) [bohr-1] (BSE)
numpw	Dimension of polarizability basis (if w_type=0) (BSE)
l_truncated_coulomb	Flag to activate truncation of Coulomb
	interaction (BSE)
truncation_radius	Radius at which Coloumb interaction
	is truncated [bohr] (BSE)

Table 2. Table 2: Input keywords for simple_ip.x . The input file must contain the string &inputsimpleip before inserting any input flag. The prefix of each input flag is simpleip_in% as reported in the table.

Input flag	Description
simpleip_in%prefix	Same as in pw.x
simpleip_in%outdir	Same as in pw.x
simpleip_in%interp_grid(1)	Interpolation k-grid ($m_1$)
simpleip_in%interp_grid(2)	Interpolation k-grid ($m_2$)
simpleip_in%interp_grid(3)	Interpolation k-grid ($m_3$)
simpleip_in%fermi_ngauss	Broadening type for Fermi energy (integer $n$)
	$n=-99$ Fermi-Dirac, $n=-1$ Marzari-Vanderbilt,
	$n\ge 0$ Methfessel-Paxton
simpleip_in%fermi_degauss	Broadening for Fermi energy [Ry]
simpleip_in%drude_degauss	Broadening for Drude plasma frequency [Ry]
simpleip_in%inter_broadening	Empirical broadening for ${\epsilon }_{\text{IP}}^{\text{inter}}(\widehat{𝐪},\omega )$ [Ry]
simpleip_in%intra_broadening	Empirical broadening for ${\epsilon }_{\text{IP}}^{\text{intra}}(\widehat{𝐪},\omega )$ [Ry]
simpleip_in%wmin	Minimum frequency for ${\epsilon }_{\text{IP}}(\widehat{𝐪},\omega )$ [Ry]
simpleip_in%wmax	Maximum frequency for ${\epsilon }_{\text{IP}}(\widehat{𝐪},\omega )$ [Ry]
simpleip_in%nw	Number of frequencies in the range [wmin, wmax]
simpleip_in%nonlocal_commutator	Flag to include the non-local commutator
simpleip_in%nonlocal_interpolation	Flag to use linear interpolation of non-local
	contributions (experimental feature)

Table 3. Table 3: Input keywords for simple_bse.x . The input file must contain the string &inputsimple before inserting any input flag. The prefix of each input flag is simple_in% as reported in the table.

Input flag	Description
simple_in%prefix	Same as in pw.x
simple_in%outdir	Same as in pw.x
simple_in%scissor	Rigid scissor energy for opening DFT gap [eV]
simple_in%spin_state	In case of collinear spin calculation selects triplet (0) or
	or singlet (1) excitons , (2) for the case of non-collinear spin
simple_in%h_level	Terms in excitonic Hamiltonian: (0) only diagonal, (1) as (0)
	plus exchange term, (2) as (1) plus direct term of bare
	Coulomb part, (3) as (2) plus direct term of correlation part
simple_in%task	If set to 0 find the lowest excitons, if set to 1 find the
	${\alpha }_i$ and ${\beta }_i$ coefficients of the Haydock recursive method
simple_in%nvec	Number of excitons to be found
simple_in%thr_evc	Threshold for defining converged excitonic eigenenergies [eV]
simple_in%diago	Algorithm for finding excitonic eigenvectors:
	(0) steepest descent (1) conjugate gradient
simple_in%max_nstep	Maximum number of steps during minimizations
simple_in%lanczos_step	Number of steps for the Haydock recursive method

Equations

$${\epsilon }_{\text{IP}}(\mathbf{q},\omega )=1-\frac{4\pi }{|\mathbf{q}{|}^2}\frac{1}{V}\sum _{\mathbf{k}}\sum _{n,n^{\prime }}\left(f_{n\mathbf{k}}-f_{n^{\prime }\mathbf{k}+\mathbf{q}}\right)\frac{|⟨{\psi }_{n^{\prime }\mathbf{k}+\mathbf{q}}|e^{i\mathbf{q}\cdot \mathbf{r}}|{\psi }_{n\mathbf{k}}⟩{|}^2}{\omega -(E_{n^{\prime }\mathbf{k}+\mathbf{q}}-E_{n\mathbf{k}})+i\eta }$$

$${\psi }_{n\mathbf{k}}(\mathbf{x})=e^{i\mathbf{k}\cdot \mathbf{r}}{u}_{n\mathbf{k}}(\mathbf{x}),$$

$${\epsilon }_{\text{BSE}}(\mathbf{q},\omega )=1-\frac{4\pi }{|\mathbf{q}{|}^2}\frac{1}{V}\sum _{\lambda }\frac{{\left|{\sum }_{v,c,\mathbf{k}}⟨{\psi }_{v\mathbf{k}+\mathbf{q}}|e^{i\mathbf{q}\cdot \mathbf{r}}|{\psi }_{c\mathbf{k}}⟩A_{\lambda ,\mathbf{q}}^{(vc\mathbf{k})}\right|}^2}{\omega -{\tilde{E}}_{\lambda }(\mathbf{q})+i\eta },$$

$$H_{vc\mathbf{k},v^{\prime }{c}^{\prime }{\mathbf{k}}^{\prime }}^{\text{exc}}(\mathbf{q})A_{\lambda ,\mathbf{q}}^{(v^{\prime }{c}^{\prime }{\mathbf{k}}^{\prime })}={\tilde{E}}_{\lambda }(\mathbf{q})A_{\lambda ,\mathbf{q}}^{(vc\mathbf{k})},$$

$$⟨\mathbf{x}⟩{\psi }_{v\mathbf{k}+\mathbf{q}}{\psi }_{c\mathbf{k}}={\psi }_{v\mathbf{k}+\mathbf{q}}\left(\mathbf{x}\right){\psi }_{c\mathbf{k}}\left(\mathbf{x}\right)$$

$${\epsilon }_{\text{IP}}(\hat{\mathbf{q}},\omega )={\epsilon }_{\text{IP}}^{\text{inter}}(\hat{\mathbf{q}},\omega )+{\epsilon }_{\text{IP}}^{\text{intra}}(\hat{\mathbf{q}},\omega ),$$

$${\epsilon }_{\text{IP}}^{\text{inter}}(\hat{\mathbf{q}},\omega )$$

$${\epsilon }_{\text{IP}}^{\text{intra}}(\hat{\mathbf{q}},\omega )$$

$${\omega }_{\text{D}}^2(\hat{\mathbf{q}})=\frac{4\pi }{V}\sum _{\mathbf{k}}\sum _n|⟨{\psi }_{n\mathbf{k}}|\hat{\mathbf{q}}\cdot \mathbf{v}|{\psi }_{n\mathbf{k}}⟩{|}^2\left(-\frac{\partial f}{\partial E}\right),$$

$${\epsilon }_{\text{BSE}}(\hat{\mathbf{q}},\omega )=1-\frac{4\pi }{V}\sum _{\lambda }\frac{1}{\omega -{\tilde{E}}_{\lambda }+i\eta }{\left|\sum _{v,c,\mathbf{k}}\frac{⟨{\psi }_{v\mathbf{k}}|\hat{\mathbf{q}}\cdot \mathbf{v}|{\psi }_{c\mathbf{k}}⟩}{E_{c\mathbf{k}}-E_{v\mathbf{k}}}{A}_{\lambda }^{(vc\mathbf{k})}\right|}^2,$$

$$H_{vc\mathbf{k},v^{\prime }{c}^{\prime }{\mathbf{k}}^{\prime }}^{\text{exc}}=\left(E_{c\mathbf{k}}-E_{v{\mathbf{k}}^{\prime }}\right){\delta }_{\mathbf{k},{\mathbf{k}}^{\prime }}{\delta }_{c,c^{\prime }}{\delta }_{v,v^{\prime }}+K_{vc\mathbf{k},v^{\prime }{c}^{\prime }{\mathbf{k}}^{\prime }}^x+K_{vc\mathbf{k},v^{\prime }{c}^{\prime }{\mathbf{k}}^{\prime }}^d.$$

$$K_{vc\mathbf{k},v^{\prime }{c}^{\prime }{\mathbf{k}}^{\prime }}^x$$

$$K_{vc\mathbf{k},v^{\prime }{c}^{\prime }{\mathbf{k}}^{\prime }}^d$$

$$u_{n\mathbf{k}}(\mathbf{x})\simeq \sum _{i=1}^{N_b}{\tilde{b}}_i^{n\mathbf{k}}{b}_i(\mathbf{x}).$$

$$H^{\text{KS}}(\mathbf{k})=\frac{(-i\nabla +\mathbf{k}{)}^2}{2}+V^{loc}+V^{nl}(\mathbf{k}).$$

$$b_i(\mathbf{x})=\sum _{\mathbf{G}}{e}^{i\mathbf{G}\cdot \mathbf{r}}{b}_i(\mathbf{G}),$$

$$H^{\text{KS}}_{ij}(\mathbf{k})=\bra{b_{i}}H^{\text{KS}}(\mathbf{k})\ket{b_{j}}=\frac{1}{2}\big{[}k^{2}\delta_{ij}+\mathbf{k}\cdot\mathbf{K}_{ij}^{(1)}+K_{ij}^{(0)}\big{]}+V_{ij}^{loc}+V_{ij}^{nl}(\mathbf{k}).$$

$${\mathbf{K}}_{ij}^{(1)}$$

$$K_{ij}^{(0)}$$

$$V_{ij}^{loc}=\int d\mathbf{r}{b}_i^{\ast }(\mathbf{x})V^{loc}(\mathbf{r})b_j(\mathbf{x}).$$

$$V_{ij}^{nl}(\mathbf{k})=\sum _{\lambda }{\beta }_{\lambda i}^{\ast }(\mathbf{k})D_{\lambda }{\beta }_{\lambda j}(\mathbf{k}),$$

$$\mathbf{v}(\mathbf{k})=-i[\mathbf{r},H^{\text{KS}}(\mathbf{k})]=-i\nabla +\mathbf{k}-i[\mathbf{r},V^{nl}(\mathbf{k})].$$

$$\begin{array}{rl}⟨u_{n\mathbf{k}}|\mathbf{v}(\mathbf{k})|u_{n^{\prime }\mathbf{k}}⟩& =\frac{1}{2}\sum _{i,j=1}^{N_b}({\tilde{b}}_i^{n\mathbf{k}}{)}^{\ast }{\tilde{b}}_j^{n^{\prime }\mathbf{k}}{\mathbf{K}}_{ij}^{(1)}+\mathbf{k}⟨u_{n\mathbf{k}}⟩u_{n^{\prime }\mathbf{k}}+\\ & +(-i)\sum _{i,j=1}^{N_b}({\tilde{b}}_i^{n\mathbf{k}}{)}^{\ast }{\tilde{b}}_j^{n^{\prime }\mathbf{k}}⟨b_i|[\mathbf{r},V^{nl}(\mathbf{k})]|b_j⟩.\end{array}$$

$${\tilde{b}}_i^{n\mathbf{k}}=⟨b_i⟩u_{n\mathbf{k}}.$$

$$b_i^{\ast }(\mathbf{x})b_j(\mathbf{x})\simeq \sum _{\alpha =1}^{N_p}{F}_{ij}^{\alpha }{B}_{\alpha }(\mathbf{x}),$$

$$F_{ij}^{\alpha }=⟨B_{\alpha }⟩b_i^{\ast }{b}_j.$$

$$|u_{n\mathbf{k}}^{\ast }{u}_{n^{\prime }\mathbf{k}}⟩\simeq \sum _{i,j=1}^{N_b}({\tilde{b}}_i^{n\mathbf{k}}{)}^{\ast }{\tilde{b}}_j^{n^{\prime }\mathbf{k}}|b_i^{\ast }{b}_j⟩\simeq \sum _{\alpha =1}^{N_p}\sum _{i,j=1}^{N_b}({\tilde{b}}_i^{n\mathbf{k}}{)}^{\ast }{\tilde{b}}_j^{n^{\prime }\mathbf{k}}{F}_{ij}^{\alpha }|B_{\alpha }⟩.$$

$$J_{\alpha ,n{n}^{\prime }}^{\mathbf{k}}=\sum _{i,j=1}^{N_b}({\tilde{b}}_i^{n\mathbf{k}}{)}^{\ast }{\tilde{b}}_j^{n^{\prime }\mathbf{k}}{F}_{ij}^{\alpha },$$

$$|u_{n\mathbf{k}}^{\ast }{u}_{n^{\prime }\mathbf{k}}⟩\simeq \sum _{\alpha =1}^{N_p}{J}_{\alpha ,n{n}^{\prime }}^{\mathbf{k}}|B_{\alpha }⟩.$$

$$K_{vc\mathbf{k},v^{\prime }{c}^{\prime }{\mathbf{k}}^{\prime }}^x$$

$$V_{\alpha \beta }^{{\mathbf{k}}^{\prime \prime }}=\int d\mathbf{x}d{\mathbf{x}}^{\prime }{B}_{\alpha }^{\ast }\left(\mathbf{x}\right)e^{i{\mathbf{k}}^{\prime \prime }\cdot \mathbf{r}}v\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)e^{-i{\mathbf{k}}^{\prime \prime }\cdot {\mathbf{r}}^{\prime }}{B}_{\beta }\left({\mathbf{x}}^{\prime }\right).$$

$$K_{vc\mathbf{k},v^{\prime }{c}^{\prime }{\mathbf{k}}^{\prime }}^{d,bare}$$

Full text

††journal: Computer Physics Communications

Program summary

Program Title: SIMPLE

Licensing provisions: GNU General Public License V 2.0

Programming language: Fortran 95

Distribution format: tar.gz

Computer: Any computer architecture

External routines: Quantum ESPRESSO distribution, BLAS, LAPACK, FFTW, MPI.

Nature of problem: First-principles calculations of the optical properties of metals and insulators.

Solution method: Shirley’s optimal basis for the complex dielectric function, at the independent-particle level for metals and by solving the Bethe-Salpeter equation for insulators.

Restrictions: Norm-conserving pseudopotentials.

1 Introduction

Density-functional theory (DFT) has largely proven itself to be a powerful and reliable tool to predict structural and vibrational properties (and, in general, any ground-state property) of complex materials. This is to be contrasted with well-known intrinsic limitations in the description of electronic excitations [1]. However, complementary first-principles methods for calculating optical properties have been developed and applied to the study of materials, giving thus birth in recent years to the field of theoretical spectroscopies. In particular, the Bethe-Salpeter equation (BSE), based on many-body perturbation theory [2, 3, 4], has become the state-of-the-art approach for the simulation of optical spectroscopies, such as absorption or reflectivity, of extended systems. Although usually at least one order of magnitude more demanding in terms of computational cost with respect to typical DFT simulations, the BSE approach and the $GW$ method [5] have proven to be the most accurate and predictive methodologies for the study of neutral and charged excitations in extended systems, respectively (we mention in passing the development of Koopmans’ compliant spectral functionals as an alternative to $GW$ [6, 7, 8, 9]). Leveraged by the computational power available nowadays and due to the large range of technologically relevant applications, such as photovoltaics, plasmonics and optoelectronics, first-principles methods for the calculation of optical properties are starting to be routinely used not just to help analyzing and interpreting experimental data but also to guide the discovery of novel materials by computational design [10]. This trend is evident from the fact that methods for theoretical spectroscopies are now being implemented and becoming available to researchers in different specialized codes, such as Yambo [11] or BerkeleyGW [12] (see Ref. [13] for a more comprehensive list), and interfaced with broadly available DFT engines (e.g. Quantum ESPRESSO [14], Abinit [15] and Siesta [16] or directly incorporated in these [17, 18, 19, 20]).

We present here the SIMPLE code, which has been implemented as a component of the Quantum ESPRESSO (QE) distribution. The main purpose of the code is to provide accurate and, at the same time, efficient first-principles simulations of the optical properties of both insulating and metallic materials. For metals, where the electron-hole interaction is screened by conduction electrons, the independent-particle (IP) approximation including both interband and intraband contributions has been implemented. For insulators instead, where the inclusion of the electron-hole interaction is often compulsory to correctly simulate optical spectra, the BSE approach has been implemented. In order to increase the efficiency of the calculation of the IP and BSE optical properties we exploit the method originally developed by E. L. Shirley in 1996 [21]. It is based on the use of an optimal basis (OB) set to represent the periodic part of the Bloch wavefunctions, and in this work we introduce a generalization of the method to represent also products of periodic functions. The most attractive feature of this approach resides in the fact that it allows a significant reduction of the computational cost of the simulations by introducing a well-controlled numerical error.

The paper is organized as follows: in Section 2 we introduce the most important theoretical concepts deployed in the code. In Section 3 we describe the main algorithms implemented, in particular the OB method and its application and extension to calculate optical properties of metals at the IP level and of insulators by solving the BSE. Then in Section 4 the structure of the SIMPLE code is presented and the main input parameters are described in detail. Finally in Section 5 and Section 6 the verification and the performance of the SIMPLE code are discussed.

2 Theoretical background

A natural physical quantity to study in first-principles simulations of optical properties is the full dielectric function, a complex (as opposite to real) dynamical (i.e. frequency-dependent) quantity. In the following we will show the basic equations that relate this macroscopic property of a material to its underling microscopic electronic structure, both within the framework of the IP approximation and of the BSE approach. Throughout the paper we use Hartree atomic units ($\hbar=m=e^{2}=1$) and Gaussian units ($4\pi\varepsilon_{0}=1$) for electromagnetism.

The general expression of the IP dielectric function111 Here we neglect the fact that the dielectric function is actually a $3\times 3$ complex tensor whose additional off-diagonal terms describe subtle magnetic phenomena beyond the scope of this work, such as the Kerr effect or the anomalous Hall effect [22, 23]. In the following we will deal exclusively with the diagonal terms which are specified by the direction of the transferred momentum $\mathbf{q}$. (valid both for metals and insulators) at a generic transferred momentum $\mathbf{q}$ and frequency $\omega$ is [24]

[TABLE]

which is written in terms of the solutions of the one-particle Schrödinger equation, $H^{\text{KS}}\ket{\psi_{n\mathbf{k}}}=E_{n\mathbf{k}}\ket{\psi_{n\mathbf{k}}}$, where $H^{\text{KS}}$ is the Kohn-Sham (KS) Hamiltonian. $f_{n\mathbf{k}}$ gives the occupation according to the Fermi-Dirac distribution of the KS Bloch state $\ket{\psi_{n\mathbf{k}}}$ with band index $n$ and wavevector $\mathbf{k}$ while $V$ is the volume of the crystal. Instead, $\eta$ is an infinitesimal broadening introduced to perform the adiabatic switching-on of the perturbation and which, in practice, is used as an empirical broadening with the purpose to account for scattering processes, always present in real materials, and/or for finite experimental resolution.

We recall here that the Bloch wavefunctions have the form

[TABLE]

where $\mathbf{x}$ is a combined position and spin coordinate, $\mathbf{x}=(\mathbf{r},\sigma)$, and the function $u_{n\mathbf{k}}(\mathbf{x})$ is periodic with the periodicity of the crystal. In the non-collinear case, where spin-orbit interaction is explicitly included in the DFT calculation, $n$ has to be understood as a spinorial band index and $\ket{\psi_{n\mathbf{k}}}$ as a two-component spinor.

On the other hand the dielectric function obtained from the solution of the BSE within the Tamm-Dancoff approximation is given by an expression [25] (valid for finite-gap systems) very similar to the IP case of Eq. 1

[TABLE]

where the broadening $\eta$ has to be interpreted as in Eq. 1. The main difference with respect to the IP expression is that Eq. 3 is written in terms of the solutions of an effective two-particle Schrödinger-like equation

[TABLE]

where $H^{\text{exc}}_{vc\mathbf{k},v^{\prime}c^{\prime}\mathbf{k}^{\prime}}(\mathbf{q})$ is the excitonic Hamiltonian matrix222In Eq. 3 it is understood that the long-range component of the bare Coulomb interaction in the excitonic Hamiltonian is not included. in the basis of the single-particle Bloch wavefunctions of electron and hole states. In this basis, $A_{\lambda,\mathbf{q}}^{(vc\mathbf{k})}=\braket{\psi_{v\mathbf{k}+\mathbf{q}}\psi_{c\mathbf{k}}}{A_{\lambda,\mathbf{q}}}$ are the excitonic amplitudes labelled by the quantum numbers $(\lambda,\mathbf{q})$ and having energy $\tilde{E}_{\lambda}(\mathbf{q})$333 We have introduced the notation used throughout the paper

$\braket{\mathbf{x}}{\psi_{v\mathbf{k}+\mathbf{q}}\psi_{c\mathbf{k}}}=\psi_{v\mathbf{k}+\mathbf{q}}\left(\mathbf{x}\right)\psi_{c\mathbf{k}}\left(\mathbf{x}\right)$

, where the indeces $v$ and $c$ stand for valence band and conduction band, respectively.. We note in passing that both for Eq. 1 and Eq. 3 $GW$ energies (at the $G_{0}W_{0}$ level) or even eigenvectors (e.g. at the self-consistent $GW$ level) could be substituted for the KS ones; this is often the practice, especially for Eq. 3, albeit at a largely increased computational cost.

We are interested in studying the optical limit, $\mathbf{q}\to 0$, for which the transferred momentum of the photons is negligible with respect to the size of the Brillouin zone (BZ). In general the dielectric function depends on the direction $\hat{\mathbf{q}}=\mathbf{q}/|\mathbf{q}|$ of the perturbing electric field and only for crystals with cubic symmetry the dielectric tensor is equivalent to a scalar. In the following we will study the optical limit of Eq. 1 and Eq. 3 and we will keep explicit the dependence on the direction $\hat{\mathbf{q}}$ for the general case of anisotropic systems.

2.1 IP

In metals, the IP dielectric function of Eq. 1 in the optical limit can be divided into two separate contributions, an intraband Drude-like term due to the conduction electrons at the Fermi surface and an interband term due to vertical transitions between occupied and unoccupied bands [26, 27, 28]

[TABLE]

where

[TABLE]

We have defined the Drude plasma frequency as

[TABLE]

where $\mathbf{v}=-i\,[\mathbf{r},H^{\text{KS}}]$ is the velocity operator and $\gamma$ is an empirical broadening representing dissipation effects of the conduction electrons, whose value can be taken from experimental data or estimated from additional first-principles calculations [29]. In the limit of zero temperature, the derivative of the Fermi-Dirac distribution in Eq. 8 is a Dirac-delta function centred at the chemical potential and the Drude plasma frequency $\omega_{\text{D}}(\hat{\mathbf{q}})$ is always zero in finite-gap systems.

2.2 BSE

The expression in Eq. 3 of the BSE dielectric function, in the optical limit $\mathbf{q}\rightarrow 0$, is given by [30]

[TABLE]

where we construct the excitonic Hamiltonian considering only vertical electron-hole transitions ($\mathbf{q}=\mathbf{0}$). This means that $\tilde{E}_{\lambda}$ and $A_{\lambda}^{(vc\mathbf{k})}$ in Eq. 9 are the eigenvalues and eigenfunctions of the Hamiltonian $H^{\text{exc}}_{vc\mathbf{k},v^{\prime}c^{\prime}\mathbf{k}^{\prime}}=H^{\text{exc}}_{vc\mathbf{k},v^{\prime}c^{\prime}\mathbf{k}^{\prime}}(\mathbf{q}=\mathbf{0})$.

The excitonic Hamiltonian at $\mathbf{q}=0$ can be divided into three different contributions [31]:

[TABLE]

The first term on the right-hand side of Eq. 10 is called the diagonal term and corresponds to the IP approximation444In the diagonal term we will consider here differences of KS energies and not of quasi-particle energies (obtained, for example, within the $GW$ approximation or from Koopmans’ compliant spectral functionals [8]). In our implementation the KS energies can be corrected by means of an user-defined scissor operator that performs a rigid shift of the bands.. Indeed it is worth noting that Eq. 9 reduces to Eq. 6 (in case of insulators) if the other remaining two terms in Eq. 10, that give rise to the electron-hole interactions, are neglected.

$K^{d}$, called direct term, describes the attraction between the electron and hole involving the screened Coulomb interaction $W$ while $K^{x}$, called exchange term, describes the repulsion between the electron and hole involving the bare Coulomb interaction $v$. Their explicit expressions are

[TABLE]

In the non-collinear case the notation $\int d\mathbf{x}$ implies that we perform both an integration over the space variable $\mathbf{r}$ and a summation over the spin variable $\sigma$ $\left(\int d\mathbf{x}\equiv\sum_{\sigma}\int d\mathbf{r}\right)$.

3 Implementation

3.1 Optimal basis (OB)

The basic idea of the optimal basis (OB) method [21] is to obtain a reduced set of basis functions, indicated with the notation $\{b_{i}\}$, to represent the periodic part of the Bloch wavefunctions at any k-point inside the BZ. The OB $\{b_{i}\}$ are constructed starting from the periodic KS states $\{u_{n\mathbf{k}}\}$ calculated on an initial grid of $N_{k}$ k-points. The band index $n$, which ranges from 1 to $N$, runs over the selected $N_{v}$ top valence and $N_{c}$ lowest conduction bands ($N=N_{v}+N_{c}$). We use a Gram-Schimdt orthonormalization algorithm with a threshold $s_{b}$ which proceeds k-point by k-point, as illustrated in Algorithm 1. With this approach we disregard basis vectors which would contribute marginally (i.e. below the threshold) to the periodic KS states; the dimension $N_{b}$ of the OB is directly dictated by the threshold $s_{b}$.

It is worth noting that the algorithm proceeds on blocks of $N$ states at each time. This permits, on one hand, the use of efficient BLAS matrix-matrix multiplication routines for the orthonormalization and, on the other hand, to avoid dealing with too large matrices, as it would be the case when working with all the periodic KS states at once.

Using this basis we can approximate a generic, periodic KS state with wavevector $\mathbf{k}$ and band index $n$ as

[TABLE]

3.2 IP

Once the OB is constructed it is possible to obtain the periodic part of the Bloch wavefunctions at a generic k-point following the interpolation procedure described in Ref. [32]. This allows to perform the fine samplings of the BZ required in Eq. 6 and in Eq. 7 for metals. For this purpose, we need to construct the matrix elements of the k-dependent KS Hamiltonian $H^{\text{KS}}(\mathbf{k})\equiv e^{-i\mathbf{k}\cdot\mathbf{r}}H^{\text{KS}}e^{i\mathbf{k}\cdot\mathbf{r}}$ in terms of the OB. By diagonalizing this matrix we can then obtain the coefficients $\tilde{b}_{i}^{n\mathbf{k}}$ and the band energies $E_{n\mathbf{k}}$ for each point $\mathbf{k}$ in the dense interpolation k-grid.

Band interpolation

The operator $H^{\text{KS}}(\mathbf{k})$ is divided into three separate contributions (kinetic energy, local potential and non-local potential):

[TABLE]

In our implementation the OB functions are expanded in a plane-wave basis as

[TABLE]

where the sum is over the reciprocal lattice vectors $\mathbf{G}$ and {$b_{i}(\mathbf{G})$} is the set of Fourier coefficients of $b_{i}(\mathbf{x})$. By using the expression above we obtain the matrix elements of the KS Hamiltonian

[TABLE]

The terms inside the square brackets in Eq. 16 refer to the matrix elements of the k-dependent kinetic energy, which have a simple quadratic polynomial dependence on $\mathbf{k}$. For convenience we define two additional quantities, $\mathbf{K}_{ij}^{(1)}$ and $K_{ij}^{(0)}$, in terms of the Fourier coefficients $b_{i}(\mathbf{G})$:

[TABLE]

The matrix elements of the local potential do not depend on $\mathbf{k}$ and are easily obtained in real space using fast-Fourier transforms

[TABLE]

For the non-local potential, in this work we only consider norm-conserving pseudopotentials so that the non-local part of the pseudopotential can be written as

[TABLE]

where the index $\lambda=(I,n,m,l)$ refers to the sites $I$ of the ions in the cell together with the associated atomic quantum numbers $(n,m,l)$ and $D_{\lambda}$ are the coefficients of the pseudopotential. In the expression above $\beta_{\lambda i}(\mathbf{k})=\bra{\beta_{\lambda}}e^{i\mathbf{k}\cdot\mathbf{r}}\ket{b_{i}}$, where $\ket{\beta_{\lambda}}$ are the pseudopotential projector functions centered on each ionic site $I$.

Since $\mathbf{K}_{ij}^{(1)}$, $K_{ij}^{(0)}$ and $V_{ij}^{loc}$ do not depend on $\mathbf{k}$, they need to be calculated only once after the OB is built. The matrix elements $V_{ij}^{nl}(\mathbf{k})$ of the non-local part of the pseudopotential instead have a non-analytic dependence over $\mathbf{k}$ and thus they should be calculated for every $\mathbf{k}$. The calculation of the matrix elements $H^{\text{KS}}_{ij}(\mathbf{k})$ and the subsequent diagonalization of the matrix for each point $\mathbf{k}$ give the coefficients $\tilde{b}_{i}^{n\mathbf{k}}$ and the band energies $E_{n\mathbf{k}}$ for all the bands $n\leq N_{b}$.

Dielectric function

Once the coefficients and the band energies at all the needed k-points are known, we have almost all the elements necessary to calculate the IP dielectric function through Eq. 6 and Eq. 7. For this purpose we also need the matrix elements of the k-dependent velocity operator

[TABLE]

In terms of the OB these are given by

[TABLE]

The local contribution to the velocity matrix elements (first two terms in the right-hand side of Eq. 22) are easily obtained once the KS Hamiltonian matrix in Eq. 16 is diagonalized. On the other hand, the matrix elements $\bra{b_{i}}[\mathbf{r},V^{nl}(\mathbf{k})]\ket{b_{j}}$ of the commutator of the non-local part of the pseudopotential need to be computed and, having a non-analytic dependence over $\mathbf{k}$, they should be calculated for each $\mathbf{k}$.

3.3 BSE

For BSE calculations in finite-gap systems we do not interpolate over the k-points, as in the IP calculation, but we work with the initial k-grid as there is no need to perform an interpolation since significantly coarser k-grids with respect to metals are usually sufficient to obtain converged results. In this case the band energies are known from the initial DFT calculation, while the coefficients $\tilde{b}_{i}^{n\mathbf{k}}$ are simply calculated through the scalar products

[TABLE]

Optimal product basis (OPB)

We then build a basis, that we call optimal product basis (OPB), for representing in real space products of the periodic parts of the Bloch wavefunctions at the same k-point: $u_{n\mathbf{k}}^{*}\left(\mathbf{x}\right)u_{n^{\prime}\mathbf{k}}\left(\mathbf{x}\right)$. We construct this optimally small basis $\left\{B_{\alpha}\right\}$ using a Gram-Schmidt orthonormalisation algorithm with a threshold starting from the products $b_{i}^{*}\left(\mathbf{x}\right)b_{j}\left(\mathbf{x}\right)$ of OB functions obtained in the first step. The threshold $s_{p}$ controls the precision of the OPB and hence its dimension, that we indicate with $N_{p}$. We outline this strategy in Algorithm 2.

Also in this case we can take advantage of the BLAS library for matrix-matrix multiplications, avoiding at the same time to deal with large matrices. Using this new basis, we approximate the product of two OB functions as

[TABLE]

with

[TABLE]

We can now express products of periodic KS states as

[TABLE]

We then define

[TABLE]

so that, similarly to Eq. 13, $u_{n\mathbf{k}}^{*}\left(\mathbf{x}\right)u_{n^{\prime}\mathbf{k}}\left(\mathbf{x}\right)$ is approximated as

[TABLE]

Exchange and direct term

The reduction of the cost of the BSE calculation is achieved by evaluating only on the product basis $\left\{B_{\alpha}\right\}$ the matrix elements of the bare ($v$) and of the screened ($W$) Coulomb interactions required for the exchange and the direct terms of the excitonic Hamiltonian. The number of such elements is much smaller then all the elements defining the $K^{x}_{vc\mathbf{k},v^{\prime}c^{\prime}\mathbf{k}^{\prime}}$ and $K^{d}_{vc\mathbf{k},v^{\prime}c^{\prime}\mathbf{k}^{\prime}}$ matrices. These are then obtained through:

[TABLE]

where we have defined

[TABLE]

We divide the screened Coulomb interaction in a bare and a correlation part: $W=v+W_{c}$. The bare part of the direct operator reads

[TABLE]

where we have introduced the matrices

[TABLE]

The correlation part of the direct operator instead reads

[TABLE]

where we have defined

[TABLE]

Treatment of the Coulomb interaction

We evaluate Eq. 29 and Eq. 31 in a plane-wave basis where the Coulomb interaction is diagonal. However, some care is needed because in reciprocal space terms like

[TABLE]

appear, which diverge for $\mathbf{k}^{\prime}-\mathbf{k}\rightarrow 0,\mathbf{\;G}\rightarrow 0$ (it is worth noting that these divergencies are not physical but come from the numerical discretization of the BZ). To avoid these divergencies we perform the following replacement

[TABLE]

i.e. we integrate numerically over the reduced volume of the Brillouin’s zone (BZR) centered in $\mathbf{k}^{\prime}-\mathbf{k}$ and of volume $\frac{\left(2\pi\right)^{3}}{\Omega N_{q}}$ , where $\Omega$ is the volume of the simulation cell and $N_{q}$ is the total number of k-points of a uniform mesh in the BZ. The integration is performed as a sum over an equally-spaced grid which does not contain the $\Gamma$ point. This is equivalent to the strategy of Ref. [33], where instead the sum is performed through Monte-Carlo sampling. In the present implementation the correlation part of the screened Coulomb interaction ($W_{c}$) is either approximated using the model described in Ref. [34] or is taken from a $GW$ calculation performed with the GWL package of the QE distribution [20].

The model screened Coulomb interaction is obtained from the following expression for the inverse microscopic dielectric function

[TABLE]

The model depends on two parameters: the macroscopic dielectric constant $\epsilon_{m}$ and the screening length $\overline{\lambda}$. From $\varepsilon^{-1}_{\mathrm{model}}(\mathbf{G},\mathbf{G}^{\prime};\mathbf{q})$ an approximated screened interaction is earned

[TABLE]

where $\mathbf{q}=\mathbf{k}^{\prime}-\mathbf{k}$ is the transferred momentum. $W_{\mathrm{model}}(\mathbf{G},\mathbf{G}^{\prime};\mathbf{q})$ is then Fourier transformed to be inserted in Eq. 31.

If $W_{c}$ is taken instead from a GWL calculation, some care is required because GWL samples the BZ at the $\Gamma$ point only, even if the long-range parts of the dielectric matrix (commonly known as “wings” of the matrix) are properly calculated on a discrete k-points mesh. In this way the convergence of $GW$ energies is achieved for not too large simulation cells. Indeed, the GWL code is optimized to work with large model structures, as described in Ref. [35]. However, denser k-point samplings are required for a subsequent BSE calculation. For this purpose, we extrapolate $W_{c}$ at any k-point $\mathbf{q}$ from the polarizability matrix $\Pi$ calculated at the $\Gamma$ point and at zero frequency

[TABLE]

where $\left\{\Phi_{\mu}\right\}$ is a generic and real basis set used to represent the polarizability operator. Although in our calculations the $\left\{\Phi_{\mu}\right\}$ are found according to the so-called optimal polarizability basis recipe [36, 35], our formalism remains unaltered if a plane-waves basis set is chosen. We stress that the expression above is formally exact only at $\mathbf{q}=\mathbf{0}$ and its use at a generic $\mathbf{q}=\mathbf{k}^{\prime}-\mathbf{k}\neq\mathbf{0}$ is justified only if we assume that KS bands and wavefunctions do not change significantly at different k-points. This flat band approximation works well in the case of large simulation cells (or supercells).

Single excitonic states

We find the lowest eigenvalues and eigenvectors of the excitonic Hamiltonian matrix by solving Eq. 4. Instead of the direct diagonalization of the entire matrix, which would be too costly, we use an iterative conjugate-gradient (CG) algorithm [37] so that only the lowest states are found. The CG algorithm proceeds optimizing the search direction: at each minimization step the minimum along the search direction is found, assuming a parabolic behaviour for the excitonic energy. The minimum is found from the energy and its line derivative at the initial point and the energy calculated at a trial point, which is defined by a factor $\mu$ and is given as an input parameter.

Dielectric function

For the calculation of the dielectric function given in Eq. 9 we neglect, in the present implementation, the contribution coming from the non-local commutator (see Eq. 21). In principle, all the excitonic eigenvectors are required so that a full diagonalization of the excitonic Hamiltonian would be necessary. However, iterative algorithms can be used for evaluating directly the frequency-dependent dielectric function. We adopt here the Haydock’s recursive method [38], so that we can calculate $\varepsilon_{\text{BSE}}(\hat{\mathbf{q}},\omega)$ without the need to diagonalize the excitonic Hamiltonian. Within this formalism Eq. 9 can be rewritten as

[TABLE]

In the case of isolated systems the components of the vectors $|P\rangle$ are

[TABLE]

for which only $\Gamma$-point sampling is meaningful, while for extended systems they read

[TABLE]

This allows to eliminate the explicit sum over $\lambda$ in Eq. 37:

[TABLE]

The Haydock’s recursive method provides us with two series of coefficients $\alpha_{i}$ and $\beta_{i}$ so that we can write the continuous fraction as

[TABLE]

The factors $\alpha_{i}$ and $\beta_{i}$ converge to the values $\alpha_{\infty}$ and $\beta_{\infty}$ for large $i$. This permits to terminate the fraction substituting the final $\beta_{n}^{2}/\cdot$ with the terminator $\beta_{n}^{2}/T(\omega)$, where

[TABLE]

4 Details of the code

The SIMPLE code is implemented as a post-processing code of the QE distribution, directly inserted within the GWL package, and its structure is schematically depicted in Fig. 1.

It is divided into three different executable: simple.x, simple_ip.x and simple_bse.x. These could be easily interfaced with different DFT codes as they are separate from the other routines of GWL. The executable simple.x builds the OB and saves to disk all the relevant matrix elements needed for the calculation of the optical properties. The actual IP or BSE calculation is then performed by simple_ip.x or simple_bse.x, respectively. simple.x uses the QE environment while simple_ip.x and simple_bse.x are independent from this environment but make use of several routines and modules of QE. Throughout the SIMPLE code, MPI parallelization is exploited and linear algebra operations are efficiently performed through calls to BLAS and LAPACK libraries. The code works only with norm-conserving pseudopotentials. Spin-orbit interactions are implemented and can be included for both IP and BSE calculations.

Two test examples, i.e. bulk silver for simple_ip.x and bulk silicon for simple_bse.x, are included within the SIMPLE code.

simple.x

simple.x relies, as a starting point, on the results of a nscf calculation of QE performed on a uniform grid of k-points without the use of symmetry (nosym=.true. and noinv=.true. in the input of pw.x). A list of the input parameters of simple.x is reported in Table 1.

As a first step it builds the OB functions as described in Section 3.1 where the threshold $s_{b}$ controlling the precision of the OB is specified by the input variable s_bands (in bohr3). The number of valence bands $N_{v}$ included, starting from the highest occupied band (HOMO), and the number of conduction bands $N_{c}$ included, starting from the lowest unoccupied band (LUMO), are specified by num_val and num_cond, respectively. It then computes and saves to disk the matrix elements needed for the subsequent BSE calculation (if calc_mode=0 in input) or for the subsequent IP calculation (if calc_mode=1 in input).

For a IP calculation the k-independent matrix elements $\mathbf{K}_{ij}^{(1)}$, $K_{ij}^{(0)}$ and $V_{ij}^{loc}$ (see Eq. 17, Eq. 18 and Eq. 19 respectively) of the KS Hamiltonian in the OB are computed. For the non-local part of the pseudopotential we store the k-dependent projectors $\beta_{\lambda j}(\mathbf{k})$ of Eq. 20, as originally proposed in Ref. [32], at each k-point $\mathbf{k}$ on a uniform $n_{1}\times n_{2}\times n_{3}$ k-grid in the BZ defined by the input parameter nkpoints. Besides, if nonlocal_commutator=.true. as it is by default, the matrix elements $\bra{b_{i}}[\mathbf{r},V^{nl}(\mathbf{k})]\ket{b_{j}}$ giving the non-local contribution to the velocity operator (see last term of Eq. 22) are also computed on the k-grid defined by nkpoints.

The calculation of the non-local commutator is the most time and memory consuming part of the simulation. By neglecting this term the computational cost is significantly reduced but the error introduced in the final optical spectra should be checked for each system and pseudopotential considered.

For a BSE calculation the OPB is first constructed following the procedure described in Section 3.3 where the threshold $s_{p}$ controlling the precision of the OPB is specified by the input parameter s_product (in bohr3). The matrix elements in the OPB defined in Eq. 29 and Eq. 31 are then computed ($\mathbf{k}^{\prime\prime}$ in Eq. 29 and Eq. 31 is constructed from the same input parameter nkpoints described above which now, differently from the IP case, has to be the same as the initial nscf k-grid).

The screened Coulomb interaction is either approximated with the model of Eq. 34 if w_type is set to 1, or is taken from a previous GWL calculation if w_type is set to 0. In the first case, the macroscopic dielectric constant (epsm) and the screening length (lambdam) must be given. In the second case, a complete $G_{0}W_{0}$ calculation is performed in order to obtain the screened Coulomb interaction. The polarizability is represented either on a common plane-waves basis or on an optimal polarizability basis set. In both cases its dimension numpw is required as an input parameter. If l_truncated_coulomb is .false. the complete long-range Coulomb interaction is considered, in the opposite case the Coulomb interaction is truncated at a distance truncation_radius. The latter case is meaningful exclusively for $\Gamma$-only sampling.

simple_ip.x

The IP calculation is performed by the executable simple_ip.x and a list of its input parameters is reported in Table 2. For a IP calculation in which a k-space interpolation is performed, the uniform k-grid specified in the starting nscf calculation should include also the seven periodic images of $\Gamma$ at the corners of the unit cube (in units of the reciprocal lattice vectors) in order to better preserve the periodicity in reciprocal space of the k-dependent Hamiltonian (see Ref. [32] for a more detailed discussion). For this purpose we supply a python script (gen_kgrid.py) that can generate a uniform grid of k-points for the QE input with or without periodic images of $\Gamma$ starting from a user-defined k-grid.

The uniform $m_{1}\times m_{2}\times m_{3}$ k-grid on which the IP dielectric function is computed is defined by interp_grid and it should be the same as the variable nkpoints of simple.x. However, in order to reduce the computational cost, a simple linear interpolation of the non-local contributions, i.e. $\beta_{\lambda j}(\mathbf{k})$ and $\bra{b_{i}}[\mathbf{r},V^{nl}(\mathbf{k})]\ket{b_{j}}$, is also implemented so that it is possible to use a denser k-grid in interp_grid with respect to nkpoints (by specifying nonlocal_interpolation = .true.). To have good precision it is not suggested to use a grid larger than the double of the grid specified by nkpoints and the reliability of the results should be always carefully checked.

The code performs the band interpolation on interp_grid and calculates the Fermi energy if it is not already specified in input by fermi_energy. The broadening type and the value of the broadening (in Ry) can be specified in input by the parameters fermi_ngauss and fermi_degauss, respectively.

The occupations of the states are then computed according to the Fermi-Dirac distribution with a broadening specified by elec_temp (in Ry) that by default is set to room temperature. The IP dielectric function is then calculated on an equispaced energy grid of nw points in the energy interval [wmin, wmax] (in Ry).

The derivative of the Fermi-Dirac distribution appearing in the expression of the Drude plasma frequency (see Eq. 8) should be in principle computed with the broadening specified by elec_temp, which however can require extremely dense k-grids. Therefore, a Gaussian broadening is used instead, defined in input by drude_degauss, because converged results can be obtained in practice using coarser k-grids and larger broadenings. The empirical interband and intraband broadenings, i.e. $\eta$ of Eq. 6 and $\gamma$ of Eq. 7, are defined by inter_broadening and intra_broadening, respectively.

The code writes in standard output the Drude plasma frequency $\omega_{\text{D}}(\hat{\mathbf{q}})$, computed in the limit $\mathbf{q}\to\mathbf{0}$ along the three Cartesian directions (i.e. with $\hat{\mathbf{q}}=\hat{\mathbf{x}},\hat{\mathbf{y}},\hat{\mathbf{z}}$), and provides two data files with the real and imaginary part of the interband contribution to the IP dielectric function, $\Re{\left[\varepsilon_{\text{IP}}^{\text{inter}}(\hat{\mathbf{q}},\omega)\right]}$ and $\Im{\left[\varepsilon_{\text{IP}}^{\text{inter}}(\hat{\mathbf{q}},\omega)\right]}$, also computed along the three Cartesian limits. The real and imaginary part of the total IP dielectric function of Eq. 5 and the electron energy loss spectrum (EELS), which is given by $-\Im{\left[\varepsilon_{\text{IP}}^{-1}(\hat{\mathbf{q}},\omega)\right]}$, are also provided.

For completeness the code also computes and writes on file the total dielectric function averaged over the three Cartesian directions

[TABLE]

together with the averaged refractive index $n(\omega)$, extinction coefficient $k(\omega)$ and reflectivity at normal incidence $R(\omega)$. The density of states (DOS) and the joint density of states (JDOS) are also computed.

simple_bse.x

The BSE calculation is performed by the executable simple_bse.x and a list of its input parameters is reported in Table 3. The number of terms in the excitonic Hamiltonian can be chosen through the parameter h_level. Considering only the diagonal term in Eq. 10 corresponds to the IP approximation, adding also the exchange term gives the time-dependent Hartree approximation while adding also the bare Coulomb part of the direct term yields the time-dependent Hartree-Fock approximation. If the correlation part of the direct term is also included, the full BSE approach is recovered. The code supports either spin collinear calculations with doubly occupied DFT states or the general case of non collinear spin (option spin_state). In the first case, exciton states have either singlet or triplet symmetry [31]. The code performs two different tasks (option task): the evaluation of the lowest eigenvectors of the excitonic Hamitonian and of the corresponding energies through the CG algorithm (task=0), and the computation of the $\alpha_{i}$ and $\beta_{i}$ coefficients of Eq. 41 through the Haydock’s recursive scheme (task=1). In the latter case, three series of coefficients relative to the three Cartesian directions are written on disk in a file named ab_coeff.dat for later calculation of the continuous fraction in order to get the components $\varepsilon_{\text{BSE}}(\hat{\mathbf{x}},\omega)$, $\varepsilon_{\text{BSE}}(\hat{\mathbf{y}},\omega)$ and $\varepsilon_{\text{BSE}}(\hat{\mathbf{z}},\omega)$ of the complex frequency-dependent dielectric function. For this purpose we include in the distribution a basic post-processing program (abcoeff_to_eps.f90).

5 Convergence and verification

In this section, first we study the convergence of the optical properties calculated with the SIMPLE code with respect to the relevant computational parameters controlling the precision of the simulations. Subsequently we verify the code by comparing our results for test systems with other well-established and publicly available codes and/or with results in the literature.

5.1 IP

We consider as test system for the IP implementation bulk silver in the FCC primitive cell with lattice parameter $a=7.869$ bohr. The DFT calculations with the pw.x code of QE are performed with the PBE approximation [39] for the exchange-correlation functional at a wavefunction cutoff of 55 Ry and using the scalar-relativistic norm-conserving pseudopotential from the SG15 library [40]. For the ground-state calculation (scf calculation in pw.x) we use a $24\times 24\times 24$ Monkhorst-Pack [41] k-grid while the periodic functions $\{u_{n\mathbf{k}}\}$ needed for the construction of the OB are obtained on a $2\times 2\times 2$ uniform k-grid including the seven periodic images of $\Gamma$ and considering 11 conduction bands555Convergence of the optical spectra with respect to the nscf k-grid has also been studied. In a primitive cell, a $2\times 2\times 2$ k-grid is more than enough to obtain converged results (see Supplementary Material) while in a supercell the $\Gamma$ point is already sufficient (in both cases including the seven periodic images of $\Gamma$). (nscf calculation in pw.x). The IP optical properties are calculated on a $44\times 44\times 44$ uniform k-grid with the inclusion of the non-local contribution to the velocity matrix elements and broadening parameters (i.e. inter_broadening, intra_broadening and drude_degauss of Table 2) of 0.1 eV. Because we are considering a FCC cubic crystal, the dielectric function is isotropic and is simply a scalar (i.e. independent from the direction $\hat{\mathbf{q}}$). The same holds for the test example of the BSE implementation, i.e. bulk silicon, which has the diamond structure.

First, we investigate the convergence of the IP dielectric function with respect to the threshold $s_{b}$ which controls the quality of the OB used for representing the periodic part of the Bloch wavefunctions. In Fig. 2 we show how $\Im{\left[\varepsilon_{\text{IP}}^{\text{inter}}(\hat{\mathbf{q}},\omega)\right]}$ for Ag is fully converged at $s_{b}=0.01$ bohr3, or equivalently with $N_{b}=103$ basis functions (to be compared with the value $N_{b}=315$ corresponding to the complete basis, i.e. to $s_{b}=0.0$ bohr3). Similar results are found for the convergence of $\Re{\left[\varepsilon_{\text{IP}}^{\text{inter}}(\hat{\mathbf{q}},\omega)\right]}$ and $\omega_{\text{D}}(\hat{\mathbf{q}})$, and are not shown here (see Supplementary Material).

Then, we verify the IP implementation of the SIMPLE code by comparing the results obtained for Ag with respect to equivalent simulations performed with the Yambo code [11] for the interband contribution to $\varepsilon_{\text{IP}}(\hat{\mathbf{q}},\omega)$, and with the BoltzWann module [42] of the Wannier90 code [43] for the intraband contribution to $\varepsilon_{\text{IP}}(\hat{\mathbf{q}},\omega)$. As both codes are interfaced with the QE distribution we can perform a straightforward comparison starting from the same DFT ground-state density used for the SIMPLE calculation using all the same relevant computational parameters (k-grid, number of bands, broadenings, etc.). To note that, in order to calculate the optical properties, Yambo does not use any interpolation method in k-space while BoltzWann, instead, uses the Wannier’s interpolation method [44]. For the SIMPLE calculation we use an OB constructed setting $s_{b}=0.01$ bohr3 that, as shown above, gives well converged spectra. Fig. 3 shows that SIMPLE and Yambo give almost identical results for $\varepsilon_{\text{IP}}^{\text{inter}}(\hat{\mathbf{q}},\omega)$. For the intraband contribution we find that the Drude plasma frequency of Ag is $\omega_{\text{D}}(\hat{\mathbf{q}})=8.93$ eV with the SIMPLE code while it is $\omega_{\text{D}}(\hat{\mathbf{q}})=8.96$ eV with the BoltzWann code. Therefore we conclude that for both the interband and intraband contributions the agreement between SIMPLE and the other codes is excellent.

It is worth mentioning here that a powerful approach for the interpolation of band-structure properties is the interpolation mentioned above based on the use of maximally localized Wannier functions (MLWF) as basis [45]. Wannier functions are physically appealing since they are localized in real space and can be used as an exact tight-binding basis. Moreover, Wannier interpolation is both an efficient and precise method for the evaluation of quantities that require a fine sampling of the BZ; in particular, in the Wannier representation the k-dependent KS Hamiltonian does not need to be explicitly constructed and the matrix elements of the velocity operator are analytic if written in terms of Wannier functions [46]. On the other hand, the method has the drawback that the construction of the MLWF’s is not an easily automatizable procedure, especially for metals, for which an additional disentanglement step of the empty bands is required [47]. However, in recent years there have been promising advances towards the development of an automatic procedure for the “Wannierization” of periodic systems [48], including also the more complex case of metals [49]. Once these new approaches will be interfaced with DFT engines, it would be interesting to systematically assess and compare performance and precision of Shirley’s and Wannier’s interpolations.

5.2 BSE

We test the BSE branch of the code considering bulk silicon as this was used, historically, as a test case for the introduction of a variety of BSE implementations [50, 3, 31, 51]. As the correlation part of the screened Coulomb interaction is taken from a GWL calculation, performed at the sole $\Gamma$-point and then extrapolated as described in Section 3.3, we study a $8$ atoms cubic supercell. We work with the local-density approximation [52] (LDA) for the DFT exchange-correlation functional, with a norm-conserving pseudopotential (Si.pz-vbc.UPF) from the QE library and with a wavefunction cutoff of $20$ Ry for defining the plane-wave basis sets for the KS wavefunctions. These parameters yield a theoretical lattice constant of $10.26$ bohr. While for the starting DFT self-consistent run we use a regular $4\times 4\times 4$ mesh k-points, for the convergence tests of the BSE implementation we use a regular $6\times 6\times 6$ mesh of k-points, shifted not to include the $\Gamma$-point. We include in the BSE simulations all the $16$ valence states, and $24$ conduction states. We focus on the complex dielectric function which is obtained using a rigid shift of the valence band manifold, with a scissor parameter of $0.6$ eV for opening the electronic gap, as done in previous studies [50, 31, 3], $1000$ iterative steps for the Haydock’s recursive algorithm and a broadening parameter $\eta$ of $0.2$ eV. We note that using only $\sim 50$ steps is enough to ensure convergence. A polarizability basis [36, 35] of $100$ vectors is enough to converge the screened Coulomb interaction (see Supplementary Material). We show in Fig. 4 (left panel) the calculated $\Im{\left[\varepsilon_{\text{BSE}}(\hat{\mathbf{q}},\omega)\right]}$ as a function of the $s_{b}$ threshold that controls the quality of the OB: a value of $0.1$ bohr3 is sufficient for converging the spectrum, so that the OB comprises only 91 vectors. This number, 8640, should be contrasted with the theoretical dimension of a complete basis (i.e. $s_{b}=0.0$ bohr3), and with the dimension of the plane-wave basis sets for representing the KS wavefunctions, that has 2205 elements. We show in Fig. 4 (right panel) how $\Im{\left[\varepsilon_{\text{BSE}}(\hat{\mathbf{q}},\omega)\right]}$ behaves if we lower the $s_{p}$ threshold controlling the quality of the OPB. For this test the OB threshold is kept fixed to $s_{b}=0.1$ bohr3 giving a basis dimension $N_{b}=91$. We see how the choice $s_{p}=1.0$ bohr3 provides a spectrum almost indistinguishable from the one obtained with $s_{p}=0.5$ bohr3, while the dimension of the OPB raises from $N_{p}=160$ to $N_{p}=443$. From our tests of the sensitivity of the dielectric function on the $s_{b}$ and $s_{p}$ thresholds we opt for the choice $s_{b}=0.1$ bohr3 and $s_{p}=1.0$ bohr3 for verifying our BSE implementation. As the spectra are quite sensitive to k-points sampling [53, 51], we take a denser k-point mesh comprising $448$ points. We compare in Fig. 5 the resulting $\Im{\left[\varepsilon_{\text{BSE}}(\hat{\mathbf{q}},\omega)\right]}$ spectrum with those taken from Refs. [50, 31, 3, 51], registering good agreement despite the different approximations the various implementations rely on. In particular, the strategy used for extrapolating the $W_{c}$ operator is verified.

6 Performance

6.1 Efficiency of the optimal basis

As the core of the present approach is the reduction of the basis sets for representing wavefunctions and their products, we are compelled to assess how such reduction improves the computational performance.

6.1.1 IP

For IP calculations we have to deal with matrices of dimensions $N_{b}\times N_{b}$ (i.e. the k-dependent Hamiltonian $H^{\text{KS}}_{ij}(\mathbf{k})$ of Eq. 16) and the computational cost is dictated by the construction (or “filling”) of this matrix through the computation of the k-dependent contributions coming from the non-local part of the pseudopotentials, rather than by its diagonalization. For this reason the wall time of a IP calculation is expected to scale quadratically with respect to $N_{b}$.

We study the computational cost of the IP approach performed with the SIMPLE code as the number of OB functions $N_{b}$ is varied, as shown in Fig. 6 for bulk silver. We set as reference time the computer time (or wall time) of the calculation with the smallest number of OB functions. As expected, timing scales with the square of $N_{b}$. In Fig. 6 we also show how the inclusion of the non-local commutator in the computation of the velocity matrix elements significantly slows down the calculations.

6.1.2 BSE

For BSE calculations we work with matrices of dimensions $N_{b}\times N_{v}$, $N_{b}\times N_{c}$ and $N_{p}\times N_{\pi}$ where $N_{\pi}$ is the dimension of the basis used to represent the polarizability operator (see Eq. 36). As $N_{\pi}$ does not depend on OB or OPB, we expect the wall time of a BSE calculation to scale linearly with respect to $N_{p}$. Besides, we also know that the number of vectors of the OPB $N_{p}$ scales linearly with $N_{b}$, instead of its square, because of the threshold $s_{p}$ used in the construction of the basis (see Section 3.3). As a consequence, the wall time of a BSE calculation is also expected to scale linearly with respect to $N_{b}$.

We show in Fig. 7 (left panel) how the total computational cost depends almost linearly on the dimension $N_{p}$ of the OPB, as expected (with the $s_{b}$ threshold kept fixed). Indeed, we register a scaling better than linear likely due to a more efficient use of the cache memory. In Fig. 7 (right panel) we show that linear scaling is also found for the variation of the $N_{b}$ size of the OB (with the $s_{p}$ threshold kept fixed). In this case the slightly worse than ideal behavior should be traced to the requirement of storing and working with larger matrices, implying significant transfers from memory to CPU and back.

6.2 Core scaling

In order to assess the performances of the SIMPLE code for its use in parallel architectures, we investigate the speedup of the code as a function of the number of single processing units, or cores, both for the IP and BSE implementation.

6.2.1 IP

For the IP part we consider Ag in a $4\times 4\times 4$ FCC supercell with 64 atoms and 1216 electrons. The total number of OB functions considered in the simulations is 2687, corresponding to $s_{b}=0.01$ bohr3. In general, for IP calculations, a significant part of the total computational cost is due to the simple.x run. Therefore the speedup of the IP section is calculated by summing the wall times of the simple.x and simple_ip.x runs performed with the same number of cores. In Fig. 8 we show the results for the case in which the number of cores corresponds to the number of MPI tasks and we compare simulations performed including and not including the computation of the non-local commutator entering in the estimation of the velocity matrix elements. In order to obtain similar wall times for a given number of cores, we use an interpolation k-grid of $8\times 8\times 8$ and $16\times 16\times 16$ when including and not including the non-local commutator, respectively. Although the scaling of simple_ip.x alone is close to linear, the total scaling of the combined simulations deviates from the ideal linear scaling, mainly because of the computations of the k-dependent matrix elements due to the non-local part of the pseudopotentials (see Eq. 20 and Eq. 22) performed by simple.x. Indeed, if the calculation of the non-local commutator of Eq. 22 is avoided, the scaling is significantly improved, as shown in Fig. 8.

6.2.2 BSE

For the BSE part, in order to illustrate the capabilities of the code, we choose the mixed organic-inorganic perovskite CH3NH3PbI3 which is widely studied for organic/inorganic perovskite solar cells [54]. It involves a primitive cell, comprising 48 inequivalent atoms, fully relativistic calculations and 200 electron. We test the evaluation of the dielectric function through the Haydock’s recursive algorithm. The code is optimized for a mixed OpenMP-MPI parallelization as it strongly relies on BLAS routines and in particular on the DGEMM matrix-matrix multiplication. First, we check the scaling with respect to the number of OpenMP threads keeping fixed the number of MPI tasks. The corresponding results are displayed in Fig. 9 (left panel), where we do not consider the cost of the preceding simple.x run as this is smaller: $3369$ s of wall time on only 128 cores, to be contrasted with $8782$ s for the simple_bse.x run with 256 cores. We note that changing from 4 to 16 the number of OpenMP threads per MPI task determines a reduction of the wall time by a factor $2.1$. Although this figure is lower than the ideal factor $4$, we register an almost linear scaling with the number of OpenMP threads. Indeed, an only slightly larger factor ($2.3$) would be expected extrapolating the scaling from 256 to 512 cores. Such linear scaling supports the parallelization of the single MPI tasks on a relatively large number of OpenMP threads. The scaling of the required wall time with respect to the number of MPI tasks keeping fixed the distribution of the single tasks on OpenMP threads is displayed in Fig. 9 (right panel). Also in this test we focus only on the simple_bse.x run as the simple.x one is much less demanding ($7840$ s of wall time on only $128$ cores versus $12251$ s for the simple_bse.x calculation on $3672$ cores). The number of OpenMP threads per MPI task is kept fixed to $68$ while the number of MPI tasks varies from $54$ to $108$. We register linear scaling and, in this case, we note also an almost ideal behavior. The favorable scaling both with respect to the number of OpenMP threads and of MPI tasks will allow the use of our BSE code on the latest high-performance computing infrastructures running on up to few tens of thousands of cores.

7 Conclusions

We have developed a new code, named SIMPLE, to calculate optical properties using the Shirley’s optimal basis method proposed in Ref. [21], both within the framework of the independent-particle approximation (for metals) and of the Bethe-Salpeter equation (for insulators).

In the independent-particle implementation, we highlight that the code computes both interband and intraband contributions to the complex dielectric function and that the matrix elements of the velocity operator are correctly calculated by including, if desired, the non-local contribution of the norm-conserving pseudopotentials. In the Bethe-Salpeter part, a generalization of the optimal basis method has been developed and implemented to deal in an efficient way also with products of periodic Kohn-Sham states. The screened Coulomb interaction is obtained either through a model for the microscopic dielectric function or starting from the results of a previous GWL calculation. To avoid the diagonalization of the excitonic Hamiltonian matrix, the lowest eigenvectors and corresponding energies can be obtained exploiting an iterative conjugate-gradient algorithm while the complex dielectric function is calculated with the Haydock’s recursive method.

The SIMPLE code has been included within the GWL package of the Quantum ESPRESSO distribution. It has been developed and optimized, as demonstrated by scalability tests, in order to run efficiently in laptops as well as in massively parallel architectures by deploying both fast linear algebra routines (BLAS, LAPACK) and parallelization strategies (MPI and OpenMP). We have verified the results of our code obtained within the framework of the independent-particle approximation and of the Bethe-Salpeter equation on bulk silver and bulk silicon, respectively. In the case of silver we have compared SIMPLE with the publicly available codes Yambo and BoltzWann and have found very good agreement. In the case of silicon, for which the Bethe-Salpeter equation is solved, SIMPLE has been compared with the literature, finding also here a good agreement for the results, considering the use of different base codes and pseudopotentials.

The systematic improvement of the optimal basis set and, as a consequence, of the precision of the results by simply increasing the number of basis functions included, provides the notable advantage that our code can be used to perform efficient computations in a straightforward way for the users: precision is controlled simply by a few input parameters, i.e. $s_{b}$ and $s_{p}$. Moreover, the simplicity of the method and the fact that all the calculations needed to obtain the optical properties from a given initial crystal structure can be run within a single open-source software, i.e. the Quantum ESPRESSO distribution, makes the automation of the sequence of computational steps required and thus the systematic evaluation of the complex dielectric function for a large number of materials much easier. Because of all these considerations, we hope that the SIMPLE code will have an impact in the field of computational materials science by helping researchers to perform efficient, systematic, and reliable first-principles studies of the optical properties of materials.

Acknowledgments

G. P. acknowledges funding from Varinor SA (CH 2800 Delémont, Switzerland). P. U. acknowledges funding from the EU-H2020 research and innovation programme under Grant Agreement No. 654360 NFFA-Europe. This research was partially supported by the NCCR MARVEL, funded by the Swiss National Science Foundation. We acknowledge PRACE (project ID 2016163963) for awarding us access to Marconi at CINECA, Italy. G. P. warmly thanks Gian-Marco Rignanese for originally suggesting the use of the optimal basis method.

Declarations of interest: none