Back to the main Index
In this notebook we discuss how to plot electron band structures and density of states (DOS) using the netcdf files produced by Abinit.
For the tutorial, we will use the netcdf files shipped with AbiPy.
The function abidata.ref_file
returns the absolute path of the reference file.
In your scripts, you have to replace data.ref_file("abipy_filename")
with a string giving
the location of your netcdf file.
Alteratively, one can use the abiopen.py
script to open the file inside the shell with the syntax:
abiopen.py out_GSR.nc
This command will start the ipython interpreter so that one can interact directly
with the GsrFile
object (named abifile
inside ipython).
To generate a jupyter notebook use:
abiopen.py out_GSR.nc -nb
For a quick visualization of the data, use the --expose
options:
abiopen.py out_GSR.nc -e
# Use this at the beginning of your script so that your code will be compatible with python3
from __future__ import division, print_function, unicode_literals
import warnings
warnings.filterwarnings("ignore") # Ignore warnings
from abipy import abilab
abilab.enable_notebook() # This line tells AbiPy we are running inside a notebook
# Import abipy reference data.
import abipy.data as abidata
# This line configures matplotlib to show figures embedded in the notebook.
# Replace `inline` with `notebook` in classic notebook
%matplotlib inline
# Option available in jupyterlab. See https://github.com/matplotlib/jupyter-matplotlib
#%matplotlib widget
The GSR
file (mnemonics: Ground-State Results) is a netcdf file with the
results produced by SCF or NSCF ground-state calculations
(band energies, forces, energies, stress tensor).
To open a GSR
file, use the abiopen
function defined in abilab
:
gsr = abilab.abiopen(abidata.ref_file("si_scf_GSR.nc"))
The gsr object has a Structure
:
print(gsr.structure)
Full Formula (Si2) Reduced Formula: Si abc : 3.866975 3.866975 3.866975 angles: 60.000000 60.000000 60.000000 Sites (2) # SP a b c cartesian_forces --- ---- ---- ---- ---- ----------------------------------------------------------- 0 Si 0 0 0 [-5.89948306e-27 -1.93366149e-27 2.91016904e-27] eV ang^-1 1 Si 0.25 0.25 0.25 [ 5.89948306e-27 1.93366149e-27 -2.91016904e-27] eV ang^-1 Abinit Spacegroup: spgid: 227, num_spatial_symmetries: 48, has_timerev: True, symmorphic: True
and an ElectronBands
object with the band energies, the occupation factors, the list of k-points:
print(gsr.ebands)
================================= Structure ================================= Full Formula (Si2) Reduced Formula: Si abc : 3.866975 3.866975 3.866975 angles: 60.000000 60.000000 60.000000 Sites (2) # SP a b c cartesian_forces --- ---- ---- ---- ---- ----------------------------------------------------------- 0 Si 0 0 0 [-5.89948306e-27 -1.93366149e-27 2.91016904e-27] eV ang^-1 1 Si 0.25 0.25 0.25 [ 5.89948306e-27 1.93366149e-27 -2.91016904e-27] eV ang^-1 Abinit Spacegroup: spgid: 227, num_spatial_symmetries: 48, has_timerev: True, symmorphic: True Number of electrons: 8.0, Fermi level: 5.598 (eV) nsppol: 1, nkpt: 29, mband: 8, nspinor: 1, nspden: 1 smearing scheme: none (occopt 1), tsmear_eV: 0.272 Direct gap: Energy: 2.532 (eV) Initial state: spin: 0, kpt: [+0.000, +0.000, +0.000], weight: 0.002, band: 3, eig: 5.598, occ: 2.000 Final state: spin: 0, kpt: [+0.000, +0.000, +0.000], weight: 0.002, band: 4, eig: 8.130, occ: 0.000 Fundamental gap: Energy: 0.562 (eV) Initial state: spin: 0, kpt: [+0.000, +0.000, +0.000], weight: 0.002, band: 3, eig: 5.598, occ: 2.000 Final state: spin: 0, kpt: [+0.375, +0.375, +0.000], weight: 0.012, band: 4, eig: 6.161, occ: 0.000 Bandwidth: 11.856 (eV) Valence maximum located at: spin: 0, kpt: [+0.000, +0.000, +0.000], weight: 0.002, band: 3, eig: 5.598, occ: 2.000 Conduction minimum located at: spin: 0, kpt: [+0.375, +0.375, +0.000], weight: 0.012, band: 4, eig: 6.161, occ: 0.000 TIP: Call set_fermie_to_vbm() to set the Fermi level to the VBM if this is a non-magnetic semiconductor
A GSR file produced by a self-consistent run, contains the values of the total energy, the forces, and the stress tensor at the end of the SCF cycle:
print("energy:", gsr.energy, "pressure:", gsr.pressure)
energy: -241.23647031327783 eV pressure: -5.211617575719521 GPa
To get a summary of the most important results:
print(gsr)
================================= File Info ================================= Name: si_scf_GSR.nc Directory: /Users/gmatteo/git_repos/abipy/abipy/data/refs/si_ebands Size: 14.83 kb Access Time: Tue Oct 29 22:38:54 2019 Modification Time: Wed Mar 20 16:53:35 2019 Change Time: Wed Mar 20 16:53:35 2019 ================================= Structure ================================= Full Formula (Si2) Reduced Formula: Si abc : 3.866975 3.866975 3.866975 angles: 60.000000 60.000000 60.000000 Sites (2) # SP a b c cartesian_forces --- ---- ---- ---- ---- ----------------------------------------------------------- 0 Si 0 0 0 [-5.89948306e-27 -1.93366149e-27 2.91016904e-27] eV ang^-1 1 Si 0.25 0.25 0.25 [ 5.89948306e-27 1.93366149e-27 -2.91016904e-27] eV ang^-1 Abinit Spacegroup: spgid: 227, num_spatial_symmetries: 48, has_timerev: True, symmorphic: True Stress tensor (Cartesian coordinates in GPa): [[5.21161758e+00 7.86452261e-11 0.00000000e+00] [7.86452261e-11 5.21161758e+00 0.00000000e+00] [0.00000000e+00 0.00000000e+00 5.21161758e+00]] Pressure: -5.212 (GPa) Energy: -241.23647031 (eV) ============================== Electronic Bands ============================== Number of electrons: 8.0, Fermi level: 5.598 (eV) nsppol: 1, nkpt: 29, mband: 8, nspinor: 1, nspden: 1 smearing scheme: none (occopt 1), tsmear_eV: 0.272 Direct gap: Energy: 2.532 (eV) Initial state: spin: 0, kpt: [+0.000, +0.000, +0.000], weight: 0.002, band: 3, eig: 5.598, occ: 2.000 Final state: spin: 0, kpt: [+0.000, +0.000, +0.000], weight: 0.002, band: 4, eig: 8.130, occ: 0.000 Fundamental gap: Energy: 0.562 (eV) Initial state: spin: 0, kpt: [+0.000, +0.000, +0.000], weight: 0.002, band: 3, eig: 5.598, occ: 2.000 Final state: spin: 0, kpt: [+0.375, +0.375, +0.000], weight: 0.012, band: 4, eig: 6.161, occ: 0.000 Bandwidth: 11.856 (eV) Valence maximum located at: spin: 0, kpt: [+0.000, +0.000, +0.000], weight: 0.002, band: 3, eig: 5.598, occ: 2.000 Conduction minimum located at: spin: 0, kpt: [+0.375, +0.375, +0.000], weight: 0.012, band: 4, eig: 6.161, occ: 0.000 TIP: Call set_fermie_to_vbm() to set the Fermi level to the VBM if this is a non-magnetic semiconductor
The different contributions to the total energy are stored in a dictionary:
print(gsr.energy_terms)
Term Value e_localpsp -68.62575673515487 eV e_eigenvalues 4.453251818689574 eV e_ewald -226.94266726068815 eV e_hartree 14.989583992487596 eV e_corepsp 2.262480265343313 eV e_corepspdc 0.0 eV e_kinetic 80.66035184382825 eV e_nonlocalpsp 52.15353088290843 eV e_entropy 0.0 eV entropy 0.0 eV e_xc -95.73399330200242 eV e_xcdc 0.0 eV e_paw 0.0 eV e_pawdc 0.0 eV e_elecfield 0.0 eV e_magfield 0.0 eV e_fermie 5.5984533251030255 eV e_sicdc 0.0 eV e_exactX 0.0 eV h0 0.0 eV e_electronpositron 0.0 eV edc_electronpositron 0.0 eV e0_electronpositron 0.0 eV e_monopole 0.0 eV
At this point, we don't need this file anymore so we close it with:
gsr.close()
Note that we don't always follow this rule inside the jupyter notebook to maintain the code readable but you should definitively close all your files, especially when writing code that may be running for hours or even more.
Let's open the GSR file produced by a NSCF calculation done on a high-symmetry k-path and extract the electronic band structure.
A warning is issued by pymatgen about the structure not being standard. Be aware that this might possibly affect the automatic labelling of the boundary k-points on the k-path. So, check carefully the k-point labels on the figures that are produced in such case. In the present case, the labelling is correct.
with abilab.abiopen(abidata.ref_file("si_nscf_GSR.nc")) as nscf_gsr:
ebands_kpath = nscf_gsr.ebands
Now we can plot the band energies with:
# The labels for the k-points are found in an internal database.
ebands_kpath.plot(with_gaps=True, title="Silicon band structure");
Alternatively you can use the optional argument klabels
to define the mapping reduced_coordinates --> name of the k-point
and pass it to the plot method
klabels = {
(0.5, 0.0, 0.0): "L",
(0.0, 0.0, 0.0): "$\Gamma$",
(0.0, 0.5, 0.5): "X"
}
# ebands_kpath.plot(title="User-defined k-labels", band_range=(0, 5), klabels=klabels);
ebands_kpath.plotly(with_gaps=True, title="Silicon band structure with plotly");