Wannier90 as a community code: new features and applications

TL;DR

Wannier90 has evolved into a community-driven software with new features for Wannier functions, property calculations, performance, usability, and modern software practices, broadening its application scope in materials science.

Contribution

The paper introduces Wannier90 v3.0 with new functionalities, performance improvements, and community-driven development practices to enhance its capabilities and user experience.

Findings

Added symmetry-adapted and selectively-localised Wannier functions

Enabled calculation of shift currents and Berry-curvature dipole

Improved performance and usability features

Abstract

Wannier90 is an open-source computer program for calculating maximally-localised Wannier functions (MLWFs) from a set of Bloch states. It is interfaced to many widely used electronic-structure codes thanks to its independence from the basis sets representing these Bloch states. In the past few years the development of Wannier90 has transitioned to a community-driven model; this has resulted in a number of new developments that have been recently released in Wannier90 v3.0. In this article we describe these new functionalities, that include the implementation of new features for wannierisation and disentanglement (symmetry-adapted Wannier functions, selectively-localised Wannier functions, selected columns of the density matrix) and the ability to calculate new properties (shift currents and Berry-curvature dipole, and a new interface to many-body perturbation theory); performance improvements, including parallelisation of the core code; enhancements in functionality (support for spinor-valued Wannier functions, more accurate methods to interpolate quantities in the Brillouin zone); improved usability (improved plotting routines, integration with high-throughput automation frameworks), as well as the implementation of modern software engineering practices (unit testing, continuous integration, and automatic source-code documentation). These new features, capabilities, and code development model aim to further sustain and expand the community uptake and range of applicability, that nowadays spans complex and accurate dielectric, electronic, magnetic, optical, topological and transport properties of materials.

Tables

Table 1. Table 1: Wall-time for some of the runs performed with the Berry module, before ( Wannier90 v2.0) and after ( Wannier90 v3.0) the optimisations, for the test system described in the main text. N c subscript 𝑁 c N_{\text{c}} indicates the number of cores used in the calculation.

Mode	$k$-grid	$N_{\text{c}}$	Time (s)
version 3.0
Berry 3D	30$\times$30$\times$30	24	6903
	30$\times$30$\times$30	48	3527
	30$\times$30$\times$30	480	441
	100$\times$100$\times$100	480	13041
	100$\times$100$\times$100	7680	957
Berry 2D	100$\times$100	24	1389
Berry 1D	10000	24	12639
version 2.0
Berry 3D	30$\times$30$\times$30	24	56497
	30$\times$30$\times$30	48	40279

Equations

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

$$|w_{n\mathbf{R}}⟩=V{\int }_{\mathrm{BZ}}\frac{\mathrm{d}\mathbf{k}}{(2\pi {)}^3}{e}^{-i\mathbf{k}\cdot \mathbf{R}}\sum _{m=1}^J|{\psi }_{m\mathbf{k}}⟩U_{mn\mathbf{k}},$$

$$\Omega =\sum _{n=1}^J\left[⟨w_{n0}⟩\mathbf{r}\cdot \mathbf{r}{w}_{n0}-{\left|⟨w_{n0}⟩\mathbf{r}{w}_{n0}\right|}^2\right].$$

$$\Omega ={\Omega }_{\mathrm{I}}+\tilde{\Omega },$$

$${\Omega }_{\mathrm{I}}=\sum _n\left[⟨w_{n0}⟩\mathbf{r}\cdot \mathbf{r}{w}_{n0}-\sum _{m\mathbf{R}}{\left|⟨w_{m\mathbf{R}}⟩\mathbf{r}{w}_{n0}\right|}^2\right]$$

$$\tilde{\Omega }=\sum _n\sum _{m\mathbf{R}\ne n0}{\left|⟨w_{m\mathbf{R}}⟩\mathbf{r}{w}_{n0}\right|}^2$$

$${\Omega }_{\mathrm{I}}=\frac{1}{N}\sum _{\mathbf{k},\mathbf{b}}{w}_b\left[J-\sum _{mn}{\left|M_{mn}^{(\mathbf{k},\mathbf{b})}\right|}^2\right]$$

$$\tilde{\Omega }$$

$$+\frac{1}{N}\sum _{\mathbf{k},\mathbf{b}}{w}_b\sum _n(-\text{Im}\ln M_{nn}^{(\mathbf{k},\mathbf{b})}-\mathbf{b}\cdot {\bar{\mathbf{r}}}_n{)}^2,$$

$$M_{mn}^{(\mathbf{k},\mathbf{b})}=⟨u_{m\mathbf{k}}⟩u_{n,\mathbf{k}+\mathbf{b}},$$

$${\bar{\mathbf{r}}}_n\equiv ⟨w_{n0}⟩\mathbf{r}{w}_{n0}=-\frac{1}{N}\sum _{\mathbf{k},\mathbf{b}}{w}_b\mathbf{b}\text{Im}\ln M_{nn}^{(\mathbf{k},\mathbf{b})}.$$

$${\mathcal{G}}_{\mathbf{k}}\equiv \frac{\mathrm{d}\Omega }{\mathrm{d}{W}_{mn\mathbf{k}}}=4\sum _{\mathbf{b}}{w}_b\left(\mathcal{A}[R_{mn}^{(\mathbf{k},\mathbf{b})}]-\mathcal{S}[T_{mn}^{(\mathbf{k},\mathbf{b})}]\right),$$

$$R_{mn}^{(\mathbf{k},\mathbf{b})}$$

$$T_{mn}^{(\mathbf{k},\mathbf{b})}$$

$$q_n^{(\mathbf{k},\mathbf{b})}$$

$$U_{\mathbf{k}}\to U_{\mathbf{k}}\text{exp}[\alpha {\mathcal{D}}_{\mathbf{k}}],$$

$$M^{(\mathbf{k},\mathbf{b})}=U_{\mathbf{k}}^{†}{M}^{(0)(\mathbf{k},\mathbf{b})}{U}_{\mathbf{k}+\mathbf{b}}^{},$$

$$M_{mn}^{(0)(\mathbf{k},\mathbf{b})}=⟨u_{m\mathbf{k}}^{(0)}⟩u_{n,\mathbf{k}+\mathbf{b}}^{(0)}$$

$$|{\tilde{\psi }}_{n\mathbf{k}}⟩=\sum _{m=1}^{{\mathcal{J}}_{\mathbf{k}}}|{\psi }_{m\mathbf{k}}⟩V_{mn\mathbf{k}},$$

$${\Omega }_{\mathrm{I}}=\frac{1}{N}\sum _{\mathbf{k},\mathbf{b}}{w}_b\mathrm{Tr}[P_{\mathbf{k}}{Q}_{\mathbf{k}+\mathbf{b}}],$$

$$\left[\sum _{\mathbf{b}}{w}_b{P}_{\mathbf{k}+\mathbf{b}}\right]|{\tilde{u}}_{n\mathbf{k}}⟩={\lambda }_{n\mathbf{k}}|{\tilde{u}}_{n\mathbf{k}}⟩.$$

$$\left[\sum _{\mathbf{b}}{w}_b{P}_{\mathbf{k}+\mathbf{b}}^{(i-1)}\right]|{\tilde{u}}_{n\mathbf{k}}^{(i)}⟩={\lambda }_{n\mathbf{k}}^{(i)}|{\tilde{u}}_{n\mathbf{k}}^{(i)}⟩,$$

$$[P_{\mathbf{k}+\mathbf{b}}^{(i)}{]}_{\mathrm{in}}=\beta P_{\mathbf{k}+\mathbf{b}}^{(i-1)}+(1-\beta )[P_{\mathbf{k}+\mathbf{b}}^{(i-1)}{]}_{\mathrm{in}},$$

$$Z_{mn\mathbf{k}}^{(i)}=⟨u_{m\mathbf{k}}^{(0)}⟩\sum _{\mathbf{b}}{w}_b[P_{\mathbf{k}+\mathbf{b}}^{(i)}{]}_{\mathrm{in}}{u}_{n\mathbf{k}}^{(0)}.$$

$$|{\varphi }_{n\mathbf{k}}⟩=\sum _{m=1}^{J\mathrm{or}{\mathcal{J}}_{\mathbf{k}}}|{\psi }_{m\mathbf{k}}⟩⟨{\psi }_{m\mathbf{k}}⟩g_n,$$

$$|{\tilde{\psi }}_{n\mathbf{k}}⟩$$

$$=$$

$$\{w_i^{(\varrho )}(\mathbf{r}-{\mathbf{q}}_a)=w_{ia}^{(\varrho )}(\mathbf{r}),i=1,\dots ,n_{\varrho }\}$$

$$w_{ia}^{(\varrho )}(\mathbf{r}-\mathbf{R})$$

$$g{\psi }_{n\mathbf{k}}(\mathbf{r})=\sum _{m=1}^J{\psi }_{m\mathcal{R}\mathbf{k}}(\mathbf{r}){\tilde{d}}_{mn\mathbf{k}}(g),g\in G.$$

$${\tilde{d}}_{mn\mathbf{k}}(g)=\int \mathrm{d}\mathbf{r}{\psi }_{m\mathcal{R}\mathbf{k}}^{\ast }(\mathbf{r}){\psi }_{n\mathbf{k}}\left(g^{-1}\mathbf{r}\right).$$

$$g{\psi }_{ia\mathbf{k}}^{(\varrho )}(\mathbf{r})$$

Full text

I Introduction

Wannier90 is an open-source code for generating Wannier functions (WFs), in particular maximally-localised Wannier functions (MLWFs), and using them to compute advanced materials properties with high efficiency and accuracy. Wannier90 is a paradigmatic example of interoperable software, achieved by ensuring that all the quantities required as input are entirely independent of the underlying electronic-structure code from which they are obtained. Most of the major and widely used electronic-structure codes have an interface to Wannier90, including Quantum ESPRESSOgiannozzi_qe_2017 , ABINITgonze-cpc180 , VASPkresse-prb47 ; kresse-cms6 ; kresse-prb54 , Siestasoler-condmatt14 , Wien2kblaha-wien2k , Fleurblugel-fleur and Octopusoctopus2015 . As a consequence, once a property is implemented within Wannier90, it can be immediately available to users of all codes that interface to it.

Over the last few years, Wannier90 has undergone a transition from a code developed by a small group of developers to a community code with a much wider developers’ base. This has been achieved in two principal ways: (i) hosting the source code and associated development efforts on a public GitHub repositoryW90repo ; and (ii) building a community of Wannier90 developers and facilitating personal interactions between individuals through community workshops, the most recent in 2016. In response, the code has grown significantly, gaining many novel features contributed by this community, as well as numerous fixes.

In this paper, we describe the most important novel contributions to the Wannier90 code, as embodied in its 3.0 release. The paper is structured as follows: In Sec. II we first summarise the background theory for the computation of MLWFs (additional details can be found in Ref. Marzari-RMP2012, ), and introduce the notation that will be used throughout the paper. In Sec. III we describe the novel features of Wannier90 that are related to the core wannierisation and disentanglement algorithms; these include symmetry-adapted WFs, selective localisation of WFs, and parallelisation using the message-passing interface (MPI). In Sec. IV we describe new functionality enhancements, including the ability to handle spinor-valued WFs and calculations with non-collinear spin that use ultrasoft pseudopotentials (within Quantum ESPRESSO); improved interpolation of the $k$-space Hamiltonian; a more flexible approach for handling and using initial projections; and the ability to plot WFs in Gaussian cube format on WF-centred grids with non-orthogonal translation vectors. In Sec. V we describe new functionalities associated with using MLWFs for computing advanced electronic-structure properties, including the calculation of shift currents, gyrotropic effects and spin Hall conductivities, as well as parallelisation improvements and the interpolation of bands originating from calculations performed with many-body perturbation theory (GW). In Sec. VI we describe the selected-columns-of-the-density-matrix (SCDM) method, which enables computation of WFs without the need for explicitly defining initial projections. In Sec. VII we describe new post-processing tools and codes, and the integration of Wannier90 with high-throughput automation and workflow management tools (specifically, the AiiDA materials’ informatics infrastructurepizzi-AiiDA ). In Sec. VIII we describe the modern software engineering practices now adopted in Wannier90, that have made it possible to improve the development lifecycle and transform Wannier90 into a community-driven code. Finally, our conclusions and outlook are presented in Sec. IX.

II Background

In the independent-particle approximation, the electronic structure of a periodic system is conventionally represented in terms of one-electron Bloch states $\psi_{n\mathbf{k}}(\mathbf{r})$, which are labelled by a band index $n$ and a crystal momentum $\mathbf{k}$ inside the first Brillouin zone (BZ), and which satisfy Bloch’s theorem:

[TABLE]

where $u_{n\mathbf{k}}(\mathbf{r})=u_{n\mathbf{k}}(\mathbf{r}+\mathbf{R})$ is a periodic function with the same periodicity of the single-particle Hamiltonian, and $\mathbf{R}$ is a Bravais lattice vector. (For the moment we ignore the spin degrees of freedom and work with spinless wave functions; spinor wave functions will be treated in Sec. IV.1.) Such a formalism is also commonly applied, via the supercell approximation, to non-periodic systems, typically used to treat point, line and planar defects in crystals, surfaces, amorphous solids, liquids and molecules.

II.1 Isolated bands

A group of bands is said to be isolated if it is separated by energy gaps from all the other lower and higher bands throughout the BZ (this isolated group of bands may still show arbitrary crossing degeneracies and hybridisations within itself). For such isolated set of $J$ bands, the electronic states can be equivalently represented by a set of $J$ WFs per cell, that are related to the Bloch states via two unitary transformations (one continuous, one discrete) wannier-pr37 :

[TABLE]

where $w_{n\mathbf{R}}(\mathbf{r})=w_{n{\bf 0}}(\mathbf{r}-\mathbf{R})$ is a periodic (but not necessarily localised) WF labelled by the quantum number $\mathbf{R}$ (the counterpart of the quasi-momentum $\mathbf{k}$ in the Bloch representation), $V$ is the cell volume and $U_{\mathbf{k}}$ are unitary matrices that mix Bloch states at a given $\mathbf{k}$ and represent the gauge freedom that exists in the definition of the Bloch states and that is inherited by the WFs.

MLWFs are obtained by choosing $U_{\mathbf{k}}$ matrices that minimise the sum of the quadratic spreads of the WFs about their centres for a reference $\mathbf{R}$ (say, $\mathbf{R}=\bm{0}$). This sum is given by the spread functional

[TABLE]

$\Omega$ may be decomposed into two positive-definite partsMV_PRB56 ,

[TABLE]

where

[TABLE]

is gauge invariant (i.e., invariant under the action of any unitary $U_{\mathbf{k}}$ on the Bloch states), and

[TABLE]

is gauge dependent. Therefore, the “wannierisation” of an isolated manifold of bands, i.e., the transformation of Bloch states into MLWFs, amounts to minimising the gauge-dependent part $\widetilde{\Omega}$ of the spread functional.

Crucially, the matrix elements of the position operator between WFs can be expressed in reciprocal space. Under the assumption that the BZ is sampled on a uniform Monkhorst–Pack mesh of $k$-points composed of $N$ points ($V\int_{\mbox{\tiny{BZ}}}\frac{\mathrm{d}\mathbf{k}}{(2\pi)^{3}}\rightarrow\frac{1}{N}\sum_{\mathbf{k}}$), the gauge-independent and gauge-dependent parts of the spread may be expressed, respectively, asMV_PRB56

[TABLE]

and

[TABLE]

where $\mathbf{b}$ are the vectors connecting a $k$-point to its neighbours, $w_{b}$ are weights associated with the finite-difference representation of $\nabla_{\mathbf{k}}$ for a given geometry, the matrix of overlaps $M^{(\mathbf{k},\mathbf{b})}$ is defined by

[TABLE]

and the centres of the WFs are given by

[TABLE]

Minimisation of the spread functional is achieved by considering infinitesimal gauge transformations $U_{mn\mathbf{k}}=\delta_{mn}+\mathrm{d}W_{mn\mathbf{k}}$, where $\mathrm{d}W$ is anti-Hermitian ($\mathrm{d}W^{\dagger}=-\mathrm{d}W$). The gradient of the spread functional with respect to such variations is given by

[TABLE]

where $\mathcal{A}$ and $\mathcal{S}$ are the super-operators $\mathcal{A}[B]=(B-B^{\dagger})/2$ and $\mathcal{S}[B]=(B+B^{\dagger})/2i$, respectively, and

[TABLE]

For the full derivation of Eq. (11) we refer to Ref. [MV_PRB56, ]. This gradient is then used to generate a search direction $\mathcal{D}_{\mathbf{k}}$ for an iterative steepest-descent or conjugate-gradient minimisation of the spread Mostofi_CPC : at each iteration the unitary matrices are updated according to

[TABLE]

where $\alpha$ is a coefficient that can either be set to a fixed value or determined at each iteration via a simple polynomial line-search, and the matrix exponential is computed in the diagonal representation of $\mathcal{D}_{\mathbf{k}}$ and then transformed back in the original representation. Once the unitary matrices have been updated, the updated set of $M^{(\mathbf{k},\mathbf{b})}$ matrices is calculated according to

[TABLE]

where

[TABLE]

is the set of initial $M^{(\mathbf{k},\mathbf{b})}$ matrices, computed once and for all, at the start of the calculation, from the original set of reference Bloch orbitals $\ket{u^{(0)}_{n\mathbf{k}}}$.

II.2 Entangled bands

It is often the case that the bands of interest are not separated from other bands in the Brillouin zone by energy gaps and are overlapping and hybridising with other bands that extend beyond the energy range of interest. In such cases, we refer to the bands as being entangled.

The difficulty in constructing MLWFs for entangled bands arises from the fact that, within a given energy window, the number of bands $\mathcal{J}_{\mathbf{k}}$ at each $k$-point $\mathbf{k}$ in the BZ is not a constant and is, in general, different from the target number $J$ of WFs: $\mathcal{J}_{\mathbf{k}}\geq J$. Even making the energy window $k$-dependent would see discontinuous inclusion and exclusion of bands as the BZ is traversed. The treatment of entangled bands requires thus a more complex approach that is typically a two-step process. In the first step, a $J$-dimensional manifold of Bloch states is selected at each $k$-point, chosen to be as smooth as possible as a function of $\mathbf{k}$. In the second step, the gauge freedom associated with the selected manifold is used to obtain MLWFs, just as described in Sec. II.1 for the case of an isolated set of bands.

Focusing on the first step, an orthonormal basis for the $J$-dimensional subspace $\mathcal{S}_{\mathbf{k}}$ at each $\mathbf{k}$ can be obtained by performing a semi-unitary transformation on the $\mathcal{J}_{\mathbf{k}}$ states at $\mathbf{k}$,

[TABLE]

where $V_{\mathbf{k}}$ is a rectangular matrix of dimension $\mathcal{J}_{\mathbf{k}}\times J$ that is semi-unitary in the sense that $V_{\mathbf{k}}^{\dagger}V_{\mathbf{k}}^{\phantom{\dagger}}=\mathbf{1}$.

To select the smoothest possible manifold, a measure of the intrinsic smoothness of the chosen subspace is needed. It turns out that such a measure is given precisely by what was the gauge-invariant part $\Omega_{\mathrm{I}}$ of the spread functional for isolated bands.SMV_PRB65 Indeed, Eq. (7) can be expressed as

[TABLE]

where ${P}_{\mathbf{k}}=\sum_{n=1}^{J}\Ket{\widetilde{u}_{n\mathbf{k}}}\Bra{\widetilde{u}_{n\mathbf{k}}}$ is the projection operator onto $\mathcal{S}_{\mathbf{k}}$, ${Q}_{\mathbf{k}}={\mathbf{1}}-{P}_{\mathbf{k}}$ is its Hilbert-space complement, and “$\mathrm{Tr}$” represents the trace over the entire Hilbert space. $\mathrm{Tr}[{P}_{\mathbf{k}}{Q}_{\mathbf{k}+\mathbf{b}}]$ measures the mismatch between the subspaces $\mathcal{S}_{\mathbf{k}}$ and $\mathcal{S}_{\mathbf{k}+\mathbf{b}}$, vanishing if they overlap identically. Hence $\Omega_{\mathrm{I}}$ measures the average mismatch of the local subspace $\mathcal{S}_{\mathbf{k}}$ across the BZ, so that an optimally-smooth subspace can be selected by minimising $\Omega_{\mathrm{I}}$. Doing this with orthonormality constraints on the Bloch-like states is equivalent to solving self-consistently the set of coupled eigenvalue equationsSMV_PRB65

[TABLE]

The solution can be achieved via an iterative procedure, whereby at the $i^{\mathrm{th}}$ iteration the algorithm traverses the entire set of $k$-points, selecting at each one the $J$-dimensional subspace $\mathcal{S}^{(i)}_{\mathbf{k}}$ that has the smallest mismatch with the subspaces $\mathcal{S}^{(i-1)}_{\mathbf{k}+\mathbf{b}}$ at the neighbouring $k$-points obtained in the previous iteration. This amounts to solving

[TABLE]

and selecting the $J$ eigenvectors with the largest eigenvaluesSMV_PRB65 . Self-consistency is reached when $\mathcal{S}^{(i)}_{\mathbf{k}}=\mathcal{S}^{(i-1)}_{\mathbf{k}}$ (to within a user-defined threshold dis_conv_tol) at all the $k$-points. To make the algorithm more robust, the projector appearing on the left-hand-side of Eq. (21) is replaced with $[{P}_{\mathbf{k}+\mathbf{b}}^{(i)}]_{\mathrm{in}}$, given by

[TABLE]

which is a linear mixture of the projector that was used as input for the previous iteration and the projector defined by the output of the previous iteration. The parameter $0<\beta\leq 1$ determines the degree of mixing, and is typically set to $\beta=0.5$; setting $\beta=1$ reverts precisely to Eq. (21), while smaller and smaller values of $\beta$ make convergence smoother (and thus more robust) but also slower.

In practice, Eq. (21) is solved by diagonalising the Hermitian operator appearing on the left-hand-side in the basis of the original $\mathcal{J}_{\mathbf{k}}$ Bloch states:

[TABLE]

Once the optimal subspace has been selected, the wannierisation procedure described in Sec. II.1 is carried out to minimise the gauge-dependent part $\widetilde{\Omega}$ of the spread functional within that optimal subspace.

II.3 Initial projections

In principle, the overlap matrix elements $M_{mn}^{(\mathbf{k,b})}$ are the only quantities required to compute and minimise the spread functional, and generate MLWFs for either isolated or entangled bands. In practice, this is generally true when dealing with an isolated set of bands, but in the case of entangled bands a good initial guess for the subspaces $\mathcal{S}_{\mathbf{k}}$ alleviates problems associated with falling into local minima of $\Omega_{\mathrm{I}}$, and/or obtaining MLWFs that cannot be chosen to be real-valued (in the case of spinless WFs). Even in the case of an isolated set of bands, a good initial guess for the WFs, whilst not usually critical, often results in faster convergence of the spread to the global minimum. (It is important to note that both for isolated and for entangled bands multiple solutions to the wannierisation or disentanglement can exist, as discussed later.)

A simple and effective procedure for selecting an initial gauge (in the case of isolated bands) or an initial subspace and initial gauge (in the case of entangled bands) is to project a set of $J$ trial orbitals $g_{n}(\mathbf{r})$ localised in real space onto the space spanned by the set of original Bloch states at each $\mathbf{k}$:

[TABLE]

where the sum runs up to either $J$ or $\mathcal{J}_{\mathbf{k}}$, depending on whether the bands are isolated or entangled, respectively, and the inner product $A_{mn\mathbf{k}}=\Braket{\psi_{m\mathbf{k}}}{g_{n}}$ is over all the Born–von Karman supercell. (In practice, the fact that the $g_{n}$ are localised greatly simplifies this calculation.) The matrices $A_{\mathbf{k}}$ are square $(J\times J)$ or rectangular $(\mathcal{J}_{k}\times J)$ in the case of isolated or entangled bands, respectively. The resulting orbitals are then orthonormalised via a Löwdin transformation lowdin1950 :

[TABLE]

where $S_{mn\mathbf{k}}={\Braket{\phi_{m\mathbf{k}}}{\phi_{n\mathbf{k}}}}=(A_{\mathbf{k}}^{\dagger}A_{\mathbf{k}}^{\phantom{\dagger}})_{mn}$, and $A_{\mathbf{k}}^{\phantom{}}S_{\mathbf{k}}^{-\nicefrac{{1}}{{2}}}$ is a unitary or semi-unitary matrix. In the case of entangled bands, once an optimally-smooth subspace has been obtained as described in Sec. II.2, the same trial orbitals $g_{n}(\mathbf{r})$ can be used to initialise the wannierisation procedure of Sec. II.1. In practice, the matrices $A_{\mathbf{k}}$ are computed once and for all at the start of the calculation, together with the overlap matrices $M^{(\mathbf{k},\mathbf{b})}$. These two operations need to be performed within the context of the electronic-structure code and basis set adopted; afterwards, all the operations of Wannier90 rely only on $A_{\mathbf{k}}$ and $M^{(\mathbf{k},\mathbf{b})}$ and not on the specific representation of $\psi_{m\mathbf{k}}$ (e.g., plane waves, linearised augmented plane waves, localised basis sets, real-space grids, …).

III New features for wannierisation and disentanglement

In this section we provide an overview of the new features associated with the core wannierisation and disentanglement algorithms in Wannier90, namely the ability to generate WFs of specific symmetry; selectively localise a subset of the WFs and/or constrain their centres to specific sites; and perform wannierisation and disentanglement more efficiently through parallelisation.

III.1 Symmetry-adapted Wannier functions

In periodic systems, atoms are usually found at sites $\mathbf{q}$ whose site-symmetry group $G_{q}$ is a subgroup of the full point group $F$ of the crystalevarestov2012site (the symmetry operations in the group $G_{q}$ are those that leave $\mathbf{q}$ fixed). The set of points $\mathbf{q}_{a}$ that are symmetry-equivalent sites to $\mathbf{q}$ is called an orbitInt_tab_cryst . These are all the points in the unit cell that can be generated from $\mathbf{q}$ by applying the symmetry operations in $G$ that do not leave $\mathbf{q}$ fixed. If $\mathbf{q}_{a}$ is a high-symmetry site then its Wyckoff position has a single orbitInt_tab_cryst ; for low-symmetry sites different orbits correspond to the same Wyckoff position. The number of points in the orbit(s) is the multiplicity $n_{q_{a}}$ of the Wyckoff position. MLWFs, however, are not bound to reside on such high-symmetry sites, and they do not necessarily possess the site symmetries of the crystal SMV_PRB65 ; Sakuma_PRB87 ; Thygesen2005 . When using MLWFs as a local orbital basis set in methods such as first-principles tight binding, DFT+U and DFT plus dynamical-mean-field theory (DMFT), which deal with beyond-DFT correlations in a local subspace such as that spanned by $3d$ orbitals (for transition metals or transition-metal oxides) or $4f$ orbitals (for rare-earth or actinide intermetallics), it is often desirable to ensure that the WFs basis possesses the local site symmetries.

SakumaSakuma_PRB87 has shown that such symmetry-adapted Wannier functions (SAWFs) can be constructed by introducing additional constraints on the unitary matrices $U_{\mathbf{k}}$ of Eq. (2) during the minimisation of the spread. SAWFs, therefore, can be fully integrated within the original maximal-localisation procedure. The SAWF approach gives the user a certain degree of control over the symmetry and centres of the Wannier functions at the expense of some localisation since the final total spread of the resulting SAWFs can only be equal to, or most often larger than, that of the corresponding MLWFs with no constraints (note that in principle some SAWFs can have a smaller individual spread than any MLWFs).

A set of SAWFs

[TABLE]

can be specified by one representative point $\mathbf{q}_{a}$ of the orbit (in the home unit cell $\mathbf{R}=\mathbf{0}$), and by the irreducible representation (irrep) $\varrho$ of the corresponding site-symmetry group $G_{a}$ (the dimension of the irrep being $n_{\varrho}$). For instance, in simple fcc crystals, like copper (space group $Fm{-}3m$), the Cu atom is placed at the origin of the unit cell $\mathbf{q}=(0,0,0)$ (Wyckoff letter $a$ with multiplicity 1, i.e., only one point in the orbit of $(0,0,0)$ in the unit cell, due to the fact that $Fm{-}3m$ is symmorphicInt_tab_cryst ), whose site-symmetry group is $m{-}3m$ (also referred to as $O_{h}$). One of the irreps of $O_{h}$ is that with character $T_{2g}$, which is 3-dimensional.

To find symmetry-adapted Wannier functions, one needs to specify the unitary transformations $U^{(\varrho)}_{mia\mathbf{k}}$ of the Bloch states, defined by

[TABLE]

Therefore, the goal is to construct basis functions of the irrep $\varrho$, $\{\psi_{ia\mathbf{k}}^{(\varrho)}(\mathbf{r})\}$ , from a linear combination of the $J$ eigenstates $\psi_{n\mathbf{k}}(\mathbf{r})$ of the Hamiltonian $H$. Since $H$ is invariant under the full space-group $G$, the representation of a given symmetry operation $g=(\mathcal{R}|\mathbf{t})\in G$ (where $\mathcal{R}$ and $\mathbf{t}$ are the rotation and fractional-translation parts of the symmetry operation, respectively) in the $\{\psi_{n\mathbf{k}}(\mathbf{r})\}$ basis must be a $J\times J$ unitary matrixevarestov2012site $\widetilde{d}_{\mathbf{k}}(g)$, i.e. $\widetilde{d}_{\mathbf{k}}(g)$ represents how the $J$ Bloch states are transformed by the symmetry operation $g$:

[TABLE]

When single-particle eigenfunctions $\psi_{n\mathbf{k}}(\mathbf{r})$ are used, as in this case, the matrix elements $\widetilde{d}_{\mathbf{k}}(g)$ can be computed directly as

[TABLE]

As for the overlap matrices $M^{(\mathbf{k},\mathbf{b})}$, the resulting $\widetilde{d}_{\mathbf{k}}(g)$ are also basis-set independent. Moreover, once computed (using the original gauge) they remain fixed during the calculation. For instance, in a plane-wave code the integrals in Eq. (30) can be easily computed in reciprocal space.

On the other hand, it can be shown that the SAWFs transform under the action of $g\in G$ with a different matrix $d^{(\varrho)}_{\mathbf{k}}(g)$ Sakuma_PRB87 ; evarestov2012site , which in turn defines how the $\psi_{ia\mathbf{k}}^{(\varrho)}(\mathbf{r})$ transform under the action of $g\in G$:

[TABLE]

where

[TABLE]

and

[TABLE]

is a lattice translation vector. It is worth to mention that ${a^{\prime}}$ in Eq. (31) is uniquely defined by specifying the symmetry operations $g\in G$; see Ref. evarestov2012site, ; Sakuma_PRB87, for details.

$D_{\mathbf{k}}(g)$ is block-diagonal in the $\varrho$ index. For a given set of Wyckoff positions, the number of blocks is given by the sum of the number of all irreps considered (if non-equivalent Wyckoff positions ($\mathbf{q}_{b}\neq\mathbf{q}_{a}$) are present then $D_{\mathbf{k}}(g)$ contains also blocks corresponding to these positions). Each block contains $n_{q_{a}}^{2}$ sub-matrices of dimension $n_{\varrho}\times n_{\varrho}$. Therefore, if only one Wyckoff position is given with multiplicity $n_{q_{a}}$, then there are $J=n_{q_{a}}n_{\varrho}$ energy bands in the representation given by the $\{\psi_{ia\mathbf{k}}^{(\varrho)}(\mathbf{r})\}$ (see Ref. Sakuma_PRB87, for full details).

To compute the $D_{\mathbf{k}}(g)$ matrices one needs to specify the centre $\mathbf{q}_{a}$ and the symmetry of the initial functions (e.g., $s,p,$ and $d$). Then, for each symmetry operation $g_{q_{a}}$ in the site-symmetry group $G_{q_{a}}$ one needs to calculate the matrix representation of the rotational part expressed in the basis of these functions.

From Eqs. (28), (29) and (31) one can show that the following relationship exists between $U_{\mathbf{k}}$ and $U_{\mathcal{R}\mathbf{k}}$

[TABLE]

Let $g_{\mathbf{k}}$ now be the symmetry operations that leave a given $\mathbf{k}$ unchanged. Then, Eq. (34) gives the condition that $U_{\mathbf{k}}$ must satisfy in this case:

[TABLE]

The initial unitary matrix $U_{\mathbf{k}}$ ($\mathbf{k}\in$ IBZ) must satisfy the constraints in Eq. (35); this can be done in an iterative fashion, as discussed in Ref. Sakuma_PRB87, . In practice, the Wannier functions are generated from a limited subspace spanned by a finite number of states inside a target “energy window”, but this does not guarantee that a $U_{\mathbf{k}}$ can be constructed for any desired irrep. In fact, if a given irrep is not compatible with the symmetry of the states within the energy window, Eq. (35) cannot be fulfilled.

For an isolated set of bands, the minimisation of $\widetilde{\Omega}$ with the constraints defined in Eq. (34) requires the gradient $\mathcal{G}^{\mbox{\scriptsize{sym}}}_{\mathbf{k}}$ of the total spread $\Omega$ with respect to a symmetry-adapted gauge variation to generate a search direction $\mathcal{D}^{\mbox{\scriptsize{sym}}}_{\mathbf{k}}$. The equation for the symmetry-adapted gradient reads

[TABLE]

where $\mathcal{G}_{\mathbf{k}}$ is the original gradient given in Eq. (11), and $n_{\mathbf{k}}$ is the number of symmetry operations in $G$ that leave $\mathbf{k}$ fixed.

The procedure described above for isolated bands has to be modified only slightly for the case of entangled bands. The main difference with respect to the unconstrained case of Sec. II is that the $J$ eigenvectors of the $J$ largest eigenvalues of the $Z_{\mathbf{k}}$ matrix in Eq. (23) do not necessarily span the same subspace spanned by the desired symmetry-adapted Wannier functions. Since direct minimisation is not bound to give symmetry-adapted WFs, SakumaSakuma_PRB87 has proposed an alternative steepest-descent approach to construct the optimal unitary matrix from a set of $\mathcal{J}_{\mathbf{k}}\geq J$ Bloch wavefunctions that also fulfil symmetry constraints. Once this step is completed and optimal symmetry-adapted Bloch functions $\widetilde{\psi}_{ia\mathbf{k}}^{(\varrho)}(\mathbf{r})$ have been computed, the algorithm proceeds as in the isolated case where one seeks the $U_{\mathbf{k}}$ that minimise $\widetilde{\Omega}$ and give the symmetry-adapted Bloch functions in terms of $\widetilde{\psi}_{ia\mathbf{k}}^{(\varrho)}(\mathbf{r})$ as in Eq. (28).

III.2 Selectively-localised Wannier functions and constrained Wannier centres

Wang et al. have proposed an alternative methodMarianetti_PRB_90 to the symmetry-adapted Wannier functions described in Section III.1. Their method permits the selective localisation of a subset of the Wannier functions, which may optionally be constrained to have specified centres. Whilst this method does not enforce or guarantee symmetry constraints, it has been observed in the cases that have been studiedMarianetti_PRB_90 that Wannier functions whose centres are constrained to a specific site typically possess the corresponding site symmetries.

For an isolated set of $J$ bands, selective localisation of a subset of $J^{\prime}\leq J$ Wannier functions is accomplished by minimising the total spread $\Omega$ with respect to only $J^{\prime}\times J^{\prime}$ degrees of freedom in the unitary matrix $U_{\mathbf{k}}$. The spread functional to minimise is then given by

[TABLE]

which reduces to the original spread functional $\Omega$ of Eq. (3) for $J^{\prime}=J$. When $J^{\prime}<J$, it is no longer possible to cast the functional $\Omega^{\prime}$ as a sum of a gauge-independent term $\Omega_{\mathrm{I}}$ and gauge-dependent one $\widetilde{\Omega}$, as done in Eq. (4) for $\Omega$. Nevertheless, the minimisation can be carried out with methods very similar to those described in Section II. In fact, for $J^{\prime}<J$, $\Omega^{\prime}$ can be written as the sum of two gauge-dependent terms, $\Omega^{\prime}=\Omega_{\mathrm{IOD}}+\Omega_{\mathrm{D}}$, where $\Omega_{\mathrm{IOD}}$ is formally given by the sum of $\Omega_{\mathrm{I}}$ and the off-diagonal term $(m\neq n),\quad m,n\leq J^{\prime}<J$ of $\widetilde{\Omega}$, and $\Omega_{\mathrm{D}}$ by the diagonal term $(m=n)$ of $\widetilde{\Omega}$. If one adopts the usual discrete representation on a uniform Monkhorst–Pack grid of $k$-points, $\Omega_{\mathrm{IOD}}$ and $\Omega_{\mathrm{D}}$ are given byMarianetti_PRB_90

[TABLE]

and

[TABLE]

With this new spread functional, we can mimic the procedure used to obtain a set of MLWFs, and derive the gradient $\mathcal{G}_{\mathbf{k}}^{\prime}$ of $\Omega^{\prime}$ which gives the search direction to be used in the minimisation. The matrix elements of $\mathcal{G}_{\mathbf{k}}^{\prime}$ read

[TABLE]

where $\mathcal{G}_{mn\mathbf{k}}$ are the matrix elements of the original gradient in Eq. (11) (see also Ref. MV_PRB56, ), and $R_{mn}^{(\mathbf{k},\mathbf{b})}$ and $T_{mn}^{(\mathbf{k},\mathbf{b})}$ are given by Eq. (12) and Eq. (13), respectively. As a result of the minimisation, we obtain a set of $J^{\prime}$ maximally-localised Wannier functions, known as selectively-localised Wannier functions (SLWFs), whose spreads are in general smaller than the corresponding MLWFs. Naturally, the remaining $J-J^{\prime}$ functions will be more delocalised than their MLWF counterparts, as they are not optimised, and the overall sum of spreads will be larger (or in the best case scenario equal).

The centres of the SLWFs may be constrained by adding a quadratic penalty function to the spread functional $\Omega^{\prime}$, defining a new functional given by

[TABLE]

where $\lambda$ is a Lagrange multiplier and $\mathbf{x}_{n}$ is the desired centre for the $n^{\rm th}$ WF. The procedure outlined above for minimising $\Omega^{\prime}$ can be also adapted to deal with $\Omega^{\prime}_{\lambda}$ (see Ref. Marianetti_PRB_90, for details), and minimising $\Omega^{\prime}_{\lambda}$ results in selectively-localised Wannier functions subject to the constraint of fixed centres (SLWF+C). As noted above, it is observed that WFs derived using the SLWF+C approach naturally possess site symmetries, and their individual spreads are usually smaller than the corresponding spreads of MLWFs, although the total spread, combination of the $J^{\prime}$ selectively optimised WFs and the $J-J^{\prime}$ unoptimised functions, is larger than the total spread of the MLWFs (see, for instance last column in Tab. 1).

In the case of entangled bands, the SLWF(+C) method implicitly assumes that a subspace selection has been performed, i.e., that a smooth $J$-dimensional manifold exists. Since for the $\Omega^{\prime}$ and $\Omega^{\prime}_{\lambda}$ functionals it is not possible to define an $\Omega_{I}$ that measures the intrinsic smoothness of the underlying manifold, the additional constraints in Eq. (37) and Eq. (III.2) can only be imposed during the wannierisation step. This means that SLWF(+C) can be seamlessly coupled with the disentanglement procedure, with no further additions to the original procedure of Sec. II.2.

SAWF and SLWF+C in GaAs

As an example of the capabilities of the SAWF and SLWF+C approaches, we show how to construct atom-centred WFs that possess the local site symmetries in gallium arsenide (GaAs). In particular, we discuss how to obtain one WF from the four valence bands of GaAs that is centred on the As atom and that transforms like the identity under the symmetry operations in $T_{d}$, the site-symmetry group of the As site (for completeness, we also show one MLWF and one SLWF without constraints). Since we only deal with the four valence bands of GaAs—an isolated manifold—no prior subspace selection is required for the wannierisation. All calculations were carried out with the plane-wave DFT code Quantum ESPRESSOgiannozzi_qe_2017 , employing PAW pseudopotentialsPAW_PRB50 ; Kresse_PRB59 from the pslibrary (v1.0)DALCORSO2014337 . For the exchange-correlation functional we use the Perdew–Burke–Ernzerhof approximationPBE_PRL77 . The energy cut-off for the plane-waves basis is set to 35.0 Ry, and a $4\times 4\times 4$ uniform grid is used to sample the Brillouin zone. The lattice parameter is set to the experimental value (5.65 Å). The overlap matrices $M_{mn}^{(\mathbf{k,b})}$ in Eq. (9), the projection matrices $A_{mn\mathbf{k}}$ in Eq. (26) and both $\widetilde{d}_{\mathbf{k}}(g)$ in Eq. (30) and $D_{\mathbf{k}}(g)$ in Eq. (32) have been computed with the pw2wannier90.x interface.

GaAs is a III-V semiconductor that crystallises in the fcc cubic structure, with a two-atom basis: the Ga cation and the As anion (space group $F{-}43m$); in our example the Ga atom is placed at the origin of the unit cell, whose Wyckoff letter is $a$ and site-symmetry group is ${-}43m$, also known as $T_{d}$. The As atom is placed at ($\nicefrac{{1}}{{4}}$,$\nicefrac{{1}}{{4}}$,$\nicefrac{{1}}{{4}}$), whose Wyckoff letter is $c$ and site-symmetry group is also $T_{d}$.

Marzari and VanderbiltMV_PRB56 have shown that the MLWFs for the 4-dimensional valence manifold are centred on the four As-Ga bonds, have $sp^{3}$ character and can be found by specifying four $s$-like orbitals on each covalent bond as initial guess (a representative is shown in Fig. 1(a)). These bond-centred functions correspond to the irreducible representation $A_{1}$ of the site-symmetry group $C_{3v}$ of the Wyckoff position e. Hence, the MLWFs can also be obtained with the SAWF approach by specifying the centres and the shapes of the initial projections, e.g. four $s$-like orbitals centred on the four As–Ga bonds, and the symmetry operations in the point group $C_{3v}$.

Using the SAWF method we can enforce the WFs to have the local site symmetries. In particular, since $T_{d}$ has 5 irreps of dimension 1, 1, 2, 3 and 3 respectively, one can form an 1+3–dimensional representation for the four SAWFs. Thus, a set of initial projections compatible with the symmetries of the valence bands is: one $s$-like orbital (1-dimensional irrep whose character is $A_{1}$) and three $p$-like orbitals (3-dimensional irrep whose character is $T_{2}$) centred on As. Fig. 1(b) shows the SAWF which corresponds to the $A_{1}$ representation and transforms like the identity under $T_{d}$.

The same result can be obtained with the SLWF+C method by selectively localising one function $J^{\prime}=1$ $(J=4)$ and constraining its centre to sit on the As site $(\nicefrac{{1}}{{4}},\nicefrac{{1}}{{4}},\nicefrac{{1}}{{4}})$. In the case of GaAs the SLWF+C method turns out to be very robust, to the point that four $s$-like orbitals randomly centred in the unit cell can be used as initial guess without affecting the result of the optimised function. Fig. 1(c) shows the resulting function using the SLWF method without constraints, while Fig. 1(d) shows the result using SLWF+C.

It is worth to note that for this particular system, it is possible to achieve this result with the maximal localisation procedure if one carefully selects the initial projections, i.e., one $s$-like and three $p$-like orbitals on the As atom. The resulting WFs will possess the local site symmetries but will not correspond to the global minimum of the spread functional $\Omega$. More precisely, they will correspond to a saddle point of $\Omega$ (unstable against small perturbations of the initial projections). In Fig. 1-(a)-(b)-(c)-(d) we show a comparison of the centre and symmetries of a MLWF, SAWF, SLWF and SLWF+C; the individual spreads and total spreads—for all four valence states—are reported in the Table below it.

III.3 Parallelisation

In Wannier90 v3.0 we have implemented an efficient parallelisation scheme for the calculation of MLWFs using the message passing interface (MPI).

Calculation of the spread. The time-consuming part in the evaluation of the spread $\Omega$ is updating the $M^{(\mathbf{k},\mathbf{b})}$ matrices according to Eq. (16), since this requires computing overlap matrix elements between all pairs of bands, and between all $k$-points $\mathbf{k}$ and their neighbours $\mathbf{k}+\mathbf{b}$. Therefore, an efficient speed up for the evaluation of the spread can be achieved by distributing over several processes the calculation of the $M^{(\mathbf{k},\mathbf{b})}$ matrices for different $k$-points. In order to compute the $M^{(\mathbf{k},\mathbf{b})}$ according to Eq. (16), the $U_{\mathbf{k}+\mathbf{b}}$ matrices are sent from process to process prior to the calculation of the overlap matrices. We stress the fact that the $U_{\mathbf{k}+\mathbf{b}}$ matrices are the only large arrays that have to be shared between processes, which limits the time spent in communication. The relatively large $M^{(\mathbf{k},\mathbf{b})}$ matrices are not sent between processes for the evaluation of Eqs. (7) and (8). Instead, it is enough to collect the contributions to the spread from the different $k$-points, i.e., a set of scalars, and then sum them up for evaluation of the total spread. This parallelisation scheme is illustrated in Fig. 2 for a $3\times 3$ mesh of $k$-points with 9 MPI processes.

Minimisation of the spread. The minimisation of the spread functional is based on an iterative steepest-descent or conjugate-gradient algorithm. In each iteration, the unitary matrices $U_{\mathbf{k}}$ are updated according to $U_{\mathbf{k}}=U_{\mathbf{k}}\exp{(\Delta W_{\mathbf{k}})}$ MV_PRB56 , where $\Delta W_{\mathbf{k}}=\alpha\mathcal{D}_{\mathbf{k}}$, see Eq. (15). Updating the $U_{\mathbf{k}}$ matrices according to this equation is by far the most time-consuming part in the iterative minimisation algorithm, as it requires a diagonalisation of the $\Delta W_{\mathbf{k}}$ matrices. A significant speed-up can be obtained, however, by distributing the diagonalisation of the different $\Delta W_{\mathbf{k}}$ matrices over several processes, and performing the calculations fully in parallel. The evaluation of $\Delta W_{\mathbf{k}}$ essentially requires the calculation of the overlap matrices $M^{(\mathbf{k},\mathbf{b})}$, as discussed above.

Disentanglement. The disentanglement procedure is concerned with finding the optimal subspace $\mathcal{S}_{\mathbf{k}}$. As the functional $\Omega_{\mathrm{I}}$ measures the global subspace dispersion across the Brillouin zone, at first sight it is not obvious that the task of minimising the spread $\Omega_{\mathrm{I}}$ can be parallelised with respect to the $k$-points. In the iterative algorithm of Eq. (21), the systematic reduction of the spread functional at the $i^{\text{th}}$ iteration is achieved by minimising the spillage of the subspace $\mathcal{S}^{(i)}_{\mathbf{k}}$ over the neighbouring subspaces from the previous iteration $\mathcal{S}^{(i-1)}_{\mathbf{k}+\mathbf{b}}$. This problem reduces to the diagonalisation of $N$ independent matrices ($N$ is the total number of $k$-points of the mesh), where an efficient speed-up of the disentanglement procedure can be achieved by distributing the diagonalisation of the $Z^{(i)}_{\mathbf{k}}$ matrices of Eq. (23) over several processes, which can be done fully in parallel. Since the construction of $Z^{(i)}_{\mathbf{k}}$ only requires the knowledge of the $U^{(i-1)}_{\mathbf{k}+\mathbf{b}}$ matrices, these must be communicated between processes, as shown in Fig. 2. This results in a similar time spent in communication for the disentanglement part of the code as for the wannierisation part.

Distribution of large matrices. The parallelisation scheme discussed above relies on the parallel evaluation of relevant matrices over $k$-points on each processor. For systems with large number of $k$-points and bands, it is also desirable to distribute the matrices themselves among the available cores so that the memory per core required to store them is reduced. For example, in the case of isolated bands, storing all the $M^{(\mathbf{k},\mathbf{b})}$ matrices requires an allocation of dimension $J\times J\times N\times N_{b}$, where $N_{b}$ is the number neighbours of each of the $N$ $k$-points of the mesh. By distributing the storage across $N_{\text{c}}$ cores, the storage requirement per core decreases accordingly by a factor of approximately $N_{\text{c}}$.

Performance. We have tested the performance of this parallelisation scheme for the calculation of the MLWFs in a $L1_{0}-$FePt(5)/Pt(18) thin film. Computational details were given in Ref. Geranton-PRB2015, . The benchmarks have been performed on the JURECA supercomputer of the Jülich Supercomputing Center. We have extracted an optimal subspace of dimension $J=414$ from a set of 580 Bloch states per $k$-point. The upper limit of the inner window was set to 5 eV above the Fermi energy, and 414 MLWFs were constructed by minimising the spread $\Omega$. The performance benchmark was based on the average wall-clock time for a single iteration of the minimisation procedure (several thousand iterations are usually needed for convergence). We first analyse the weak scaling of our implementation, i.e., how the computation time varies with the number of cores $N_{\text{c}}$ for a fixed number of $k$-points per process. We show in Fig. 3(a) the time per iteration for the disentanglement and wannierisation parts of the minimisation, always using one $k$-point per process. As we vary the number of $k$-points $N$ from 4 to 144, the computation time increases only by a factor of 1.3 and 1.8 for disentanglement and wannierisation, respectively. We then demonstrate the strong scaling of our parallelisation scheme in Fig. 3(b), i.e., how the computation time varies with the number of cores $N_{\text{c}}$ for a fixed number $N=64$ of $k$-points. When varying the number of cores from 4 to 64, we observe a decrease of the computation time per iteration by a factor of 12.6 and 9.5 for disentanglement and wannierisation, respectively. The deviation from ideal scaling is mostly explained by the time spent in inter-core communication of the $U_{\mathbf{k}+\mathbf{b}}$ matrices.

IV Enhancements in functionality

In this section we describe a number of enhancements to the functionality of the core Wannier90 code, namely: the ability to compute and visualise spinor-valued WFs, including developments to the interface with the Quantum ESPRESSO package to cover also the case of non-collinear spin calculations performed with ultrasoft pseudopotentials (previously not implemented); an improvement to the method for interpolating the $k$-space Hamiltonian; the ability to select a subset from a larger set of projections of localised trial orbitals onto the Bloch states for initialising the WFs; and new functionality for plotting WFs in Gaussian cube format on WF-centred grids with non-orthogonal translation vectors.

IV.1 Spinor-valued Wannier functions with ultrasoft and projector-augmented-wave pseudopotentials

The calculation of the overlap matrix in Eq. (17) within the ultrasoft-pseudopotential formalism proceeds via the inclusion of so-called augmentation functions,Ferretti2007

[TABLE]

where $\ket{\psi^{\mathrm{ps}}_{m\mathbf{k}}}$ is the pseudo-wavefunction,

[TABLE]

is the Fourier transform of the augmentation charge, and $B^{(\mathbf{k},\mathbf{b})}_{Iij}=\ket{\beta^{\mathbf{k}}_{Ii}}\bra{\beta^{\mathbf{k}+\mathbf{b}}_{Ij}}$, where $\ket{\beta^{\mathbf{k}}_{Ii}}$ denotes the $i^{\rm th}$ projector of the pseudopotential on the $I^{\rm th}$ atom in the unit cell. We refer to Appendix B of Ref. Ferretti2007, for detailed expressions.

When spin-orbit coupling is included, the Bloch functions become two-component spinors $(\psi^{\uparrow}_{n\mathbf{k}}(\mathbf{r}),\psi^{\downarrow}_{n\mathbf{k}}(\mathbf{r}))^{\textrm{T}}$, where $\psi^{\sigma}_{n\mathbf{k}}(\mathbf{r})$ is the spin-up (for $\sigma=\uparrow$) or spin-down (for $\sigma=\downarrow$) component with respect to the chosen spin quantisation axis. Accordingly, $Q^{I}_{ij}(\mathbf{b})$ becomes $Q_{ij}^{I\sigma\sigma^{\prime}}(\mathbf{b})$ (see Eq. (18) in Ref. DalCorso2005, ) and Eq. (42) becomes

[TABLE]

The above expressions, together with the corresponding ones for the matrix elements of the spin operator, have been implemented in the pw2wannier90.x interface between Quantum ESPRESSO and Wannier90.

The plotting routines of Wannier90 have also been adapted to work with the complex-valued spinor WFs obtained from calculations with spin-orbit coupling. It then becomes necessary to decide how to represent graphically the information contained in the two spinor components.

One option is to only plot the norm $|\psi_{n\mathbf{k}}(\mathbf{r})|=\sqrt{|\psi^{\uparrow}_{n\mathbf{k}}(\mathbf{r})|^{2}+|\psi^{\downarrow}_{n\mathbf{k}}(\mathbf{r})|^{2}}$ of spinor WFs, which is reminiscent of the total charge density in the case of a 2$\times$2 density matrix in non-collinear DFT. Another possibility is to plot independently the up- and down-spin components of the spinor WF. Since each of them is in general complex-valued, two options are provided in the code: (i) to plot only the magnitudes $|\psi^{\uparrow}_{n\mathbf{k}}(\mathbf{r})|$ and $|\psi^{\downarrow}_{n\mathbf{k}}(\mathbf{r})|$ of the two components; or (ii) to encode the phase information by outputting $|\psi^{\uparrow}_{n\mathbf{k}}(\mathbf{r})|\text{sgn}(\text{Re}\{\psi^{\uparrow}_{n\mathbf{k}}(\mathbf{r})\})$ and $|\psi^{\downarrow}_{n\mathbf{k}}(\mathbf{r})|\text{sgn}(\text{Re}\{\psi^{\downarrow}_{n\mathbf{k}}(\mathbf{r})\})$, where sgn is the sign function. Which of these various options is adopted by the Wannier90 code is controlled by two input parameters, wannier_plot_spinor_mode and wannier_plot_spinor_phase.

Finally we note that, for WFs constructed from ultrasoft pseudopotentials or within the projector-augmented-wave (PAW) method, only pseudo-wavefunctions represented on the soft FFT grid are considered in plotting WFs within the present scheme, that is, the WFs are not normalised.

IV.2 Improved Wannier interpolation by minimal-distance replica selection

The interpolation of band structures (and many other quantities) based on Wannier functions is an extremely powerful toolLee_PRL05 ; wang-prb06 ; Yates_PRB07 . In many respects it resembles Fourier interpolation, which uses discrete Fourier transforms to reconstruct faithfully continuous signals from a discrete sampling, provided that the signal has a finite bandwidth and that the sampling rate is at least twice the bandwidth (the so-called Nyquist–Shannon condition).

In the context of Wannier interpolation, the “sampled signal” is the set of matrix elements

[TABLE]

of a lattice-periodic operator such as the Hamiltonian, defined on the same uniform grid $\{\mathbf{k}_{j}\}$ that was used to minimise the Wannier spread functional (see Sec. II.1). The states $\ket{\chi_{n\mathbf{k}_{j}}}$ are the Bloch sums of the WFs, related to ab initio Bloch eigenstates by $\ket{\chi_{n\mathbf{k}_{j}}}=\sum_{m}\ket{\psi_{m\mathbf{k}_{j}}}U_{mn\mathbf{k}_{j}}$.

To reconstruct the “continuous signal” $H_{nm\mathbf{k}}$ at arbitrary $\mathbf{k}$, the matrix elements of Eq. (45) are first mapped onto real space using the discrete Fourier transform

[TABLE]

where $N=N_{1}\times N_{2}\times N_{3}$ is the grid size (which is also the number of $k$-points in Wannier90). The matrices $H_{mn\mathbf{k}_{j}}$ are then interpolated onto an arbitrary $\mathbf{k}$ using an inverse discrete Fourier transform,

[TABLE]

where the sum is over $N$ lattice vectors $\mathbf{R}^{\prime}$, and the interpolated energy eigenvalues are obtained by diagonalising $H_{\mathbf{k}}$. In the limit of an infinitely dense grid of $k$-points the procedure is exact and the sum in Eq. (47) becomes an infinite series. Owing to the real-space localisation of the Wannier functions, the matrix elements $\widetilde{H}_{mn\mathbf{R}}$ become vanishingly small when the distance between the Wannier centres exceeds a critical value $L$ (the “bandwidth” of the Wannier Hamiltonian), so that actually only a finite number of terms contributes significantly to the sum in Eq. (47). This means that, even with a finite $N_{1}\times N_{2}\times N_{3}$ grid, the interpolation is still accurate provided that – by analogy with the Nyquist–Shannon condition – the “sampling rate” $N_{i}$ along each cell vector $\mathbf{a}_{i}$ is sufficiently large to ensure that $N_{i}|\mathbf{a}_{i}|>2L$.

Still, the result of the interpolation crucially depends on the choice of the $N$ lattice vectors to be summed over in Eq. (47). Indeed, when using a finite grid, there is a considerable freedom in choosing the set $\{\mathbf{R}^{\prime}\}$ as $\widetilde{H}_{mn\mathbf{R}}$ is invariant under $\mathbf{R}\rightarrow\mathbf{R}+\mathbf{T}$ for any vector $\mathbf{T}$ of the Born–von Karman superlattice generated by $\{\mathbf{A}_{i}=N_{i}\mathbf{a}_{i}\}$. The phase factor in Eq. (47) is also invariant when $\mathbf{k}\in\{\mathbf{k}_{j}\}$, but not for arbitrary $\mathbf{k}$. Hence we need to choose, among the infinite set of “replicas” $\mathbf{R}^{\prime}=\mathbf{R}+\mathbf{T}$ of $\mathbf{R}$, which one to include in Eq. (47). We take the original vectors $\mathbf{R}$ to lie within the Wigner–Seitz supercell centred at the origin. If some of them fall on its boundary then their total number exceeds $N$ and weight factors must be introduced in Eq. (47). For each combination of $m$, $n$ and $\mathbf{R}$, the optimal choice of $\mathbf{T}$ is the one that minimises the distance

[TABLE]

between the two Wannier centres. With this choice, the spurious effects arising from the artificial supercell periodicity are minimised.

Earlier versions of Wannier90 implemented a simplified procedure whereby the vectors $\mathbf{R}^{\prime}$ in Eq. (47) were chosen to coincide with the unshifted vectors $\mathbf{R}$ that are closer to the origin than to any other point $\mathbf{T}$ on the superlattice, irrespective of the WF pair $(m,n)$. As illustrated in Fig. 4, this procedure does not always lead to the shortest distance between the pair of WFs, especially when some of the $N_{i}$ are small and the Wannier centres are far from the origin of the cell.

Wannier90 now implements an improved algorithm that enforces the minimal-distance condition of Eq. (48), yielding a more accurate Fourier interpolation. The algorithm is the following:

1. (a)

For each term in Eq. (47) pick, among all the replicas $\mathbf{R}^{\prime}=\mathbf{R}+\mathbf{T}$ of $\mathbf{R}$, the one that minimises the distance between Wannier centres (Eq. (48)). 2. (b)

If there are $\mathcal{N}_{mn\mathbf{R}}$ different vectors $\mathbf{T}$ for which the distance of Eq. (48) is minimal, then include all of them in Eq. (47) with a weight factor $1/\mathcal{N}_{mn\mathbf{R}}$.

An equivalent way to describe these steps is that (a) we choose $\mathbf{T}$ such that $\mathbf{r}_{n}+\mathbf{R}+\mathbf{T}$ falls inside the Wigner–Seitz supercell centred at $\mathbf{r}_{m}$ (see Fig. 4), and that (b) if it falls on a face, edge or vertex of the Wigner–Seitz supercell, we keep all the equivalent replicas with an appropriate weight factor. In practice the condition in step (b) is enforced within a certain tolerance, to account for the numerical imprecision in the values of the Wannier centres and in the definition of the unit cell vectors. Although step (b) is much less important than (a) for obtaining a good Fourier interpolation, it helps ensuring that the interpolated bands respect the symmetries of the system; if step (b) is skipped, small artificial band splittings may occur at high-symmetry points, lines, or planes in the BZ.

The procedure outlined above amounts to replacing Eq. (47) with

[TABLE]

where $\{\mathbf{T}^{(j)}_{mn\mathbf{R}}\}$ are the $\mathcal{N}_{mn\mathbf{R}}$ vectors $\mathbf{T}$ that minimise the distance of Eq. (48) for a given combination of $m$, $n$ and $\mathbf{R}$; $\mathbf{R}$ lies within the Wigner–Seitz supercell centred on the origin.

The benefits of this modified interpolation scheme are most evident when considering a large unit cell sampled at the $\Gamma$ point only. In this case $N=1$ so that Eq. (47) with $\{\mathbf{R}^{\prime}\}=\{\mathbf{R}\}=\{\mathbf{0}\}$ would reduce to $H_{mn\mathbf{k}}=\widetilde{H}_{mn\mathbf{0}}$, yielding interpolated bands that do not disperse with $\mathbf{k}$. This is nonetheless an artefact of the choice $\{\mathbf{R}^{\prime}\}=\{\mathbf{0}\}$ (of earlier versions of Wannier90) and not an intrinsic limitation of Wannier interpolation, as first demonstrated in Ref. Lee_PRL05, for one-dimensional systems. Indeed, equation (49), which in a sense extends Ref. Lee_PRL05, to any spatial dimension, becomes in this case

[TABLE]

which can produce dispersive bands. This is illustrated in Fig. 5(a) for the case of a one-dimensional chain of carbon atoms: the interpolated bands obtained from Eq. (47) with $\{\mathbf{R}^{\prime}\}=\{\mathbf{R}\}=\{\mathbf{0}\}$ (earlier version of Wannier90) are flat, while those obtained from Eq. (49) (new versions of Wannier90) are in much better agreement with the dispersive ab initio bands up to a few eV above the Fermi energy.

Clear improvements in the interpolated bands are also obtained for bulk solids, as shown in Fig. 5(b) for the case of silicon. The earlier implementation breaks the two-fold degeneracy along the X$-$W line, with one of the two bands becoming flat. The new procedure recovers the correct degeneracies, and reproduces more closely the ab initio band structure (the remaining small deviations are due to the use of a coarse $k$-point mesh that does not satisfy the Nyquist–Shannon condition, and would disappear for denser $k$-grids together with the differences between the two interpolation procedures).

IV.3 Selection of projections

In many cases, and particularly for entangled bands, it is necessary to have a good initial guess for the MLWFs in order to properly converge the spread to the global minimum. Determining a good initial guess often involves a trial and error approach, using different combinations of orbital types, orientations and positions. While for small systems performing many computations of the projection matrices is relatively cheap, for large systems there is a cost associated with storing and reading the wavefunctions to compute new projection matrices for each new attempt at a better initial guess. Previously, the number of projections that could be specified had to be equal to the number $J$ of WFs to be constructed. The latest version of the code lifts this restriction, making it possible to define in the pre-processing step a larger number $J_{+}>J$ of projection functions to consider as initial guesses. In this way, the computationally expensive and potentially I/O-heavy construction of the projection matrices $A_{\mathbf{k}}$ is performed only once for all possible projections that a user would like to consider.

Once the $A_{\mathbf{k}}$ matrices (of dimension $J\times J_{+}$ at each $\mathbf{k}$) have been obtained, one proceeds with constructing the MLWFs by simply selecting, via a new input parameter (select_projections) of the Wannier90 code, which $J$ columns to use among the $J_{+}$ that were computed by the interface code. Experimenting with different trial orbitals can thus be achieved by simply selecting a different set of projections within the Wannier90 input file, without the need to perform the pre-processing step again.

Similarly, another use case for this new option is the construction of WFs for the same material but for different groups of bands. Typically one would have to modify the Wannier90 input file and run the interface code multiple times, while now the interface code may compute $A_{\mathbf{k}}$ for a superset of trial orbitals just once, and then different subsets may be chosen by simple modification of a single input parameter. As a demonstration, we have adapted example11 of the Wannier90 distribution (silicon band structure), that considers two band groups: (a) the valence bands only, described by four bond-centred $s$ orbitals, and (b) the four valence and the four lowest-lying conduction bands together, described by atom-centred $sp^{3}$ orbitals. In the example, projections onto all 12 trial orbitals are provided, and the different cases are covered by specifying in the Wannier90 input file which subset of projections is required.

IV.4 Plotting cube files with non-orthogonal vectors

In Wannier90 v3.0 it is possible to plot the MLWFs in real-space in Gaussian cube format, including the case of non-orthogonal cell lattice vectors. Many modern visualisation programs such as Vestavesta are capable of handling non-orthogonal cube files and the cube file format can be read by many computational chemistry programs. Wannier90’s representation of MLWFs in cube format can be significantly more compact than using the alternative xsf format. With the latter, MLWFs are calculated (albeit with a coarse sampling) on a supercell of the computational cell that can be potentially large (the extent of the supercell is controlled by an input parameter wannier_plot_supercell). Whereas, with the cube format, each Wannier function is represented on a grid that is centred on the Wannier function itself and has a user-defined extent, which is the smallest parallelepiped (whose sides are aligned with the cell vectors) that can enclose a sphere with a user-defined radius wannier_plot_radius. Because MLWFs are strongly localised in real space, relatively small cut-offs are all that is required, significantly smaller than the length-scale over which the MLWFs themselves are periodic. As a result, the cube format is particularly useful when a more memory-efficient representation is needed. The cube format can be activated by setting the input parameter wannier_plot_mode to cube, and the code can handle both isolated molecular systems (treated within the supercell approximation) as well as periodic crystals by setting wannier_plot_mode to either molecule or crystal, respectively.

V New post-processing features

Once the electronic bands of interest have been disentangled and wannierised to obtain well-localised WFs, the Wannier90 software package includes a number of modules and utilities that use these WFs to calculate various electronic-structure properties. Much of this functionality exists within postw90.x, an MPI-parallel code that forms an integral part of the Wannier90 package. In v2.x of Wannier90, postw90.x included functionality for computing densities of states and partial densities of states, energy bands and Berry curvature along specified lines and planes in $k$-space, anomalous Hall conductivity, orbital magnetisation and optical conductivity, Boltzmann transport coefficients within the relaxation time approximation, and band energies and derivatives on a generic user-defined list of $k$-points. Some further functionality exists in a set of utilities that are provided as part of the Wannier90 package, including a code (w90pov.F90) to plot WFs rendered using the Persistence of Vision Raytracer (POV-Ray)povray code and to compute van der Waals interactions with WFs (w90vdw.F90).

In addition, there are a number of external packages for computing advanced properties based on WFs and which interface to Wannier90. These include codes to generate tight-binding models such as pythTB pythtb and tbmodels tbmodels , quantum transport codes such as sisl zerothi_sisl , gollum gollum , omen omen and nanoTCAD-ViDES nanotcadvides , the EPW epw code for calculating properties related to electron-phonon interactions and WannierTools WannierTools for the investigation of novel topological materials.

Below we describe some of the new post-processing features of Wannier90 that have been introduced in the latest version of the code, v3.0.

V.1 postw90.x: Shift Current

The photogalvanic effect (PGE) is a nonlinear optical response that consists in the generation of a direct current (DC) when light is absorbed.belinicher-spu80 ; sturman-book92 ; ivchenko-book97 It can be divided phenomenologically into linear (LPGE) and circular (CPGE) effects, which have different symmetry requirements within the acentric crystal classes. The CPGE requires elliptically-polarised light, and occurs in gyrotropic crystals (see next subsection). The LPGE occurs with linearly or unpolarised light as well; it is present in piezoelectric crystals and is given by

[TABLE]

where ${\bf J}(0)$ is the induced DC photocurrent density, ${\bf E}(\omega)=\bf E^{*}(-\omega)$ is the amplitude of the optical electric field, and $\sigma_{abc}=\sigma_{acb}=\sigma_{abc}^{*}$ is a nonlinear photoconductivity tensor.

The shift current is the part of the LPGE photocurrent generated by interband light absorption.tan-cm16 Intuitively, it arises from a coordinate shift accompanying the photoexcitation of electrons from one band to another. Like the intrinsic anomalous Hall effect nagaosa-rmp10 , the shift current involves off-diagonal velocity matrix elements between occupied and empty bands, depending not only on their magnitudes but also on their phases baltz-prb81 ; belinicher-jetp82 ; kristoffel-zpb82 ; sipe-prb00 .

The shift current along direction $a$ induced by light that is linearly polarised along $b$ is described by the following photoconductivity tensor:sipe-prb00 ; fregoso-prb17

[TABLE]

Here, $f_{nm\mathbf{k}}=f_{n\mathbf{k}}-f_{m\mathbf{k}}$ is the difference between occupation factors, $\hbar\omega_{mn\mathbf{k}}=\epsilon_{m\mathbf{k}}-\epsilon_{n\mathbf{k}}$ is the difference between energy eigenvalues of the Bloch bands, $r^{b}_{nm\mathbf{k}}$ is the $b^{\rm th}$ Cartesian component of the interband dipole matrix (the off-diagonal part of the Berry connection matrix ${\bf A}_{nm\mathbf{k}}=i\langle u_{n\mathbf{k}}|\partial_{\mathbf{k}}u_{m\mathbf{k}}\rangle$), and

[TABLE]

is the shift vector (not to be confused with the lattice vector ${\bf R}$, or with the matrix $R^{({\bf k},{\bf b})}$ defined in Eq. (12)). The shift vector has units of length, and it describes the real-space shift of wavepackets under photoexcitation.

The numerical evaluation of Eq. (53) is tricky because the individual terms therein are gauge-dependent, and only their sum is unique. Different strategies were discussed in the early literature in the context of model calculationskristoffel-zpb82 ; presting-pssb82 and more recently for ab initio calculations. The ab initio implementation of Young and Rappeyoung-prl12 employed a gauge-invariant $k$-space discretisation of Eq. (53), inspired by the discretised Berry-phase formula for electric polarisation.king-smith-prb93

The implementation in Wannier90 is based instead on the formulation of Sipe and co-workers.aversa-prb95 ; sipe-prb00 In this formulation, the shift (interband) contribution to the LPGE tensor in Eq. (51) is expressed as

[TABLE]

where

[TABLE]

is the generalised derivative of the interband dipole. When $b=c$, Eq. (54) becomes equivalent to Eq. (52).sipe-prb00

The generalised derivative $r^{b;a}_{nm\mathbf{k}}$ is a well-behaved (covariant) quantity under gauge transformation but – as in the case of the shift vector – this is not the case for the individual terms in Eq. (55), leading to numerical instabilities. To circumvent this problem, Sipe and co-workers used $\mathbf{k}\cdot\mathbf{p}$ perturbation theory to recast Eq. (55) as a summation over intermediate virtual states where the individual terms are gauge-covariant.aversa-prb95 ; sipe-prb00 That strategy has been successfully employed to evaluate the shift-current spectrum from first principles.nastos-prb06 ; rangel-prl17

As it is well known, similar “sum rule” expressions can be written for other quantities involving $\mathbf{k}$ derivatives, such as the inverse effective-mass tensor and the Berry curvature tensor. When evaluating such expressions, a sufficiently large number of virtual states must be included to achieve convergence. Alternatively, one can work with a basis spanning a finite number of bands, such as a tight-binding or Wannier basis, and carefully reformulate $\mathbf{k}\cdot\mathbf{p}$ perturbation theory within that incomplete basis to avoid truncation errors. This was done first for the inverse effective-mass tensorgraf-prb95 ; boykin-prb95 and later for the Berry curvature,wang-prb06 and is at the heart of the Wannier-interpolation technique for calculating the intrinsic anomalous Hall conductivity (AHC).wang-prb06

A truncation-free tight-binding sum rule for the generalised derivative of Eq. (55) was given in Ref. cook-nc17, . Contrary to the inverse effective-mass tensor, which only depends on the Hamiltonian matrix elements,graf-prb95 ; boykin-prb95 the generalised derivative also depends – in some cases rather strongly – on the intracell coordinates of the basis orbitals. cook-nc17 Building on that formulation, Wannier interpolation schemes for calculating the shift current were recently introducedwang-prb17 ; azpiroz-prb18 (the implementation in Wannier90 follows Ref. azpiroz-prb18, ). In addition to the Hamiltonian matrix elements and Wannier centres, the shift current was found to depend sensitively on the off-diagonal position matrix elements.

The generalised derivative can be used to evaluate other nonlinear optical responses, such as second-harmonic generation.sipe-prb00 ; wang-prb17 While these are not currently implemented in Wannier90, it should be straightforward to adapt the shift-current routines for that purpose.

V.2 postw90.x: Gyrotropic module

The spontaneous magnetisation of ferromagnets endows their linear conductivity tensor $\sigma_{ab}(\omega)$ with an antisymmetric part. At $\omega=0$ that part describes the anomalous Hall conductivity (AHC), and at finite frequencies it gives rise to magneto-optical effects such as Faraday rotation in transmission and magnetic circular dichroism in absorption. In paramagnets, those effects appear under an external magnetic field.

Interestingly, an antisymmetric conductivity can be induced in certain nonmagnetic (semi)conductors by purely electrical means, namely, by passing a current through the sample baranova-oc77 ; ivchenko-jetp78 . Symmetry arguments indicate that this is allowed in the gyrotropic crystal classes, a subset of the acentric crystal classes that includes those that are chiral, polar, or optically active.belinicher-spu80 The first experimental demonstration consisted in the measurement of a current-induced change in the rotatory power of $p$-doped trigonal tellurium.vorobev-jetp79 ; shalygin-pss12 The DC or transport limit of this current-induced Faraday effect is the current-induced anomalous Hall effect (AHE), which has become known in the recent literature as the nonlinear AHE sodemann-prl15 ; zhang-prb18 ; zhang-2dmater18 ; you-prb18 ; ma-nat19 . Like the linear (spontaneous) AHE in ferromagnetic metals, the nonlinear (current-induced) AHE in gyrotropic conductors has an intrinsic contribution associated with the Berry curvature in momentum space.sodemann-prl15

Along with nonlinear magneto-optical and anomalous Hall effects, the flow of electrical current in a gyrotropic conducting medium also generates a net magnetisation. This kinetic magnetoelectric effect was originally proposed for bulk chiral conductors, ivchenko-jetp78 ; levitov-jetp85 and later for two-dimensional (2D) inversion layers with an out-of-plane polar axis edelstein-ssc90 ; aronov-jetp89 , where it has been studied intensively ganichev-book12 . The kinetic magnetoelectric effect in 2D – also known as the Edelstein effect – is a purely spin effect, whereas in bulk crystals an orbital contribution is also present. levitov-jetp85 The orbital kinetic magnetoelectric effect was recently formulated in terms of the intrinsic orbital moment of the Bloch electrons,yoda-sr15 ; zhong-prl16 a quantity closely related to the Berry curvature.

Another phenomenon characteristic of gyrotropic crystals is the circular photogalvanic effect (CPGE) that was mentioned briefly in Sec. V.1. This nonlinear optical effect consists in the generation of a photocurrent that reverses sign with the helicity of light ivchenko-jetp78 ; asnin-jetp78 ; belinicher-spu80 ; sturman-book92 ; ivchenko-book97 , and it occurs when light is absorbed via interband or intraband scattering processes. The intraband contribution to the CPGE is closely related to the nonlinear AHE, as both arise from the Berry curvature of the conduction electrons deyo-arxiv09 ; moore-prl10 ; sodemann-prl15 .

The gyrotropic effects listed above are being very actively investigated in connection with novel materials ranging from topological semimetals juan-natcomms17 ; zhang-prb18 ; flicker-prb18 to monolayer and bilayer transition-metal dichalcogenides zhang-2dmater18 ; you-prb18 ; ma-nat19 . The sensitivity of both the Berry curvature and the intrinsic orbital moment to the details of the electronic structure, together with the need to sample them on a dense mesh of $k$-points, calls for the development of accurate and efficient ab initio methodologies for this class of problems.

Building on existing Wannier-interpolation schemes for calculating the spontaneous intrinsic AHC and orbital magnetisation,wang-prb06 ; lopez-prb12 the corresponding methodology for gyrotropic effects was presented in Ref. tsirkin-prb18, , where it was applied to $p$-doped trigonal tellurium (in that work, only the intraband contribution to the CPGE was considered). The resulting computer code has been incorporated in Wannier90 as the gyrotropic.F90 module.

The central task of that module is to evaluate response tensors such as the “Berry-curvature dipole” sodemann-prl15

[TABLE]

where $f_{0}$ is the equilibrium occupation factor and ${\bm{\Omega}}_{n\mathbf{k}}=\bm{\nabla}_{\mathbf{k}}\times{\bf A}_{nn\mathbf{k}}$ is the Berry curvature of the Bloch bands (the curl of the band-diagonal Berry connection ${\bf A}_{nn\mathbf{k}}$ introduced in Sec. V.1). Also of interest is the tensor $K_{ab}$, obtained by replacing the Berry curvature in Eq. (56) with the intrinsic magnetic moment of the Bloch states.yoda-sr15 ; zhong-prl16 ; tsirkin-prb18

The two tensors $D_{ab}$ and $K_{ab}$ describe several of the aforementioned gyrotropic effects as follows:

• •

intraband CPGE: $j_{a}\propto\frac{\omega\tau^{2}D_{ab}}{1+\omega^{2}\tau^{2}}{\rm Im}\left[{\bf E}(\omega)\times{\bf E}^{*}(\omega)\right]_{b}$,

• •

nonlinear AHE: $j_{a}\propto\tau\varepsilon_{adc}D_{bd}E_{b}E_{c}$,

• •

kinetic magnetoelectric effect: $M_{a}\propto\tau K_{ba}E_{b}$,

where ${\bf j}$ is the induced current density, ${\bf M}$ is the induced magnetisation, ${\bf E}$ or ${\bf E}(\omega)$ is the amplitude of the static or optical electric field, $\tau$ is the relaxation time of the conduction electrons, and $\varepsilon_{abc}$ is the alternating tensor. The reader is referred to Ref. tsirkin-prb18, for more details such as the prefactors in the expressions above, as well as the formulas for the current-induced Faraday effect and natural optical activity, both of which are also implemented in the gyrotropic.F90 module.

V.3 postw90.x: Spin Hall conductivity

The spin Hall effect (SHE) is a phenomenon in which a spin current is generated by applying an electric field. The current is often transverse to the field (Hall-like), but this is not always the case.wimmer-prb15 The SHE is characterised by the spin Hall conductivity (SHC) tensor $\sigma_{ab}^{\text{spin},c}$ as follows:

[TABLE]

where $J_{a}^{\text{spin},c}$ is the spin-current density along direction $a$ with its spin pointing along $c$, and $E_{b}$ is the external electric field of frequency $\omega$ applied along $b$. In non-magnetic materials the equal number of up- and down-spin electrons forces the AHE to vanish, resulting in a pure spin current.

Like the AHC, the SHC contains both intrinsic and extrinsic contributions.RevModPhys.87.1213 The intrinsic contribution to the SHC can be calculated from the following Kubo formula,PhysRevB.98.214402

[TABLE]

where ${s}_{c}$, ${v}_{a}$ and ${j}_{a}^{\text{spin},c}=\frac{1}{2}\{{s}_{c},{v}_{a}\}$ are the spin, velocity and spin current operators, respectively; $V$ is the cell volume, and $N$ is the total number of $k$-points used to sample the BZ. Equations (58) are very similar to the Kubo formula for the AHC, except for the replacement of a velocity matrix element by a spin-current matrix element. As mentioned in the previous two subsections, Wannier-interpolation techniques are very efficient at calculating such quantities.

A Wannier-interpolation method scheme for evaluating the intrinsic SHC was developed in Ref. PhysRevB.98.214402, (see also Ref. ryoo-prb19, for a related but independent work). The required quantities from the underlying ab initio calculation are the spin matrix elements $S_{mn\mathbf{k},a}^{(0)}=\langle\psi_{m\mathbf{k}}^{(0)}|{s}_{a}|\psi_{n\mathbf{k}}^{(0)}\rangle$, the Hamiltonian matrix elements $H_{mn\mathbf{k}}^{(0)}=\langle\psi_{m\mathbf{k}}^{(0)}|{H}|\psi_{n\mathbf{k}}^{(0)}\rangle=\epsilon^{(0)}_{m\bf k}\delta_{mn}$, and the overlap matrix elements of Eq. (17). Since the calculation of all these quantities has been previously implemented in pw2wannier90.x (the interface code between pwscf and Wannier90), this advantageous interpolation scheme can be readily used while keeping to a minimum the interaction between the ab initio code and Wannier90.

The application of the method to fcc Pt is illustrated in Fig. 6. Panel (a) shows the calculated SHC as a function of the Fermi-level position, and panel (b) depicts the “spin Berry curvature” of Eq. (58b) that gives the contribution from each band state to the SHC. The aforementioned functionalities have been incorporated in the berry.F90, kpath.F90 and kslice.F90 modules of postw90.x.

V.4 postw90.x: Parallelisation improvements

The original implementation of the berry.F90 module in postw90.x (for computing Berry-phase properties such as orbital magnetisation and anomalous Hall conductivitylopez-prb12 ), introduced in Wannier90 v2.0, was written with code readability in mind and had not been optimised for computational speed. In Wannier90 v3.0, all parts of the berry.F90 module have been parallelised while keeping the code readable; moreover, its scalability has been improved, accelerating its performance by several orders of magnitude.whitepaper

To illustrate the improvements in performance we present calculations on a 128-atom supercell of GaAs interstitially doped with Mn. We use a lattice constant of the elementary cell of 5.65 Å. We use norm-conserving relativistic pseudopotentials with the PBE exchange-correlation functional. The energy cut-off for the plane waves is set to 40 Ry, and the Brillouin-zone sampling of the supercell is $3\times 3\times 3$. We use a Gaussian metallic smearing with a broadening of 0.015 Ry. For the non-self-consistent step of the calculation, 600 bands are computed and used to construct 517 Wannier functions. The initial projections are chosen as a set of $sp^{3}$ orbitals centred on each Ga and As atom, and a set of $d$ orbitals on Mn. The calculations were performed on the Prometheus supercomputer of PL-GRID (in Poland).

The Berry-phase calculations can be performed in three distinct ways: (i) 3D quantities in $k$-space (routine berry_main), (ii) the same quantities resolved on 2D planes (routine kslice.F90), and (iii) 1D paths (routine kpath.F90) in the Brillouin zone. In the benchmarks, we will refer to these three cases as “Berry 3D”, “Berry 2D”, and “Berry 1D”, respectively.

The first optimisation target was the function utility_rotate in the module utility.F90, which calculates a matrix product of the form $B=R^{\dagger}AR$ using Fortran’s built-in matmul function. The new routine utility$\_$rotate$\_$new uses instead BLAS and performs about 5.7 times better than the original one, giving a total speedup for berry$\_$main of about 55$\%$.

A second performance-critical section of code was identified in the routine get_imfgh_k_list, which took more than 50$\%$ of the total run-time of berry_main. This routine computes three quantities: $F_{\alpha\beta}$, $G_{\alpha\beta}$ and $H_{\alpha\beta}$, which are defined in Eqs. (51), (66) and (56) of Ref. lopez-prb12, . By some algebraic transformations, it was possible to reduce 25 calls to matmul, carried out in the innermost runtime-critical loop, to only 5 calls. After replacement of matmul with the Basic Linear Algebra Subprogram (BLAS), the speed up of this routine exceeds a factor of 11, and the total time spent in berry_main is 2.5 times shorter (including the speed-up from the first optimisation).

In the third step, a bottleneck was eliminated in the initialisation phase, where mpi_bcast was waiting more than two minutes for the master rank to broadcast the parameters. The majority of this time was spent in loops computing matrix products of the form $S=(V_{1})^{\dagger}S_{0}V_{2}$. Again, we replaced this with two calls to the BLAS gemm routine. This resulted in a speed-up of a factor of 610 for the calculation of this matrix product in our test case, and the total initialisation time dropped to less than 15 seconds. In total, the berry_main routine runs about 5 times faster than it did originally.

Finally, the routines kslice.F90 and kpath.F90 were parallelised. The scalability results of berry_main, kslice.F90 and kpath.F90 are presented in Fig. 7, and a comparison with the scalability of the previous version of berry_main is also given. Absolute times for some of the calculations are reported in Table 1.

V.5 GW bands interpolation

While density-functional theory (DFT) is the method of choice for most applications in materials modelling, it is well known that DFT is not meant to provide spectral properties such as band structures, band gaps and optical spectra. Green’s function formulation of many-body perturbation theory (MBPT) fetter overcomes this limitation, and allows the excitation spectrum to be obtained from the knowledge of the Green’s function. Within MBPT the interacting electronic Green’s function $G(\mathbf{r},\mathbf{r^{\prime}},\omega)$ may be expressed in terms of the non-interacting Green’s function $G^{0}(\mathbf{r},\mathbf{r^{\prime}},\omega)$ and the so-called self-energy $\Sigma(\mathbf{r},\mathbf{r^{\prime}},\omega)$, where several accurate approximations for $\Sigma$ have been developed and implemented into first-principles codesmartin_reining_ceperley_2016 . While maximally-localised Wannier functions for self-consistent GW quasiparticles have been discussed in Ref. hamann_gw_09, , here we focus on the protocol to perform bands interpolation at the one-shot G0W0 level. For solids, the G0W0 approximation has proven to be an excellent compromise between accuracy and computational cost and it has become the most popular MBPT technique in computational materials science lucia_review_18 . In the standard one-shot G0W0 approach, $\Sigma$ is written in terms of the Kohn–Sham (KS) Green’s function and the RPA dielectric matrix, both obtained from the knowledge of DFT-KS orbitals and eigenenergies. Quasi-particle (QP) energies are obtained from:

[TABLE]

where $\psi_{n\mathbf{k}}$ and $\epsilon_{n\mathbf{k}}$ are the KS orbitals and eigenenergies, $Z_{n\mathbf{k}}$ is the so-called renormalisation factor and $V_{\rm xc}$ is the DFT exchange-correlation potential. In addition, in the standard G0W0 approximation the QP orbitals are approximated by the KS orbitals. At variance with DFT, QP corrections for a given $k$-point require knowledge of the KS orbitals and eigenenergies at all points $(\mathbf{k}+\mathbf{q})$ in reciprocal space. In practice, codes such as Yambo yambo_paper_09 compute QP corrections on a regular grid and rely on interpolation schemes to obtain the full band structure along high-symmetry lines. Wannier90 supports the use of G0W0 QP corrections through the general interface gw2wannier90.py distributed with Wannier90, while dedicated tools for Quantum ESPRESSO and Yambo allow an efficient use of symmetries. Here we briefly outline the procedure for performing Wannier interpolation at the G0W0 level using Quantum ESPRESSO and Yambo. The procedure starts by obtaining MLWFs at the DFT level. While Wannier90 works with uniform coarse meshes on the full BZ (FBZ), Yambo uses symmetries to compute quantities on the irreducible BZ (IBZ). In addition, the G0W0 self-energy may require finer $k$-point grids to achieve convergence compared to those required for the charge density in DFT or the Wannier interpolation itself. The k-mapper.py utility allows the user to quickly select only the symmetry-inequivalent $k$-points in the IBZ that belong to the grid used by Wannier90. At this point, Yambo computes the QP corrections on these selected $k$-points only. After that, a post-processing code of Yambo (ypp) unfolds the QP corrections onto the full BZ as required by Wannier90. Using the unfolded QP corrections, the utility gw2wannier90.py corrects and reorders in energy both the KS eigenvalues and the input matrices. After reading these eigenvalues and matrices, Wannier90 can proceed as usual to interpolate the desired quantities such as the band structure, but now at the G0W0 level.

VI Automatic Wannier functions: the SCDM method

An alternative method for generating localised Wannier functions, known as the selected columns of the density matrix (SCDM) algorithm, has been proposed by Damle, Lin and Ying DL_2015_SCDM ; DL_2018_SIAM . At its core the scheme exploits the information stored in the real-space representation of the single-particle density matrix, a gauge-invariant quantity. Localisation of the resulting functions is a direct consequence of the well-known nearsightedness principleKohn_PRB7 ; Prodan_PNAS_2005 of electronic structure in extended systems with a gapped Hamiltonian, i.e., insulators and semiconductors. In these cases, the density matrix is exponentially localised along the off-diagonal direction in its real-space representation $\rho(\mathbf{r},\mathbf{r}^{\prime})$ and it is generally accepted that Wannier functions with an exponential decay also exist; numerical studies have confirmed this claim for a number of materials, and there exist formal proofs for multiband time-reversal-invariant insulatorsBrouder_PRL_2007 ; He_PRL_86 ; Fiorenza2016 . Since the SCDM method does not minimise a given gauge-dependent localisation measure via a minimisation procedure, it is free from any issue regarding the dependence on initial conditions, i.e., it does not require a good initial guess of localised orbitals. It also avoids other problems associated with a minimisation procedure, such as getting stuck in local minima. More generally, the localised Wannier functions provided by the SCDM method can be used as starting points for the MLWF minimisation procedure, by using them to generate the $A_{\mathbf{k}}$ projection matrices needed by Wannier90.

For extended insulating systems, the density matrix is given by

[TABLE]

As shown in Sec. II, the $P_{\mathbf{k}}$ are the spectral projectors associated with the crystal Hamiltonian operator $H_{\mathbf{k}}$ onto the valence space $\mathcal{S}_{\mathbf{k}}$, hence their rank is $N_{\rm e}$ (number of valence electrons). Moreover, they are analytic functions of $\mathbf{k}$ and also manifestly gauge invariantNenciu_RMP_63 ; Panati_CMP_2013 . As mentioned above, the nearsightedness principleProdan_PNAS_2005 guarantees that the columns of the kernels $P_{\mathbf{k}}(\mathbf{r},\mathbf{r}^{\prime})=\braket{\mathbf{r}}{P_{\mathbf{k}}}{\mathbf{r}^{\prime}}$ are localised along the off-diagonal direction and therefore they may be used to construct a localised basis. If we consider a discretisation of the $J$ Bloch states at each $\mathbf{k}$ on a real-space grid of $N_{\rm g}$ points, we can arrange the wavefunctions into the columns of a unitary $N_{\rm g}\times J$ $k$-dependent matrix $\Psi_{\mathbf{k}}$

[TABLE]

such that $P_{\mathbf{k},ij}=\left(\Psi_{\mathbf{k}}\Psi_{\mathbf{k}}^{\dagger}\right)\phantom{}_{ij}$ is a $N_{\rm g}{\times}N_{\rm g}$ matrix. In this representation, it is straightforward to see that the columns of $P_{\mathbf{k}}(\mathbf{r}_{i},\mathbf{r}_{j})$ are projections of extremely localised functions (i.e., Dirac-delta functions localised on the grid points) onto the valence eigenspace. As a result, selecting any linearly-independent subset of $J$ of them will yield a localised basis for the span of $P(\mathbf{r},\mathbf{r}^{\prime})$. However, randomly selecting $J$ columns does not guarantee that a well-conditioned basis will be obtained. For instance, there could be too much overlap between the selected columns. Conceptually, the most well conditioned columns may be found via a QR factorisation with column pivoting (QRCP) applied to $P(\mathbf{r},\mathbf{r}^{\prime})$, in the form $P\Pi=QR$, with $\Pi$ being a matrix permuting the columns of $P$, $Q$ a unitary matrix and $R$ an upper-triangular matrix (not to be confused with the lattice vector ${\bf R}$, or with the matrix $R^{({\bf k},{\bf b})}$ defined in Eq. (12), or with the shift vector of Eq. (53)), and where $\Pi$ is chosen so that $|R_{11}|\geq|R_{22}|\geq\cdots\geq|R_{nn}|$. Then the $J$ columns forming a localised basis set are chosen to be the first $J$ of the matrix with permuted columns $P\Pi$.

The SCDM-$k$DL_2018_SIAM method suggests that it is sufficient to apply the QRCP factorisation at $\mathbf{k}=\bm{0}$ ($\Gamma$ point) only, and use the same selection of columns at all $k$-points. However, this is still often impractical since $P_{\bm{\Gamma}}$ is prohibitively expensive to construct and store in memory. Therefore an alternative procedure is proposed, for which the columns can be computed via the QRCP of the (smaller) matrix $\Psi_{\bm{\Gamma}}^{\dagger}$ instead:

[TABLE]

i.e., the same $\Pi$ matrix is obtained by computing a QRCP on $\Psi^{\dagger}$ only. Once the set of columns has been obtained, we need to impose the orthonormality constraint on the chosen columns without destroying their locality in real space. This can be achieved by a Löwdin orthogonalisation, similarly to Eq. (26). In particular, the selection of columns of $\Psi_{\bm{\Gamma}}$ can be used to select the columns of all $\Psi_{\mathbf{k}}$, which in turn define the $A_{mn\mathbf{k}}$ matrices needed as input by Wannier90 to start the MLWF minimisation procedure, by defining $A_{mn\mathbf{k}}=\psi^{\ast}_{m\mathbf{k}}(\mathbf{r}_{\Pi(n)})$, where the $\Pi(n)$ is the index of the $n^{\rm th}$ column of $P$ after permutation with $\Pi$. In fact, we can write the $n^{\rm th}$ column of $P$ after permutation, $P_{\mathbf{k}}(\mathbf{r},\mathbf{r}_{\Pi(n)})$, as

[TABLE]

The unitary matrix $U_{\mathbf{k}}$ sought for is then constructed via Löwdin orthogonalisation

[TABLE]

We can also extend the SCDM-$k$ method to the case where the Bloch states are represented as two-component spinor wavefunctions $\psi_{n\mathbf{k}}(\mathbf{r},\alpha)$, e.g., when including spin-orbit interaction in the Hamiltonian. Here, $\alpha=\uparrow,\downarrow$ is the spinor index. In this case, we include the spin index as well as the position index to perform QRCP. First, we define the $2N_{g}\times J$ matrix $\Psi_{\mathbf{k}}$

[TABLE]

Next, as in the spinless case, the QRCP of $\Psi_{\mathbf{\Gamma}}^{\dagger}$ is computed, and the first $J$ columns of the $\Pi$ matrix are selected. Now, $\Pi(n)$, the index of the $n^{\rm th}$ column of $P$ after permutation with $\Pi$, determines both the position index $\mathbf{r}_{\Pi(n)}$ and the spin index $\alpha_{\Pi(n)}$. We define $A_{mn\mathbf{k}}=\psi^{*}_{m\mathbf{k}}(\mathbf{r}_{\Pi(n)},\alpha_{\Pi(n)})$ and perform Löwdin orthogonalisation to obtain the unitary matrix $U_{\mathbf{k}}$.

In the case of entangled bands, we need to introduce a so-called quasi-density matrix defined as

[TABLE]

where $f(\epsilon_{n\mathbf{k}})\in[0,1]$ is a generalisation of the Fermi-Dirac probability for the occupied states. Also in this case we only use the information at $\Gamma$ to generate the permutation matrix. Depending on what kind of entangled manifold one is interested in, $f(\epsilon)$ can be modelled with various functional forms. In particular, the authors of Ref. DL_2018_SIAM, suggest the following three forms:

Isolated manifold, e.g., the valence bands of an insulator or a semiconductor: $f(\epsilon)$ is a step function, with the step inside the energy gap $\Delta\epsilon_{\rm g}=\epsilon_{\rm c}-\epsilon_{\rm v}$, where $\epsilon_{\rm c(v)}$ represents the minimum (maximum) of the conduction (valence) band:

[TABLE]

Both $\Delta\epsilon_{\rm g}$ and $\epsilon_{\rm v}$ are not free parameters, as they may be obtained directly from the ab initio calculation. 2. 2.

Entangled manifold (case I), e.g., the valence bands and low-lying conduction bands in a semiconductor: $f(\epsilon)$ is a complementary error function:

[TABLE]

where $\mu$ is used to shift the mid-value of the complementary error function, so that states with energy equal to $\mu$ have a weight of $f(\mu)=1/2$. The parameter $\sigma$ is used to gauge the “broadness” of the distribution function. 3. 3.

Entangled manifold (case II), e.g., the $d$ bands in a transition metal: $f(\epsilon)$ is a Gaussian function

[TABLE]

The procedure then follows as in the previous case, by computing a QRCP factorisation on the quasi-density matrix. It is worth to note that in the case of an entangled manifold, the SCDM method requires the selection of two real numbers: $\mu$ and $\sigma$, as well as the number of Wannier functions to disentangle $J$. These parameters play a crucial role in the selection of the columns of the density matrix. While the selection of these parameters requires some care, as a rule of thumb (e.g., in entangled case I) $\sigma$ is of the order $2-5$ eV (which is the energy range of a typical bandwidth), while $\mu$ can often be set around the Fermi energy (but the exact value depends on various factors, including the number $J$ of bands chosen and the specific properties of the bands of interest). It is worth to mention that since the SCDM-$k$ method is employed as an alternative way of specifying a set of initial projections and hence to compute the $A_{\mathbf{k}}$ matrices in Eq. (26), the disentanglement procedure can be used in exactly the same way as described in Sec. II.2. However, in the case of entangled bands the column selection is done on a quasi-density matrix, which implicitly defines a working subspace larger than the target subspace of dimension $J$. We find that for well-known systems SCDM-$k$ is typically already capable of selecting a smooth manifold and no further subspace selection is needed.

This method is now implemented as part of the pw2wannier90.x interface code to Quantum ESPRESSO. We have decided to implement the algorithm in the interface code(s) rather than in Wannier90 itself, because the SCDM method requires knowledge of the wavefunctions $\psi_{n\mathbf{k}}$, which are only available in the ab initio code.

In Wannier90 only a single new input parameter auto_projections is required. This disables the check on the number of projections specified in the input file (as we rely on SCDM to provide us with the initial guesses) and adds a new entry to the <seedname>.nnkp file (which is read by pw2wannier90.x in order to compute the quantities required by Wannier90) that specifies the number of Wannier functions required. The remaining control parameters for the SCDM method are specified in the input file for the pw2wannier90.x code, including whether to use the SCDM method, the functional form of the $f(\epsilon)$ function in Eq. (67) and, optionally, the values of $\mu$ and $\sigma$ in the definition of $f(\epsilon)$.

VII Automation and workflows: AiiDA-Wannier90 plugin

AiiDA pizzi-AiiDA (Automated Interactive Infrastructure and Database for Computational Science) is an informatics infrastructure that helps researchers in managing, automating, storing and sharing their computations and results. AiiDA automatically tracks the entire provenance of every calculation to ensure full reproducibility, which is also stored in a tailored database for efficient querying of previous results. Moreover, it provides a workflow engine, allowing researchers to implement high-level workflows to automate sequences of tedious or complex calculation steps. AiiDA supports simulation codes via a plugin interface, and over 30 different plugins are available to date aiida-registry .

Among these, the AiiDA-Wannier90 plugin provides support for the Wannier90 code. Users interact with the code (to submit calculations and retrieve the results) via the high-level python interface provided by AiiDA rather than directly creating the Wannier90 input files. AiiDA will then handle automatically the various steps involved in submitting calculations to a cluster computer, retrieving and storing the results, and parsing them into a database. Furthermore, using the AiiDA workflow system users can chain pre-processing and post-processing calculations automatically (e.g., the preliminary electronic structure calculation with an ab initio code). These scientific workflows, moreover, can encode in a reproducible form the scientific knowledge of expert computational researchers in the field on how to run the simulations, choose the numerical parameters and recover from potential errors. In turn, their availability reduces the training time of new researchers, eliminates sources of error and enables large-scale high-throughput simulations.

The AiiDA-Wannier90 plugin expects that each calculation takes a few well-defined input parameters. Among the most important ones, a Wannier90 calculation run via AiiDA requires that the following input nodes are provided: an input crystal structure, a node of parameters with a dictionary of input flags for Wannier90, a node with the list of kpoints, a node representing the atomic projections, and a local_input_folder or remote_input_folder node containing the necessary input files (.amn, .mmn, .nnkp, .eig, .dmn) for the Wannier90 calculation as generated by an ab initio code.

All of these parameters, with the exception of projections, are generic to AiiDA to facilitate their reuse with different simulation codes. More detailed information on all inputs can be found in the AiiDA-Wannier90 package documentation aiidawan-docs .

After the Wannier90 execution is completed, the AiiDA-Wannier90 plugin provides parsers that are able to detect whether the convergence was successful and retrieve key parameters including the centres of the Wannier functions and their spread, as well as the different components of the spread ($\Omega_{\mathrm{I}}$, $\Omega_{\mathrm{D}}$, $\Omega_{\mathrm{OD}}$ and $\Omega$), and (if computed) the maximum imaginary/real ratio of the Wannier functions and the interpolated band structure.

The whole simulation is stored in the form of a graph, representing explicitly the provenance of the data generated including all inputs and outputs of the codes used in the workflow. An example of a provenance graph, automatically generated by AiiDA when running a Quantum ESPRESSO calculation followed by a Wannier90 calculation, is shown in Fig. 8.

We emphasise that the availability of a platform to run Wannier90 in a fully-automated high-throughput way via the AiiDA-Wannier90 plugin has already proved to be beneficial for the Wannier90 code itself. Indeed, it has pushed the development of additional features or improvements now part of Wannier90 v3.0, including additional output files to facilitate output parsing and improvements in some of the algorithms and their default parameters to increase robustness.

VIII Modern software engineering practices

In this section, we describe a number of modern software engineering practices that are now part of the development cycle of the Wannier90 code. In particular, Wannier90 includes a number of tests that are run at every commit via a continuous integration approach, as well as nightly in a dedicated test farm. Version control is handled using git and the code is hosted on the GitHub platform GitHub . We follow the fork and pull-request model, in which users can duplicate (fork) the project into their own private repository, make their own changes, and make a pull request (i.e., request that their changes be incorporated back into the main repository). When a pull request is made, a series of tests are automatically performed: the test suite is run both in serial and parallel using the Travis continuous integration platform travis , and code coverage is checked using codecov codecov . If these tests are successful then the changes are reviewed by members of the Wannier90 developers group and, if the code meets the published coding guidelines, it can be merged into the development branch.

In addition, while interaction with end users happens via a mailing-list forum, discussion among developers is now tracked using GitHub issues. This facilitates the maintenance of independent conversation threads for each different code issue, new feature proposal or bug. These can easily reference code lines as well as be referenced in code commit messages. Moreover, for every new bug report a new issue is opened, and pull requests that close the issue clearly refer to it. This approach facilitates tracking back the reasoning behind the changes in case a similar problem resurfaces.

In the remainder of this section we describe more in detail some of these modern software engineering practices.

VIII.1 Code documentation (FORD)

The initial release of Wannier90 came with extensive documentation in the form of a User Guide describing the methodology, input flags to the program and format of the input and output files. This document was aimed at the end users running the software. Documentation of the code itself was done via standard code comments. In order to foster not only a community of users but also of code contributors to Wannier90, we have now created an additional documentation of the internal structure of the code. This makes the code more approachable, particularly for new contributors. To create this code documentation in a fully automated fashion, we use the FORD (FORtran Documenter)ford-docs documentation generator. We have chosen this over other existing documentation solutions because of FORD’s specific support for Fortran. This tool parses the Fortran source, and generates a hyperlinked (HTML) index of source files, modules, procedures, types and programs defined in the code. Furthermore, it constructs graphs showing the dependencies between different modules and subroutines. Additional information can be provided in the form of special in-code comments (marked with double exclamation marks) describing in more detail variables, modules or subroutines. By tightly coupling the code to its documentation using in-code comments, the documentation maintenance efforts are greatly reduced, decreasing the risk of having outdated documentation. The compiled version of the documentation for the most recent code version is made available on the Wannier90 website aiidawan-ford-docs .

VIII.2 Testing infrastructure and continuous integration

With the recent opening to the community of the Wannier90 development, it has become crucial to create a non-regression test suite to ensure that new developments do not break existing functionalities of the code. Its availability facilitates the maintenance of the code and ensures its long-term stability.

The Wannier90 test suite relies on a modified version of James Spencer’s python testcode.py Spencer-testcode . This provides the functionality to run tests and compare selected quantities parsed from the output files against benchmarked values.

At present, the Wannier90 test suite includes over 50 tests which are run both in serial and parallel and cover over 60% of the source code (with many modules exceeding 80% coverage). The code coverage is measured with the codecov software codecov . Developers are now required to add tests when adding new features to the code to ensure that their additions work as expected. This also ensures that future changes to the code will never break that functionality. Two different test approaches are implemented, serving different purposes.

First, the Wannier90 repository is now linked with the Travis continuous integration platform travis to prevent introducing errors and bugs into the main code branch. Upon any commit to the GitHub repository, the test suite is run both in serial and in parallel. Any test failure is reported back to the GitHub webpage. Additionally, for tests run against pull requests, any failed test results in the pull request being blocked and not permitted to merge. Contributors will first need to change their code to fix the problems highlighted in the tests; pull requests are able to be merged only after all tests pass successfully.

Second, nightly automatic tests are run on a Buildbot test-farm. The test-farm compiles and runs the code with a combination of compilers and libraries (current compilers include GFortran v6.4.0 and v7.3.0, Intel Fortran Compiler v17 and v18, and PGI compiler v18.05; current MPI libraries include Open MPI v1.10.7 and v3.1.3, Intel MPI v17 and MVAPICH v2.3b). This ensures that the code runs correctly on various high-performance computer (HPC) architectures. More information on the test-farm can be found on the Wannier90 GitHub wiki websitegithub-wiki .

In addition to these tests, we have implemented git pre-commit hooks to help keep the same code style in all source files. The current pre-commit hooks run Patrick Seewald’s Fortran source code formatter fprettify fprettify to remove trailing whitespaces at the end of a line and to enforce a consistent indentation style. These precommit hooks, besides validating the code, can reformat it automatically. Developers may simply run the formatter code to convert the source to a valid format. If a developer installs the pre-commit hooks, these will be run automatically before every commit. Even if this is not the case, these tests are also run on Travis; therefore, a pull request that does not conform to the standard code style cannot be merged before the style is fixed.

VIII.3 Command-line interface and dry-run

The command-line interface of the code has been improved. Just running wannier90.x without parameters shows a short explanation of the available command line options. In addition, a -v flag has been added to print the version of the code, as well as a new -d dry-run mode, that just parses the input file to perform all needed checks of the inputs without running the actual calculation. The latter functionality is particularly useful to be used in input validators for Wannier90 or to precalculate quantities computed by the code at the beginning of the simulation (such as nearest-neighbour shells, $b$-vectors or expected memory usage) and use this information to validate the run or optimise it (e.g., to decide the parallelisation strategy within automated AiiDA workflows).

VIII.4 Library mode

Wannier90 also comes with a library mode, where the core code functionality can be compiled into a library that can then be linked by external programs. This library mode is used as the default interaction protocol by some interface codes. The library mode provides only support for a subset of the full functionality, in particular at the moment it only supports serial execution. We have now added and improved support for the use of excluded bands also within the library mode. Moreover, beside supporting the generation of a statically-linked library, we now also support the generation of dynamically-linked versions. Finally, we have added a minimal test code, run together with all other tests in the test suite, that serves both to verify that the library functionality works as expected, and as an example of the interface of the library mode.

IX Conclusions and outlook

Wannier90 v2.0 was released in October 2013 with a small update for v2.1 in January 2017. The results and developments of the past years, presented in this work, were released in Wannier90 v3.0 in February 2019. Thanks to the transition of Wannier90 to a community code, Wannier90 includes now a large number of new functionalities and improvements that make it very robust, efficient and rich with features. These include the implementation of new methods for the calculation of WFs and for the generation of the initial projections; parallelisation and optimisations; interfaces with new codes, methods and infrastructures; new user functionality; improved documentation; and various bug fixes. The effect of enlarging the community of developers is not only visible in the large number of contributions to the code, but also in the modern software engineering practices that we have put in place, that help improve the robustness and reliability of the code and facilitate its maintenance by the core Wannier90 developers group and its long-term sustainability.

The next major improvement that we are planning is the implementation of a more robust and general library mode. The features that we envision are: (1) the possibility to call the code from C or Fortran codes without the need to store files but by passing all variables from memory; (2) a more general library interface that is easily extensible in the future when new functionality is added; and (3) the possibility to run Wannier90 from a parallel MPI code, both by running each instance in parallel and by allowing massively-parallel codes to call, in parallel, various instances of Wannier90 on various structures or with different parameters. This improvement will demand a significant restructuring of most of the codebase and requires a good design of the new interface. Currently we are drafting the new library interface, by collecting feedback and use cases from the various contributors and users of the code, to ensure that the new library mode can be beneficial to all different possible use cases.

Acknowledgements.

We acknowledge code contributions by Daniel Aberg (w90pov code), Lampros Andrinopoulos (w90vdw code), Pablo Aguado Puente (gyrotropic.F90 module), Raffaello Bianco ($k$-slice plotting), Marco Buongiorno Nardelli (dosqc v1.0 subroutines upon which some of transport.F90 is based), Stefano de Gironcoli (pw2wannier90.x interface to Quantum ESPRESSO), Pablo Garcia Fernandez (matrix elements of the position operator), Nicholas Hine (w90vdw code), Young-Su Lee (specialised $\Gamma$-point routines and transport), Antoine Levitt (preconditioning), Graham Lopez (extension of pw2wannier90.x to add terms needed for orbital magnetisation), Radu Miron (constrained centres), Nicolas Poilvert (transport routines), Michel Posternak (original plotting routines), Rei Sakuma (symmetry-adapted Wannier functions), Gabriele Sclauzero (disentanglement in spheres in $k$-space), Matthew Shelley (transport routines), Christian Stieger (routine to print the $U$ matrices), David Strubbe (various bug fixes and improvements), Timo Thonhauser (extension of pw2wannier90.x to add terms needed for orbital magnetisation), as well as the participants of the first Wannier90 developers meeting in San Sebastián (Spain) in 2016 for useful discussions (Daniel Fritsch, Victor Garcia Suarez, Pablo Garcia Fernandez, Jan-Philipp Hanke, Ji Hoon Ryoo, Jürg Hutter, Javier Junquera, Liang Liang, Michael Obermeyer, Gianluca Prandini, Christian Stieger, Paolo Umari). The WDG acknowledges financial support from the NCCR MARVEL of the Swiss National Science Foundation, the European Union’s Centre of Excellence E-CAM (grant no. 676531), and the Thomas Young Centre for Theory and Simulation of Materials for support (grant no. TYC-101).