View source for Lattice constant, DOS, and band structure of Si

==Lattice constant of Si==

For bulk materials it is a common task in DFT calculations, given an atomic structure, to find the equilibrium volume of the unit cell. Here we study Silicon which has diamond lattice structure shown in the Figure below:


{| cellpadding="2" style="border:0px solid darkgray;"
|- border=0
| [[Image:si.png|thumb|200px|Diamond cubic crystal structure of silicon.]]
|}

To determine the equilibrium volume and thereby the lattice constant <math> a </math>, a series of calculations is made where <math> a </math> is swept over a range around some initial guess (below we start from <math> a = 5.43 </math> Angstrom and sweep the lattice parameter <math> \pm 5\% </math> in five steps).

*si_lattice.py
<pre>
import numpy as np
from ase import Atoms
from ase.lattice import bulk
from ase.io.trajectory import PickleTrajectory
from ase.optimize.bfgs import BFGS
from gpaw import GPAW
from gpaw import PW
from gpaw import FermiDirac

silicon = bulk('Si', 'diamond', a=5.459)

cell = silicon.get_cell()
traj = PickleTrajectory('silicon.traj', 'w')

calc=GPAW(mode=PW(200),                # Energycutoff for planewaves [eV]
          h=0.2,                       # The distance between gridpoints AA^-1
          xc='PBE',                    # xc-functional
          nbands=8,                    # number of bands
          kpts=(2,2,2),                # number of k-points
          occupations=FermiDirac(0.1), # Fermi temperature [eV]
          txt='sibulk_a.txt')
          
silicon.set_calculator(calc)

for x in np.linspace(0.95, 1.05, 7):
silicon.set_cell(cell * x, scale_atoms=True)
traj.write(silicon)
</pre>

After the calculation has finished, you may analyze the results using the following script:

*plot_eos.py
<pre>
from ase.io import read
from ase.units import kJ
from ase.utils.eos import EquationOfState

configs = read('silicon.traj@0:7') # read 7 configurations

# Extract volumes and energies:
volumes = [atoms.get_volume() for atoms in configs]
energies = [atoms.get_potential_energy() for atoms in configs]
eos = EquationOfState(volumes , energies)
v0, e0, B = eos.fit()

print v0, 'AA', B / kJ * 1.0e24, 'GPa'
eos.plot('Silicon-eos.png')
</pre>

The information is stored in '''silicon.traj''' file which we will use for post processing. Open the trajectory file with:

 ase-gui silicon.traj 

and choose ’Tools’ ! ‘Bulk modulus’. This will produce a fit of the energy vs. unit cell volume for the five calculations, the bulk modulus B and equilibrium volume V0 can be read from the plot title. Be sure to calculate the equilibrium lattice constant a from V0 as we will use it for all subsequent calculations on Silicon.

You can redo the above calculations by increasing the plane wave cutoff in steps of 100 eV until you reach 500 eV, and by increasing the k-point mesh in steps of two until you reach 12x12x12. This should allow you to check how the lattice constant and Bulk modulus converge with respect to the number of  plane waves and k-points.

==Density of States==

To calculate the [https://wiki.fysik.dtu.dk/ase/ase/dft/dos.html Density of States] (DOS),  we start from our previous calculation and increase the number
of k-points in a second calculation as the DOS converges slowly with respect to the k-point sampling.

*si_dos.py
<pre>
import numpy as np
from ase import Atoms
from ase.lattice import bulk
from ase.io.trajectory import PickleTrajectory
from ase.optimize.bfgs import BFGS
from gpaw import GPAW
from gpaw import PW
from gpaw import FermiDirac

silicon = bulk('Si', 'diamond', a=5.459)

calc=GPAW(mode=PW(400),                # Energycutoff for planewaves [eV]
          h=0.2,                       # The distance between gridpoints AA^-1
          xc='PBE',                    # xc-functional
          nbands=20,                   # number of bands
          kpts=(20,20,20),             # number of k-points
          occupations=FermiDirac(0.1), # Fermi temperature [eV]
          txt='sibulk_dos.txt')

silicon.set_calculator(calc)
silicon.get_potential_energy()
calc.write('silicon.gpw', mode='all')
</pre>

To plot the DOS, you can use ASE utility function:

*plot_dos.py
<pre>
from ase import *
from ase.dft.dos import DOS
from gpaw import GPAW , restart
import pylab as p

slab , calc = restart('silicon.gpw')
e, dos = calc.get_dos(spin=0, npts=2001, width=0.1)
e_f = calc.get_fermi_level()

p.subplot(211)
p.plot(e-e_f, dos)
p.grid(True)
p.axis([-15, 10, None , None])
p.ylabel('DOS')

p.subplot(212)
p.plot(e-e_f, dos)
p.grid(True)
p.axis([-1, 1, None , None])
p.ylabel('DOS')
p.show()
</pre>

You can repeat the same calculation using xc='LDA'. Note that neither the LDA nor the GGA are able to calculate the badgap with any high accuracy (there are
even examples were semiconductors becomes metallic when using local and semi-local exchange-correlation functionals).

==Electronic band structure==

To calculate the electronic band structure of Si, we also use the constant charge density mode in DFT and calculate the energy spectra along a number of high-symmetry lines in the Brillouin-zone using the following script:

*si_bands.py
<pre>
import pickle
import numpy as np
from ase.dft.kpoints import ibz_points , get_bandpath
from gpaw import GPAW

points = ibz_points['fcc']

G = points['Gamma']
X = points['X']
W = points['W']
K = points['K']
L = points['L']

calc = GPAW('silicon.gpw',
             txt='bands_out.txt',
             parallel={'domain': 1},
             fixdensity=True,
             usesymm=None,
             maxiter=150,
             convergence={'bands': 12})

kpts, x, X = get_bandpath([W, L, G, X, W, K], calc.atoms.cell)
calc.set(kpts=kpts)
calc.get_potential_energy()

e_kn = np.array([calc.get_eigenvalues(k) for k in range(len(kpts))])
pickle.dump((x, X, e_kn), open('eigenvalues.pckl', 'w'))
</pre>

The result can be post-processed to produce a plot using:

*plot_bands.py
<pre>
import pickle
import matplotlib.pyplot as plt

point_names = ['W', 'L', '\Gamma', 'X', 'W', 'K']
x, X, e_kn = pickle.load(open('eigenvalues.pckl'))
e_kn -= e_kn[:, :8].max()
emin = e_kn.min() - 1
emax = 1

plt.figure(figsize=(6, 4))
for n in range(8):
    plt.plot(x, e_kn[:, n])
for p in X:
    plt.plot([p, p], [emin , emax], 'k-')
plt.xticks(X, ['$%s$' % n for n in point_names])
plt.axis(xmin=0, xmax=X[-1], ymin=emin , ymax=emax)
plt.xlabel('k-vector')
plt.ylabel('Energy (eV)')
plt.title('PBE bandstructure of Silicon')
plt.savefig('bs-PBE.png')
</pre>