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$ Angstrom 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='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)

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

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 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.

import numpy as np
from ase import Atoms
from ase.structure 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.txt')

silicon.set_calculator(calc)
silicon.get_potential_energy()
calc.write('silicon.gpw', mode='all')