Band structure of bulk graphene

Graphene using grid

• graphene_grid.py:

from gpaw import GPAW, FermiDirac
from ase import Atoms
from ase.io import read, write
from gpaw import GPAW, PoissonSolver, Mixer
from ase.structure import bulk 


# -------------------------------------------------------------
# Bulk configuration
# -------------------------------------------------------------

gnr = bulk('C', 'hcp', a=2.4612, c=6.709) # hexagonal close-packed (hpc)
gnr.positions=([[ 0. ,  0.,  0.], [ 1.23060,  0.7104873,  0. ]])

write('gnr.traj', gnr)

# Make self-consistent calculation and save results
calc = GPAW(h=0.18,
            mode='fd', #finite difference(fd)
            xc='PBE',
            basis='szp(dzp)',
            kpts=(5,5,1),
            occupations=FermiDirac(width=0.05, maxiter=2000),
            mixer=Mixer(beta=0.010, nmaxold=8, weight=100.0),
            poissonsolver=PoissonSolver(eps=1e-12),
            txt='band_sc.txt')

gnr.set_calculator(calc)
gnr.get_potential_energy()
calc.write('band_sc.gpw')



from ase.dft.kpoints import ibz_points, get_bandpath
points = ibz_points['hexagonal']
G = points['Gamma']
K = points['K']
M = points['M']
kpts, x, X = get_bandpath([K, G, M, K], gnr.cell, 60)

calc = GPAW('band_sc.gpw',
            mode='fd',
            kpts=kpts,
            txt='band_harris.txt',
            fixdensity=True,
            parallel={'domain': 1},
            eigensolver='cg', # 'cg' is allowed for grid method only
            usesymm=None,
            convergence={'bands': 'all'})

if calc.input_parameters['mode'] == 'lcao':
    calc.scf.reset()


calc.get_potential_energy()
ef = calc.get_fermi_level()
calc.write('band_harris.gpw')

# Extract eigevanlues in a file for plotting with some external package

import numpy as np

calc = GPAW('band_harris', txt=None)
eps_skn = np.array([[calc.get_eigenvalues(k,s)
                     for k in range(40)]
                    for s in range(1)]) - ef

f = open('bands.dat', 'w')
for n in range(8):
    for k in range(40):
        print >>f, k, eps_skn[0, k, n]
    print >>f

Graphene using LCAO

• graphene_lcao.py:

from gpaw import GPAW, FermiDirac
from ase import Atoms
from ase.io import read, write
from gpaw import GPAW, PoissonSolver, Mixer
from ase.structure import bulk 

# -------------------------------------------------------------
# Bulk configuration
# -------------------------------------------------------------

gnr = bulk('C', 'hcp', a=2.4612, c=6.709)

gnr.positions=([[ 0. ,  0.,  0.], [ 1.23060,  0.7104873,  0. ]])

#gnr.center()
write('gnr.traj', gnr)

# Make self-consistent calculation and save results
calc = GPAW(h=0.18,
            mode='lcao',
            xc='PBE',
            basis='szp(dzp)',
            kpts=(5,5,1),
            occupations=FermiDirac(width=0.05, maxiter=2000),
            mixer=Mixer(beta=0.010, nmaxold=8, weight=100.0),
            poissonsolver=PoissonSolver(eps=1e-12),
            txt='band_sc.txt')

gnr.set_calculator(calc)
gnr.get_potential_energy()
calc.write('band_sc.gpw')

from ase.dft.kpoints import ibz_points, get_bandpath
points = ibz_points['hexagonal']
G = points['Gamma']
K = points['K']
M = points['M']
kpts, x, X = get_bandpath([K, G, M, K], gnr.cell, 40)

calc = GPAW('band_sc.gpw',
            mode='lcao',
            xc='PBE',
            basis='szp(dzp)',
            kpts=kpts,
            txt='band_harris.txt',
            fixdensity=True,
            parallel={'domain': 1},
            usesymm=None,
            convergence={'bands': 'all'})

if calc.input_parameters['mode'] == 'lcao':
    calc.scf.reset()

calc.get_potential_energy()
calc.write('band_harris.gpw')

• find_band_dat.py

from ase.lattice import bulk
from gpaw import GPAW

calc = GPAW('band_harris', txt=None)
import numpy as np
eps_skn = np.array([[calc.get_eigenvalues(k,s)
                     for k in range(40)]
                    for s in range(1)]) + 0.53129

for n in range(8):
    for k in range(40):
        print k, eps_skn[0, k, n]
    print