Partly occupied Wannier Functions

from gpaw import GPAW

from ase.build import molecule

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

calc = GPAW(h=.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')

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)

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

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

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)) + [1e-5, 0, 0]
calc = GPAW(h=.21, xc='PBE', kpts=kpts, nbands=12, txt='poly.txt',
            eigensolver='cg', convergence={'bands': 9})

CC = 1.38
CH = 1.094
a = 2.45
x = a / 2.
y = np.sqrt(CC**2 - x**2)
atoms = Atoms('C2H2', pbc=(True, False, False), cell=(a, 8., 6.),
              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')

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(-.5, .5, 100):
        ham = wan.get_hamiltonian_kpoint([k, 0, 0])
        for eps in np.linalg.eigvalsh(ham).real:
            print(k, eps - ef, file=fd)

import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure(1, dpi=80, figsize=(4.2, 6))
fig.subplots_adjust(left=.16, right=.97, top=.97, bottom=.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([-.5, .5], [1, 1], 'k:', label='_nolegend_')
plt.text(-.5, 1, 'fixedenergy', ha='left', va='bottom')
plt.axis('tight')
plt.xticks([-.5, -.25, 0, .25, .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()