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:

Diamond cubic crystal structure of silicon.

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

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='out_task6_Sibulk.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)

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

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