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

Integration with the materials project database

AbiPy, pymatgen and fireworks have been used by Petretto et al to compute the vibrational properties of more than 1500 compounds with Abinit. The results are available on the materials project website. The results for the rocksalt phase of MgO are available at

To fetch the DDB file from the materials project database and build a DdbFile object, use:

ddb = abilab.DdbFile.from_mpid("mp-1009129")

Remember to set the PMG_MAPI_KEY in your ~/.pmgrc.yaml as described here.

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 [58]:
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 as abidata
import os

# 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
#%matplotlib widget  

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)

A DdbFile has a structure object:

In [3]:
print(ddb.structure)  # Lengths 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]:
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]:

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 do not need to specify this parameter when calling anaddb:

In [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 [9]:
print("Contains macroscopic dielectric tensor:", ddb.has_epsinf_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 [10]:
from pprint import pprint

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

In [11]:
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 [12]:

# 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: Tue Oct 15 00:49:59 2019
Modification Time: Wed Mar 20 16:53:35 2019
Change Time: Wed Mar 20 16:53:35 2019

================================= 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)
ecut = 44.000000, ecutsm = 0.000000, nkpt = 405, nsym = 12, usepaw = 0
nsppol 1, nspinor 1, nspden 1, ixc = -116133, occopt = 1, tsmear = 0.010000

Has total energy: False, Has forces: False
Has stress tensor: False

Has (at least one) atomic pertubation: True
Has (at least one diagonal) electric-field perturbation: True
Has (at least one) Born effective charge: True
Has (all) strain terms: False
Has (all) internal strain terms: False
Has (all) piezoelectric terms: False

If you are a terminal aficionado, remember that one can use the 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

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$ limit 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 [13]:
# 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 [14]:

and plot the phonon bands along this path with:

In [15]:

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 [16]:

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

In [17]:

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 [18]:
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 [19]:
# 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 [20]:
phbands.plot_fatbands(phdos_file=phdosnc, colormap="rainbow", alpha=0.4, units="cm-1");

We can also plot the PJDOS summed over directions and atomic types without fatbands with the command:

In [21]:

The netcdf file contains the individual contributions to the total DOS for each atomic site and the three cartesian directions. So there are several quantities we can plot to understand the vibrational spectrum of our systems. For example, we can decide to sum over all atoms of the same type while keeping the dependence on the Cartesian direction:

In [22]:

This analysis tells us that the peak at ~16 Thz mainly consists of oxygen vibrations along z. We could now extend this analysis by looking at the contributions arising from the different sites with:

In [23]:
#phdosnc.plot_pjdos_cartdirs_site(view="inequivalent", units="eV", stacked=True);

but we prefer to stop here and discuss other tools that can be used to analyze individual phonon modes.

Visualizing atomic displacements

[back to top]

In you need to visualize the lattice vibrations in 3D to gain a better insight about the nature of the phonon modes you may want to use the phononwebsite. One can either convert the produced by anaddb into json format and upload it to the phononwebsite server or, alternatively, open the terminal and execute the AbiPy script: ddb out_DDB --phononwebsite

to automate the entire process (replace out_DDB with the name of your DDB file).

Note that there are other AbiPy methods that are quite handy if we need to investigate the nature of the phonon modes at a particular q-point without a 3D visualization tool. For example, it is possible to analyze the contribution given by the different types of atoms to the phonon displacements at a given q-point with:

In [24]:
phbands.plot_phdispl(qpoint=(0, 0, 0), tight_layout=True);

As expected the first three (acoustic) modes have zero frequency and the two atoms oscillate with the same amplitude. These modes indeed correspond to a rigid translation of the crystal hence the amplitude (and the direction) does not depend on the atomic site. The other three optical modes are (almost) degenerate, while LO-TO splitting should be present, but this is due to the fact that we are using the frequencies at the $\Gamma$ point without the inclusion of the non-analytical term.

In [25]:
#phbands.plot_phdispl((0, 0, 0.1), is_non_analytical_direction=True, tight_layout=True);

To project the phonon displacements along the three cartesian directions, use:

In [26]:
phbands.plot_phdispl_cartdirs(qpoint=(0, 0, 0.0), tight_layout=True);

This plot confirms that the first three modes correspond to a rigid translation along x, y and z, respectively.

Analyzing the breaking of the acoustic sum rule

[back to top]

Due to the invariance of the system under an infinitesimal rigid translation, the frequency of the lowest three modes at $\Gamma$ should be zero. Unfortunately, all the terms that are evaluated on the real-space FFT mesh (e.g. $V_{xc}$, non-linear core-correction) break this kind of translational invariance. The error depends on several factors: the density of the FFT mesh, pseudopotentials with hard model core charges, XC functional, etc.) Note that it is not always possible to reduce the error to zero by just increasing the convergence parameters but fortunately it is possible to restore the acoustic sum rule via the asr input variable.

One can easily compare the phonons bands obtained with different values of asr with:

In [27]:
asr_plotter = ddb.anacompare_asr()

This method invokes anaddb with different values of asr and returns a plotter object we can call to compare the phonon band structures:

In [28]:

Now we can perform a similar test for the treatment of the non-analytical term in the $q \rightarrow 0$ limit. We compute the phonon band dispersion for dipdip in [0, 1] and compare the results on the same figure with the commands:

In [29]:
dipdip_plotter = ddb.anacompare_dipdip(nqsmall=0)
In [30]:

The figure above shows that the (Fourier interpolated) bands obtained with dipdip = 0 have unphysical oscillations around the $\Gamma$ point. These oscillation are due to the long-range behavior in real space of the interatomic force constants in heteropolar semiconductors. The correct description of this long-range term without dipdip = 1 would require using an extremely dense q-point mesh in the DFPT calculation.

With dipdip = 1, on the other hand, we can model this long-range behavior in terms of a dipole-dipole interaction involving the Born effective charges and the macroscopic dielectric tensor. This allows us to decompose the full dynamical matrix into:

\begin{equation} D(q) = D^{sr}(q) + D^{dip-dip}(q) \end{equation}

The analytical part of the dynamical matrix, $D^{sr}(q)$, is short-ranged and can be Fourier-interpolated with a relatively coarse q-mesh. Then the model for the non-analytical term, $D^{dip-dip}(q)$, is added back to the interpolated matrix to get the full dynamical matrix. This procedure solves the problem with the unphysical oscillations and is required to describe correctly the non-analytical behavior of the optical modes for $q \rightarrow 0$.

Computing DOS with different q-meshes

[back to top]

Phonon DOS and derived quantites (e.g. thermodynamic properties) are sensitive to the BZ sampling and dense meshes may be required to converge the final results.

The method anacompare_phdos provides a simple interface to compare phonon DOS computed with different $q$-meshes. We only need to provide a list of integers (nqsmalls). Each integer defines the number of divisions used to sample the smallest reciprocal lattice vector while the other two vectors are sampled such that proportions are preserved.

To calculate three phonon DOS with increasing number of q-points use:

In [31]:
res = ddb.anacompare_phdos(nqsmalls=[8, 12, 24])

The return value is a named tuple with the phonon DOSes in res.phdoses while res.plotter is PhononDosPlotter. We can easily compare our results with:

In [32]:

Thermodynamic properties in the harmonic approximation

[back to top]

The thermodynamic properties of an ensemble of non-interacting phonons can be expressed in terms of integrals of the phonon DOS $g(\omega)$ using:

\begin{equation} %\label{eq:helmholtz} \Delta F = 3nNk_BT\int_{0}^{\omega_L}\text{ln}\left(2\text{sinh}\frac{\hbar\omega}{2k_BT}\right)g(\omega)d\omega \end{equation}\begin{equation} %\label{eq:free_en} \Delta E = 3nN\frac{\hbar}{2}\int_{0}^{\omega_L}\omega\text{coth}\left(\frac{\hbar\omega}{2k_BT}\right)g(\omega)d\omega \end{equation}\begin{equation} %\label{eq:c_v} C_v = 3nNk_B\int_{0}^{\omega_L}\left(\frac{\hbar\omega}{2k_BT}\right)^2\text{csch}^2\left(\frac{\hbar\omega}{2k_BT}\right)g(\omega)d\omega \end{equation}\begin{equation} %\label{eq:entropy} S = 3nNk_B\int_{0}^{\omega_L}\left(\frac{\hbar\omega}{2k_BT}\text{coth}\left(\frac{\hbar\omega}{2k_BT}\right) - \text{ln}\left(2\text{sinh}\frac{\hbar\omega}{2k_BT}\right)\right)g(\omega)d\omega, \end{equation}

where $k_B$ is the Boltzmann constant. This should represent a reasonable approximation especially in the low-temperature regime in which anharmonic effects can be neglected.

Let's plot the vibrational contributions to the thermodynamic properties as function of temperature $T$:

In [33]:

and the zero-point energy $\dfrac{1}{2 N_q} \sum_{q\nu} \omega_{q\nu}$

In [34]:
zpe = phdos.zero_point_energy
print("Zero-point energy", zpe,"Ha"))
Zero-point energy 0.13173567745342335 eV 0.00484119689285712 Ha

To get the free energy for a range of temperatures in Kelvin degrees, use:

In [35]:
f = phdos.get_free_energy(tstart=10, tstop=100)

Macroscopic dielectric tensor and Born effective charges

[back to top]

Let us call anaddb to compute the electronic contribution to the macroscopic dielectric tensor, $\epsilon^{\infty}_{\alpha\beta}$, and the Born effective charges $Z^*_{\kappa,\alpha\beta}$:

In [38]:
emacro, becs = ddb.anaget_epsinf_and_becs(chneut=1)

Note the use of chneut to enforce charge neutrality.

In [39]:
x y z
x 3.438497 0.000000 0.000000
y 0.000000 3.438497 0.000000
z 0.000000 0.000000 3.499416
In [40]:
element site_index frac_coords cart_coords wyckoff xx yy zz yz xz xy
0 Mg 0 [0.0, 0.0, 0.0] [0.0, 0.0, 0.0] 1a 1.98721 1.98721 2.10659 0.0 0.0 0.0
1 O 1 [0.33333, 0.66667, 0.5] [1.45432, 0.83965, 1.32842] 1d -1.98721 -1.98721 -2.10659 0.0 0.0 0.0

In the computation of the low-frequency (infrared) dielectric tensor, $\epsilon^{0}_{\alpha\beta}$, one has to include the response of the ions, whose motion will be triggered by the electric field. One can show that:

\begin{equation} %\label{eq:dielectric} \epsilon^{0}_{\alpha\beta}(\omega) = \epsilon^{\infty}_{\alpha\beta} + 4\pi\sum_m\frac{S_{m,\alpha\beta}}{(\omega^{\Gamma}_m - \omega)^2}. \end{equation}

where $\omega^{\Gamma}_m$ are the phonon frequencies at the center of the BZ and $S_{m,\alpha\beta}$ is the so-called mode-oscillator strengh tensor that depends on the phonon displacement and the Born effective charges.

To compute the dielectric tensor using the data stored in the DDB file, use:

In [59]:
dtgen = ddb.anaget_dielectric_tensor_generator()
In [51]:
# call print to get useful information
================================= 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

============================ Oscillator strength ============================
Real part in Cartesian coordinates. a.u. units; 1 a.u. = 253.2638413 m3/s2. Set to zero below 1.00e-08.
                          xx                      yy                     zz   yz   xz   xy
0                        0.0                     0.0                    0.0  0.0  0.0  0.0
1                        0.0                     0.0                    0.0  0.0  0.0  0.0
2                        0.0                     0.0                    0.0  0.0  0.0  0.0
3     0.00022453328281351966                     0.0                    0.0  0.0  0.0  0.0
4                        0.0  0.00022453328281351966                    0.0  0.0  0.0  0.0
5                        0.0                     0.0  0.0002523207013218457  0.0  0.0  0.0

============================= Dielectric Tensors =============================
Electronic dielectric tensor (eps_inf) in Cartesian coordinates. Set to zero below 1.00e-03.
          x         y         z
x  3.438497  0.000000  0.000000
y  0.000000  3.438497  0.000000
z  0.000000  0.000000  3.499416

Zero-frequency dielectric tensor (eps_zero) in Cartesian coordinates. Set to zero below 1.00e-03.
           x          y          z
x  13.101919   0.000000   0.000000
y   0.000000  13.101919   0.000000
z   0.000000   0.000000  13.779005

To compute the tensor for a given frequency:

In [45]:
e0 = dtgen.tensor_at_frequency(0)
[[ 1.31018619e+01+0.02381961j  2.54893071e-16+0.j
   1.47502725e-25+0.j        ]
 [ 0.00000000e+00+0.j          1.31018619e+01+0.02381961j
  -2.55482215e-25+0.j        ]
 [-3.68756813e-25+0.j         -1.27741107e-25+0.j

To plot the frequency dependence with a damping factor of 1e-4 eV (phonon linewidth):

In [53]:
dtgen.plot(w_max=None, gamma_ev=1e-4, component='diag', units='eV');

Using DdbRobot to perform convergence studies

[back to top]

A DdbRobot receives a list of DDB files and provides methods to construct pandas dataframes and analyze the results of multiple calculations. DdbRobots, in particular, are extremely useful to study the convergence of the phonon frequencies with respect to some computational parameters e.g. the number of k-points and the electronic smearing in metallic systems.

In this example, we are interested in the effect of the k-point sampling and the smearing parameter on the vibrational properties of magnesium diboride. $MgB_2$ is a metallic system with a critical temperature of 39 K that is the highest among conventional (phonon-mediated) superconductors. We use precomputed DDB files obtained by running GS+DFPT calculations with different values of nkpt and tsmear.

Let's build our DdbRobot object with:

In [54]:
import os
paths = [

paths = [os.path.join(abidata.dirpath, "refs", "mgb2_phonons_nkpt_tsmear", f) for f in paths]

robot = abilab.DdbRobot.from_files(paths)
  1. ../../abipy/abipy/data/refs/mgb2_phonons_nkpt_tsmear/mgb2_888k_0.01tsmear_DDB
  2. ../../abipy/abipy/data/refs/mgb2_phonons_nkpt_tsmear/mgb2_888k_0.02tsmear_DDB
  3. ../../abipy/abipy/data/refs/mgb2_phonons_nkpt_tsmear/mgb2_888k_0.04tsmear_DDB
  4. ../../abipy/abipy/data/refs/mgb2_phonons_nkpt_tsmear/mgb2_121212k_0.01tsmear_DDB
  5. ../../abipy/abipy/data/refs/mgb2_phonons_nkpt_tsmear/mgb2_121212k_0.02tsmear_DDB
  6. ../../abipy/abipy/data/refs/mgb2_phonons_nkpt_tsmear/mgb2_121212k_0.04tsmear_DDB

The script provides a command line interface to build robots from a list of files/directories given as arguments

The DDB files are now stored in the robot with a label constructed from the file path. These labels, however, are not very informative. In principle we would like to have a label that reflects the value of (nkpt, tsmear) also because these labels will be used to generate the labels in our plots.

Let's fix it with a function that recomputes the labels from the metadata available in ddb.header:

In [55]:
function = lambda ddb: "nkpt: %s, tsmear: %.2f" % (ddb.header["nkpt"], ddb.header["tsmear"])
  1. nkpt: 256, tsmear: 0.01
  2. nkpt: 256, tsmear: 0.02
  3. nkpt: 256, tsmear: 0.04
  4. nkpt: 864, tsmear: 0.01
  5. nkpt: 864, tsmear: 0.02
  6. nkpt: 864, tsmear: 0.04

We are usually interested in the convergence behavior with respect to one or two parameters of the calculations. Let's build a pandas dataframe with the most important parameters extracted from the DDB header:

In [56]:
nkpt nsppol ecut tsmear occopt ixc nband usepaw
nkpt: 256, tsmear: 0.01 256 1 35.0 0.01 4 1 8 0
nkpt: 256, tsmear: 0.02 256 1 35.0 0.02 4 1 8 0
nkpt: 256, tsmear: 0.04 256 1 35.0 0.04 4 1 8 0
nkpt: 864, tsmear: 0.01 864 1 35.0 0.01 4 1 8 0
nkpt: 864, tsmear: 0.02 864 1 35.0 0.02 4 1 8 0
nkpt: 864, tsmear: 0.04 864 1 35.0 0.04 4 1 8 0

Now we tell the robot to invoke anaddb to compute the phonon bands for all DDB files. Since we are not interested in the phonon DOS, nqsmall is set to 0

In [60]:
r = robot.anaget_phonon_plotters(nqsmall=0)

Now we can plot all the phonon band structures on the same figure with:

In [61]:

The plot is a bit crowded. Still, it is clear that there are portions of the vibration spectrum that are quite sensitive to the values of (nkpt, tsmear).

In metals, it's common to analyze the convergence of the physical properties by plotting the results as function of the k-point sampling for fixed value of tsmear. Let's do something similar for the phonon band structures with the command:

In [62]:
r.phbands_plotter.gridplot_with_hue("tsmear", units="Thz");