Skip to main content
`Ctrl`+`K`

# Partly occupied Wannier Functions

This tutorial walks through building partly occupied Wannier functions with the `ase.dft.wannier` module and the GPAW electronic structure code. For more information on the details of the method and the implementation, see

> K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen
> Partly occupied Wannier functions: Construction and applications
> Phys. Rev. B 72, 125119, (2005)

## Benzene molecule

### Step 1 – Ground-state calculation

Run the script below to obtain the ground-state density and the Kohn–Sham (KS) orbitals. The result is stored in `benzene.gpw`.

    from gpaw import GPAW

    from ase.build import molecule

    atoms = molecule('C6H6')
    atoms.center(vacuum=3.5)

    calc = GPAW(mode='fd', h=0.21, xc='PBE', txt='benzene.txt', nbands=18)
    atoms.calc = calc
    atoms.get_potential_energy()

    calc = calc.fixed_density(
        txt='benzene-harris.txt',
        nbands=40,
        eigensolver='cg',
        convergence={'bands': 35},
    )
    atoms.get_potential_energy()

    calc.write('benzene.gpw', mode='all')

### Step 2 – Maximally localized WFs for the occupied subspace (15 WFs)

There are 15 occupied bands in the benzene molecule. We construct one Wannier function per occupied band by setting `nwannier = 15`. By calling `wan.localize()`, the code attempts to minimize the spread functional using a gradient-descent algorithm. The resulting WFs are written to .cube files, which allows them to be inspected using e.g. VESTA.

    from gpaw import restart

    from ase.dft.wannier import Wannier

    atoms, calc = restart('benzene.gpw', txt=None)

    # Make wannier functions of occupied space only
    wan = Wannier(nwannier=15, calc=calc)
    wan.localize()
    for i in range(wan.nwannier):
        wan.write_cube(i, 'benzene15_%i.cube' % i)

### Step 3 – Adding three extra degrees of freedom (18 WFs)

To improve localization we augment the basis with three extra Wannier functions - so-called extra degrees of freedom (`nwannier = 18`, `fixedstates = 15`). This will allow the Wannierization procedure to use the unoccupied states to minimize spread functional.

    from gpaw import restart

    from ase.dft.wannier import Wannier

    atoms, calc = restart('benzene.gpw', txt=None)

    # Make wannier functions using (three) extra degrees of freedom.
    wan = Wannier(nwannier=18, calc=calc, fixedstates=15)
    wan.localize()
    wan.save('wan18.json')
    for i in range(wan.nwannier):
        wan.write_cube(i, 'benzene18_%i.cube' % i)

### Step 4 – Spectral-weight analysis

The script below projects the WFs on the KS eigenstates. You should see the 15 lowest bands perfectly reconstructed (weight ≃ 1.0) while higher bands are only partially represented.

    import matplotlib.pyplot as plt
    import numpy as np
    from gpaw import restart

    from ase.dft.wannier import Wannier

    atoms, calc = restart('benzene.gpw', txt=None)
    wan = Wannier(nwannier=18, calc=calc, fixedstates=15, file='wan18.json')

    weight_n = np.sum(abs(wan.V_knw[0]) ** 2, 1)
    N = len(weight_n)
    F = wan.fixedstates_k[0]
    plt.figure(1, figsize=(12, 4))
    plt.bar(
        range(1, N + 1),
        weight_n,
        width=0.65,
        bottom=0,
        color='k',
        edgecolor='k',
        linewidth=None,
        align='center',
        orientation='vertical',
    )
    plt.plot([F + 0.5, F + 0.5], [0, 1], 'k--')
    plt.axis(xmin=0.32, xmax=N + 1.33, ymin=0, ymax=1)
    plt.xlabel('Eigenstate')
    plt.ylabel('Projection of wannier functions')
    plt.savefig('spectral_weight.png')
    plt.show()

## Polyacetylene chain (1-D periodic)

We now want to construct partially occupied Wannier functions to describe a polyacetylene chain.

### Step 1 – Structure & ground-state calculation

Polyacetylene is modelled as an infinite chain; we therefore enable periodic boundary conditions along x.

    import numpy as np
    from gpaw import GPAW

    from ase import Atoms
    from ase.dft.kpoints import monkhorst_pack

    kpts = monkhorst_pack((13, 1, 1))
    calc = GPAW(
        mode='fd',
        h=0.21,
        xc='PBE',
        kpts=kpts,
        nbands=12,
        txt='poly.txt',
        eigensolver='cg',
        convergence={'bands': 9},
        symmetry='off',
    )

    CC = 1.38
    CH = 1.094
    a = 2.45
    x = a / 2.0
    y = np.sqrt(CC**2 - x**2)
    atoms = Atoms(
        'C2H2',
        pbc=(True, False, False),
        cell=(a, 8.0, 6.0),
        calculator=calc,
        positions=[[0, 0, 0], [x, y, 0], [x, y + CH, 0], [0, -CH, 0]],
    )
    atoms.center()
    atoms.get_potential_energy()
    calc.write('poly.gpw', mode='all')

### Step 2 – Wannierization

We repeat the localization procedure, keeping the five lowest bands fixed and adding one extra degree of freedom to aid localization.

    import numpy as np
    from gpaw import restart

    from ase.dft.wannier import Wannier

    atoms, calc = restart('poly.gpw', txt=None)

    # Make wannier functions using (one) extra degree of freedom
    wan = Wannier(
        nwannier=6,
        calc=calc,
        fixedenergy=1.5,
        initialwannier='orbitals',
        functional='var',
    )
    wan.localize()
    wan.save('poly.json')
    wan.translate_all_to_cell((2, 0, 0))
    for i in range(wan.nwannier):
        wan.write_cube(i, 'polyacetylene_%i.cube' % i)

    # Print Kohn-Sham bandstructure
    ef = calc.get_fermi_level()
    with open('KSbands.txt', 'w') as fd:
        for k, kpt_c in enumerate(calc.get_ibz_k_points()):
            for eps in calc.get_eigenvalues(kpt=k):
                print(kpt_c[0], eps - ef, file=fd)

    # Print Wannier bandstructure
    with open('WANbands.txt', 'w') as fd:
        for k in np.linspace(-0.5, 0.5, 100):
            ham = wan.get_hamiltonian_kpoint([k, 0, 0])
            for eps in np.linalg.eigvalsh(ham).real:
                print(k, eps - ef, file=fd)

### Step 3 – High-resolution band structure

Using the Wannier Hamiltonian we can interpolate the band structure on a fine 100-point k mesh and compare it to the original DFT result.

    import matplotlib.pyplot as plt
    import numpy as np

    fig = plt.figure(dpi=80, figsize=(4.2, 6))
    fig.subplots_adjust(left=0.16, right=0.97, top=0.97, bottom=0.05)

    # Plot KS bands
    k, eps = np.loadtxt('KSbands.txt', unpack=True)
    plt.plot(k, eps, 'ro', label='DFT', ms=9)

    # Plot Wannier bands
    k, eps = np.loadtxt('WANbands.txt', unpack=True)
    plt.plot(k, eps, 'k.', label='Wannier')

    plt.plot([-0.5, 0.5], [1, 1], 'k:', label='_nolegend_')
    plt.text(-0.5, 1, 'fixedenergy', ha='left', va='bottom')
    plt.axis('tight')
    plt.xticks(
        [-0.5, -0.25, 0, 0.25, 0.5],
        [r'$X$', r'$\Delta$', r'$\Gamma$', r'$\Delta$', r'$X$'],
        size=16,
    )
    plt.ylabel(r'$E - E_F\  \rm{(eV)}$', size=16)
    plt.legend()
    plt.savefig('bands.png', dpi=80)
    plt.show()

Within the fixed-energy window—that is, for energies below the fixed-energy line—the Wannier-interpolated bands coincide perfectly with the DFT reference (red circles). Above this window the match is lost, because the degrees of freedom deliberately mix several Kohn–Sham states to achieve maximal real-space localisation.