Back to the main Index

Post-processing DFPT calculations with the DDB file

This notebook explains how to use AbiPy and the DDB file produced by Abinit to analyze:

  • Phonon band structures including the LO-TO splitting in heteropolar semiconductors
  • Phonon fatbands, phonon DOS and projected DOS
  • Born effectives charges $Z^*_{\kappa,\alpha\beta}$ and the dielectric tensors $\epsilon^{\infty}_{\alpha\beta}$, $\epsilon^{0}_{\alpha\beta}$
  • Thermodynamic properties in the harmonic approximation

In the last part, we discuss how to use the DdbRobot to analyze multiple DDB files and perform typical convergence studies.

Table of Contents

Suggested references

How to create a DdbFile object

[back to top]

Let us start by importing the basic AbiPy modules we have already used in the other examples:

In [1]:
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.data as abidata

# This line configures matplotlib to show figures embedded in the notebook, 
# instead of popping up a new window. 
%matplotlib notebook

To open a DDB file, use the high-level interface provided by abiopen:

In [2]:
ddb_filepath = abidata.ref_file("mp-1009129-9x9x10q_ebecs_DDB")
ddb = abilab.abiopen(ddb_filepath)

# To get the DDB from the materials project database, use:
#ddb = abilab.DdbFile.from_mpid("mp-1009129")

A DdbFile has a structure object:

In [3]:
print(ddb.structure)  # Lengths are in Angstrom.
Full Formula (Mg1 O1)
Reduced Formula: MgO
abc   :   2.908638   2.908638   2.656848
angles:  90.000000  90.000000 120.000000
Sites (2)
  #  SP           a         b    c
---  ----  --------  --------  ---
  0  Mg    0         0         0
  1  O     0.333333  0.666667  0.5

Abinit Spacegroup: spgid: 0, num_spatial_symmetries: 12, has_timerev: True, symmorphic: False

and a list of q-points associated to the dynamical matrix $D(q)$:

In [4]:
ddb.qpoints
Out[4]:
0) [+0.000, +0.000, +0.000]
1) [+0.111, +0.000, +0.000]
2) [+0.222, +0.000, +0.000]
3) [+0.333, +0.000, +0.000]
4) [+0.444, +0.000, +0.000]
5) [+0.111, +0.111, +0.000]
6) [+0.222, +0.111, +0.000]
7) [+0.333, +0.111, +0.000]
8) [+0.444, +0.111, +0.000]
9) [+0.222, +0.222, +0.000]
10) [+0.333, +0.222, +0.000]
11) [+0.333, +0.333, +0.000]
12) [+0.000, +0.000, +0.100]
13) [+0.111, +0.000, +0.100]
14) [+0.222, +0.000, +0.100]
15) [+0.333, +0.000, +0.100]
16) [+0.444, +0.000, +0.100]
17) [+0.111, +0.111, +0.100]
18) [+0.222, +0.111, +0.100]
19) [+0.333, +0.111, +0.100]
20) [+0.444, +0.111, +0.100]
21) [+0.222, +0.222, +0.100]
22) [+0.000, +0.000, +0.200]
23) [+0.333, +0.222, +0.100]
24) [+0.333, +0.333, +0.100]
25) [+0.111, +0.000, +0.200]
26) [+0.222, +0.000, +0.200]
27) [+0.333, +0.000, +0.200]
28) [+0.444, +0.000, +0.200]
29) [+0.111, +0.111, +0.200]
30) [+0.222, +0.111, +0.200]
31) [+0.333, +0.111, +0.200]
32) [+0.444, +0.111, +0.200]
33) [+0.222, +0.222, +0.200]
34) [+0.000, +0.000, +0.300]
35) [+0.333, +0.222, +0.200]
36) [+0.333, +0.333, +0.200]
37) [+0.111, +0.000, +0.300]
38) [+0.222, +0.000, +0.300]
39) [+0.333, +0.000, +0.300]
40) [+0.444, +0.000, +0.300]
41) [+0.111, +0.111, +0.300]
42) [+0.222, +0.111, +0.300]
43) [+0.333, +0.111, +0.300]
44) [+0.444, +0.111, +0.300]
45) [+0.222, +0.222, +0.300]
46) [+0.333, +0.222, +0.300]
47) [+0.000, +0.000, +0.400]
48) [+0.333, +0.333, +0.300]
49) [+0.111, +0.000, +0.400]
50) [+0.222, +0.000, +0.400]
51) [+0.333, +0.000, +0.400]
52) [+0.444, +0.000, +0.400]
53) [+0.111, +0.111, +0.400]
54) [+0.222, +0.111, +0.400]
55) [+0.333, +0.111, +0.400]
56) [+0.444, +0.111, +0.400]
57) [+0.222, +0.222, +0.400]
58) [+0.333, +0.222, +0.400]
59) [+0.000, +0.000, +0.500]
60) [+0.333, +0.333, +0.400]
61) [+0.111, +0.000, +0.500]
62) [+0.222, +0.000, +0.500]
63) [+0.333, +0.000, +0.500]
64) [+0.444, +0.000, +0.500]
65) [+0.111, +0.111, +0.500]
66) [+0.222, +0.111, +0.500]
67) [+0.333, +0.111, +0.500]
68) [+0.444, +0.111, +0.500]
69) [+0.222, +0.222, +0.500]
70) [+0.333, +0.222, +0.500]
71) [+0.333, +0.333, +0.500]

At this point, it is worth mentioning that Abinit takes advantage of symmetries to reduce the number of q-points as well as the number of perturbations that must be computed explicitly within DFPT.

The set of q-points in the DDB file (usually) does not form a homogeneous sampling of the Brillouin zone (BZ). Actually they correspond to the sampling of the irreducible wedge (IBZ), and this sampling is obtained from an initial q-mesh specified in terms of divisions along the three reduced directions (ngqpt).

In [5]:
ddb.qpoints.plot();

Note that the DDB file does not contain any information about the value of ngqpt because one can merge an arbitrary list of q-points in the same DDB. The algorithms implemented in anaddb, however, need to know the divisions to compute integrals in the full BZ (this is indeed one of the variables that must be provided by the user in the anaddb input file).

AbiPy uses a heuristic method to guess the q-mesh from this scattered list of q-points so that you don not need to specify this parameter when calling anaddb:

In [6]:
ddb.guessed_ngqpt
Out[6]:
array([ 9,  9, 10])

To test whether the DDB file contains the entries associated to $\epsilon^{\infty}_{\alpha\beta}$ and $Z^*_{\kappa,\alpha\beta}$, use:

In [7]:
print("Contains macroscopic dielectric tensor:", ddb.has_emacro_terms())
print("Contains Born effective charges:", ddb.has_bec_terms())
Contains macroscopic dielectric tensor: True
Contains Born effective charges: True

Metadata are stored in the header of the DDB file. AbiPy parses this initial section and stores the values in a dict-like object. Let us print a sorted list with all the keys in the header:

In [8]:
from pprint import pprint
pprint(sorted(list(ddb.header.keys())))
['acell',
 'amu',
 'dfpt_sciss',
 'dilatmx',
 'ecut',
 'ecutsm',
 'intxc',
 'iscf',
 'ixc',
 'kpt',
 'kptnrm',
 'lines',
 'natom',
 'nband',
 'ngfft',
 'nkpt',
 'nspden',
 'nspinor',
 'nsppol',
 'nsym',
 'ntypat',
 'occ',
 'occopt',
 'rprim',
 'spinat',
 'symafm',
 'symrel',
 'tnons',
 'tolwfr',
 'tphysel',
 'tsmear',
 'typat',
 'usepaw',
 'version',
 'wtk',
 'xred',
 'zion',
 'znucl']

and use the standard syntax for dictionaries to access the keys:

In [9]:
print("This DDB has been generated with ecut", ddb.header["ecut"], "Ha and nsym:", ddb.header["nsym"])
This DDB has been generated with ecut 44.0 Ha and nsym: 12

We can also print the DDB object to get a summary of the most important parameters and dimensions:

In [10]:
print(ddb)   

# If more info is needed use: 
# print(ddb.to_string(verbose=1)
================================= File Info =================================
Name: mp-1009129-9x9x10q_ebecs_DDB
Directory: /Users/gmatteo/git_repos/abipy/abipy/data
Size: 218.73 kb
Access Time: Wed Feb  7 12:14:51 2018
Modification Time: Thu Jan 18 15:17:04 2018
Change Time: Thu Jan 18 15:17:04 2018

================================= Structure =================================
Full Formula (Mg1 O1)
Reduced Formula: MgO
abc   :   2.908638   2.908638   2.656848
angles:  90.000000  90.000000 120.000000
Sites (2)
  #  SP           a         b    c
---  ----  --------  --------  ---
  0  Mg    0         0         0
  1  O     0.333333  0.666667  0.5

Abinit Spacegroup: spgid: 0, num_spatial_symmetries: 12, has_timerev: True, symmorphic: False

================================== DDB Info ==================================

Number of q-points in DDB: 72
guessed_ngqpt: [ 9  9 10] (guess for the q-mesh divisions made by AbiPy)
Has electric-field perturbation: True
Has Born effective charges: True

If you are a terminal aficionado, remember that one can use the abiopen.py script to open the DDB file directly from the shell and generate a jupyter notebook with the -nb option. For a quick visualization script use abiview.py.

Invoking Anaddb from the DdbFile object

[back to top]

The DdbFile object provides specialized methods to invoke anaddb and compute important physical properties such as the phonon band structure, the phonon density of states etc. All these methods have a name that begins with the ana* prefix followed by a verb (anaget, anacompare). These specialized methods

  • build the anaddb input file
  • run anaddb
  • parse the netcdf files produced by the Fortran code
  • build and return AbiPy objects that can be used to plot/analyze the data.

Note that in order to run anaddb from AbiPy, you need a manager.yml with configuration options. For further details, please consult the TaskManager documentation.

The python API is flexible and exposes several anaddb input variables. The majority of the arguments have default values covering the most common cases so that you will need to specify these arguments explicitly only if the default behavior does not suit your needs. The most important parameters to remember are:

  • ndivsm: Number of divisions used for the smallest segment of the high-symmetry q-path
  • nqsmall: Defines the q-mesh for the phonon DOS in terms of the number of divisions used to sample the smallest reciprocal lattice vector. 0 to disable DOS computation
  • lo_to_splitting: Activate the computation of the frequencies in the $q\rightarrow 0$ with the inclusion of the non-analytical term (requires dipdip 1 and DDB with $Z^*_{\kappa,\alpha\beta}$ and $\epsilon^{\infty}_{\alpha\beta}$).

The high-symmetry q-path is automatically selected assuming the structure fulfills the conventions introduced by Setyawan and Curtarolo but you can also specify your own q-path if needed.

Plotting phonon bands and DOS

[back to top]

To compute phonon bands and DOS, use:

In [11]:
# Call anaddb to compute phonon bands and DOS. Return PHBST and PHDOS netcdf files.
phbstnc, phdosnc = ddb.anaget_phbst_and_phdos_files(
    ndivsm=20, nqsmall=20, lo_to_splitting=True, asr=2, chneut=1, dipdip=1, dos_method="tetra")

# Extract phbands and phdos from the netcdf object.
phbands = phbstnc.phbands
phdos = phdosnc.phdos

Let us have a look at the high symmetry q-path automatically selected by AbiPy with:

In [12]:
phbands.qpoints.plot();

and plot the phonon bands along this path with:

In [13]:
phbands.plot();

Note the discontinuity of the optical modes when we cross the $\Gamma$ point. In heteropolar semiconductors, indeed, the dynamical matrix is non-analytical for $q \rightarrow 0$. Since lo_to_splitting was set to True, AbiPy has activated the calculation of the phonon frequencies for all the $q \rightarrow \Gamma$ directions present in the path.

There are several band crossings and anti-crossings hence it's not easy to understand how the branches should be connected. Fortunately, there is a heuristic method to estimate band connection from the overlap of the eigenvectors at adjacent q-points. To connect the modes and plot the phonon branches with different colors, use:

In [14]:
phbands.plot_colored_matched();

To plot the DOS, $g(\omega)$, and the integrated $IDOS(\omega) = \int^{\omega}_{-\infty} g(\omega')\,d\omega'$, use:

In [15]:
phdos.plot();

Note how the phonon DOS integrates to $3 * N_{atom} = 6$

To plot the phonon bands and the DOS on the same figure use:

In [16]:
phbands.plot_with_phdos(phdos, units="meV");

Fatbands and projected DOS

[back to top]

The phbands object stores the phonon displacements, $\vec{d}_{q\nu}$ and the eigenvectors, $\vec{\epsilon}_{q\nu}$ obtained by diagonalizing the dynamical matrix $D(q)$.

\begin{equation} D(q) \vec{\epsilon}_{q\nu} = \omega_{q\nu}^2 \vec{\epsilon}_{q\nu} \end{equation}

We can therefore use the eigenvectors (or the displacements) to associate a width to the different bands (a.k.a. fatbands). This width gives us a qualitative understanding of the vibrational mode: what are the atomic types involved in the vibrations at a given energy, their direction of oscillation and the amplitude (related to the displacement).

In [17]:
# NB: LO-TO is not available in fatbands
phbands.plot_fatbands(use_eigvec=True, units="Thz");

To plot the fatbands with the type-projected DOS stored in phdocsnc use:

In [18]:
phbands.plot_fatbands(phdos_file=phdosnc, colormap="rainbow", alpha=0.4, units="cm-1");