Tutorial 4 Localized Wannier Orbitals

Introduction

This tutorial describes how to use the ASE Wannier class for making Wannier orbitals. The Wannier orbitals are maximally localized orbitals constructed from a set of extended bloch eigenstates.

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

At the moment localized Wannier orbitals can be made from dacapo and GPAW.

For a full description of the Wannier class including a description of how a calculator can interface the Wannier class, see the Wannier class section in the ASE manual.

Gamma-points calculations

The first examples are for molecular systems which can be treated in a large supercell where a gamma-point sampling of the Brillouin zone is sufficient. The Wannier scheme is applied to an ethylene molecule and a chain of Pt atoms.

The ethylene molecule

In this tutorial we will make localized Wannier orbitals for ethylene (C2H4).

The ethylene molecule is defined below:

from ASE import ListOfAtoms, Atom
a = 6.0  # Size of unit cell (Angstrom)

atoms = ListOfAtoms([
                       Atom('H', (-1.235, 0.936 , 0 ),tag=0),
                       Atom('H', ( 1.235,-0.936 , 0 ),tag=1),
                       Atom('H', ( 1.235, 0.936 , 0 ),tag=1),
                       Atom('H', (-1.235,-0.936 , 0 ),tag=1),
                       Atom('C', ( 0.660, 0.000 , 0 ),tag=1),
                       Atom('C', (-0.660, 0.000 , 0 ),tag=1)],
                       cell=(a,a,a), periodic=True)

Now a calculator must be defined, in this tutorial we will make localized Wannier orbitals from the dacapo planewave code.

The calculator for the dacapo code is defined as follows:

from Dacapo import Dacapo
calc = Dacapo(planewavecutoff=400,nbands=8, xc='PBE',out='ethylene.nc')
atoms.SetCalculator(calc)
energy = atoms.GetPotentialEnergy()

We now assume that you have a ethylene.nc file from dacapo. You can now find the 6 localized Wannier orbitals from the 6 occupied eigenstates:

from ASE.Utilities.Wannier.Wannier import Wannier
atoms = Calculator.ReadAtoms('ethylene.nc')  # Dacapo.ReadAtoms(..) for dacapo
calc = atoms.GetCalculator()
wannier = Wannier(numberofwannier=6,calculator=calc)
wannier.Localize()

The centers and radii of the Wannier orbitals can be obtained in the following way:

centers=wannier.GetCenters()
for center in centers:
    print center

One can also get the centers for the Wannier orbitals as a ListOfAtoms and plot them, togeher with the atoms, using:

from ASE.Visualization.VMD import VMD
centers = wannier.GetCentersAsAtoms()
VMD(atoms,centers)

or if you want to use RasMol:

from ASE.Visualization.RasMol import RasMol
centers = wannier.GetCentersAsAtoms()
centers.extend(atoms)
rasmol = RasMol(centers)

The resulting plot is shown below (the red spheres indicate the Wannier centers):

frontview sideview

For a 3D isosurface plot of the first Wannier orbital one can use:

from ASE.Visualization.VTK import VTKPlotElectronicState,VTKPlotAtoms
state = wannier.GetElectronicState(0)
plot = VTKPlotElectronicState(state)
plot.SetRepresentationToIsoSurface2([0.1])
atomsplot = VTKPlotAtoms(atoms,parent=plot)
plot.Update()

The isosurface plots for one of the four C-H wannier functions and one of the two C==C wannier functions are shown below:

one of four C-H bond One of two double bonds
wannier1 wannier2

The above example constructs 6 Wannier functions from the 6 occupied eigenstates. This results in double-bond Wannier functions mixing pi and sigma character. A set of Wannier functions which separates pi and sigma character is obtained by choosing 7 Wannier functions instead of 6:

from ASE.Utilities.Wannier.Wannier import Wannier
atoms = Calculator.ReadAtoms('ethylene.nc')  # Dacapo.ReadAtoms(..) for dacapo
calc = atoms.GetCalculator()
wannier = Wannier(numberofwannier=7,calculator=calc)
wannier.Localize()

In this case the double bond Wannier functions are replaced by two p-like orbitals centered at each of the C atoms and a single sigma orbital centered in the C-C bond (see figures below):

The single C-C bond One of two p-like orbitals
ethylene1 ethylene2

A Platinum wire

In this example we construct Wannier functions for an infinite chain of Pt atoms. We use a supercell containing 4 Pt atoms and sample the Brillouin zone at the gamma point only. The wire and a calculator are defined below:

# Pt wire
lat = 2.41
atoms = ListOfAtoms([
      Atom('Pt', ([0*lat,2*lat,2*lat]),tag=0),
      Atom('Pt', ([1*lat,2*lat,2*lat]),tag=0),
      Atom('Pt', ([2*lat,2*lat,2*lat]),tag=0),
      Atom('Pt', ([3*lat,2*lat,2*lat]),tag=0)],
      cell=([4*lat,4*lat,4*lat]), periodic=True)

      # Dacapo calculator:
      calc = Dacapo(planewavecutoff=350, nbands=30, xc='PBE',out='pt4.nc')
      atoms.SetCalculator(calc)
      energy = atoms.GetPotentialEnergy()

Each Pt has 10 valence electrons so the unitcell contains a total of 40 electrons accomodated in the lowest 20 eigenstates. The valence bands of an infinite Pt wire consists of five d-bands and a single s-band. To describe these bands we need minimum of 6 Wannier functions per atom, and we therefore calculate for 4*(5+1)=24 Wannier orbitals. We furthermore would like the Wannier functions to reproduce all eigenstates with an energy below 2 eV above the Fermi level exactly. The Wannier class is initialized like this:

from ASE.Utilities.Wannier.Wannier import Wannier
atoms = Dacapo.ReadAtoms('pt4.nc')
calc = atoms.GetCalculator()
wannier = Wannier(numberofwannier=24,calculator=calc,occupationenergy=2.0)
wannier.SaveZIBlochMatrix('zibloch.pickle')
wannier.Localize()

The following 3 lines show how to save the localization matrix (needed to carry out the localization), run the actual localization procedure, and finally save the resulting unitary transformation that produces the Wannier functions from the original eigenstates:

wannier.SaveZIBlochMatrix('zibloch.pickle')
wannier.Localize()
wannier.SaveRotation('Ptwire')

In a succeeding run we can then (after having prepared a Wannier object), read the localization matrix and the unitary transformation:

wannier.ReadZIBlochMatrix('zibloch.pickle')
wannier.ReadRotation('Ptwire')

For each Pt atom we obtain 5 d-like orbitals centered on each Pt atom, and 1 sigma-orbital centered between each Pt-Pt pair. The wire and the 24 Wannier centers are shown below:

wire

The Platinum atoms are here shown in transparent red and the Wannier centers are shown in green. This plot is constructed using:

from ASE.Visualization.VMD import VMD
centers = wannier.GetCentersAsAtoms()
VMD(atoms,centers)

Below is shown an example of a d-orbital and a sigma-orbital:

A d orbital centered on a Pt atom A sigma-orbital centered between a Pt-Pt pair
dlike slike

One on these plot is constructed like this:

from ASE.Visualization.VTK import VTKPlotElectronicState,VTKPlotAtoms
state = wannier.GetElectronicState(0)
plot = VTKPlotElectronicState(state)
plot.SetRepresentationToIsoSurface2([0.1])
atomsplot = VTKPlotAtoms(atoms,parent=plot)
plot.Update()

K-point calculations

This section describes how the Wannier class is used to generate Wannier functions in the general case where a gamma-point description is not appropriate.

Introduction

In the first example we re-visit the the infinite Pt wire which was treated within the gamma point approximation in the last section. Here we take advantage of the periodicity of the system and construct the Wannier functions from a k-point sampled calculation of the electronic structure. This also allows us to calculate a band structure of the wire based on the Wannier functions. As a second example we consider ferromagnetic Iron. This example illustrates the method for a system periodic in three directions, and shows how to obtain Wannier functions for spin polarized systems.

NOTE: In general the k-point grid must be uniform and include the Gamma-point, i.e. a Monkhorst-Pack grid with an odd number of points in each direction. It is important that the DFT program does not use any symmetry operations to reduce k-point grid.

The Platinum wire re-visited

We use a unit cell containg a single Pt atom. We include 8 bands in the calculation (this is enough to include the 6 valence bands) and sample the Brillouin zone by 15 k-points along the wire axis. The calculator is defined below:

# Pt wire
from ASE import Atom,ListOfAtoms
from Dacapo import Dacapo
lat = 2.41
atoms = ListOfAtoms([
      Atom('Pt', ([0*lat,2*lat,2*lat]),tag=0)],
      cell=([lat,4*lat,4*lat]), periodic=True)

# Dacapo calculator:
calc = Dacapo(planewavecutoff=350, nbands=8,kpts=(15,1,1),xc='PBE',out='pt1_15kpt.nc')
atoms.SetCalculator(calc)

# Displace kpoints sligthly, so that the symmetry program in dacapo does
# not use inversion symmetry to reduce kpoints.
kpoints = calc.GetBZKPoints()
kpoints[[,0]] += 2e-5
calc.SetBZKPoints(kpoints)

energy = atoms.GetPotentialEnergy()

To get a minimal set of Wannier functions describing the 6 Pt valence bands we choose to calculate 6 Wannier functions (per unit cell). Also we would like to include all eigenstates below 2 eV above the Fermi level in the construction of the Wannier functions and thus we set the parameter 'occupationenergy'=2:

from ASE.Utilities.Wannier.Wannier import Wannier
atoms = Dacapo.ReadAtoms('pt1_15kpt.nc')
calc = atoms.GetCalculator()
wannier = Wannier(numberofwannier=6,calculator=calc,occupationenergy=2.0)
wannier.Localize()

Since the start guess at the Wannier functions is a random rotation of the eigenstates, the resulting Wannier functions will be localized in random unit cells between unit cell 0 and unit cell 14 (defined by the number of k-points). To translate all 6 Wannier functions to the unit cell next to the one at the origin we write:

wannier.TranslateAllWannierFunctionsToCell([1,0,0])

Here the vector [n,m,k] spefifies the target unit cell in scaled coordinates, i.e. n,m,k are integers. The resulting Wannier functions can now be written to a cube file. To get a plot of Wannier function number 3 in a box containing 5 unit cells along the wire we write:

wannier.WriteCube(3,'Ptwire.cube',repeat=(5,1,1))

The obtained Wannier functions are similar to those obtained in the gamma-point calculation in the preceding section, i.e. 5 d-like orbitals centered at the Pt and a single sigma orbital centered at the Pt-Pt bond.

Instead of using a random start guess (which is default), one can specify a start guess in terms of atomic-like orbitals. The same calculation as above using a start guess of 5 d-orbitals centered at atom number 0 (this is the only possibility since there is one atom in the unit cell), and 1 s-orbital centered at the Pt-Pt bond center, reads like this:

initialwannier=[
[[0],2,0.4],
[[0.5,0.5,0.5],0,0.4]
]
wannier = Wannier(numberofwannier=6,calculator=calc,occupationenergy=2.0,initialwannier=initialwannier)
wannier.Localize()

In this case the resulting Wannier functions will be automatically centered in the unit cell at the origin.

To access the hopping parameters (Hamiltonian matrix elements) between the various Wannier functions or to calculate band diagrams based on the Wannier functions we initialize an instance of the HoppingParameter class:

from ASE.Utilities.Wannier.HoppingParameters import HoppingParameters
hop=HoppingParameters(wannier,couplingradius=7.0)

Here 'couplingradius' denotes the distance (in Angstrom) beyond which the coupling matrix element between two given orbitals is set to zero. To get the Hamiltonian matrix element between Wannier function 0 in cell [0,0,0] and Wannier function number 2 in the neighboring cell [1,0,0] we write:

hop.GetHoppingParameter([1,0,0],0,2)

Note, that the matrix element depends only on the relative distance between the two orbitals. Thus the expression above generally gives the matrix element between Wannier function 0 in cell [n,0,0] and orbital 2 in cell [n+1,0,0].

To construct a band diagram using our minimal basis of 6 Wannier functions per Pt atom, and plot them using Gnuplot we write:

k,band=hop.GetBandDiagram(50,[0,0,0],[0.5,0,0])
import Gnuplot
plot=Gnuplot.Gnuplot()
plot('set data style lines')
plot.plot(Gnuplot.Data(k,b[0]),Gnuplot.Data(k,b[1]),Gnuplot.Data(k,b[2]),Gnuplot.Data(k,b[3]),Gnuplot.Data(k,b[4]),Gnupl

System Message: WARNING/2 (<string>, line 394)

Literal block ends without a blank line; unexpected unindent.

ot.Data(k,b[5]))

Here 'k' contains the k-points in units of per Angstrom, and 'band' contain the energies in units of eV. To directly save the band structure as an NetCDF file write:

hop.WriteBandDiagramToNetCDFFile('band.nc',50,[0,0,0],[0.5,0,0])

The result is seen in the figure below:

bandstructure plot

Band structure for an infinite chain of Pt atoms. Red lines are the bands calculated using a minimal basis set consisting of 6 Wannier functions per Pt atom. The circles are the exact Dacapo eigenvalues.

We notice that there is some deviation between the exact and Wannier function bands. This is a result of the finite 'cutoffradius' which here was set to 7.0 Angstroms. By increasing this parameter one can in principle converge the Wannier function bands to the exact values for energies below the 'occupationenergy', which in this case was set to 2 eV above the Fermi level. Above the 'occupationenergy' the Wannier function bands may deviate from the exact bands.

Bandstructure of ferromagnetic Iron

In this example we construct the Wannier functions of ferromagnetic Iron. The example demonstrates how the Wannier function method is applied to spin polarized systems. The band structure obtained from the Wannier function basis set is compared with the exact bandstructure.

The python code below defines the BCC structure and a calculator, using a 11x11x11 Monkhorst-Pack kpoint set:

from Dacapo import Dacapo
from ASE import Atom,ListOfAtoms

lat = 2.87
atoms = ListOfAtoms([Atom('Fe', ([0,0,0]),magmom=2.32)])

cell = num.array([[-1/2.,1/2.,1/2.],[1/2.,-1/2.,1/2.],[1/2.,1/2.,-1/2.]])
cell *= lat
atoms.SetUnitCell(cell)

# Dacapo calculator:
calc = Dacapo(planewavecutoff=400,nbands=16,kpts=(11,11,11),xc='PBE',
              out='fe-bcc.nc', spinpol=True,usesymm=False)

atoms.SetCalculator(calc)

# displace kpoints sligthly, so that the symmetry program in dacapo does
# not use inversion symmetry to remove kpoints.
kpoints = calc.GetBZKPoints()
kpoints[[,0]] += 2e-5
calc.SetBZKPoints(kpoints)

energy = atoms.GetPotentialEnergy()

Note, that as usual the k-points are displaced so that they are not reduced by the symmetry program in Dacapo.

Now the eigenstates are contained in 'fe-bcc.nc' and we can proceed to generate the set of localized Wannier functions.

To describe the 6 valence bands (5 d- and 1 s-band) of iron we will construct 6 Wannier orbitals per Fe. We would like the Wannier functions to reproduce the exact bandstructure up to the Fermi level, and thus we set 'occupationenergy' to zero. If, for a given k-point in the Brillouin zone, the number of states below the Fermi level, is less than the number of Wannier orbitals (6 in this case), an optimal linear combination for the extra degress of freedom is selected.

The localized Wannier orbitals are calculated using:

from Dacapo import Dacapo
from ASE.Utilities.Wannier.Wannier import Wannier
import Numeric as num

atoms = Dacapo.ReadAtoms('fe-bcc.nc')
calc = atoms.GetCalculator()

initialwannier = [[[0],2,0.4]]
spin = 0
nwannier = 6
wannier = Wannier(numberofwannier=nwannier,
                  calculator=calc,
                  occupationenergy=0.0,
                  initialwannier=initialwannier,
                  spin=spin)


wannier.Localize(tolerance=1.0e-8)
print 'Localization ',wannier.GetFunctionalValue()

wannier.TranslateAllWannierFunctionsToCell((3,3,3))

centers = wannier.GetCenters()
for n in range(len(centers)):
   print n,centers[n]

centers = wannier.GetCentersAsAtoms()
traj = NetCDFTrajectory('centers-febcc-'+str(spin)+'.traj',centers)
traj.Update()
traj.Close()


for n in range(nwannier):
   wannier.WriteCube(n,'fe_'+str(n)+'_spin_'+str(spin)+'.cube',
                     repeat=(7,7,7))

This script will produce the following output:

0 {'radius': 0.104295, 'pos': array([ 4.3007,  4.3100,  4.3660])}
1 {'radius': 0.094142, 'pos': array([ 4.3062,  4.3050,  4.2761])}
2 {'radius': 0.081718, 'pos': array([ 4.3049,  4.2933,  4.3394])}
3 {'radius': 0.116333, 'pos': array([ 4.3105,  4.3367,  4.3053])}
4 {'radius': 0.105191, 'pos': array([ 4.3012,  4.3046,  4.2404])}
5 {'radius': 0.311760, 'pos': array([ 4.2728,  5.5450,  3.5339])}

From the radius and positions of the Wannier orbitals, Wannier orbital 5 is s-like, in an interstitial positions, while orbitals 0-4 are d-like, centered at the Fe atom.

A Cube formatted file is written for each Wannier orbital. The Cube file is generated by repeating the original unit cell 5 times, in each direction. The default repetition will here be 11 times, equal to the number of kpoints. This will correspond to the maximum periodicity of the Wannier orbitals, however due to the localized nature of the Wannier orbitals and for saving disk space we chose a smaller repetition. The Wannier orbitals are initially placed in cell (3,3,3). This should place the Wannier orbitals well within the 7x7x7 cell.

The BCC Brillouin zone is shown in the figure below:

BCC Brillouin zone

The high symmetry points H,P and N (see figure above) are defined, in the basis of the reciprocal unit vectors:

G = (0,0,0)
P = (0.25,0.25,0.25)
N = (0.0,0.0,0.5)
H = (-0.5,0.5,0.5)

The bandstructure can now be constructed, using the HoppingParameters class:

fermilevel = calc.GetFermilevel()
G =  (0,0,0)
P =  (0.25,0.25,0.25)
N =  (0.0,0.00,0.5)
H =  (-0.5,0.5,0.5)

from ASE.Utilities.Wannier import HoppingParameters

npoints = 80
cutoff = 12.0

hop=HoppingParameters.HoppingParameters(wannier,cutoff)
hop.WriteBandDiagramToNetCDFFile('GH.nc',npoints,G,H)
hop.WriteBandDiagramToNetCDFFile('HP.nc',npoints,H,P)
hop.WriteBandDiagramToNetCDFFile('PN.nc',npoints,P,N)
hop.WriteBandDiagramToNetCDFFile('NG.nc',npoints,N,G)
hop.WriteBandDiagramToNetCDFFile('GP.nc',npoints,G,P)

In order to describe completely localized Wannier orbitals, we must truncate the hopping matrix elements beyond a certain radius specified by the keyword 'cutoff' (in Angstroms). This is necessary since the Wannier orbitals by construction are periodic in the repeated cell (remember the repeated cell is constructed by repeating the unit cell basis vector i (i=1,2,3) by the number of k-points in direction i). Consequently, to avoid unphysical coupling between repeated images, the 'cutoff' parameter should be chosen < 0.5*(length of smallest basis vector of the repeated cell). In this case we use 11 k-points in all three directions and the length of the unit cell basis vectors is sqrt(3)*2.87. Thus we should have: 'cutoff'<0.5*(11*sqrt(3)*2.87)=13.67 Angstroms.

The figure below show the resulting band structure plot. The exact band structure from dacapo is made using the same method as in the dacapo example calculating a band diagram.

bandstructure plot

Majority spin bandstructure of ferromagnetic Iron. Black line are calculated using the set of localized Wannier orbitals. Red dots are the exact dacapo bandstructure.

Notice that the Wannier function band structure follows the exact bands perfectly below the 'occupationenergy' which we chose to be the Fermi level in this example. Above this energy there are deviations, indicating that the localization algorithm has selected non/trivial mixtures of Bloch states as the extra degrees of freedom.

The figure below show Wannier orbital 0 (red) and Wannier orbital 5 (blue), plottet together with the atoms (orange):

one d orbital centered on the Fe atom s like orbital in an interstitial position
fe-dlike fe-slike

The complete script can be found in the dacapo examples page.

System Message: WARNING/2 (<string>, line 648)

Explicit markup ends without a blank line; unexpected unindent.

http://arxiv.org/pdf/cond-mat/0506487