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