Crystals and band structure
In this tutorial we calculate properties of crystals.
Setting up bulk structures
ASE provides three frameworks for setting up bulk structures:
ase.build.bulk(). Knows lattice types and lattice constants for elemental bulk structures and a few compounds, but with limited customization.
ase.spacegroup.crystal(). Creates atoms from typical crystallographic information such as spacegroup, lattice parameters, and basis.
ase.lattice. Creates atoms explicitly from lattice and basis.
Let’s run a simple bulk calculation.
Exercise
Use ase.build.bulk() to get a primitive cell of silver,
then visualize it.
Silver is known to form an FCC structure, so presumably the function returned a primitive FCC cell. But it’s always nice to be sure what we have in front of us. Can you recognize it as FCC? You can e.g. use the ASE GUI to repeat the structure and recognize the A-B-C stacking. ASE should also be able to verify that it really is a primitive FCC cell and tell us what lattice constant was chosen:
print(atoms.cell.get_bravais_lattice())
Periodic structures in ASE are represented using atoms.cell and atoms.pbc.
The cell is a Cell object which represents
the crystal lattice with three vectors.
pbc is an array of three booleans indicating whether the
system is periodic in each direction.
Bulk calculation
For periodic DFT calculations we should generally use a number of
k-points which properly samples the Brillouin zone.
Many calculators including GPAW and Aims
accept the kpts keyword which can be a tuple such as
(4, 4, 4). In GPAW, the planewave mode
is very well suited for smaller periodic systems.
Using the planewave mode, we should also set a planewave cutoff (in eV):
from gpaw import GPAW, PW
calc = GPAW(mode=PW(600), kpts=(8, 8, 8),
setups={'Ag': '11'}, ...)
Here we have used the setups keyword to specify that we want the
11-electron PAW dataset instead of the default which has 17 electrons,
making the calculation faster.
(In principle, we should be sure to converge both kpoint sampling and planewave cutoff – I.e., write a loop and try different samplings so we know both are good enough to accurately describe the quantity we want.)
Exercise
Run a single-point calculation of bulk silver with GPAW. Save the
ground-state in GPAW’s own format using calc.write('Ag.gpw').
Density of states
Having saved the ground-state, we can reload it for ASE to extract the density of states:
import matplotlib.pyplot as plt
from ase.dft.dos import DOS
from gpaw import GPAW
calc = GPAW('groundstate.gpw')
dos = DOS(calc, npts=500, width=0)
energies = dos.get_energies()
weights = dos.get_dos()
plt.plot(energies, weights)
plt.show()
Calling the DOS class with width=0 means ASE well calculate the DOS using
the linear tetrahedron interpolation method, which takes time but gives a nicer
representation. We could also have given it
a nonzero width such as the default value of 0.1 (eV). In that case it
would have used a simple Gaussian smearing with that width, but we would need
more k-points to get a plot of the same quality.
Note that the zero point of the energy axis is the Fermi energy.
Exercise
Plot the DOS.
You probably recall that an Ag atom has 10 d electrons and one s electron.
Which parts of the spectrum do you think originate (mostly) from s electrons? And which parts (mostly) from d electrons?
Time for analysis. As we probably know, the d-orbitals in a transition metal atom are localized close to the nucleus while the s-electron is much more delocalized.
In bulk systems, the s-states overlap a lot and therefore split into a very broad band over a wide energy range. d-states overlap much less and therefore also split less: They form a narrow band with a very high DOS. Very high indeed because there are 10 times as many d electrons as there are s electrons.
So to answer the question, the d-band accounts for most of the states forming the big, narrow chunk between -6.2 eV to -2.6 eV. Anything outside that interval is due to the much broader s band.
The DOS above the Fermi level may not be correct, since the SCF convergence criterion (in this calculation) only tracks the convergenece of occupied states. Hence, the energies over the Fermi level 0 are probably wrong.
What characterises the noble metals Cu, Ag, and Au, is that the d-band is fully occupied. I.e.: The whole d-band lies below the Fermi level (energy=0). If we had calculated any other transition metal, the Fermi level would lie somewhere within the d-band.
Note
We could calculate the s, p, and d-projected DOS to see more conclusively which states have what character. In that case we should look up the GPAW documentation, or other calculator-specific documentation. So let’s not do that now.
Band structure
Let’s calculate the band structure of silver.
First we need to set up a band path. Our favourite image search engine can show us some reference graphs. We might find band structures from both Exciting and GPAW with Brillouin-zone path \(\mathrm{W L \Gamma X W K}\). Luckily ASE knows these letters and can also help us visualize the reciprocal cell:
lat = atoms.cell.get_bravais_lattice()
print(lat.description())
lat.plot_bz(show=True)
In general, the ase.lattice module provides
BravaisLattice classes used to represent each
of the 14 + 5 Bravais lattices in 3D and 2D, respectively.
These classes know
about the high-symmetry k-points and standard Brillouin-zone paths
(using the AFlow conventions).
Exercise
Build a band path for \(\mathrm{W L \Gamma X W K}\).
You can use path = atoms.cell.bandpath(...)
— see the Cell documentation for which
parameters to supply.
This gives us a BandPath object.
You can print() the band path object to see some basic information about it,
or use its write() method to save the
band path to a json file such as path.json. Then visualize it
using the command:
$ ase reciprocal path.json
Once we are sure we have a good path with a reasonable number of k-points, we can run the band structure calculation. How to trigger a band structure calculation depends on which calculator we are using, so we would typically consult the documentation for that calculator (ASE will one day provide shortcuts to make this easier with common calculators):
calc = GPAW('groundstate.gpw')
atoms = calc.get_atoms()
path = atoms.cell.bandpath(<...>)
calc = calc.fixed_density(kpts=path, symmetry='off')
We have here told GPAW to use our bandpath for k-points, not to perform symmetry-reduction of the k-points, and to fix the electron density.
Then we trigger a new calculation, which will be non-selfconsistent, and extract and save the band structure:
atoms.get_potential_energy()
bs = calc.band_structure()
bs.write('bs.json')
Again, the ASE command-line tool offers a helpful command to plot the band structure from a file:
$ ase band-structure bs.json
Exercise
Calculate, save, and plot the band structure of silver for the path \(\mathrm{W L \Gamma X W K}\).
You may need to zoom around a bit to see the whole thing at once. The plot will show the Fermi level as a dotted line (but does not define it as zero like the DOS plot before). Looking at the band structure, we see the complex tangle of what must be mostly d-states from before, as well as the few states with lower energy (at the \(\Gamma\) point) and higher energy (crossing the Fermi level) attributed to s.
Equation of state
We can find the optimal lattice parameter and calculate the bulk modulus by doing an equation-of-state calculation. This means sampling the energy and lattice constant over a range of values to get the minimum as well as the curvature, which gives us the bulk modulus.
The online ASE docs already provide a tutorial on how to do this, using the empirical EMT potential: https://wiki.fysik.dtu.dk/ase/tutorials/eos/eos.html
Exercise
Run the EOS tutorial.
Complex crystals and cell optimisation
(If time is scarce, please consider skipping ahead to do the remaining exercises before returning here.)
For the simple FCC structure we only have a single parameter, a, and the EOS fit tells us everything there is to know.
For more complex structures we first need a more advanced framework
to build the atoms, such as the
ase.spacegroup.crystal() function.
The documentation helpfully tells us how to build a rutile structure, saving us the trouble of looking up the atomic basis and other crystallographic information. Rutile is a common mineral form of TiO2
Exercise
Build and visualize a rutile structure.
Let’s uptimise the structure. In addition to the positions, we must optimise the unit cell which, being tetragonal, is characterised by the two lengths a and c.
Optimising the cell requires the energy derivatives with respect to the
cell parameters accessible through the stress tensor.
atoms.get_stress() calculates and returns the stress
as a vector of the 6 unique components (Voigt form). Using it requires
that the attached calculator supports the stress tensor. GPAW’s planewave
mode does this.
The ase.constraints.ExpCellFilter allows us to optimise
cell and positions simultaneously. It does this by exposing the degrees
of freedom to the optimiser as if they were additional positions —
hence acting as a kind of filter.
We use it by wrapping it around the atoms:
from ase.optimize import BFGS
from ase.constraints import ExpCellFilter
opt = BFGS(ExpCellFilter(atoms), ...)
opt.run(fmax=0.05)
Exercise
Use GPAW’s planewave mode to optimize the rutile unit cell. You will probably need a planewave cutoff of at least 500 eV. What are the optimised lattice constants a and c?
Exercise
Calculate the band structure of rutile. Does it agree with your favourite internet search engine?
Solutions
Ag ground state:
from gpaw import GPAW, PW
from ase.build import bulk
atoms = bulk('Ag')
calc = GPAW(mode=PW(350), kpts=[8, 8, 8], txt='gpaw.bulk.Ag.txt',
setups={'Ag': '11'})
atoms.calc = calc
atoms.get_potential_energy()
calc.write('bulk.Ag.gpw')
Ag DOS:
import matplotlib.pyplot as plt
from gpaw import GPAW
from ase.dft.dos import DOS
calc = GPAW('bulk.Ag.gpw')
dos = DOS(calc, npts=800, width=0)
energies = dos.get_energies()
weights = dos.get_dos()
ax = plt.gca()
ax.plot(energies, weights)
ax.set_xlabel('Energy [eV]')
ax.set_ylabel('DOS [1/eV]')
plt.savefig('dos.png')
plt.show()
Ag band structure:
from gpaw import GPAW
calc = GPAW('bulk.Ag.gpw')
atoms = calc.get_atoms()
path = atoms.cell.bandpath('WLGXWK', density=10)
path.write('path.json')
calc = calc.fixed_density(kpts=path, symmetry='off')
bs = calc.band_structure()
bs.write('bs.json')
Rutile cell optimisation:
from gpaw import GPAW, PW
from ase.constraints import ExpCellFilter
from ase.io import write
from ase.optimize import BFGS
from ase.spacegroup import crystal
a = 4.6
c = 2.95
# Rutile TiO2:
atoms = crystal(['Ti', 'O'], basis=[(0, 0, 0), (0.3, 0.3, 0.0)],
spacegroup=136, cellpar=[a, a, c, 90, 90, 90])
write('rutile.traj', atoms)
calc = GPAW(mode=PW(800), kpts=[2, 2, 3],
txt='gpaw.rutile.txt')
atoms.calc = calc
opt = BFGS(ExpCellFilter(atoms), trajectory='opt.rutile.traj')
opt.run(fmax=0.05)
calc.write('groundstate.rutile.gpw')
print('Final lattice:')
print(atoms.cell.get_bravais_lattice())
Rutile band structure:
from gpaw import GPAW
calc = GPAW('groundstate.rutile.gpw')
atoms = calc.get_atoms()
path = atoms.cell.bandpath(density=7)
path.write('path.rutile.json')
calc = calc.fixed_density(kpts=path, symmetry='off')
bs = calc.band_structure()
bs.write('bs.rutile.json')