Maximally localized Wannier functions

This page describes how to construct the Wannier orbitals using the class Wannier. The page is organized as follows:

Introduction: A short summary of the basic theory.
The Wannier class : A description of how the Wannier class is used, and the methods defined within.

Introduction

The point of Wannier functions is the transform the extended Bloch eigenstates of a DFT calculation, into a smaller set of states designed to facilitate the analysis of e.g. chemical bonding. This is achieved by designing the Wannier functions to be localized in real space instead of energy (which would be the eigen states).

The standard Wannier transformation is a unitary rotation of the Bloch states. This implies that the Wannier functions (WF) span the same Hilbert space as the Bloch states, i.e. they have the same eigenvalue spectrum, and the original Bloch states can all be exactly reproduced from a linear combination of the WF. For maximally localized Wannier functions (MLWF), the unitary transformation is chosen such that the spread of the resulting WF is minimized.

The standard choice is to make a unitary transformation of the occupied bands only, thus resulting in as many WF as there are occupied bands. If you make a rotation using more bands, the localization will be improved, but the number of wannier functions increase, thus making orbital based analysis harder.

The class defined here allows for construction of partly occupied MLWF. In this scheme the transformation is still a unitary rotation for the lowest states (the fixed space), but it uses a dynamically optimized linear combination of the remaining orbitals (the active space) to improve localization. This implies that e.g. the eigenvalues of the Bloch states contained in the fixed space can be exactly reproduced by the resulting WF, whereas the largest eigenvalues of the WF will not necessarily correspond to any “real” eigenvalues (this is irrelevant, as the fixed space is usually chosen large enough, i.e. high enough above the fermilevel, that the remaining DFT eigenvalues are meaningless anyway).

The theory behind this method can be found in “Partly Occupied Wannier Functions” Thygesen, Hansen, and Jacobsen, Phys. Rev. Lett, Vol. 94, 26405 (2005), doi: 10.1103/PhysRevLett.94.026405.

An improved localization method based on penalizing the spread functional proportionally to the variance of spread distributions has since been published as “Spread-balanced Wannier functions: Robust and automatable orbital localization” Fontana, Larsen, Olsen, and Thygesen, Phys. Rev. B 104, 125140 (2021), doi: 10.1103/PhysRevB.104.125140.

The Wannier class

Usual invocation:

from ase.dft import Wannier
wan = Wannier(nwannier=18, calc=GPAW('save.gpw'), fixedstates=15)
wan.localize() # Optimize rotation to give maximal localization
wan.save('wannier.json') # Save localization and rotation matrix

# Re-load using saved wannier data
wan = Wannier(nwannier=18, calc=calc, fixedstates=15, file='wannier.json')

# Write a cube file
wan.write_cube(index=5, fname='wannierfunction5.cube')

For examples of how to use the Wannier class, see the Partly occupied Wannier Functions tutorial.

To enable the improved localization method from the 2021 paper, use Wannier(functional='var', ...). Computations from the 2021 paper were performed using Wannier(functional='var', initialwannier='orbitals', ...).

class ase.dft.wannier.Wannier(nwannier, calc, file=None, nbands=None, fixedenergy=None, fixedstates=None, spin=0, initialwannier='orbitals', functional='std', rng=<module 'numpy.random' from '/scratch/jensj/ase-docs/venv/lib/python3.11/site-packages/numpy/random/__init__.py'>, log=<function silent>)[source]

Partly occupied Wannier functions

Find the set of partly occupied Wannier functions according to Thygesen, Hansen and Jacobsen PRB v72 i12 p125119 2005.

Required arguments:

nwannier: The number of Wannier functions you wish to construct.
This must be at least half the number of electrons in the system and at most equal to the number of bands in the calculation. It can also be set to ‘auto’ in order to automatically choose the minimum number of needed Wannier function based on the fixedenergy / fixedstates parameter.

calc: A converged DFT calculator class.
If file arg. is not provided, the calculator must provide the method get_wannier_localization_matrix, and contain the wavefunctions (save files with only the density is not enough). If the localization matrix is read from file, this is not needed, unless get_function or write_cube is called.

Optional arguments:

nbands: Bands to include in localization.
The number of bands considered by Wannier can be smaller than the number of bands in the calculator. This is useful if the highest bands of the DFT calculation are not well converged.

spin: The spin channel to be considered.
The Wannier code treats each spin channel independently.

fixedenergy / fixedstates: Fixed part of Hilbert space.
Determine the fixed part of Hilbert space by either a maximal energy or a number of bands (possibly a list for multiple k-points). Default is None meaning that the number of fixed states is equated to nwannier. The maximal energy is relative to the CBM if there is a finite bandgap or to the Fermi level if there is none.

file: Read localization and rotation matrices from this file.

initialwannier: Initial guess for Wannier rotation matrix.
Can be ‘bloch’ to start from the Bloch states, ‘random’ to be randomized, ‘orbitals’ to start from atom-centered d-orbitals and randomly placed gaussian centers (see init_orbitals()), ‘scdm’ to start from localized state selected with SCDM or a list passed to calc.get_initial_wannier.

functional: The functional used to measure the localization.
Can be ‘std’ for the standard quadratic functional from the PRB paper, ‘var’ to add a variance minimizing term.

rng: Random number generator for initialwannier.

log: Function which logs, such as print().

distances(R)[source]

Relative distances between centers.

Returns a matrix with the distances between different Wannier centers. R = [n1, n2, n3] is in units of the basis vectors of the small cell and allows one to measure the distance with centers moved to a different small cell. The dimension of the matrix is [Nw, Nw].

get_centers(scaled=False)[source]

Calculate the Wannier centers

pos =  L / 2pi * phase(diag(Z))

get_function(index, repeat=None)[source]

Get Wannier function on grid.

Returns an array with the funcion values of the indicated Wannier function on a grid with the size of the repeated unit cell.

For a calculation using k-points the relevant unit cell for eg. visualization of the Wannier orbitals is not the original unit cell, but rather a larger unit cell defined by repeating the original unit cell by the number of k-points in each direction. Note that for a \(\Gamma\)-point calculation the large unit cell coinsides with the original unit cell. The large unitcell also defines the periodicity of the Wannier orbitals.

index can be either a single WF or a coordinate vector in terms of the WFs.

get_functional_value()[source]

Calculate the value of the spread functional.

Tr[|ZI|^2]=sum(I)sum(n) w_i|Z_(i)_nn|^2,

where w_i are weights.

If the functional is set to ‘var’ it subtracts a variance term

Nw * var(sum(n) w_i|Z_(i)_nn|^2),

where Nw is the number of WFs nwannier.

get_hamiltonian(k)[source]

Get Hamiltonian at existing k-vector of index k

        dag
H(k) = V    diag(eps )  V
        k           k    k

get_hamiltonian_kpoint(kpt_c)[source]

Get Hamiltonian at some new arbitrary k-vector

        _   ik.R
H(k) = >_  e     H(R)
        R

Warning: This method moves all Wannier functions to cell (0, 0, 0)

get_hopping(R)[source]

Returns the matrix H(R)_nm=<0,n|H|R,m>.

                      1   _   -ik.R
H(R) = <0,n|H|R,m> = --- >_  e      H(k)
                      Nk  k

where R is the cell-distance (in units of the basis vectors of the small cell) and n,m are indices of the Wannier functions.

get_optimal_nwannier(nwrange=5, random_reps=5, tolerance=1e-06)[source]

The optimal value for ‘nwannier’, maybe.

The optimal value is the one that gives the lowest average value for the spread of the most delocalized Wannier function in the set.

nwrange: number of different values to try for ‘nwannier’, the values will span a symmetric range around nwannier if possible.

random_reps: number of repetitions with random seed, the value is then an average over these repetitions.

tolerance: tolerance for the gradient descent algorithm, can be useful to increase the speed, with a cost in accuracy.

get_pdos(w, energies, width)[source]

Projected density of states (PDOS).

Returns the (PDOS) for Wannier function w. The calculation is performed over the energy grid specified in energies. The PDOS is produced as a sum of Gaussians centered at the points of the energy grid and with the specified width.

get_radii()[source]

Calculate the spread of the Wannier functions.

              --  /  L  \ 2       2
radius**2 = - >   | --- |   ln |Z|
              --d \ 2pi /

Note that this function can fail with some Bravais lattices, see \(get_spreads()\) for a more robust alternative.

get_spreads()[source]

Calculate the spread of the Wannier functions in Å². The definition is based on eq. 13 in PRB61-15 by Berghold and Mundy.

           / 1  \ 2  --                2
spread = - |----|    >   W_d  ln |Z_dw|
           \2 pi/    --d

initialize(file=None, initialwannier='random', rng=<module 'numpy.random' from '/scratch/jensj/ase-docs/venv/lib/python3.11/site-packages/numpy/random/__init__.py'>)[source]

Re-initialize current rotation matrix.

Keywords are identical to those of the constructor.

localize(step=0.25, tolerance=1e-08, updaterot=True, updatecoeff=True)[source]: Optimize rotation to give maximal localization

save(file)[source]: Save information on localization and rotation matrices to file.

translate(w, R)[source]

Translate the w’th Wannier function

The distance vector R = [n1, n2, n3], is in units of the basis vectors of the small cell.

translate_all_to_cell(cell=(0, 0, 0))[source]

Translate all Wannier functions to specified cell.

Move all Wannier orbitals to a specific unit cell. There exists an arbitrariness in the positions of the Wannier orbitals relative to the unit cell. This method can move all orbitals to the unit cell specified by cell. For a \(\Gamma\)-point calculation, this has no effect. For a k-point calculation the periodicity of the orbitals are given by the large unit cell defined by repeating the original unitcell by the number of k-points in each direction. In this case it is useful to move the orbitals away from the boundaries of the large cell before plotting them. For a bulk calculation with, say 10x10x10 k points, one could move the orbitals to the cell [2,2,2]. In this way the pbc boundary conditions will not be noticed.

translate_to_cell(w, cell)[source]: Translate the w’th Wannier function to specified cell

write_cube(index, fname, repeat=None, angle=False)[source]

Dump specified Wannier function to a cube file.

Parameters:

index – Integer, index of the Wannier function to save.
repeat – Array of integer, repeat supercell and Wannier function.
fname – Name of the cube file.
angle – If False, save the absolute value. If True, save the complex phase of the Wannier function.

Note

For calculations using k-points, make sure that the \(\Gamma\)-point is included in the k-point grid. The Wannier module does not support k-point reduction by symmetry, so you must use the usesymm=False keyword in the calc, and shift all k-points by a small amount (but not less than 2e-5 in) in e.g. the x direction, before performing the calculation. If this is not done the symmetry program will still use time-reversal symmetry to reduce the number of k-points by a factor 2. The shift can be performed like this:

from ase.dft.kpoints import monkhorst_pack
kpts = monkhorst_pack((15, 9, 9)) + [2e-5, 0, 0]