index | modules | gitlab | page source

# Partly occupied Wannier Functions

    from ase.build import molecule
    from gpaw import GPAW

    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.set(fixdensity=True, 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 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.pickle')
    for i in range(wan.nwannier):
        wan.write_cube(i, 'benzene18_%i.cube' % i)

    import numpy as np
    import matplotlib.pyplot as plt
    from ase.dft import Wannier
    from gpaw import restart

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

    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 ase import Atoms
    from ase.dft.kpoints import monkhorst_pack
    from gpaw import GPAW

    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 ase.dft import Wannier
    from gpaw import restart

    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)
    wan.localize()
    wan.save('poly.pickle')
    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 numpy as np
    import matplotlib.pyplot as plt

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