Advanced capabilities for materials modelling with Quantum ESPRESSO

I Introduction

Numerical simulations based on density-functional theory (DFT) Hohenberg and Kohn (1964); Kohn and Sham (1965) have become a powerful and widely used tool for the study of materials properties. Many of such simulations are based upon the “plane-wave pseudopotential method”, often using ultrasoft pseudopotentials Vanderbilt (1990) or the projector augmented wave method (PAW) Blöchl (1994) (in the following, all of these modern developments will be referred to under the generic name of “pseudopotentials”). An important role in the diffusion of DFT-based techniques has been played by the availability of robust and efficient software implementations Lejaeghere et al. (2016), as is the case for Quantum ESPRESSO, which is an open-source software distribution—i.e., an integrated suite of codes—for electronic-structure calculations based on DFT or many-body perturbation theory, and using plane-wave basis sets and pseudopotentials Giannozzi et al. (2009).

The core philosophy of Quantum ESPRESSO can be summarized in four keywords: openness, modularity, efficiency, and innovation. The distribution is based on two core packages, PWscf and CP, performing self-consistent and molecular-dynamics calculations respectively, and on additional packages for more advanced calculations. Among these we quote in particular: PHonon, for linear-response calculations of vibrational properties; PostProc, for data analysis and postprocessing; atomic, for pseudopotential generation; XSpectra, for the calculation of X-ray absorption spectra; GIPAW, for nuclear magnetic resonance and electron paramagnetic resonance calculations.

In this paper we describe and document the novel or improved capabilities of Quantum ESPRESSO up to and including version 6.2. We do not cover features already present in v.4.1 and described in Ref. Giannozzi et al., 2009, to which we refer for further details. The list of enhancements includes theoretical and methodological extensions but also performance enhancements for current parallel machines and modularization and extended interoperability with other software.

Among the theoretical and methodological extensions, we mention in particular:

Fast implementations of exact (Fock) exchange for hybrid functionals Lin (2016); Wu et al. (2009); DiStasio Jr. et al. (2014); Ko et al. ; implementation of non-local van der Waals functionals Berland et al. (2015) and of explicit corrections for van der Waals interactions Grimme (2006); Tkatchenko and Scheffler (2009); Becke and Johnson (2007); Johnson (2017); improvement and extensions of Hubbard-corrected functionals Sclauzero and Dal Corso (2013); Himmetoglu et al. (2011).

Excited-state calculations within time-dependent density-functional and many-body perturbation theories.

Relativistic extension of the PAW formalism, including spin-orbit interactions in density-functional theoryDal Corso (2010, 2012).

Continuum embedding environments (dielectric solvation models, electronic enthalpy, electronic surface tension, periodic boundary corrections) via the Environ module Andreussi et al. (2012); Andreussi and Marzari (2014) and its time-dependent generalization Timrov et al. (2015a).

Several new packages, implementing the calculation of new properties, have been added to Quantum ESPRESSO. We quote in particular:

turboTDDFT Walker et al. (2006); Rocca et al. (2008); Malcioğlu et al. (2011); Ge et al. (2014) and turboEELS Timrov et al. (2013, 2015b), for excited-state calculations within time-dependent DFT (TDDFT), without computing virtual orbitals, also interfaced with the Environ module (see above).

QE-GIPAW, replacing the old GIPAW package, for nuclear magnetic resonance and electron paramagnetic resonance calculations.

EPW, for electron-phonon calculations using Wannier-function interpolation Poncé et al. (2016).

GWL and SternheimerGW for quasi-particle and excited-state calculations within many-body perturbation theory, without computing any virtual orbitals, using the Lanczos bi-orthogonalization Umari et al. (2009, 2010) and multi-shift conjugate-gradient methods Schlipf et al. (2017), respectively.

thermo_pw, for computing thermodynamical properties in the quasi-harmonic approximation, also featuring an advanced master-slave distributed computing scheme, applicable to generic high-throughput calculations Dal Corso (a).

d3q and thermal2, for the calculation of anharmonic 3-body interatomic force constants, phonon-phonon interaction and thermal transport Paulatto et al. (2013); Fugallo et al. (2013).

Improved parallelization is crucial to enhance performance and to fully exploit the power of modern parallel architectures. A careful removal of memory bottlenecks and of scalar sections of code is a pre-requisite for better and extending scaling. Significant improvements have been achieved, in particular for hybrid functionals Varini et al. (2013); Barnes et al. (2017).

Complementary to this, a complete pseudopotential library, pslibrary, including fully-relativistic pseudopotentials, has been generated Dal Corso (b, 2015). A curation effort Castelli et al. on all the pseudopotential libraries available for Quantum ESPRESSO has led to the identification of optimal pseudopotentials for efficiency or for accuracy in the calculations, the latter delivering an agreement comparable to any of the best all-electron codes Lejaeghere et al. (2016). Finally, a significant effort has been dedicated to modularization and to enhanced interoperability with other software. The structure of the distribution has been revised, the code base has been re-organized, the format of data files re-designed in line with modern standards. As notable examples of interoperability with other software, we mention in particular the interfaces with the LAMMPS molecular dynamics (MD) code Plimpton (1995) used as molecular-mechanics “engine” in the Quantum ESPRESSO implementation of the QM-MM methodology Ma et al. (2015), and with the i-PI MD driver Ceriotti et al. (2014), also featuring path-integral MD.

All advances and extensions that have not been documented elsewhere are described in the next sections. For more details on new packages we refer to the respective references.

The paper is organized as follows. Sec. II contains a description of new theoretical and methodological developments and of new packages distributed together with Quantum ESPRESSO. Sec. III contains a description of improvements of parallelization, updated information on the philosophy and general organization of Quantum ESPRESSO, notably in the field of modularization and interoperability. Sec. IV contains an outlook of future directions and our conclusions.

II Theoretical, algorithmic, and methodological extensions

In the following, CGS units are used, unless noted otherwise.

Hybrid functionals are already the de facto standard in quantum chemistry and are quickly gaining popularity in the condensed-matter physics and computational materials science communities. Hybrid functionals reduce the self-interaction error that plagues lower-rung exchange-correlation functionals, thus achieving more accurate and reliable predictive capabilities. This is of particular importance in the calculation of orbital energies, which are an essential ingredient in the treatment of band alignment and charge transfer in heterogeneous systems, as well as the input for higher-level electronic-structure calculations based on many-body perturbation theory. However, the widespread use of hybrid functionals is hampered by the often prohibitive computational requirements of the exact-exchange (Fock) contribution, especially when working with a plane-wave basis set. The basic ingredient here is the action $(\hat{V}_{x}\phi_{i})({\bf r})$ of the Fock operator $\hat{V}_{x}$ onto a (single-particle) electronic state $\phi_{i}$, requiring a sum over all occupied Kohn-Sham (KS) states $\{\psi_{j}\}$. For spin-unpolarized systems, one has:

where $-e$ is the charge of the electron. In the original algorithm Giannozzi et al. (2009) implemented in PWscf, self-consistency is achieved via a double loop: in the inner one the $\psi$’s entering the definition of the Fock operator in Eq. (1) are kept fixed, while the outer one cycles until the Fock operator converges to within a given threshold. In the inner loop, the integrals appearing in Eq. (1):

are computed by solving the Poisson equation in reciprocal space using fast Fourier transforms (FFT). This algorithm is straightforward but slow, requiring {\cal O}\bigl{(}(N_{b}N_{k})^{2}\bigr{)} FFTs, where $N_{b}$ is the number of electronic states (“bands” in solid-state parlance) and $N_{k}$ the number of ${\bf k}$ points in the Brillouin zone (BZ). While feasible in relatively small cells, this unfavorable scaling with the system size makes calculations with hybrid functionals challenging if the unit cell contains more than a few dozen atoms.

To enable exact-exchange calculations in the condensed phase, various ideas have been conceived and implemented in recent Quantum ESPRESSO versions. Code improvements aimed at either optimizing or better parallelizing the standard algorithm are described in Sec. III.1. In this section we describe two important algorithmic developments in Quantum ESPRESSO, both entailing a significant reduction in the computational effort: the adaptively compressed exchange (ACE) concept Lin (2016) and a linear-scaling (${\cal O}(N_{b})$) framework for performing hybrid-functional ab initio molecular dynamics using maximally localized Wannier functions (MLWF) Wu et al. (2009); DiStasio Jr. et al. (2014); Ko et al. .

The simple formal derivation of ACE allows for a robust implementation, which applies straightforwardly both to isolated or aperiodic systems ($\Gamma-$only sampling of the BZ, that is, ${\bf k}=0$) and to periodic ones (requiring sums over a grid of ${\bf k}$ points in the BZ); to norm conserving and ultrasoft pseudopotentials or PAW; to spin-unpolarized or polarized cases or to non-collinear magnetization. Furthermore, ACE is compatible with, and takes advantage of, all available parallelization levels implemented in Quantum ESPRESSO: over plane waves, over ${\bf k}$ points, and over bands.

With ACE, the action of the exchange operator is rewritten as

where $|\xi_{i}\rangle=\hat{V}_{x}|\psi_{i}\rangle$ and $M_{jm}=\langle\psi_{j}|\xi_{m}\rangle$. At self-consistency, ACE becomes exact for $\phi_{i}$’s in the occupied manifold of KS states. It is straightforward to implement ACE in the double-loop structure of PWscf. The new algorithm is significantly faster while not introducing any loss of accuracy at convergence. Benchmark tests on a single processor show a $3\times$ to $4\times$ speedup for typical calculations in molecules, up to $6\times$ in extended systems Carnimeo et al. .

An additional speedup may be achieved by using a reduced FFT cutoff in the solution of Poisson equations. In Eq. (1), the exact FFT algorithm requires a FFT grid containing G-vectors up to a modulus $G_{max}=2G_{c}$, where $G_{c}$ is the largest modulus of G-vectors in the plane-wave basis used to expand $\psi_{i}$ and $\phi_{j}$, or, in terms of kinetic energy cutoff, up to a cutoff $E_{x}=4E_{c}$, where $E_{c}$ is the plane-wave cutoff. The presence of a $1/G^{2}$ factor in the reciprocal space expression suggests, and experience confirms, that this condition can be relaxed to $E_{x}\sim 2E_{c}$ with little loss of precision, down to $E_{x}=E_{c}$ at the price of increasing somewhat this loss Marsili and Umari (2013). The kinetic-energy cutoff for Fock-exchange computations can be tuned by specifying the keyword ecutfock in input.

Hybrid functionals have also been extended to the case of ultrasoft pseudopotentials and to PAW, following the method of Ref. Paier et al., 2005. A large number of integrals involving augmentation charges $q_{lm}$ are needed in this case, thus offsetting the advantage of a smaller plane-wave basis set. Better performances are obtained by exploiting the localization of the $q_{lm}$ and computing the related terms in real space, at the price of small aliasing errors.

These improvements allow to significantly speed up a calculation, or to execute it on a larger number of processors, thus extending the reach of calculations with hybrid functionals. The bottleneck represented by the sum over bands and by the FFT in Eq. (1) is however still present: ACE just reduces the number of such expensive calculations, but doesn’t eliminate them. In order to achieve a real breakthrough, one has to get rid of delocalized bands and FFT’s, moving to a representation of the electronic structure in terms of localized orbitals. Work along this line using the selected column density matrix localization scheme Anil et al. (2015, 2017) is ongoing. In the next section we describe a different approach, implemented in the CP code, based on maximally localized Wannier functions (MLWF).

The CP code can now perform highly efficient hybrid-functional ab initio MD using MLWFs Marzari and Vanderbilt (1997) $\{\overline{\varphi}_{i}\}$ to represent the occupied space, instead of the canonical KS orbitals $\{\psi_{i}\}$, which are typically delocalized over the entire simulation cell. The MLWF localization procedure can be written as a unitary transformation, $\overline{\varphi}_{i}({\bf r})=\sum_{j}U_{ij}\psi_{j}({\bf r})$, where $U_{ij}$ is computed at each MD time step by minimizing the total spread of the orbitals via a second-order damped dynamics scheme, starting with the converged $U_{ij}$ from the previous time step as initial guesses Sharma et al. (2003).

The natural sparsity of the exchange interaction provided by a localized representation of the occupied orbitals (at least in systems with a finite band gap) is efficiently exploited during the evaluation of exact-exchange based applications (e.g., hybrid DFT functionals). This is accomplished by computing each of the required pair-exchange potentials $\overline{v}_{ij}({\bf r})$ (corresponding to a given localized pair-density $\overline{\rho}_{ij}({\bf r})$) through the numerical solution of the Poisson equation:

using finite differences on the real-space grid. Discretizing the Laplacian operator ($\nabla^{2}$) using a 19-point central-difference stencil (with an associated $\mathcal{O}(h^{6})$ accuracy in the grid spacing $h$), the resulting sparse linear system of equations is solved using the conjugate-gradient technique subject to the boundary conditions imposed by a multipolar expansion of $\overline{v}_{ij}({\bf r})$:

in which the $Q_{lm}$ are the multipoles describing $\overline{\rho}_{ij}({\bf r})$ Wu et al. (2009); DiStasio Jr. et al. (2014); Ko et al. .

Since $\overline{v}_{ij}({\bf r})$ only needs to be evaluated for overlapping pairs of MLWFs, the number of Poisson equations that need to be solved is substantially decreased from ${\cal O}(N_{b}^{2})$ to ${\cal O}(N_{b})$. In addition, $\overline{v}_{ij}({\bf r})$ only needs to be solved on a subset of the real-space grid (that is in general of fixed size) that encompasses the overlap between a given pair of MLWFs. This further reduces the overall computational effort required to evaluate exact-exchange related quantities and results in a linear-scaling (${\cal O}(N_{b})$) algorithm. As such, this framework for performing exact-exchange calculations is most efficient for non-metallic systems (i.e., systems with a finite band gap) in which the occupied KS orbitals can be efficiently localized.

The MLWF representation not only yields the exact-exchange energy $E_{\rm xx}$,

at a significantly reduced computational cost, but it also provides an amenable way of computing the exact-exchange contributions to the (MLWF) wavefunction forces, $\overline{D}_{xx}^{i}({\bf r})=e^{2}\sum_{j}\overline{v}_{ij}({\bf r})\overline{\varphi}_{j}({\bf r})$, which serve as the central quantities in Car-Parrinello MD simulations Car and Parrinello (1985). Moreover, the exact-exchange contributions to the stress tensor are readily available, thereby providing a general code base which enables hybrid DFT based simulations in the NVE, NVT, and NPT ensembles for simulation cells of any shape Ko et al. . We note in passing that applications of the current implementation of this MLWF-based exact-exchange algorithm are limited to $\Gamma$–point calculations employing norm-conserving pseudo-potentials.

The MLWF-based exact-exchange algorithm in CP employs a hybrid MPI/OpenMP parallelization strategy that has been extensively optimized for use on large-scale massively-parallel (super-) computer architectures. The required set of Poisson equations—each one treated as an independent task—are distributed across a large number of MPI ranks/processes using a task distribution scheme designed to minimize the communication and to balance computational workload. Performance profiling demonstrates excellent scaling up to 30,720 cores (for the $\alpha$-glycine molecular crystal, see Fig. 1) and up to 65,536 cores (for (H2O)256, see Ref. DiStasio Jr. et al., 2014) on Mira (BG/Q) with extremely promising efficiency. In fact, this algorithm has already been successfully applied to the study of long-time MD simulations of large-scale condensed-phase systems such as (H2O)128 DiStasio Jr. et al. (2014); Santra et al. (2015). For more details on the performance and implementation of this exact-exchange algorithm, we refer the reader to Ref. Ko et al., .

II.1.2 Dispersion interactions

Dispersion, or van der Waals, interactions arise from dynamical correlations among charge fluctuations occurring in widely separated regions of space. The resulting attraction is a non-local correlation effect that cannot be reliably captured by any local (such as local density approximation, LDA) or semi-local (generalized gradient approximation, GGA) functional of the electron density French et al. (2010). Such interactions can be either accounted for by a truly non-local exchange-correlation (XC) functional, or modeled by effective interactions amongst atoms, whose parameters are either computed from first principles or estimated semi-empirically. In Quantum ESPRESSO both approaches are implemented. Non-local XC functionals are activated by selecting them in the input_dft variable, while explicit interactions are turned on with the vdw_corr option. From the latter group, DFT-D2 Grimme (2006), Tkatchenko-Scheffler Tkatchenko and Scheffler (2009), and exchange-hole dipole moment models Becke and Johnson (2007); Johnson (2017) are currently implemented (DFT-D3 Grimme et al. (2010) and the many-body dispersion (MBD) Tkatchenko et al. (2012); Ambrosetti et al. (2014); Blood-Forsythe et al. (2016) approaches are already available in a development version).

A fully non-local correlation functional able to account for van der Waals interactions for general geometries was first developed in 2004 and named vdW-DF Dion et al. (2004). Its development is firmly rooted in many-body theory, where the so-called adiabatic connection fluctuation-dissipation theorem (ACFD) Langreth and Perdew (1977) provides a formally exact expression for the XC energy through a coupling constant integration over the response function—see Sec. II.1.4. A detailed review of the vdW-DF formalism is provided in Ref. Berland et al., 2015. The overall XC energy given by the ACFD theorem—as a functional of the electron density $n$—is then split in vdW-DF into a GGA-type XC part $E_{\rm xc}^{0}[n]$ and a truly non-local correlation part $E_{\rm c}^{\rm nl}[n]$, i.e.

where the non-local part is responsible for the van der Waals forces. Through a second-order expansion in the plasmon-response expression used to approximate the response function, the non-local part turns into a computationally tractable form involving a universal kernel $\Phi({\bf r},{\bf r}^{\prime})$,

The kernel $\Phi({\bf r},{\bf r}^{\prime})$ depends on ${\bf r}$ and ${\bf r}^{\prime}$ only through $q_{0}({\bf r})|{\bf r}-{\bf r}^{\prime}|$ and $q_{0}({\bf r}^{\prime})|{\bf r}-{\bf r}^{\prime}|$, where $q_{0}({\bf r})$ is a function of $n({\bf r})$ and $\nabla n({\bf r})$. As such, the kernel can be pre-calculated, tabulated, and stored in some external file. To make the scheme self-consistent, the XC potential $V_{\rm c}^{\rm nl}({\bf r})=\delta E_{\rm c}^{\rm nl}[n]/\delta n({\bf r})$ also needs to be computed Thonhauser et al. (2007). The evaluation of $E_{\rm c}^{\rm nl}[n]$ in Eq. (8) is computationally expensive. In addition, the evaluation of the corresponding potential $V_{\rm c}^{\rm nl}({\bf r})$ requires one spatial integral for each point ${\bf r}$. A significant speedup can be achieved by writing the kernel in terms of splines Román-Pérez and Soler (2009)

where $q_{\alpha}$ are fixed values and $p_{\alpha}$ are cubic splines. Equation (8) then becomes a convolution that can be simplified to

Here \theta_{\alpha}({\bf r})=n({\bf r})p_{\alpha}\big{(}q_{0}({\bf r})\big{)} and $\theta_{\alpha}({\bf k})$ is its Fourier transform. Accordingly, $\Phi_{\alpha\beta}(k)$ is the Fourier transform of the original kernel $\Phi_{\alpha\beta}(r)=\Phi(q_{\alpha},q_{\beta},|{\bf r}-{\bf r}^{\prime}|)$. Thus, two spatial integrals are replaced by one integral over Fourier transformed quantities, resulting in a considerable speedup. This approach also provides a convenient evaluation for $V_{\rm c}^{\rm nl}({\bf r})$.

The vdW-DF functional was implemented in Quantum ESPRESSO version 4.3, following Eq. (10). As a result, in large systems, compute times in vdW-DF calculations are only insignificantly longer than for standard GGA functionals. The implementation uses a tabulation of the Fourier transformed kernel $\Phi_{\alpha\beta}(k)$ from Eq. (10) that is computed by an auxiliary code, generate_vdW_kernel_table.x, and stored in the external file vdW_kernel_table. The file then has to be placed either in the directory where the calculation is run or in the directory where the corresponding pseudopotentials reside. The formalism for vdW-DF stress was derived and implemented in Ref. Sabatini et al., 2012. The proper spin extension of vdW-DF, termed svdW-DF Thonhauser et al. (2015), was implemented in Quantum ESPRESSO version 5.2.1.

Although the ACFD theorem provides guidelines for the total XC functional in Eq. (7), in practice $E_{\rm xc}^{0}[n]$ is approximated by simple GGA-type functional forms. This has been used to improve vdW-DF—and correct the often too large binding separations found in its original form—by optimizing the exchange contribution to $E_{\rm xc}^{0}[n]$. The naming convention for the resulting variants is that the extension should describe the exchange functional used. In this context, the functionals vdW-DF-C09 Cooper (2010), vdW-DF-obk8 Klimeš et al. (2010), vdW-DF-ob86 Klimeš et al. (2011), and vdW-DF-cx Berland and Hyldgaard (2014) have been developed and implemented in Quantum ESPRESSO. While all of these variants use the same kernel to evaluate $E_{\rm c}^{\rm nl}[n]$, advances have also been made in slightly adjusting the kernel form, which is referred to and implemented as vdW-DF2 Lee et al. (2010). A corresponding variant, i.e., vdW-DF2-b86r Hamada and Otani (2010), is also implemented. Note that vdW-DF2 uses the same kernel file as vdW-DF.

The functional VV10 Vydrov and Van Voorhis (2010) is related to vdW-DF, but adheres to fewer exact constraints and follows a very different design philosophy. It is implemented in Quantum ESPRESSO in a form called rVV10 Sabatini et al. (2013) and uses a different kernel and kernel file that can be generated by running the auxiliary code generate_rVV10_kernel_table.x.

An alternative approach to accounting for dispersion forces is to add to the XC energy $E^{0}_{\rm xc}$ a dispersion energy, $E_{\rm disp}$, written as a damped asymptotic pairwise expression:

where $I$ and $J$ run over atoms, $R_{IJ}=|\mathbf{R}_{I}-\mathbf{R}_{J}|$ is the interatomic distance between atoms $I$ and $J$, and $f_{n}(R)$ is a suitable damping function. The interatomic dispersion coefficients $C^{(n)}_{IJ}$ can be derived from fits, as in DFT-D2 Grimme (2006), or calculated non-empirically, as in the Tkatchenko-Scheffler (TS-vdW) Tkatchenko and Scheffler (2009) and exchange-hole dipole moment (XDM) models Becke and Johnson (2007); Johnson (2017).

In XDM, the $C^{(n)}_{IJ}$ coefficients are calculated assuming that dispersion interactions arise from the electrostatic attraction between the electron-plus-exchange-hole distributions on different atoms Becke and Johnson (2007); Johnson (2017). In this way, XDM retains the simplicity of a pairwise dispersion correction, like in DFT-D2, but derives the $C^{(n)}_{IJ}$ coefficients from the electronic properties of the system under study. The damping functions $f_{n}$ in Eq. (11) suppress the dispersion interaction at short distances, and serve the purpose of making the link between the short-range correlation (provided by the XC functional) and the long-range dispersion energy, as well as mitigating erroneous behavior from the exchange functional in the representation of intermolecular repulsion Johnson (2017). The damping functions contain two adjustable parameters, available online sch for a number of popular density functionals. Although any functional for which damping parameters are available can be used, the functionals showing best performance when combined with XDM appear to be B86bPBE Becke (1986); Perdew et al. (1996) and PW86PBE Perdew and Yue (1986); Perdew et al. (1996), thanks to their accurate modeling of Pauli repulsion Johnson (2017). Both functionals have been implemented in Quantum ESPRESSO since version 5.0.

In the canonical XDM implementation, recently included in Quantum ESPRESSO and described in detail elsewhere Otero-de-la-Roza and Johnson (2012), the dispersion coefficients are calculated from the electron density, its derivatives, and the kinetic energy density, and assigned to the different atoms in the system using a Hirshfeld atomic partition scheme Hirshfeld (1977). This means that XDM is effectively a meta-GGA functional of the dispersion energy whose evaluation cost is small relative to the rest of the self-consistent calculation. Despite the conceptual and computational simplicity of XDM, and because the dispersion coefficients depend upon the atomic environment in a physically meaningful way, the XDM dispersion correction offers good performance in the calculation of diverse properties, such as lattice energies, crystal geometries, and surface adsorption energies. XDM is especially good for modeling organic crystals and organic/inorganic interfaces. For a recent review, see Ref. 13.

The XDM dispersion calculation is turned on by specifying vdw_corr=’xdm’ and optionally selecting appropriate damping function parameters (with the xdm_a1 and xdm_a2 keywords). Because the reconstructed all-electron densities are required during self-consistency, XDM can be used only in combination with a PAW approach. The XDM contribution to forces and stress is not entirely consistent with the energies because the current implementation neglects the change in the dispersion coefficients. Work is ongoing to remove this limitation, as well as to make XDM available for Car-Parrinello MD, in future Quantum ESPRESSO releases.

In the TS-vdW approach (vdw_corr=’ts-vdw’), all vdW parameters (which include the atomic dipole polarizabilities, $\alpha_{I}$, vdW radii, $R^{0}_{I}$, and interatomic $C^{(6)}_{IJ}$ dispersion coefficients) are functionals of the electron density and computed using the Hirshfeld partitioning scheme Hirshfeld (1977) to account for the unique chemical environment surrounding each atom. This approach is firmly based on a fluctuating quantum harmonic oscillator (QHO) model and results in highly accurate $C^{(6)}_{IJ}$ coefficients with an associated error of approximately 5.5% Tkatchenko and Scheffler (2009). The TS-vdW approach requires a single empirical range-separation parameter based on the underlying XC functional and is recommended in conjunction with non-empirical DFT functionals such as PBE and PBE0. For a recent review of the TS-vdW approach and several other vdW/dispersion corrections, please see Ref. 79.

The implementation of the density-dependent TS-vdW correction in Quantum ESPRESSO is fully self-consistent Ferri et al. (2015) and currently available for use with norm-conserving pseudo-potentials. An efficient linear-scaling implementation of the TS-vdW contribution to the ionic forces and stress tensor allows for Born-Oppenheimer and Car-Parrinello MD simulations at the DFT+TS-vdW level of theory; this approach has already been successfully employed in long-time MD simulations of large-scale condensed-phase systems such as (H2O)128 DiStasio Jr. et al. (2014); Santra et al. (2015). We note in passing that the Quantum ESPRESSO implementation of the TS-vdW correction also includes analytical derivatives of the Hirshfeld weights, thereby completely reflecting the change in all TS-vdW parameters during geometry/cell optimizations and MD simulations.

II.1.3 Hubbard-corrected functionals: DFT+U

Most approximate XC functionals used in modern DFT codes fail quite spectacularly on systems with atoms whose ground-state electronic structure features partially occupied, strongly localized orbitals (typically of the $d$ or $f$ kind), that suffer from strong self-interaction effects and a poor description of electronic correlations. In these circumstances, DFT+U is often, although not always, an efficient remedy. This method is based on the addition to the DFT energy functional $E_{\text{DFT}}$ of a correction $E_{U}$, shaped on a Hubbard model Hamiltonian: $E_{\text{DFT}+U}=E_{\text{DFT}}+E_{U}$. The original implementation in Quantum ESPRESSO, extensively described in Refs. Cococcioni and de Gironcoli, 2005; Himmetoglu et al., 2014, is based on the simplest rotationally invariant formulation of $E_{U}$, due to Dudarev and coworkers Dudarev et al. (1998):

$|\psi^{\sigma}_{{\bf k}\nu}\rangle$ is the valence electronic wave function for state ${\bf k}\nu$ of spin $\sigma$, $f^{\sigma}_{{\bf k}\nu}$ the corresponding occupation number, $|\phi^{I}_{m}\rangle$ is the chosen Hubbard manifold of atomic orbitals, centered on atomic site $I$, that may be orthogonalized or not. The presence of the Hubbard functional results in extra terms in energy derivatives such as forces, stresses, elastic constants, or force-constant (dynamical) matrices. For instance, the additional term in forces is

where $R_{I\alpha}$ is the $\alpha$ component of position for atom $I$ in the unit cell,

As a correction to the total energy, the Hubbard functional naturally contributes an extra term to the total potential that enters the KS equations. An alternative formulation Sclauzero and Dal Corso (2013) of the DFT+U method, recently introduced and implemented in Quantum ESPRESSO for transport calculations, eliminates the need of extra terms in the potential by incorporating the Hubbard correction directly into the (PAW) pseudopotentials through a renormalization of the coefficients of their non-local terms.

A simple extension to the Dudarev functional, DFT+U+J0, was proposed in Ref. 15 and used to capture the insulating ground state of CuO. In CuO the localization of holes on the $d$ states of Cu and the consequent on-set of a magnetic ground state can only be stabilized against a competing tendency to hybridize with oxygen $p$ states when on-site exchange interactions are precisely accounted for. A simplified functional, depending upon the on-site (screened) Coulomb interaction $U$ and the Hund’s coupling $J$, can be obtained from the full second-quantization formulation of the electronic interaction potential by keeping only on-site terms that describe the interaction between up to two orbitals and by approximating on-site effective interactions with the (orbital-independent) atomic averages of Coulomb and exchange terms:

The on-site exchange coupling $J^{I}$ not only reduces the effective Coulomb repulsion between like-spin electrons as in the simpler Dudarev functional (first term of the right-hand side), but also contributes a second term that acts as an extra penalty for the simultaneous presence of anti-aligned spins on the same atomic site and further stabilizes ferromagnetic ground states.

The fully rotationally invariant scheme of Liechtenstein et al. Liechtenstein et al. (1995), generalized to non-collinear magnetism and two-component spinor wave-functions, is also implemented in the current version of Quantum ESPRESSO. The corrective energy term for each correlated atom can be quite generally written as:

where the average is taken over the DFT Slater determinant, $U_{\alpha\beta\gamma\delta}$ are Coulomb integrals, and some set of orthonormal spin-space atomic functions, $\{\alpha\}$, is used to calculate the occupation matrix, $n_{\alpha\beta}$, via Eq. (13). These basis functions could be spinor wave functions of total angular momentum $j=l\pm 1/2$, originated from spherical harmonics of orbital momentum $l$, which is a natural choice in the presence of spin-orbit coupling. Another choice, adopted in our implementation, is to use the standard basis of separable atomic functions, $R_{l}(r)Y_{lm}(\theta,\phi)\chi(\sigma)$, where $\chi(\sigma)$ are spin up/down projectors and the radial function, $R_{l}(r)$, is an eigenfunction of the pseudo-atom. In the presence of spin-orbit coupling, the radial function is constructed by averaging between the two radial functions $R_{l\pm 1/2}$. These radial functions are read from the file containing the pseudopotential, in this case a fully-relativistic one. In this conventional basis, the corrective functional takes the form:

where $\{ijkl\}$ run over azimuthal quantum number $m$. The second term contains a spin-flip contribution if $\sigma^{\prime}\neq\sigma$. For collinear magnetism, when $n_{ij}^{\sigma\sigma^{\prime}}=\delta_{\sigma\sigma^{\prime}}n_{ij}^{\sigma}$, the present formulation reduces to the scheme Liechtenstein et al. (1995) of Liechtenstein et al. All Coulomb integrals, $U_{ijkl}$, can be parameterized by few input parameters such as $U$ ($s$-shell); $U$ and $J$ ($p$-shell); $U,J$ and $B$ ($d$-shell); $U,J,E_{2}$, and $E_{3}$ ($f$-shell), and so on. We note that if all parameters but $U$ are set to zero, the Dudarev functional is recovered.

The Hubbard corrective functional $E_{U}$ depends linearly upon the effective on-site interactions, $U^{I}$. Therefore, using a proper value for these interaction parameters is crucial to obtain quantitatively reliable results from DFT+U calculations. The Quantum ESPRESSO implementation of DFT+U has also been the basis to develop a method for the calculation of $U$ Cococcioni and de Gironcoli (2005), based on linear-response theory. This method is completely ab initio and provides values of the effective interactions that are consistent with the system and with the ground state that the Hubbard functional aims at correcting. A comparative analysis of this method with other approaches proposed in the literature to compute the Hubbard interactions has been initiated in Ref. 82 and will be further refined in forthcoming publications by the same authors.

Within linear-response theory, the Hubbard interactions are the elements of an effective interaction matrix, computed as the difference between bare and screened inverse susceptibilities Cococcioni and de Gironcoli (2005):

In Eq. (19) the susceptibilities $\chi$ and $\chi_{0}$ measure the response of atomic occupations to shifts in the potential acting on the states of single atoms in the system. In particular, $\chi$ is defined as $\chi_{IJ}=\sum_{m\sigma}\left(dn^{I\sigma}_{mm}/d\alpha^{J}\right)$ and is evaluated at self consistency, while $\chi_{0}$ has a similar definition but is computed before the self-consistent re-adjustment of the Hartree and XC potentials. In the current implementation these susceptibilities are computed from a series of self-consistent total energy calculations (varying the strength $\alpha$ of the perturbing potential over a range of values) performed on supercells of sufficient size for the perturbations to be isolated from their periodic replicas. While easy to implement, this approach is quite cumbersome to use, requiring multiple calculations, expensive convergence tests of the resulting parameters and complex post-processing tools.

These difficulties can be overcome by using density-functional perturbation theory (DFpT) to automatize the calculation of the Hubbard parameters. The basic idea is to recast the entries of the susceptibility matrices into sums over the BZ:

where $I\equiv(l,s)$ and $J\equiv(l^{\prime},s^{\prime})$, $l$ and $l^{\prime}$ label unit cells, $s$ and $s^{\prime}$ label atoms in the unit cell, $\mathbf{R}_{l}$ and $\mathbf{R}_{l^{\prime}}$ are Bravais lattice vectors, and $\Delta_{\bf q}^{s^{\prime}}n^{s\,\sigma}_{mm^{\prime}}$ represent the (lattice-periodic) response of atomic occupations to monochromatic perturbations constructed by modulating the shift to the potential of all the periodic replica of a given atom by a wave-vector ${\bf q}$. This quantity is evaluated within DFpT (see Sec. II.2), using linear-response routines contained in LR_Modules (see Sec. III.4.3). This approach eliminates the need for supercell calculations in periodic systems (along with the cubic scaling of their computational cost) and automatizes complex post-processing operations needed to extract $U$ from the output of calculations. The use of DFpT also offers the perspective to directly evaluate inverse susceptibilities, thus avoiding the matrix inversions of Eq. (19), and to calculate the Hubbard parameters for closed-shell systems, a notorious problem for schemes based on perturbations to the potential. Full details about this implementation will be provided in a forthcoming publication Timrov et al. and the corresponding code will be made available in one of the next Quantum ESPRESSO releases.

II.1.4 Adiabatic-connection fluctuation-dissipation theory

In the quest for better approximations for the unknown XC energy functional in KS-DFT, the approach based on the adiabatic connection fluctuation-dissipation (ACFD) theorem Langreth and Perdew (1977) has received considerable interest in recent years. This is largely due to some attractive features: (i) a formally exact expression for the XC energy in term of density linear response functions can be derived providing a promising way for a systematic improvement of the XC functional; (ii) the method treats the exchange energy exactly, thus canceling out the spurious self-interaction error present in the Hartree energy; (iii) the correlation energy is fully non local and automatically includes long-range van der Waals interactions (see Sec. II.1.2).

Within the ACFD framework a formally exact expression for the XC energy $E_{\rm xc}$ of an electronic system can be derived:

where $\hbar=h/2\pi$ and $h$ is the Planck constant, $\chi_{\lambda}({\bf r},{\bf r}^{\prime};iu)$ is the density response function at imaginary frequency $iu$ of a system whose electrons interact via a scaled Coulomb interaction, i.e., $\lambda e^{2}/|{\bf r}-{\bf r}^{\prime}|$, and are subject to a local potential such that the electronic density $n({\bf r})$ is independent of $\lambda$, and is thus equal to the ground-state density of the fully interacting system $(\lambda=1)$. The XC energy, Eq. (21), can be further separated into a KS exact-exchange energy $E_{\rm xx}$, Eq. (6), and a correlation energy $E_{\rm c}$. The former is routinely evaluated as in any hybrid functional calculation (see Sec. II.1.1). Using a matrix notation, the latter can be expressed in a compact form in terms of the Coulomb interaction, $v_{c}=e^{2}/|{\bf r}-{\bf r}^{\prime}|$, and of the density response functions:

For $\lambda>0$, $\chi_{\lambda}$ can be related to the noninteracting density response function $\chi_{0}$ via a Dyson equation obtained from TDDFT:

The exact expression of the XC kernel $f_{\rm xc}$ is unknown, and in practical applications one needs to approximate it. In the ACFDT package, the random phase approximation (RPA), obtained by setting $f_{\rm xc}^{\lambda}=0$, and the RPA plus exact-exchange kernel (RPAx), obtained by setting $f_{\rm xc}^{\lambda}=\lambda f_{\rm x}$, are implemented. The evaluation of the RPA and RPAx correlation energies is based on an eigenvalue decomposition of the non-interacting response functions and of its first-order correction in the limit of vanishing electron-electron interaction Wilson et al. (2008); Nguyen and de Gironcoli (2009); Colonna et al. (2014). Since only a small number of these eigenvalues are relevant for the calculation of the correlation energy, an efficient iterative scheme can be used to compute the low-lying modes of the RPA/RPAx density response functions.

The basic operation required for the eigenvalue decomposition is a number of loosely coupled DFpT calculations for different imaginary frequencies and trial potentials. Although the global scaling of the iterative approach is the same as for implementations based on the evaluation of the full response matrices ($N^{4}$), the number of operation involved is 100 to 1000 times smaller Nguyen and de Gironcoli (2009), thus largely reducing the global scaling pre-factor. Moreover, the calculation can be parallelized very efficiently by distributing different trial potentials on different processors or groups of processors.

In addition, the local EXX and RPA-correlation potentials can be computed through an optimized effective potential (OEP) scheme fully compatible with the eigenvalue decomposition strategy adopted for the evaluation of the EXX/RPA energy. Iterating the energy and the OEP calculations and using an effective mixing scheme to update the KS potential, a self-consistent minimization of the EXX/RPA functional can be achieved Nguyen et al. (2014).

II.2 Linear response and excited states without virtual orbitals

One of the key features of modern DFT implementations is that they do not require the calculation of virtual (unoccupied) orbitals. This idea, first pioneered by Car and Parrinello in their landmark 1985 paper Car and Parrinello (1985) and later adopted by many groups world-wide, found its way in the computation of excited-state properties with the advent of density-functional perturbation theory (DFpT) Baroni et al. (1987); Giannozzi et al. (1991); Gonze (1995); Baroni et al. (2001). DFpT is designed to deal with static perturbations and its use is therefore restricted to those excitations that can be described in the Born-Oppenheimer approximation, such as lattice vibrations. The main idea underlying DFpT is to represent the linear response of KS orbitals to an external perturbation as generic orbitals satisfying an orthogonality constraint with respect to the occupied-state manifold and a self-consistent Sternheimer equation Sternheimer (1954); Mahan (1980), rather than as linear combinations of virtual orbitals (which would require the computation of all, or a large number, of them).

Substantial progress has been made over the past decade, allowing one to extend DFpT to the dynamical regime, and thus simulate sizable portions of the optical and loss spectra of complex molecular and extended systems, without making any explicit reference to their virtual states. Although the Sternheimer approach can be easily extended to time-dependent perturbations Schwartz and Tiemann (1959); Zernik (1964); Baroni and Quattropani (1985), its use is hampered in practice by the fact that a different Sternheimer equation has to be solved for each different value of the frequency of the perturbation. When the perturbation acting on the system vanishes, the frequency-dependent Sternheimer equation becomes a non-Hermitian eigenvalue equation, whose eigenvalues are the excitation energies of the system. In the TDDFT community, this equation is known as the Casida equation Casida (1996); Jamorski et al. (1996), which is the immediate translation to the DFT parlance of the time-dependent Hartree-Fock equation McLachlan and Ball (1964). This approach to excited states is optimal in those cases where one is interested in a few excitations only, but can hardly be extended to continuous spectra, such as those arising in extended systems or above the ionization threshold of even finite ones. In those cases where extended portions of a continuous spectrum is required, a new method has been developed, based on the Lanczos (bi-) orthogonalization algorithm, and dubbed the Liouville-Lanczos approach to time-dependent density-functional perturbation theory (TDDFpT). This method allows one to reuse intermediate products of an iterative process, essentially identical to that used for static perturbations, to build dynamical response functions from which spectral properties can be computed for a whole wide spectral range at once Walker et al. (2006); Rocca et al. (2008). A similar approach to linear optical spectroscopy was proposed later, based on the multi-shift conjugate gradient algorithm Hübener and Giustino (2014), instead of Lanczos. This powerful idea has been generalized to the solution of the Bethe-Salpeter equation, which is formally very similar to the eigenvalue equations arising in TDDFpT Rocca et al. (2010, 2012); Marsili et al. (2017), and to the computation of the polarization propagator and self-energy operator appearing in the $GW$ equations Umari et al. (2009, 2010); Govoni and Galli (2015). It is presently exploited in several components of the Quantum ESPRESSO distribution, as well as in other advanced implementations of many-body perturbation theory Govoni and Galli (2015).

The computation of vibrational properties in extended systems is one of the traditional fields of application of DFpT, as thoroughly described, e.g., in Ref. Baroni et al., 2001. The latest releases of Quantum ESPRESSO feature the linear-response implementation of several new functionals in the van der Waals and DFT+U families. Explicit expressions of the XC kernel, implementation details, and a thorough benchmark are reported in Ref. Sabatini et al., 2016 for the first case. As for the latter, DFpT+U has been implemented for both the Dudarev Dudarev et al. (1998) and DFT+U+J0 functionals Himmetoglu et al. (2011), allowing one to account for electronic localization effects acting selectively on specific phonon modes at arbitrary wave-vectors, thus substantially improving the description of the vibrational properties of strongly correlated systems with respect to “standard” LDA/GGA functionals. The current implementation allows for both norm-conserving and ultrasoft pseudopotentials, insulators and metals alike, also including the spin-polarized case. The implementation of DFpT+U requires two main additional ingredients with respect to standard DFpT Floris et al. (2011). First, the dynamical matrix contains an additional term coming from the second derivative of the Hubbard term $E_{U}$ with respect to the atomic positions (denoted $\lambda$ or $\mu$), namely:

where the notations are the same as in Eq. (12). The symbols $\partial$ and $\Delta$ indicate, respectively, a bare derivative (leaving the KS wavefunctions unperturbed) and a total derivative (including also linear-response contributions). Second, in order to obtain a consistent electronic density response to the atomic displacements from the DFT+U ground state, the perturbed KS potential $\Delta V_{SCF}$ in the Sternheimer equation is augmented with the Hubbard perturbed potential $\Delta^{\lambda}V_{U}$:

where the notations are the same as in Eq. (13). The unperturbed Hamiltonian in the Sternheimer equation is the DFT+U Hamiltonian (including the Hubbard potential $V_{U}$). More implementation details will be given in a forthcoming publication Floris et al. .

Applications of DFpT+U include the calculation of the vibrational spectra of transition-metal monoxides like MnO and NiO Floris et al. (2011), investigations of properties of materials of geophysical interest like goethite Blanchard et al. (2014a, b), iron-bearing Shukla et al. (2015); Shukla and Wentzcovitch (2016) and aluminum-bearing bridgmanite Shukla et al. (2016). These results feature a significantly better agreement with experiment of the predictions of various lattice-dynamical properties, including the LO-TO and magnetically-induced TO splittings, as compared with standard LDA/GGA calculations.

II.2.2 Dynamic perturbations: optical, electron energy loss, and magnetic spectroscopies

where $\hat{Q}$ is the projector on the unoccupied-state manifold. The response orbitals $x_{{\bf k}\nu}$ and $y_{{\bf k}\nu}$ can be collected in so-called batches $X=\{x_{{\bf k}\nu}\}$ and $Y=\{y_{{\bf k}\nu}\}$, which uniquely determine the response density matrix. In a similar way, the perturbing potential $\hat{V}^{\prime}$ can be represented by the batch $Z=\{z_{{\bf k}\nu}\}=\{\hat{Q}\hat{V}^{\prime}\psi_{{\bf k}\nu}\}$. Using these definitions, the linear-response equations of TDDFpT take the simple form:

where the super-operators $\hat{D}$ and $\hat{K}$, which enter the definition of the Liouvillian super-operator, $\hat{\mathcal{L}}$, are defined in terms of the unperturbed Hamiltonian and of the perturbed Hartree-plus-XC potential Walker et al. (2006); Rocca et al. (2008); Malcioğlu et al. (2011); Ge et al. (2014); Baroni and Gebauer (390). This implies that a Liouvillian build costs roughly twice as much as a single iteration in time-independent DFpT. It is important to note that $\hat{D}$ and $\hat{K}$, and therefore $\hat{\mathcal{L}}$, do not depend on the frequency $\omega$. For this reason, when in Eq. (28) the vector on the right-hand side, $(0,Z)^{\top}$, is set to zero, a linear eigenvalue equation is obtained (Casida’s equation).

The quantum Liouville equation (28) can be seen as the equation for the response density matrix operator $\hat{\rho}^{\prime}(\omega)$, namely $(\hbar\omega-\hat{\mathcal{L}})\cdot\hat{\rho}^{\prime}(\omega)=[\hat{V}^{\prime},\hat{\rho}^{\circ}]$, where $[\cdot,\cdot]$ is the commutator and $\hat{\rho}^{\circ}$ is the ground-state density matrix operator. With this at hand, we can define a generalized susceptibility $\chi_{AV}(\omega)$, which characterizes the response of an arbitrary one-electron Hermitian operator $\hat{A}$ to the external perturbation $\hat{V}^{\prime}$ as

where $\langle\cdot|\cdot\rangle$ denotes a scalar product in operator space. For instance, when both $\hat{A}$ and $\hat{V}^{\prime}$ are one of the three Cartesian components of the dipole (position) operator, Eq. (29) gives the dipole polarizability of the system, describing optical absorption spectroscopy Walker et al. (2006); Rocca et al. (2008); setting $\hat{A}$ and $\hat{V}^{\prime}$ to one of the space Fourier components of the electron charge-density operator would correspond to the simulation of electron energy loss or inelastic X-ray scattering spectroscopies, giving access to plasmon and exciton excitations in extended systems Timrov et al. (2013, 2015b); two different Cartesian components of the Fourier transform of the spin polarization would give access to spectroscopies probing magnetic excitations (e.g. inelastic neutron scattering or spin-polarized electron energy loss) Gorni et al. , and so on. When dealing with macroscopic electric fields, the dipole operator in periodic boundary conditions is treated using the standard DFpT prescription, as explained in Refs. Baroni and Resta (1986); Tobik and Dal Corso (2004).

The Quantum ESPRESSO distribution contains several codes to solve the Casida’s equation or to directly compute generalized susceptibilities according to Eq. (29) and by solving Eq. (28) using different approaches for different pairs of $\hat{A}$/$\hat{V}^{\prime}$, corresponding to different spectroscopies. In particular, Eq. (28) can be solved iteratively using the Lanczos recursion algorithm, which allows one to avoid computationally expensive inversion of the Liouvillian. The basic principle of how matrix elements of the resolvent of an operator can be calculated using a Lanczos recursion chain has been worked out by Haydock, Heine, and Kelly Haydock et al. (1972, 1975) for the case of Hermitian operators and diagonal matrix elements. The quantity of interst can be written as

where $\beta_{n+1}$ is given by the condition $\langle q_{n+1}|q_{n+1}\rangle=1$. The vectors $|q_{n}\rangle$ created by this recursive chain are orthonormal. Furthermore, the operator $\hat{L}$, written in the basis of these vectors, is tridiagonal. If one limits the chain to the $M$ first vectors $|q_{0}\rangle,|q_{1}\rangle,\cdots,|q_{M}\rangle$, then the resulting representation of $\hat{L}$ is a $M\times M$ square matrix $T_{M}$ which reads

Using such a truncated representation of $\hat{L}$, the resolvent in Eq. (30) can be approximated as

Thanks to the tridiagonal form of $T_{M}$, the approximate resolvent can finally be written as a continued fraction

Note that performing the recursion (31) – (34) is the computational bottleneck of this algorithm, while evaluating the continued fraction in Eq. (37) is very fast. The recursion being independent of the frequency $\omega$, a single recursion chain yields information about any desired number of frequencies, at negligible additional computational cost. It is also important to note that at any stage of the recursion chain, only three vectors need to be kept in memory, namely $|q_{n-1}\rangle$, $|q_{n}\rangle$, and $|q_{n-1}\rangle$. This is a considerable advantage with respect to the direct calculation of $N$ eigenvectors of $\hat{L}$ where all $N$ vectors need to be kept in memory in order to enforce orthogonality.

The Liouvillian $\hat{\cal L}$ in Eq. (28) is not a Hermitian operator. For this reason, the Lanczos algorithm presented above cannot be directly applied to the calculation of the generalized susceptibility (29). There are two distinct algorithms that can be applied in the non-Hermitian case. The non-Hermitian Lanczos biorthogonalization algorithm Rocca et al. (2008); Malcioğlu et al. (2011) amounts to recursively applying the operator $\hat{\cal L}$ and it Hermitian conjugate $\hat{\cal L}^{{\dagger}}$ to two previous Lanczos vectors $|v_{n}\rangle$ and $|w_{n}\rangle$. In this way, a pair of bi-orthogonal basis sets is created. The operator $\hat{\cal L}$ can then be represented in this basis as a tridiagonal matrix, similarly to the Hermitian case, Eq. (35). The Liouvillian $\hat{\cal L}$ of TDDFT belongs to a special class of non-Hermitian operators which are called pseudo-Hermitian Grüning et al. (2011); Ge et al. (2014). For such operators, there exists a second recursive algorithm to compute the resolvent – pseudo-Hermitian Lanczos algorithm – which is numerically more stable and requires only half the numbers of operations per Lanczos step Grüning et al. (2011); Ge et al. (2014). Both algorithms have been implemented in Quantum ESPRESSO. Because of its speed and numerical stability, the use of the pseudo-Hermitian method is recommended.

This methodology has also been extended—presently only in the case of absorption spectroscopy—to employ hybrid functionals Rocca et al. (2010, 2012); Ge et al. (2014) (see Sec. II.1.1). In this case the calculation requires the evaluation of the linear response of the non-local Fock potential, which is readily available from the response density matrix, represented by the batches of response orbitals. The corresponding hybrid-functional Liouvillian features additional terms with respect to the definition in Eq. (28), but presents a similar structure and similar mathematical properties. Accordingly, semi-local and hybrid-functional TDDFpT employ the same numerical algorithms in practical calculations.

The turbo_lanczos.x Malcioğlu et al. (2011); Ge et al. (2014) and turbo_davidson.x Ge et al. (2014) codes are designed to simulate the optical response of molecules and clusters. turbo_lanczos.x computes the dynamical dipole polarizability [see Eq. (29)] of finite systems over extended frequency ranges without ever computing any eigenpairs of the Casida equation. This goal is achieved by applying a recursive non-Hermitian or pseudo-Hermitian Lanczos algorithm. The two flavours of the Lanczos algorithm implemented in turbo_lanczos.x are particularly suited in those cases where one is interested in the spectrum over a wide frequency range comprising a large number of individual excitations. In turbo_davidson.x a Davidson-like algorithm Davidson (1975) is used to selectively compute a few eigenvalues and eigenvectors of ${\hat{\mathcal{L}}}$. This is useful when very few low-lying excited states are needed and/or when the excitation eigenvector is explicitly needed, e.g., to compute ionic forces on excited potential energy surfaces, a feature that will be implemented in one of the forthcoming releases. Both turbo_lanczos.x and turbo_davidson.x are interfaced with the Environ module Andreussi et al. (2012), to simulate the absorption spectra of complex molecules in solution using the self-consistent continuum solvation model Timrov et al. (2015a) (see Sec. II.5.1).

The turbo_eels.x code Timrov et al. (2015b) computes the response of extended systems to an incoming beam of electrons or X rays, aimed at simulating electron energy loss (EEL) or inelastic X-ray scattering (IXS) spectroscopies, sensitive to collective charge-fluctuation excitations, such as plasmons. Similarly to the description of vibrational modes in a lattice by the PHonon package, here the perturbation can be represented as a sum of monochromatic components corresponding to different momenta, ${\bf q}$, and energy transferred from the incoming electrons to electrons of the sample. The quantum Liouville equation (28) in the batch representation can be formulated for individual ${\bf q}$ components of the perturbation, which can be solved independently Timrov et al. (2013). The recursive Lanczos algorithm is used to solve iteratively the quantum Liouville equation, much like in the case of absorption spectroscopy. The entire EEL/IXS spectrum is obtained in an arbitrarily wide energy range (up to the core-loss region) with only one Lanczos chain. Such a numerical algorithm allows one to compute directly the diagonal element of the charge-density susceptibility, see Eq. (29), by avoiding computationally expensive matrix inversions and multiplications characteristic of standard methods based on the solution of the Dyson equation Onida et al. (2002). The current version of turbo_eels.x allows to explicitly account for spin-orbit coupling effects Timrov et al. (2017).

II.2.3 Many-body perturbation theory

Many-body perturbation theory refers to a set of computational methods, based on quantum field theory, that are designed to calculate electronic and optical excitations beyond standard DFT Onida et al. (2002). The most popular among such methods are the $GW$ approximation and the Bethe-Salpeter equation (BSE) approach. The former is intended to calculate accurate quasiparticle excitations, e.g., ionization energies and electron affinities in molecules, band structures in solids, and accurate band gaps in semiconductor and insulators. The latter is employed to study optical excitations by including electron-hole interactions.

In the $GW$ method the XC potential of DFT is corrected using a many-body self-energy consisting of the product of the electron Green’s function $G$ and the screened Coulomb interaction $W$ Hedin (1965); Hybertsen and Louie (1985), which represents the lowest-order term in the diagrammatic expansion of the exact electron self-energy. In the vast majority of $GW$ implementations, the evaluation of $G$ and $W$ requires the calculation of both occupied and unoccupied KS eigenstates. The convergence of the resulting self-energy correction with respect to the number of unoccupied states is rather slow, and in many cases it constitutes the main bottleneck in the calculations. During the past decade there has been a growing interest in alternative techniques which only require the calculation of occupied electronic states, and several computational strategies have been developed Reining et al. (1997); Wilson et al. (2009); Umari et al. (2010); Giustino et al. (2010). The common denominator to all these strategies is that they rely on linear-response DFpT and the Sternheimer equation, as in the PHonon package.

In Quantum ESPRESSO the $GW$ approximation is realized based on a DFpT representation of response and self-energy operators, thus avoiding any explicit reference to unoccupied states. There are two different implementations: the GWL ($GW$-Lanczos) package Umari et al. (2009, 2010) and the SternheimerGW package Schlipf et al. (2017). The former focuses on efficient $GW$ calculations for large systems (including disordered solids, liquids, and interfaces), and also supports the calculations of optical spectra via the Bethe-Salpeter approach Marsili et al. (2017). The latter focuses on high-accuracy calculations of band structures, frequency-dependent self-energies, and quasi-particle spectral functions for crystalline solids. In addition to these, the WEST code Govoni and Galli (2015), not part of the Quantum ESPRESSO distribution, relies on Quantum ESPRESSO as an external library to perform similar tasks and to achieve similar goals.

The GWL package consists of four different codes. The pw4gww.x code reads the KS wave-functions and charge density previously calculated by PWscf and prepares a set of data which are used by code gww.x to perform the actual $GW$ calculation. While pw4gww.x uses the plane-wave representation of orbitals and charges and the same Quantum ESPRESSO environment as all other linear response codes, gww.x does not rely on any specific representation of the orbitals. Its parallelization strategy is based on the distribution of frequencies. Only a few basic routines, such as the MPI drivers, are common with the rest of Quantum ESPRESSO.

GWL supports three different basis sets for representing polarisability operators: i) plane wave-basis set, defined by an energy cutoff; ii) the basis set formed by the most important eigenvectors (i.e., corresponding to the highest eigenvalues) of the actual irreducible polarisability operator at zero frequency calculated through linear response; iii) the basis set formed by the most important eigenvectors of an approximated polarisability operator. The last choice permits the control of the balance between accuracy and dimension of the basis. The $GW$ scheme requires the calculation of products in real space of KS orbitals with vectors of the polarisability basis. These are represented in GWL through dedicated additional basis sets of reduced dimensions.

GWL supports only the $\Gamma-$point sampling of the BZ and considers only real wave-functions. However, ordinary ${\bf k}$-point sampling of the BZ can be used for the long-range part of the (symmetrized) dielectric matrix. These terms are calculated by the head.x code. In this way reliable calculations for extended materials can be performed using quite small simulation cells (with cell edges of the order of 20 Bohr). Self-consistency is implemented in GWL, although limited to the quasi-particle energies; the so-called vertex term, arising in the diagrammatic expansion of the self-energy, is not yet implemented.

Usually ordinary $GW$ calculations for transition elements require the explicit inclusion of semicore orbitals in the valence manifold, resulting in a significantly higher computational cost. To cope with this issue, an approximate treatment of semicore orbitals has been introduced in GWL as described in Ref. Umari and Fabris, 2012. In addition to collinear spin polarization, GWL provides a fully relativistic non collinear implementation relying on the scalar relativistic calculation of the screened Coulomb interactions Umari et al. (2014).

The bse.x code of the GWL package performs BSE calculations and permits to evaluate either the entire frequency-dependent complex dielectric function through the Lanczos algorithm or a discrete set of excited states and their energies through a conjugate gradient minimization. In contrast to ordinary implementations, bse.x scales as $N^{3}$ instead of $N^{4}$ with respect to the system size $N$ (e.g., the number of atoms) thanks to the use of maximally localized Wannier functions for representing the valence manifold Marsili et al. (2017). The bse.x code, apart from reading the screened Coulomb interaction at zero frequency from a gww.x calculation, works as a separate code and uses the Quantum ESPRESSO environment. Therefore it could be easily be interfaced with other $GW$ codes.

The SternheimerGW package calculates the frequency-dependent $GW$ self-energy and the corresponding quasiparticle corrections at arbitrary ${\bf k}$-points in the BZ. This feature enables accurate calculations of band structures and effective masses without resorting to interpolation. The availability of the complete $GW$ self-energy (as opposed to the quasiparticle shifts) makes it possible to calculate spectral functions, for example including plasmon satellites Caruso et al. (2015). The spectral function can be directly compared to angle-resolved photoelectron spectroscopy (ARPES) experiments. In SternheimerGW the screened Coulomb interaction $W$ is evaluated for wave-vectors in the irreducible BZ by exploiting crystal symmetries. Calculations of $G$ and $W$ for multiple frequencies $\omega$ rely on the use of multishift linear system solvers that construct solutions for all frequencies from the solution of a single linear system Giustino et al. (2010); Lambert and Giustino (2013). This method is closely related to the Lanczos approach. The convolution in the frequency domain required to obtain the self energy from $G$ and $W$ can be performed either on the real axis or the imaginary axis. Padé functions are employed to perform approximate analytic continuations from the imaginary to the real frequency axis; the standard Godby-Needs plasmon pole model is also available to compare with literature results. The stability and portability of the SternheimerGW package are verified via a test-suite and a Buildbot test-farm (see Sec. III.6).

II.3 Other spectroscopies

The QE-GIPAW package allows for the calculation of various physical parameters measured in nuclear magnetic resonance (NMR) and electron paramagnetic resonance (EPR) spectroscopies. These encompass (i) NMR chemical shift tensors and magnetic susceptibility, (ii) electric field gradient (EFG) tensors, (iii) EPR g-tensor, and (iv) hyperfine coupling tensor.

In QE-GIPAW, the NMR and EPR parameters are obtained from the orbital linear response to an external uniform magnetic field. The response depends critically upon the exact shape of the electronic wavefunctions near the nuclei. Thus, all-electron wavefunctions are reconstructed from the pseudo-wavefunctions in a gauge- and translationally invariant way using the gauge-including projector augmented-wave (GIPAW) method Pickard and Mauri (2001). The description of a uniform magnetic field within periodic boundary conditions is achieved by the long-wavelength limit ($q\ll 1$) of a sinusoidally modulated field in real space. In practice, for each ${\bf k}$ point, we calculate the first order change of the wavefunctions at ${\bf k}+{\bf q}$, where ${\bf q}$ runs over a star of 6 points. The magnetic susceptibility and the induced orbital currents are then evaluated by finite differences, in the limit of small $q$. The induced magnetic field at the nucleus, which is the central quantity in NMR, is obtained from the induced current via the Biot-Savart law. In QE-GIPAW, the NMR orbital chemical shifts and magnetic susceptibility can be calculated both for insulators Varini et al. (2013) and for metals d’Avezac et al. (2007) (the additional contribution for metals coming from the spin-polarization of valence electrons, namely the Knight shift, can also be computed but it is not yet ready for production at the time of writing). Similarly to the NMR chemical shift, the EPR g-tensor is calculated as the cross product of the induced current with the spin-orbit operator Pickard and Mauri (2002).

For the quantities defined in zero magnetic field, namely the EFG, Mössbauer and relativistic hyperfine tensors, the usual PAW reconstruction of the wavefunctions is sufficient and these are computed as described in Refs. Petrilli et al. (1998); Zwanziger (2009). The hyperfine Fermi contact term, proportional to the spin density evaluated at the nuclear coordinates, however requires the relaxation of the core electrons in response to the magnetization of valence electrons. We implemented the core relaxation in perturbation theory, according to Ref. Bahramy et al., 2007. Basically we compute the spherically averaged PAW spin density around each atom. Then we compute the change of the XC potential, $\Delta V_{\text{XC}}$, on a radial grid, and compute in perturbation theory the core radial wavefunction, both for spin up and spin down. This provides an extra contribution to the Fermi contact, in most cases opposite in sign to and as large as that of valence electrons.

By combining the quadrupole coupling constants derived from EFG tensors and hyperfine splittings, electron nuclear double resonance (ENDOR) frequencies can be calculated. Applications highlighting all these features of the QE-GIPAW package can be found in Ref. von Bardeleben et al., 2014. These quantities are also needed to compute NMR shifts in paramagnetic systems, like novel cathode materials for Li batteries Pigliapochi et al. (2017). Previously restricted to norm-conserving pseudopotentials only, all features are now applicable using any kind of pseudization scheme and to PAW, following the theory described in Yates et al. (2007). The use of smooth pseudopotentials allows for the calculation of chemical shifts in systems with several hundreds of atoms Küçükbenli et al. (2012).

The starting point of all QE-GIPAW calculations is a previous calculation of KS orbitals via PWscf. Hence, much like other linear response routines, the QE-GIPAW code uses many subroutines of PWscf and of the linear response module. As usually done in linear response methods, the response of the unoccupied states is calculated using the completeness relation between occupied and unoccupied manifolds de Gironcoli (1995). As a result, for insulating as well as metallic systems, the linear response of the system is efficiently obtained without the need to include virtual orbitals.

As an alternative to linear response method, the theory of orbital magnetization via Berry curvature Xiao et al. (2005); Thonhauser et al. (2005) can be used to calculate the NMR Thonhauser et al. (2009) and EPR parameters Ceresoli et al. (2010). Specifically, it can be shown that the variation of the orbital magnetization $M^{orb}$ with respect to spin flip is directly related to the g-tensor: $g_{\mu\nu}=g_{e}-\frac{2}{\alpha S}\mathbf{e}_{\mu}\cdot\mathbf{M}^{\rm orb}(\mathbf{e}_{\nu})$, where $g_{e}=2.002319$, $\alpha$ is the fine structure constant, $S$ is the total spin, $\mathbf{e}$ are Cartesian unit vectors, provided that the spin-orbit interaction is explicitly considered in the Hamiltonian. This converse method of calculating the g-tensor has been implemented in an older version of QE-GIPAW. It is especially useful in critical cases where linear response is not appropriate, e.g., systems with quasi-degenerate HOMO-LUMO levels. A demonstration of this method applied to delocalized conduction band electrons can be found in Ref. 151.

The converse method will be shortly ported into the current QE-GIPAW and we will explore the possibility of computing in a converse way the Knight shift as the response to a small nuclear magnetic dipole. The present version of the code allows for parameter-free calculations of g-tensors, hyperfine splittings, and ENDOR frequencies also for systems with total spin $S>1/2$. Such triplet or even higher-spin states give rise to additional spin-spin interactions, that can be calculated within the magnetic dipole-dipole interaction approximation. This interaction results in a fine structure which can be measured in zero magnetic field. This so-called zero-field splitting is being implemented following the methodology described in Ref. Bodrog and Gali (2014).

II.3.2 XSpectra: L2,3 X-ray absorption edges

The XSpectra code Gougoussis et al. (2009a, b) has been extended to the calculation of X-ray absorption spectra at the L2,3-edges Bunău and Calandra (2013). The XSpectra code uses the self-consistent charge density produced by PWscf and acts as a post processing tool Gougoussis et al. (2009a, b); Taillefumier et al. (2002). The spectra are calculated for the L2 edge, while the L3 edge is obtained by multiplying by two (single-particle statistical branching ratio) the L2 edge spectrum and by shifting it by the value of the spin-orbit splitting of the $2p_{1/2}$ core levels of the absorbing atom. The latter can be taken either from a DFT relativistic all-electron calculation on the isolated atom, or from experiments.

In practice, the L3 edge is obtained from the L2 with the spectra_correction.x tool. Such tool contains a table of experimental $2p$ spin-orbit splittings for all the elements. In addition to computing L3 edges, spectra_correction.x allows one to remove states from the spectrum below a certain energy, and to convolute the calculated spectrum with more elaborate broadenings. These operations can be applied to any edge.

To evaluate the X-ray absorption spectrum for a system containing various atoms of the same species but in different chemical environments, one has to sum the contribution by each atom. This could be the case, for example of an organic molecule containing various C atoms in inequivalent sites. Such individual contributions can be computed separately by XSpectra, and the tool molecularnexafs.x allows one to perform their weighted sum taking into account the proper energy reference (initial and final state effects) Fratesi et al. (2013, 2014). One should in fact notice that the reference for initial state effects will depend upon the environment (e.g., the vacuum level for gas phase molecules, or the Fermi level for molecules adsorbed on a metal).

II.4 Other lattice-dynamical and thermal properties

thermo_pw Dal Corso (a) is a collection of codes aimed at computing various thermodynamical quantities in the quasi-harmonic approximation. The key ingredient is the vibrational contribution, $F_{ph}$, to the Helmholtz free energy at temperature $T$:

where $\omega_{{\bf q},\nu}$ are phonon frequencies at wave-vector ${\bf q}$, $k_{B}$ is the Boltzmann constant. thermo_pw works by calling Quantum ESPRESSO routines from PWscf and PHonon, that perform one of the following tasks: i) compute the DFT total energy and possibly the stress for a given crystal structure; ii) compute for the same system the electronic band structure along a specified path; iii) compute for the same system phonon frequencies at specified wave-vectors. Using such quantities, thermo_pw can calculate numerically the derivatives of the free energy with respect to the external parameters (e.g., different volumes). Several calls to such routines, with slightly different geometries, are typically needed in a run. All such tasks can be independently performed on different groups of processors (called images).

When the tasks carried out by different images require approximately the same time, or when the amount of numerical work needed to accomplish each task is easy to estimate a priori, it would be possible to statically assign tasks to images at the beginning of the run so that images do not need to communicate during the run. However, such conditions are seldom met in thermo_pw and therefore it would be impossible to obtain a good load balancing between images. thermo_pw takes advantage of an engine that controls these different tasks in an asynchronous way, dynamically assigning tasks to the images at run time.

At the core of thermo_pw there is a module mp_asyn, based on MPI routines, that allows for asynchronous communication between different images. One of the images is the “master” and assigns tasks to the other images (the “slaves”) as soon as they communicate that they have accomplished the previously assigned task. The master image also accomplishes some tasks but once in a while, with negligible overhead, it checks if there is an image available to do some work; if so, it assigns to it the next task to do. The code stops when the master recognizes that all the tasks have been done and communicates this information to the slaves. The routines of this communication module are quite independent of the thermo_pw variables and in principle can be used in conjunction with other codes to set up complex workflows to be executed in a massively parallel environment. It is assumed that each processor of each image reads the same input and that the only information that the image needs to synchronize with the other images is which tasks to do. The design of thermo_pw makes it easily extensible to the calculation of new properties in an incremental way.

II.4.2 thermal2: phonon-phonon interaction and thermal transport

Phonon-phonon interaction (ph-ph) plays a role in different physical phenomena: phonon lifetime (and its inverse, the linewidth), phonon-driven thermal transport in insulators or semi-metals, thermal expansion of materials. Ph-ph is possible because the harmonic Hamiltonian of ionic motion, of which phonons are stationary states, is only approximate. At first order in perturbation theory we have the third derivative of the total energy with respect to three phonon perturbations, which we compute ab-initio. This calculation is performed by the d3q code via the $2n+1$ theorem Paulatto et al. (2013); Lazzeri and de Gironcoli (2002). The d3q code is an extension of the old D3 code, which only allowed the calculation of zone-centered phonon lifetimes and of thermal expansion. The current version can compute the three-phonon matrix element of arbitrary wave-vectors $D^{(3)}({\bf q}_{1},{\bf q}_{2},{\bf q}_{3})=\partial^{3}E/\partial u_{{\bf q}_{1}}\partial u_{{\bf q}_{2}}\partial u_{{\bf q}_{3}}$, where $u$ are the phonon displacement patterns, the momentum conservation rule imposes ${\bf q}_{1}+{\bf q}_{2}+{\bf q}_{3}=0$. The current version of the code can treat any kind of crystal geometry, metals and insulators, both local density and gradient-corrected functionals, and multi-projector norm-conserving pseudopotentials. Ultrasoft pseudopotentials, PAW, spin polarization and non-collinear magnetization are not implemented. Higher order derivative of effective chargesDeinzer et al. (2003) are not implemented.

The ph-ph matrix elements, computed from linear response, can be transformed, via a generalized Fourier transform, to the real-space three-body force constants which could be computed in a supercell by finite difference derivation:

where $F^{3}(\boldsymbol{0},\mathbf{R}^{\prime},\mathbf{R}^{\prime\prime})=\partial^{3}E/\partial\tau\partial(\tau^{\prime}+\mathbf{R}^{\prime})\partial(\tau^{\prime\prime}+\mathbf{R}^{\prime\prime})$ is the derivative of the total energy w.r.t. nuclear positions of ions with basis $\tau$, $\tau^{\prime}$, $\tau^{\prime\prime}$ from the unit cells identified by direct lattice vectors $\boldsymbol{0}$ (the origin), $\mathbf{R}^{\prime}$ and $\mathbf{R}^{\prime\prime}$. The sum over $\mathbf{R}^{\prime}$ and $\mathbf{R}^{\prime\prime}$ runs, in principle, over all unit cells, however the terms of the sum quickly decay as the size of the triangle $\boldsymbol{0}-\mathbf{R}^{\prime}-\mathbf{R}^{\prime\prime}$ increases. The real-space finite-difference calculation, performed by some external softwaresLi et al. (2014), has some advantages: it is easier to implement and it can readily include all the capabilities of the self-consistent code; on the other hand it is much more computationally expensive than the linear-response method we use, its cost scaling with the cube of the supercell volume, or the 9th power of the number of side units of an isotropic system. We use the real-space formalism to apply the sum rule corresponding to translational invariance to the matrix elements. This is done with an iterative method that alternatively enforces the sum rule on the first matrix index and restores the invariance on the order of the derivations. We also use the real-space force constants to Fourier-interpolate the ph-ph matrices on a finer grid, assuming that the contribution from triangles $\boldsymbol{0}-\mathbf{R}^{\prime}-\mathbf{R}^{\prime\prime}$ which we have not computed is zero; it is important in this stage to consider the periodicity of the system.

From many-body theory we get the first-order phonon linewidthCalandra et al. (2007) ($\gamma_{\nu}$) of mode $\nu$ at ${\bf q}$, which is a sum over all the possible $N_{\bf q}$’s final and initial states (${\bf q}^{\prime}$,$\nu^{\prime}$,$\nu^{\prime\prime}$) with conservation of energy ($\hbar\omega$) and momentum (${\bf q}^{\prime\prime}=-{\bf q}-{\bf q}^{\prime}$), Bose-Einstein occupations ($n({\bf q},\nu)=(\exp(\hbar\omega_{{\bf q},\nu}/k_{B}T)-1)^{-1}$) and an amplitude $V^{(3)}$, proportional to the $D^{(3)}$ matrix element but renormalized with phonon energies and ion masses:

This sum is computed in the thermal2 suite, which is bundled with d3q. A similar expression can be written for the phonon scattering probability which appears in the Boltzmann transport equation. In order to properly converge the integral of the Dirac delta function, we express it as finite-width Gaussian function and use an interpolation grid. This equation can be solved either exactly or in the single mode approximation (SMA) Callaway (1959). The SMA is a good tool at temperatures comparable to or larger than the Debye temperature, but is known to be inadequate at low temperatures Markov et al. (2016, ) or in the case of 2D materials Ward et al. (2009); Fugallo et al. (2014); Cepellotti et al. (2015). The exact solution is computed using a variational form, minimized via a preconditioned conjugate gradient algorithm, which is guaranteed to converge, usually in less than 10 iterations Fugallo et al. (2013).

On top of intrinsic ph-ph events, the thermal2 codes can also treat isotopic disorder and substitution effects and finite transverse dimension using the Casimir formalism. In addition to using our force constants from DFpT, the code supports importing 3-body force constants computed via finite differences with the thirdorder.py code Li et al. (2014). Parallelization is implemented with both MPI (with great scalability up to thousand of CPUs) and OpenMP (optimal for memory reduction).

II.4.3 EPW: Electron-phonon coefficients from Wannier interpolation

The electron-phonon-Wannier (EPW) package is designed to calculate electron-phonon coupling using an ultra-fine sampling of the BZ by means of Wannier interpolation. EPW employs the relation between the electron-phonon matrix elements in the Bloch representation $g_{mn\nu}({\bf k},{\bf q})$, and in the Wannier representation, $g_{ij\kappa\alpha}(\mathbf{R},\mathbf{R}^{\prime})$ Giustino (2017),

in order to interpolate from coarse ${\bf k}$-point and ${\bf q}$-point grids into dense meshes. In the above expression ${\bf k}$ and ${\bf q}$ represent the electron and phonon wave-vector, respectively, the indices $m,n$ and $i,j$ refer to Bloch states and Wannier states, respectively, and $\mathbf{R},\mathbf{R}^{\prime}$ are direct lattice vectors. The matrices $U_{mi{\bf k}}$ are unitary transformations and the vector $u_{\kappa\alpha,{\bf q}\nu}$ is the displacements of the atom $\kappa$ along the Cartesian direction $\alpha$ for the phonon of wavevector ${\bf q}$ and branch $\nu$. The interpolation is performed with ab initio accuracy by relying on the localization of maximally-localized Wannier functions Marzari et al. (2012). During its execution EPW invokes the Wannier90 software Mostofi et al. (2014) in library mode in order to determine the matrices $U_{mi{\bf k}}$ on the coarse ${\bf k}$-point grid.

EPW can be used to compute the following physical properties: the electron and phonon linewidths arising from electron-phonon interactions; the scattering rates of electrons by phonons; the total, averaged electron-phonon coupling strength; the electrical resistivity of metals, see Fig. 2(b); the critical temperature of electron-phonon superconductors; the anisotropic superconducting gap within the Eliashberg theory, see Fig. 2(c); the Eliashberg spectral function, transport spectral function, see Fig. 2(d) and the nesting function. The calculation of carrier mobilities using the Boltzmann transport equation in semiconductors is under development.

The epw.x code exploits crystal symmetry operations (including time reversal) in order to limit the number of phonon calculations to be performed using PHonon to the irreducible wedge of the BZ. The code supports calculations of electron-phonon couplings in the presence of spin-orbit coupling. The current version does not support spin-polarized calculations, ultrasoft pseudopotentials nor the PAW method. As shown in Fig. 2(a), epw.x scales reasonably up to 2,000 cores using MPI. A test farm (see Sec. III.6) was set up to ensure portability of the code on many architecture and compilers. Detailed information about the EPW package can be found in Ref. Poncé et al., 2016.

II.4.4 Non-perturbative approaches to vibrational spectroscopies

Although DFpT is in many ways the state of the art in the simulation of vibrational spectroscopies in extended systems, and in fact one of defining features of Quantum ESPRESSO, it is sometimes convenient to compute lattice-dynamical properties, the response to macroscopic electric fields, or combinations thereof (such as e.g., the infrared or Raman activities), using non-perturbative methods. This is so because DFpT requires the design of dedicated codes, which have to be updated and maintained separately, and which therefore not always follow the pace of the implementation of new features, methods, and functionals (such as e.g., DFT+U, vdW-DF, hybrid functionals, or ACBN0 Agapito et al. (2015)) in their ground-state counterparts. Such a non-perturbative approach is followed in the FD package, which implements the “frozen-phonon” method for the computation of phonons and vibrational spectra: the interatomic Force Constants (IFCs) and electronic dielectric constant are computed as finite differences of forces and polarizations, with respect to finite atomic displacements or external electric fields, respectively Calzolari and Buongiorno Nardelli (2013); Umari and Pasquarello (2002). IFC’s are computed in two steps: first, code fd.x generates the symmetry-independent displacements in an appropriate supercell; after the calculations for the various displacements are completed, code fd_ifc.x reads the forces and generates IFC’s. These are further processed in matdyn.x, where non-analytical long-ranged dipolar terms are subtracted out from the IFCs following the recipe of Ref. 175. The calculation of dielectric tensor and of the Born effective charges proceeds from the evaluation of the electronic susceptibility following the method proposed by Umari and Pasquarello Umari and Pasquarello (2002), where the introduction of a non local energy functional $E^{\mathbfcal{E}}_{tot}[{\psi}]=E^{0}[{\psi}]-\mathbfcal{E}\cdot(\mathbf{P}^{ion}+\mathbf{P}^{el}[{\psi}])$ allows to compute the electronic structure for periodic systems under finite homogeneous electric fields. $E^{0}$ is the ground state total energy in the absence of external electric fields; $\mathbf{P}^{ion}$ is the usual ionic polarization, and $\mathbf{P}^{el}$ is given as a Berry phase of the manifold of the occupied bands King-Smith and Vanderbilt (1993). The high-frequency dielectric tensor $\epsilon^{\infty}$ is then computed as $\epsilon^{\infty}_{ij}=\delta_{i,j}+4\pi\chi_{ij}$, while Born effective-charge tensors $Z^{*}_{I,ij}$ are obtained as the polarization induced along the direction $i$ by a unit displacement of the $I$-th atom in the $j$ direction; alternatively, as the force induced on atom $I$ by an applied electric field, $\mathbfcal{E}$.

The calculation of the Raman spectra proceeds along similar lines. Within the finite-field approach, the Raman tensor is evaluated in terms of finite differences of atomic forces in the presence of two electric fields Umari et al. (2003). In practice, the tensor $\chi^{(1)}_{ijIk}$ is obtained from a set of calculations combining finite electric fields along different Cartesian directions. $\chi^{(1)}_{ijIk}$ is then symmetrized to recover the full symmetry of the structure under study.

II.5 Multi-scale modeling

Continuum models are among the most popular multiscale approaches to treat solvation effects in the quantum-chemistry community Tomasi et al. (2005). In this class of models, the degrees of freedom of solvent molecules are effectively integrated out and their statistically-averaged effects on the solute are mimicked by those of a continuous medium surrounding a cavity in which the solute is thought to dwell. The most important interaction usually handled with continuum models is the electrostatic one, in which the solvent is described as a dielectric continuum characterized by its experimental dielectric permittivity.

Following the original work of Fattebert and Gygi Fattebert and Gygi (2002) , a new class of continuum models was designed, in which a smooth transition from the QM-solute region to the continuum-environment region of space is introduced and defined in terms of the electronic density of the solute. The corresponding free energy functional is optimized using a fully variational approach, leading to a generalized Poisson equation that is solved via a multi-grid solverFattebert and Gygi (2002). This approach, ideally suited for plane-wave basis sets and tailored for MD simulations, has been featured in the Quantum ESPRESSO distribution since v. 4.1. This approach was recently revisedAndreussi et al. (2012), by defining an optimally smooth QM/continuum transition, reformulated in terms of iterative solversFisicaro et al. (2016) and extended to handle in a compact and effective way non-electrostatic interactions Andreussi et al. (2012). The resulting self-consistent continuum solvation (SCCS) model, based on a very limited number of physically justified parameters, allows one to reproduce experimental solvation energies for aqueous solutions of neutral Andreussi et al. (2012) and chargedDupont et al. (2013) species with accuracies comparable to or higher than state-of-the-art quantum-chemistry packages.

The SCCS model involves different embedding terms, each representing a specific interaction with an external continuum environment and contributing to the total energy, KS potential, and interatomic forces of the embedded QM system. Every contribution may depend explicitly on the ionic (rigid schemes) and/or electronic (self-consistent or soft schemes) degrees of freedom of the embedded system. All the different terms are collected in the stand-alone Environ module Andreussi et al. (2016). The present discussion refers to release 0.2 of Environ, which is compatible with Quantum ESPRESSO starting from versions 5.1. The module requires a separate input file with the specifications of the environment interactions to be included and of the numerical parameters required to compute their effects. Fully parameterized and tuned SCCS environments, e.g., corresponding to water solutions for neutral and charged species, are directly available to the users. Otherwise individual embedding terms can be switched on and tuned to the specific physical conditions of the required environment. Namely, the following terms are currently featured in Environ:

Smooth continuum dielectric, with the associated generalized Poisson problem solved via a direct iterative approach or through a preconditioned conjugate gradient algorithm Fisicaro et al. (2016).

Electronic enthalpy functional, introducing an energy term proportional to the quantum-volume of the system and able to describe finite systems under the effect of an applied external pressureCococcioni et al. (2005).

Electronic cavitation functional, introducing an energy term proportional to the quantum-surface able to describe free energies of cavitation and other surface-related interaction termsScherlis et al. (2006).

Parabolic corrections for periodic boundary conditions in aperiodic and partially periodic (slab) systems Dabo et al. (2008); Andreussi and Marzari (2014).

Fixed dielectric regions, allowing for the modelling of complex inhomogenous dielectric environments.

Fixed Gaussian-smoothed distributions of charges, allowing for a simplified modelling of countercharge distributions, e.g., in electrochemical setups.

Different packages of the Quantum ESPRESSO distribution have been interfaced with the Environ module, including PWscf, CP, PWneb, and turboTDDFT, the latter featuring a linear-response implementation of the SCCS model (see Sec.II.2.2). Moreover, continuum environment effects are fully compatible with the main features of Quantum ESPRESSO, and in particular, with reciprocal space integration and smearing for metallic systems, with both norm-conserving and ultrasoft pseudopotentials and PAW, with all XC functionals.

II.5.2 QM-MM

QM-MM was implemented in v.5.0.2 using the method documented in Ref. Ma et al., 2015. Such methodology accounts for both mechanical and electrostatic coupling between the QM (quantum-mechanical) and MM (molecular-mechanics) regions, but not for bonding interactions (i.e., bonds between the QM and MM regions). In practice, we need to run two different codes, Quantum ESPRESSO for the QM region and a classical force-field code for the MM region, that communicate atomic positions, forces, electrostatic potentials.

LAMMPS Plimpton (1995) is the software chosen to deal with the classical (MM) degrees of freedom. This is a well-known and well-maintained package, released under an open-source license that allows redistribution together with Quantum ESPRESSO. The communications between the QM and MM regions use a “shared memory” approach: the MM code runs on a master node, communicates directly via the memory with the QM code, which is typically running on a massively parallel machine. Such approach has some advantages: the MM part is typically much faster than the QM one and can be run in serial execution, wasting no time on the HPC machine; there is a clear and neat separation between the two codes, and very small code changes in either codes are needed. It has however also a few drawbacks, namely: the serial computation of the MM part may become a bottleneck if the MM region contains many atoms; direct access to memory is often restricted for security reasons on HPC machines.

An alternative approach has been implemented in v.5.4. A single (parallel) executable runs both the MM and the QM codes. The two codes exchange data and communicate via MPI. This approach is less elegant than the previous one and requires a little bit more coding, but its implementation is quite straightforward thanks also to the changes in the logic of parallelization mentioned in Sec. III.4. The coupling of the two codes has required some modifications also to the qmmm library inside LAMMPS and to the related fix qmmm (a “fix” in LAMMPS is any operation that is applied to the system during the MD run). In particular, electrostatic coupling has been introduced, following the approach described in Ref. Laio et al., 2002. The great advantage of this approach is that its performance on HPC machines is as good as the separate performances of the QM and MM codes. Since LAMMPS is very well parallelized, this is a significant advantage if the MM region contains many atoms. Moreover, it can be run without restrictions on any parallel machine. This new QM-MM implementation is an integral part of the Quantum ESPRESSO distribution and will soon be included into LAMMPS as well (the “fix” is currently under testing) and it is straightforward to compile and execute it.

II.6 Miscellaneous feature enhancements and additions

Here $T_{D}$ is the Dirac kinetic energy:

II.6.2 Electronic and structural properties in field-effect configuration

Since Quantum ESPRESSO v.6.0 it is possible to compute the electronic structure under a field-effect transistor (FET) setup in periodic boundary conditions Brumme et al. (2014). In physical FETs, a voltage is applied to a gate electrode, accumulating charges at the interface between the gate dielectric and a semiconducting system (see Fig. 3). The gate electrode is simulated with a charged plate, henceforth referred to as the gate. Since the interaction of this charged plate with its periodic image generates a spurious nonphysical electric field, a dipolar correction, equivalent to two planes of opposite charge, is added Bengtsson (1999), canceling out the field on the left side of the gate. In order to prevent electrons from spilling towards the gate for large electron doping Topsakal and Ciraci (2012), a potential barrier can be added to the electrostatic potential, mimicking the effect of the gate dielectric.

II.6.3 Cold restart of Car-Parrinello molecular dynamics

In the standard Lagrangian formulation of ab initio molecular dynamics Car and Parrinello (1985), the coefficients of KS molecular orbitals over a given basis set (i.e.. their Fourier coefficients, in the case of plane waves) are treated as classical degrees of freedom obeying Newton’s equations of motion that derive from a suitably defined extended Lagrangian. This Lagrangian is obtained from the Born-Oppenheimer total energy by augmenting it with a fictitious electronic kinetic-energy term and relaxing the constraint that the molecular orbitals stay at each instant of the trajectory in their instantaneous KS ground state. The idea is that, by choosing a suitably small fictitious electronic mass, the thermalization time of the electronic degrees of freedom can be made much longer than the typical simulation times, so that if the system is prepared in its electronic KS ground state at the start of the simulation, the electronic dynamics would follow almost adiabatically the nuclear one all over the simulation, thus effectively mimicking a bona fide Born-Oppenheimer dynamics.

While in Car-Parrinello MD both the physical nuclear and fictitious electronic velocities are determined by the equations of motion on a par, the question still remains as to how choose them at the start of the simulation. Initial nuclear velocities are dictated by physical considerations (e.g., thermal equilibrium) or may be taken from a previously interrupted MD run. Electronic velocities (i.e., the time derivatives of the KS molecular orbitals), instead, are not available when the simulation is started from scratch. and are not independent of the physical nuclear ones, but are determined by the adiabatic time evolution of the system. Moreover, the projection over the occupied-state manifold of the electronic velocities, $\dot{\psi}_{v}^{\parallel}\doteq\hat{P}\dot{\psi}_{v}$ is ill-defined because the KS ground-state solution is defined modulo a unitary transformation within this manifold. This means that the starting electronic velocities may not be simply obtained as finite differences of KS orbitals at times $t=0$ and $t=\Delta t$. Here and in the following $\hat{P}$ indicates the projector over the occupied-state manifold, and $\hat{Q}=1-\hat{P}$ its complement (i.e. the projector over the virtual-orbital manifold).

The component of the electronic velocities over the virtual-state manifold, $\dot{\psi}_{v}^{\perp}\doteq\hat{Q}\dot{\psi}_{v}$, is instead well defined and can be formally written using standard first-order perturbation theory:

While this could in principle be done using density-functional perturbation theory Baroni et al. (1987, 2001), it is more convenient to compute them numerically, following the procedure described below. At $t=0$ the KS molecular orbitals are initialized from a ground-state computation, performed with whatever method is available or preferred (standard SCF calculation or global optimization, such as e.g., with conjugate gradients Štich et al. (1989)). The KS molecular orbitals that would result from a perfectly adiabatic propagation at $t=\Delta t$ are then determined from a second ground-state computation, performed after half a “velocity-Verlet” MD step, i.e., at nuclear positions $\mathbf{R}(\Delta t)=\mathbf{R}(0)+\dot{\mathbf{R}}(0)\Delta t$. The initial velocities are then obtained from the relation:

which is obtained by simply differentiating the definition of occupied-state projector, $\hat{P}\psi_{v}=\psi_{v}$. The right-hand side of Eq. (48) is finally easily computed by subtracting from each KS orbital at time $t=0$, its component over the occupied-state manifold at $t=\Delta t$ and dividing by $\Delta t$.

II.6.4 Optimized tetrahedron method

The integration over ${\bf k}$-points in the BZ is a crucial step in the calculation of the electronic structure of a periodic system, affecting not only the ground state but linear response as well. This is especially true for metallic systems where the integrand is discontinuous at the Fermi level. Even more problematic is the integration of Dirac delta functions, such as those appearing in the density of states (DOS), partial DOS and in the electron-phonon coupling constant.

Quantum ESPRESSO has always implemented a variety of “smearing” methods, in which the delta function is replaced by a function of finite width (e.g., a Gaussian function, or more sophisticated choices). It has also always implemented the linear tetrahedron method Jepsen and Andersen (1971) with the correction proposed by Blöchl Blöchl et al. (1994), in which the BZ is divided into tetrahedra and the integration is performed analytically by linear interpolation of KS eigenvalues in each tetrahedron. Such method is however limited in its convenience and range of applicability: in fact the linear interpolation systematically overestimates convex functions, thus making the convergence against the number of ${\bf k}$-points slow. The linear tetrahedron method was thus mostly restricted to the calculation of DOS and partial DOS.

Since Quantum ESPRESSO v.6.1, the optimized tetrahedron method Kawamura et al. (2014) is implemented. Such method overcomes the drawback of the linear tetrahedron method using an interpolation that accounts for the curvature of the interpolated function. The optimized tetrahedron method has better convergence properties and an extended range of applicability: in addition to the calculation of the ground-state charge density, DOS and partial DOS, it can be used in linear-response calculation of phonons and of the electron-phonon coupling constant.

II.6.5 Wyckoff positions

In Quantum ESPRESSO the crystal geometry is traditionally specified by a Bravais lattice index (called ibrav), by the crystal parameters (celldm, or a, b, c, cosab, cosac, cosbc) describing the unit cell, and by the positions of all atoms in the unit cell, in crystal or Cartesian axis.

Since v.5.1.1, it is possible to specify the crystal geometry in crystallographic styleZadra and Dal Corso , according to the notations of the International Tables of Crystallography (ITA)Hahn (2005). A complete description of the crystal structure is obtained by specifying the space-group number according to the ITA and the positions of symmetry-inequivalent atoms only in the unit cell. The latter can be provided either in the crystal axis of the conventional cell, or as Wyckoff positions: a set of special positions, listed in the ITA for each space group, that can be fully specified by a number of parameters, none to three depending upon the site symmetry. Table 1 reports a few examples of accepted syntax.

The code generates the symmetry operations for the specified space group and applies them to inequivalent atoms, thus finding all atoms in the unit cell.

For some crystal systems there are alternate descriptions in the ITA, so additional input parameters may be needed to select the desired one. For the monoclinic system the “c-unique” orientation is the default and bunique=.TRUE. must be specified in input if the “b-unique” orientation is desired. For some space groups there are two possible choices of the origin. The origin appearing first in the ITA is chosen by default, unless origin_choice=2 is specified in input. Finally, for trigonal space groups the atomic coordinates can be referred to the rhombohedral or to the hexagonal Bravais lattices. The default is the rhombohedral lattice, so rhombohedral=.FALSE. must be specified in input to use the hexagonal lattice.

A final comment for centered Bravais lattices: in the crystallographic literature, the conventional unit cell is usually assumed. Quantum ESPRESSO however assumes the primitive unit cell, having a smaller volume and a smaller number of atoms, and discards atoms outside the primitive cell. Auxiliary code supercell.x, available in thermo_pw (see Sec.II.4.1), prints all atoms in the conventional cell when necessary.

III Parallelization, modularization, interoperability and stability

The basic modules of Quantum ESPRESSO are characterized by a hierarchy of parallelization levels, described in Ref.Giannozzi et al., 2009. Processors are divided into groups, labeled by a MPI communicator. Each group of processors distributes a specific subset of computations. The growing diffusion of HPC machines based on nodes with many cores (32 and more) makes however pure MPI parallelization not always ideal: running one MPI process per core has a high overhead, limiting performances. It is often convenient to use mixed MPI-OpenMP parallelization, in which a small number of MPI processes per node use OpenMP threads, either explicitly (i.e., with compiler directives) or implicitly (i.e., via calls to OpenMP-aware library). Explicit OpenMP parallelization, originally confined to computationally intensive FFT’s, has been extended to many more parts of the code.

One of the challenges presented by massively parallel machine is to get rid of both memory and CPU time bottlenecks, caused respectively by arrays that are not distributed across processors and by non-parallelized sections of code. It is especially important to distribute all arrays and parallelize all computations whose size/complexity increases with the dimensions of the unit cell or of the basis set. Non-parallelized computations hamper “weak” scalability, that is, parallel performance while increasing both the system size and the amount of computational resources, while non-distributed arrays may become an unavoidable RAM bottleneck with increasing problem size. “Strong” scalability (that is, at fixed problem size and increasing number of CPUs) is even more elusive than weak scalability in electronic-structure calculations, requiring, in addition to systematic distribution of computations, to keep to the minimum the ratio between time spent in communications and in computation, and to have a nearly perfect load balancing. In order to achieve strong scalability, the key is to add more parallelization levels and to use algorithms that permit to overlap communications and computations.

For what concerns memory, notable offenders are arrays of scalar products between KS states $\psi_{i}$: $O_{ij}=\langle\psi_{i}|\widehat{O}|\psi_{j}\rangle$, where $\widehat{O}$ can be either the Hamiltonian or an overlap matrix; and scalar products between KS states and pseudopotential projectors $\beta$, $p_{ij}=\langle\psi_{i}|\beta_{j}\rangle$. The size of such arrays grows as the square of the size of the cell. Almost all of them are now distributed across processors of the “linear-algebra group”, that is, the group of processors taking care of linear-algebra operations on matrices. The most expensive of such operations are subspace diagonalization (used in PWscf in the iterative diagonalization) and iterative orthonormalization (used by CP). In both cases, a parallel dense-matrix diagonalization on distributed matrix is needed. In addition to ScaLAPACK, Quantum ESPRESSO can now take advantage of newer ELPA libraries (Ref. Marek et al., 2014), leading to significant performance improvements.

The array containing the plane-wave representation, $c_{{\bf k},n}({\bf G})$, of KS orbitals is typically the largest array, or one of the largest. While plane waves are already distributed across processors of the “plane-wave group” as defined in Ref. Giannozzi et al., 2009, it is now possible to distribute KS orbitals as well. Such a parallelization level is located between the ${\bf k}$-point and the plane-wave parallelization levels. The corresponding MPI communicator defines a subgroup of the “${\bf k}$-point group” of processors and is called “band group communicator”. In the CP package, band parallelization is implemented for almost all available calculations. Its usefulness is better appreciated in simulations of large cells — several hundreds of atoms and more — where the number of processors required by memory distribution would be too large to get good scalability from plane-wave parallelization only.

In PWscf, band parallelization is implemented for calculations using hybrid functionals. The standard algorithm to compute Hartree-Fock exchange in a plane-wave basis set (see Sec. II.1.1) contains a double loop on bands that is by far the heaviest part of computation. A first form of parallelization, described in Ref. Varini et al., 2013, was implemented in v.5.0. In the latest version, this has been superseded by parallelization of pairs of bands, Ref. Barnes et al., 2017. Such algorithm is compatible with the “task-group” parallelization level (that is: over KS states in the calculation of $V\psi_{i}$ products) described in Ref. Giannozzi et al., 2009.

In addition to the above-mentioned groups, that are globally defined and in principle usable in all routines, there are a few additional parallelization levels that are local to specific routines. Their goal is to reduce the amount of non-parallel computations that may become significant for many-atom systems. An example is the calculation of DFT+U (Sec. II.1.3) terms in energy and forces, Eqs. (12) and (14) respectively. In all these expressions, the calculation of the scalar products between valence and atomic wave functions is in principle the most expensive step: for $N_{b}$ bands and $N_{pw}$ plane waves, ${\cal O}(N_{pw}N_{b})$ floating-point operations are required (typically, $N_{pw}\gg N_{b}$). The calculation of these terms is however easily and effectively parallelized, using standard matrix-matrix multiplication routines and summing over MPI processes with a mpi_reduce operation on the plane-wave group. The sum over ${\bf k}$-points can be parallelized on the ${\bf k}$-point group. The remaining sums over band indices $\nu$ and Hubbard orbitals $I,m$ may however require a significant amount of non-parallelized computation if the number of atoms with a Hubbard $U$ term is not small. The sum over band indices is thus parallelized by simply distributing bands over the plane-wave group. This is a convenient choice because all processors of the plane-wave group are available once the scalar products are calculated. The addition of band parallelization speeds up the computation of such terms by a significant factor. This is especially important for Car-Parrinello dynamics, requiring the calculation of forces at each time step, when a sizable number of Hubbard manifolds is present.

III.2 Aspects of interoperability

One of the original goals of Quantum ESPRESSO was to assemble different pieces of rather similar software into an integrated software suite. The choice was made to focus on the following four aspects: input data formats, output data files, installation mechanism, and a common base of code. While work on the first three aspects is basically completed, it is still ongoing on the fourth. It was however realized that a different form of integration — interoperability, i.e., the possibility to run Quantum ESPRESSO along with other software — was more useful to the community of users than tight integration. There are several reasons for this, all rooted in new or recent trends in computational materials science. We mention in particular the usefulness of interoperability for

excited-states calculations using many-body perturbation theory, at various levels of sophistication: $GW$, TDDFT, BSE (e.g., yambo Marini et al. (2009), SaX Martin-Samos and Bussi (2009), or BerkeleyGW Deslippe et al. (2012));

calculations using quantum Monte Carlo methods;

configuration-space sampling, using such algorithms as nudged elastic band (NEB), genetic/evolutionary algorithms, meta-dynamics;

inclusion of quantum effects on nuclei via path-integral Monte Carlo;

multi-scale simulations, requiring different theoretical approaches, each valid in a given range of time and length scale, to be used together;

high-throughput, or “exhaustive”, calculations (e.g., AiiDA Pizzi et al. (2016); Mounet et al. (2016) and AFLOW$\pi$ Supka et al. (2017)) requiring automated submission, analysis, retrieval of a large number of jobs;

“steering”, i.e., controlling the computation in real time using either a graphical user interface (GUI) or an interface in a high-level interpreted language (e.g., python).

It is in principle possible, and done in some cases, to implement all of the above into Quantum ESPRESSO, but this is not always the best practice. A better option is to use Quantum ESPRESSO in conjunction with external software performing other tasks.

Cases 1 and 2 mentioned above typically use as starting step the self-consistent solution of KS equations, so that what is needed is the possibility for external software to read data files produced by the main Quantum ESPRESSO codes, notably the self-consistent code PWscf and the molecular dynamics code CP.

Cases 3 and 4 typically require many self-consistent calculations at different atomic configurations, so that what is needed is the possibility to use the main Quantum ESPRESSO codes as “computational engine”, i.e., to call PWscf or CP from an external software, using atomic configurations supplied by the calling code.

The paradigmatic case 5 is QM-MM (Sec.II.5.2), requiring an exchange of data, notably: atomic positions, forces, and some information on the electrostatic potential, between Quantum ESPRESSO and the MM code – typically a classical MD code.

Case 6 requires easy access to output data from one simulation, and easy on-the-fly generation of input data files as well. This is also needed for case 7, which however may also require a finer-grained control over computations performed by Quantum ESPRESSO routines: in the most sophisticated scenario, the GUI or python interface should be able to perform specific operations “on the fly”, not just running an entire self-consistent calculation. This scenario relies upon the existence of a set of application programming interfaces (API’s) for calls to basic computational tasks.

III.3 Input/Output and data file format

On modern machines, characterized by fast CPU’s and large RAM’s, disk input/output (I/O) may become a bottleneck and should be kept to a strict minimum. Since v.5.3 both PWscf and CP do not perform by default any I/O at run time, except for the ordinary text output (printout), for checkpointing if required or needed, and for saving data at the end of the run. The same is being gradually extended to all codes. In the following, we discuss the case of the final data writing.

The original organization of output data files (or more exactly, of the output data directory) was based on a formatted “head” file, with a XML-like syntax, containing general information on the run, and on binary data files containing the KS orbitals and the charge density. We consider the basic idea of such approach still valid, but some improvements were needed. On one hand, the original head file format had a number of small issues—inconsistencies, missing pieces of relevant information—and used a non-standard syntax, lacking a XML “schema” for validation. On the other hand, data files suffered from the lack of portability of Fortran binary files and had to be transformed into text files, sometimes very large ones, in order to become usable on a different machine.

Since v.6.0, the “head” file is a true XML file using a consistent syntax, described by a XML schema, that can be easily parsed with standard XML tools. It also contains complete information on the run, including all data needed to reproduce the results, and on the correct execution and exit status. This aspect is very useful for high-throughput applications, for databasing of results and for verification and validation purposes.

The XML file contains an input section and can thus be used as input file, alternative to the still existing text-based input. It supersedes the previous XML-based input, introduced several years ago, that had a non-standard syntax, different from and incompatible with the one of the original head file. Implementing a different input is made easy by the clear separation existing between the reading and initialization phases: input data is read, stored in a separate module, copied to internal variables.

The current XML file can be easily parsed and generated using standard XML tools and is especially valuable in conjunction with GUI’s. The schema can be found at the URL: http://www.quantum-espresso.org/ns/qes/qes-1.0.xsd.

III.3.2 Large-record data file format

Although not as I/O-bound as other kinds of calculations, electronic-structure simulations may produce a sizable amount of data, either intermediate or needed for further processing. The largest array typically contains the plane-wave representation of KS orbitals; other sizable arrays contain the charge and spin density, either in reciprocal or in real space. In parallel execution using MPI, large arrays are distributed across processors, so one has two possibilities: let each MPI process write its own slice of the data (“distributed” I/O), or collect the entire array on a single processor before writing it (“collected” I/O). In distributed I/O, coding is straightforward and efficient, minimizing file size and achieving some sort of I/O parallelization. A global file system, accessible to all MPI processes, is needed. The data is spread into many files that are directly usable only by a code using exactly the same distribution of arrays, that is, exactly the same kind of parallelization. In collected I/O, the coding is less straightforward. In order to ensure portability, some reference ordering, independent upon the number of processors and the details of the parallelization, must be provided. For large simulations, memory usage and communication pattern must be carefully optimized when a distributed array is collected into a large array on a single processor.

In the original I/O format, KS orbitals were saved in reciprocal space, in either distributed or collected format. For the latter, a reproducible ordering of plane waves (including the ordering within shells of plane waves with the same module), independent upon parallelization details and machine-independent, ensures data portability. Charge and spin density were instead saved in real space and in collected format. In the new I/O scheme, available since v.6.0, the output directory is simplified, containing only the XML data file, one file per ${\bf k}$-point with KS orbitals, one file for the charge and spin density. Both files are in collected format and both quantities are stored in reciprocal space. In addition to Fortran binary, it is possible to write data files in HDF5 formatThe HDF Group (2010). HDF5 offers the possibility to write structured record and portability across architectures, without significant loss in performances; it has an excellent support and is the standard for I/O in other fields of scientific computing. Distributed I/O is kept only for checkpointing or as a last-resort alternative.

In spite of its advantages, such a solution has still a bottleneck in large-scale computations on massively parallel machines: a single processor must read and write large files. Only in the case of parallelization over ${\bf k}$-points is I/O parallelized in a straightforward way. More general solutions to implement parallel I/O using parallel extensions of HDF5 are currently under examination in view of enabling Quantum ESPRESSO towards “exascale” computing (that is: towards ${\cal O}(10^{18})$ floating-point operations per second).

III.4 Organization of the distribution

Codes contained in Quantum ESPRESSO have evolved from a small set of original codes, born with rather restricted goals, into a much larger distribution via continuous additions and extensions. Such a process - presumably common to most if not all scientific software projects - can easily lead to uncoordinated growth and to bad decisions that negatively affect maintainability.

In order to make the distribution easier to maintain, extend and debug, the distribution has been split into

base distribution, containing common libraries, tools and utilities, core packages PWscf, CP, PostProc, plus some commonly used additional packages, currently: atomic, PWgui, PWneb, PHonon, XSpectra, turboTDDFT, turboEELS, GWL, EPW;

external packages such as SaX Martin-Samos and Bussi (2009), yambo Marini et al. (2009), Wannier90 Mostofi et al. (2014), WanT Calzolari et al. (2004); Ferretti et al. (2007), that are automatically downloaded and installed on demand.

The directory structure now explicitly reflects the structure of Quantum ESPRESSO as a “federation” of packages rather than a monolithic one: a common base distribution plus additional packages, each of which fully contained into a subdirectory.

In the reorganization process, the implementation of the NEB algorithm was completely rewritten, following the paradigm sketched in Sec. III.2. PWneb is now a separate package that implements the NEB algorithm, using PWscf as the computational engine. The separation between the NEB algorithm and the self-consistency algorithm is quite complete: PWneb could be adapted to work in conjunction with a different computational engine with a minor effort.

The implementation of meta-dynamics has also been re-considered. Given the existence of a very sophisticated and well-maintained packageBonomi et al. (2009) Plumed for all kinds of meta-dynamics calculations, the PWscf and CP packages have been adapted to work in conjunction with Plumed v.1.x, removing the old internal meta-dynamics code. In order to activate meta-dynamics, a patching process is needed, in which a few specific “hook” routines are modified so that they call routines from Plumed.

III.4.2 Modular parallelism

The logic of parallelism has also evolved towards a more modular approach. It is now possible to have all Quantum ESPRESSO routines working inside a MPI communicator, passed as argument to an initialization routine. This allows in particular the calling code to have its own parallelization level, invisible to Quantum ESPRESSO routines; the latter can thus perform independent calculations, to be subsequently processed by the calling code. For instance: the “image” parallelization level, used by NEB calculations, is now entirely managed by PWneb and no longer in the called PWscf routines. Such a feature is very useful for coupling external codes to Quantum ESPRESSO routines. To this end, a general-purpose library for calling PWscf or CP from external codes (either Fortran or C/C++ using the Fortran 2003 ISO C binding standard) is provided in the directory COUPLE/.

III.4.3 Reorganization of linear-response codes

All linear-response codes described in Secs. II.2 and II.1.4 share as basic computational step the self-consistent solution of linear systems $Ax=b$ for different perturbations $b$, where the operator $A$ is derived from the KS Hamiltonian $H$ and the linear-response potential. Both the perturbations and the methods of solution differ by subtle details, leading to a plethora of routines, customized to solve slightly different versions of the same problem. Ideally, one should be able to solve any linear-response problem by using a suitable library of existing code. To this end, a major restructuring of linear-response codes has been started. Several routines have been unified, generalized and extended. They have been collected into the same subdirectory, LR_Modules, that will be the container of “generic” linear-response routines. Linear-response-related packages now contain only code that is specific to a given perturbation or property calculation.

III.5 Quantum ESPRESSO and scripting languages

A desirable feature of electronic-structure codes is the ability to be called from a high-level interpreted scripting language. Among the various alternatives, python has emerged in the last years due to its simple and powerful syntax and to the availability of numerical extensions (NumPy). Since v.6.0, an interface between PWscf and the path integral MD driver i-PI Ceriotti et al. (2014) is available and distributed together with Quantum ESPRESSO. Various implementations of an interface between Quantum ESPRESSO codes and the atomic simulation environment (ASE) Bahn and Jacobsen (2002) are also available. In the following we briefly highlight the integration of Quantum ESPRESSO with AiiDA, the pwtk toolkit for PWscf, and the QE-emacs-modes package for user-friendly editing of Quantum ESPRESSO with the Emacs editor Moon et al. (2017).

AiiDA Pizzi et al. (2016) is a comprehensive python infrastructure aimed at accelerating, simplifying, and organizing major efforts in computational science, and in particular computational materials science, with a close integration with the Quantum ESPRESSO distribution. AiiDA is structured around the four pillars of the ADES model (Automation, Data, Environment, and Sharing, Ref. 204)), and provides a practical and efficient implementation of all four. In particular, it aims at relieving the work of a computational scientist from the tedious and error-prone tasks of running, overseeing, and storing hundreds or more of calculations daily (Automation pillar), while ensuring that strict protocols are in place to store these calculations in an appropriately structured database that preserves the provenance of all computational steps (Data pillar). This way, the effort of a computational scientist can become focused on developing, curating, or exploiting complex workflows (Environment pillar) that calculate in a robust manner e.g. the desired materials properties of a given input structure, recording expertise in reproducible sequences that can be progressively perfected, while being able to share freely both the workflows and the data generated with public or private common repositories (Sharing). AiiDA is built using an agnostic structure that allows to interface it with any given code — through plugins and a plugin repository — or with different queuing systems, transports to remote HPC resources, and property calculators. In addition, it allows to use arbitrary object-relational mappers (ORMs) as backends (currently, Django and SQLAlchemy are supported). These ORMs map the AiiDA objects (“Codes”, “Calculations” and “Data”) onto python classes, and lead to the representation of calculations through Directed Acyclic Graphs (DAGs) connecting all objects with directional arrows; this ensures both provenance and reproducibility of a calculation. As an example, in Fig. 4 we present a simple DAG representing a PWscf calculation on BaTiO3.

III.5.2 Pwtk: a toolkit for PWscf

The pwtk, standing for PWscf ToolKit, is a Tcl scripting interface for PWscf set of programs contained in the Quantum ESPRESSO distribution. It aims at providing a flexible and productive framework. The basic philosophy of pwtk is to lower the learning curve by using syntax that closely matches the input syntax of Quantum ESPRESSO. Pwtk features include: (i) assignment of default values of input variables on a project basis, (ii) reassignment of input variables on the fly, (iii) stacking of input data, (iv) math-parser, (v) extensible and hierarchical configuration (global, project-based, local), (vi) data retrieval functions (i.e., either loading the data from pre-existing input files or retrieving the data from output files), and (vii) a few predefined higher-level tasks, that consist of several seamlessly integrated calculations. Pwtk allows to easily automate large number of calculations and to glue together different computational tasks, where output data of preceding calculations serve as input for subsequent calculations. Pwtk and related documentation can be downloaded from http://pwtk.quantum-espresso.org.

III.5.3 QE-emacs-modes

The QE-emacs-modes package is an open-source collection of Emacs major-modes for making the editing of Quantum ESPRESSO input files easier and more comfortable with Emacs. The package provides syntax highlighting (see Fig. 5a), auto-indentation, auto-completion, and a few utility commands, such as M-x prog-insert-template that inserts a respective input file template for the prog program (e.g., pw, neb, pp, projwfc, dos, bands). The QE-emacs-modes are aware of all namelists, variables, cards, and options that are explicitly documented in the INPUT_PROG.html files, which describe the respective input file syntax (see Fig. 5b), where PROG stands for the uppercase name of a given program of Quantum ESPRESSO. The reason for this is that both INPUT_PROG.html files and QE-emacs-modes are automatically generated by the internal helpdoc utility of Quantum ESPRESSO.

III.6 Continuous Integration and testing

The modularization of Quantum ESPRESSO reduces the extent of code duplication, thus improving code maintainability, but it also creates interdependencies between the modules so that changes to one part of the code may impact other parts. In order to monitor and mitigate these side effects we developed a test-suite for non-regression testing. Its purpose is to increase code stability by identifying and correcting those changes that break established functionalities. The test-suite relies on a modified version of python script testcode Spencer (2017).

The layout of the test-suite is illustrated in Fig. 6. The suite is invoked via a Makefile that accepts several options to run sequential or parallel tests or to test one particular feature of the code. The test-suite runs the various executables of Quantum ESPRESSO, extracts the numerical data of interest, compares them to reference data, and decides whether the test is successful using specified thresholds. At the moment, the test-suite contains 181 tests for PW, 14 for PH, 17 for CP, 43 for EPW, and 6 for TDDFpT covering 43%, 30%, 29%, 63% and 25% of the blocks, respectively. Moreover, 60%, 44%, 47%, 76% and 32% of the subroutines in each of these codes are tested, respectively.

The test-suite also contains the logic to automatically create reference data by running the relevant executables and storing the output in a benchmark file. These benchmarks are updated only when new tests are added or bugfixes modify the previous behavior.

The test-suite enables automatic testing of the code using several Buildbot test farms. The test farms monitor the code repository continuously, and trigger daily builds in the night after every new commit. Several compilers (Intel, GFortran, PGI) are tested both in serial and in parallel (openmpi, mpich, Intel mpi and mvapich2) execution with different mathematical libraries (LAPACK, BLAS, ScaLAPACK, FFTW3, MKL, OpenBlas). More information can be found at test-farm.quantum-espresso.org.

The official mirror of the development version of Quantum ESPRESSO (https://github.com/QEF/q-e) employs a subset of the test-suite to run Travis CI. This tool rapidly identifies erroneous commits and can be used to assist code review during a pull request.

IV Outlook and conclusions

This paper describes the core methodological developments and extensions of Quantum ESPRESSO that have become available, or are about to be released, after Ref. 6 appeared. The main goal of Quantum ESPRESSO to provide an efficient and extensible framework to perform simulations with well-established approaches and to develop new methods remains firm, and it has nurtured an ever growing community of developers and contributors.

Achieving such goal, however, becomes increasingly challenging. On one hand, computational methods become ever more complex and sophisticated, making it harder not only to implement them on a computer but also to verify the correctness of the implementation (for a much needed initial effort on verification of electronic-structure codes based on DFT, see Ref. Lejaeghere et al., 2016). On the other hand, exploiting the current technological innovations in computer hardware can requires massive changes to software and even algorithms. This is especially true for the case of “accelerated” architectures (GPUs and the like), whose exceptional performance can translate to actual calculations only after heavy restructuring and optimization. The complexity of existing codes makes a rewrite for new architectures a challenging choice, and a risky one given the fast evolution of computer architectures.

We think that the main directions followed until now in the development of Quantum ESPRESSO are still valid, not only for new methodologies, but also for adapting to new computer architectures and future “exascale” machines. Namely, we will continue pushing towards code reusability, by removing duplicated code and/or replacing it with routines performing well-defined tasks, by identifying the time-intensive sections of the code for machine-dependent optimization, by having documented APIs with a predictable behavior and with limited dependency upon global variables, and we will continue to optimize performance and reliability. Finally, we will push towards extended interoperability with other software, also in view of its usefulness for data exchange and for cross-verification, or to satisfy the needs of high-throughput calculations.

Still, the investment in the development and maintenance of state-of-the-art scientific software has historically lagged behind compared to the investment in the applications that use such software, and one wonders is this the correct or even forward-looking approach given the strategic importance of such tools, their impact, their powerful contribution to open science, and their full and complete availability to the entire community. In all of this, the future of materials simulations appear ever more brightMarzari (2016), and the usefulness and relevance of such tools to accelerating invention and discovery in science and technology is reflected in its massive uptake by the community at large.

Acknowledgments. This work has been partially funded by the European Union through the MaX Centre of Excellence (Grant No. 676598) and by the Quantum ESPRESSO Foundation. SdG acknowledges support from the EU Centre of Excellence E-CAM (Grant No. 676531). OA, MC, NC, NM, NLN, and IT acknowledge support from the SNSF National Centre of Competence in Research MARVEL, and from the PASC Platform for Advanced Scientific Computing. TT acknowledges support from NSF Grant No. DMR-1145968. SP, MS, and FG are supported by the Leverhulme Trust (Grant RL-2012-001). MBN acknowledges support by DOD-ONR (N00014-13-1-0635, N00014-11-1-0136, N00014-15-1-2863) and the Texas Advanced Computing Center at the University of Texas, Austin. RD acknowledges partial support from Cornell University through start-up funding and the Cornell Center for Materials Research (CCMR) with funding from the NSF MRSEC program (DMR-1120296). This research used resources of the Argonne Leadership Computing Facility at Argonne National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-06CH11357. This research also used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

References