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).

For the theory behind this method see the paper “Partly Occupied Wannier Functions” Thygesen, Hansen and Jacobsen, Phys. Rev. Lett, Vol. 94, 26405 (2005).

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('file.pickle') # Save localization and rotation matrix

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

# 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.

class ase.dft.wannier.Wannier(nwannier, calc, file=None, nbands=None, fixedenergy=None, fixedstates=None, spin=0, initialwannier='random', rng=<module 'numpy.random' from '/usr/lib/python3/dist-packages/numpy/random/__init__.py'>, verbose=False)

Maximally localized Wannier Functions

Find the set of maximally localized Wannier functions using the spread functional of Marzari and Vanderbilt (PRB 56, 1997 page 12847).

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.

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 Heilbert 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.

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, or a list passed to calc.get_initial_wannier.

rng: Random number generator for initialwannier.

verbose: True / False level of verbosity.

distances(R)

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)

Calculate the Wannier centers

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

get_function(index, repeat=None)

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()

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.

get_hamiltonian(k=0)

Get Hamiltonian at existing k-vector of index k

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

get_hamiltonian_kpoint(kpt_c)

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)

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_pdos(w, energies, width)

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()

Calculate the spread of the Wannier functions.

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

initialize(file=None, initialwannier='random', rng=<module 'numpy.random' from '/usr/lib/python3/dist-packages/numpy/random/__init__.py'>)

Re-initialize current rotation matrix.

Keywords are identical to those of the constructor.

localize(step=0.25, tolerance=1e-08, updaterot=True, updatecoeff=True)

Optimize rotation to give maximal localization

save(file)

Save information on localization and rotation matrices to file.

translate(w, R)

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])

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)

Translate the w’th Wannier function to specified cell

write_cube(index, fname, repeat=None, real=True)

Dump specified Wannier function to a cube file

In Dacapo, the inialwannier keyword can be a list as described below:

Setup an initial set of Wannier orbitals. initialwannier can set up a starting guess for the Wannier functions. This is important to speed up convergence in particular for large systems For transition elements with d electrons you will always find 5 highly localized d-orbitals centered at the atom. Placing 5 d-like orbitals with a radius of 0.4 Angstroms and center at atom no. 7, and 3 p-like orbitals with a radius of 0.4 Angstroms and center at atom no. 27 looks like this:

initialwannier = [[[7],2,0.4],[[27],1,0.4]]

Placing only the l=2, m=-2 and m=-1 orbitals at atom no. 7 looks like this:

initialwannier = [[[7],2,-2,0.4],[[7],2,-1,0.4]]

I.e. if you do not specify the m quantum number all allowed values are used. Instead of placing an orbital at an atom, you can place it at a specified position. For example the following:

initialwannier = [[[0.5,0.5,0.5],0,0.5]]

places an s orbital with radius 0.5 Angstroms at the position (0.5, 0.5, 0.5) in scaled coordinates of the unit cell.

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]