Back to the main Index

This lesson discusses how to compute phonon band structures, DOS and Born effective charges with Abinit and AbiPy. The discussion closely follows the second lesson on DFPT available on the Abinit web site. More specifically, we will discuss how to

- Perform a convergence study for the phonon frequencies at $\Gamma$ as function of
`ecut`

- Compute the full phonon band structure of
`AlAs`

with the inclusion of LO-TO splitting - Obtain thermodynamic properties within the harmonic approximation

We assume that you have read the references mentioned in the first Abinit lesson on DFPT. You might find additional material, related to the present section, in the following references:

- Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory
- Phonons and related crystal properties from density-functional perturbation theory

If you are already familiar with python and AbiPy-Abinit are already installed and configured, you may want to use directly the command line interface. See the README.md file in the directory of this lesson explaining how to analyze the data from the shell using ipython and matplotlib.

In [1]:

```
# Use this at the beginning of your script so that your code will be compatible with python3
from __future__ import print_function, division, unicode_literals
import numpy as np
import warnings
warnings.filterwarnings("ignore") # to get rid of deprecation warnings
from abipy import abilab
abilab.enable_notebook() # This line tells AbiPy we are running inside a notebook
import abipy.flowtk as flowtk
# 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
```

and an useful function from the `lesson_dfpt`

module that will be used to generate our DFPT flows:

In [2]:

```
from lesson_dfpt import make_scf_input
abilab.print_source(make_scf_input)
```

Out[2]:

```
def make_scf_input(ecut=2, ngkpt=(4, 4, 4)):
"""
This function constructs an `AbinitInput` to perform a GS-SCF calculation in crystalline AlAs.
Args:
ecut: cutoff energy in Ha.
ngkpt: 3 integers specifying the k-mesh for the electrons.
Return:
`AbinitInput` object
"""
# Initialize the AlAs structure from an internal database. Use the pseudos shipped with AbiPy.
gs_inp = abilab.AbinitInput(structure=abidata.structure_from_ucell("AlAs"),
pseudos=abidata.pseudos("13al.981214.fhi", "33as.pspnc"))
# Set the value of the Abinit variables needed for GS runs.
gs_inp.set_vars(
nband=4,
ecut=ecut,
ngkpt=ngkpt,
nshiftk=4,
shiftk=[0.0, 0.0, 0.5, # This gives the usual fcc Monkhorst-Pack grid
0.0, 0.5, 0.0,
0.5, 0.0, 0.0,
0.5, 0.5, 0.5],
ixc=1,
nstep=25,
diemac=9.0,
tolvrs=1.0e-10,
#iomode=3,
)
return gs_inp
```

In [3]:

```
scf_input = make_scf_input()
scf_input
```

Out[3]:

############################################################################################

# SECTION: basic

############################################################################################

nband 4

ecut 2

ngkpt 4 4 4

nshiftk 4

shiftk

0.0 0.0 0.5

0.0 0.5 0.0

0.5 0.0 0.0

0.5 0.5 0.5

ixc 1

nstep 25

tolvrs 1e-10

############################################################################################

# SECTION: gstate

############################################################################################

diemac 9.0

############################################################################################

# STRUCTURE

############################################################################################

natom 2

ntypat 2

typat 1 2

znucl 13 33

xred

0.0000000000 0.0000000000 0.0000000000

0.2500000000 0.2500000000 0.2500000000

acell 1.0 1.0 1.0

rprim

0.0000000000 5.3050000000 5.3050000000

5.3050000000 0.0000000000 5.3050000000

5.3050000000 5.3050000000 0.0000000000

# SECTION: basic

############################################################################################

nband 4

ecut 2

ngkpt 4 4 4

nshiftk 4

shiftk

0.0 0.0 0.5

0.0 0.5 0.0

0.5 0.0 0.0

0.5 0.5 0.5

ixc 1

nstep 25

tolvrs 1e-10

############################################################################################

# SECTION: gstate

############################################################################################

diemac 9.0

############################################################################################

# STRUCTURE

############################################################################################

natom 2

ntypat 2

typat 1 2

znucl 13 33

xred

0.0000000000 0.0000000000 0.0000000000

0.2500000000 0.2500000000 0.2500000000

acell 1.0 1.0 1.0

rprim

0.0000000000 5.3050000000 5.3050000000

5.3050000000 0.0000000000 5.3050000000

5.3050000000 5.3050000000 0.0000000000

In [4]:

```
print(scf_input.structure)
```

In [5]:

```
scf_input.structure.plot();
```

`ixc`

in the input file.

In [6]:

```
for pseudo in scf_input.pseudos:
print(pseudo, "\n")
```

As you can see, we essentialy have a standard input to perform a GS calculation. This object will represent
the **building block** for our DFPT calculation with AbiPy.

It this is not your first time you use the DFPT part of Abinit, you already know that phonon calculations
require an initial GS run to produce the `WFK`

file
followed by a DFPT run that reads the `WFK`

file and solves the Sternheimer equations for $N_{\text{irred}}(q)$
atomic perturbations where $N_{\text{irred}}$ is the number of independent atomic displacements (assuming $q$ belongs to the k-mesh).

If you try to do a convergence study wrt `ecut`

**without multi-datasets**, you will likely start from an initial GS input file with a given value of `ecut`

, use it as a template to generate the DFPT input files, create symbolic
links to the `WFK`

file produced in the first GS step and then instruct Abinit to read this file with `irdwfk`

.
Once you have a set of input files that work for a particular `ecut`

, one can simply replicate the set of
directories and files and use a script to change the value of `ecut`

in the input files.
Then, of course, one has to run the calculations manually, collect the results and produce nice plots to understand
what is happening.

This approach is obviously boring and error-prone if you are a human being but it is easy to implement in an algorithm
and machines do not complain if they have a lot of repetive work to do!
There are also several **technical advantages** in using this **task-based approach vs multi-datasets** but we discuss this point in more details afterwards.

If the machine could speak, it will tell you: give me an object that represents an input for GS calculations,
give me the list of q-points you want to compute as well as the parameters that must be changed in the initial input
and I will generate a `Flow`

for DFPT calculations.
This logic appears so frequenty that we decided to encapsulate it in the `flowtk.phonon_conv_flow`

factory function:

In [7]:

```
from lesson_dfpt import build_flow_alas_ecut_conv
abilab.print_source(build_flow_alas_ecut_conv)
```

Out[7]:

```
def build_flow_alas_ecut_conv(options):
"""
Build and return Flow for convergence study of phonon frequencies at Gamma as function of ecut.
"""
scf_input = make_scf_input()
return flowtk.phonon_conv_flow("flow_alas_ecut_conv", scf_input, qpoints=(0, 0, 0),
params=["ecut", [4, 6, 8]])
```

Let's call the function to build our flow:

In [8]:

```
flow = build_flow_alas_ecut_conv(options=None)
flow.show_info()
```

and call the `get_graphviz`

method to visualize the connection among the `Tasks`

:

In [9]:

```
flow.get_graphviz()
```

Out[9]:

In [10]:

```
# matplotlib version based on networkx
#flow.plot_networkx(with_edge_labels=True);
```

`ecut`

specified in `params`

.

In [11]:

```
for work in flow:
for task in work:
print(task.pos_str, "uses ecut:", task.input["ecut"])
```

In [12]:

```
flow.get_vars_dataframe("ecut")
```

Out[12]:

ecut | class | |
---|---|---|

w0_t0 | 4 | ScfTask |

w1_t0 | 4 | PhononTask |

w1_t1 | 4 | PhononTask |

w2_t0 | 6 | ScfTask |

w3_t0 | 6 | PhononTask |

w3_t1 | 6 | PhononTask |

w4_t0 | 8 | ScfTask |

w5_t0 | 8 | PhononTask |

w5_t1 | 8 | PhononTask |

Each group represents a `Workflow`

and consists of one `ScfTask`

(red circle) that solves the `KS`

equations self-consistently producing a `WFK`

file that will be used by the two children (`PhononTasks`

- blue circles)
to compute the first-order change of the wavefunctions due to one of the *irreducible* atomic pertubations.

Note that `phonon_conv_flow`

invokes Abinit under the hood to get the list of irreducible perturbations
and uses this information to build the flow.
This explains why we have two `PhononTasks`

per $q$-point instead of the total number of phonon modes that
equals $3*N_{atom}=6$.

Perhaps a table with the values of the input variables associated to the DFPT perturbation will help.
`None`

means that the variable is not defined in that particular input.

In [13]:

```
flow.get_vars_dataframe("rfphon", "rfatpol", "rfdir", "qpt", "kptopt")
```

Out[13]:

rfphon | rfatpol | rfdir | qpt | kptopt | class | |
---|---|---|---|---|---|---|

w0_t0 | None | None | None | None | None | ScfTask |

w1_t0 | 1 | [1, 1] | [1, 0, 0] | [0.0, 0.0, 0.0] | 2 | PhononTask |

w1_t1 | 1 | [2, 2] | [1, 0, 0] | [0.0, 0.0, 0.0] | 2 | PhononTask |

w2_t0 | None | None | None | None | None | ScfTask |

w3_t0 | 1 | [1, 1] | [1, 0, 0] | [0.0, 0.0, 0.0] | 2 | PhononTask |

w3_t1 | 1 | [2, 2] | [1, 0, 0] | [0.0, 0.0, 0.0] | 2 | PhononTask |

w4_t0 | None | None | None | None | None | ScfTask |

w5_t0 | 1 | [1, 1] | [1, 0, 0] | [0.0, 0.0, 0.0] | 2 | PhononTask |

w5_t1 | 1 | [2, 2] | [1, 0, 0] | [0.0, 0.0, 0.0] | 2 | PhononTask |

In [14]:

```
abilab.docvar("rfatpol")
```

Out[14]:

[1, 1]

Control the range of atoms for which displacements will be considered in
phonon calculations (atomic polarizations), using the 2n+1 theorem.
These values are only relevant to phonon response function calculations.
May take values from 1 to **natom**, with **rfatpol**(1)<=**rfatpol**(2).
The atoms to be moved will be defined by the
do-loop variable iatpol:

- do iatpol=
**rfatpol**(1),**rfatpol**(2)

For the calculation of a full dynamical matrix, use **rfatpol**(1)=1 and
**rfatpol**(2)=**natom**, together with **rfdir** 1 1 1. For selected
elements of the dynamical matrix, use different values of **rfatpol** and/or
**rfdir**. The name 'iatpol' is used for the part of the internal variable
ipert when it runs from 1 to **natom**. The internal variable ipert can also
assume values larger than **natom**, denoting perturbations of electric field
or stress type (see [help:respfn|the response function help file]).

Now we can generate the `flow_alas_ecut`

directory with the input files by executing
the `lesson_dfpt.py`

script.
Then use the `abirun.py`

script to launch the entire calculation with:

```
abirun.py flow_alas_ecut_conv scheduler
```

You will see that all `PhononTasks`

will be executed in parallel on your machine.

Please make sure that AbiPy is properly configured by running abicheck --with flow

If you prefer to skip this part, you may want to jump to the next section, that presents the post-processing of the results. Note that the output files are already available in the repository so it is also possible to try the AbiPy post-processing tools without having to run the flow.

There are several output files located inside the `outdata`

directories:

In [15]:

```
!find flow_alas_ecut_conv/ -name "*_DDB"
```

Remember that our goal is to analyze the convergence of the phonon frequencies at $\Gamma$
as function of `ecut`

.
So we are mainly interested in the DDB files located in the `outdata`

directories
of the `PhononWorks`

(`w0/outdata`

, `w1/outdata`

, `w2/outdata`

).
These are indeed the DDB files with all the information needed to reconstruct the
dynamical matrix at $\Gamma$ and to compute the phonon frequencies (AbiPy calls `mrgddb`

to merge the DDB files when all the perturbations in the `PhononWork`

have been computed).

The code below tells our robot that we would like to analyze all the DDB files located in the output directories of the works:

In [16]:

```
robot = abilab.DdbRobot.from_dir_glob("./flow_alas_ecut_conv/w*/outdata/")
robot
```

Out[16]:

- flow_alas_ecut_conv/w1/outdata/out_DDB
- flow_alas_ecut_conv/w3/outdata/out_DDB
- flow_alas_ecut_conv/w5/outdata/out_DDB

`anaddb`

to compute the phonon frequencies at $\Gamma$ for all DDBs
and return a pandas `DataFrame`

:

In [17]:

```
data_gamma = robot.get_dataframe_at_qpoint((0, 0, 0))
```

The `DataFrame`

is a dict-like object whose keys are the name of the colums in the table

In [18]:

```
print(data_gamma.keys())
```

where `mode-i`

is the frequency in eV of the i-th phonon mode.

We are mainly interested in the convergence of the phonon frequencies versus `ecut`

so we filter these columns with:

In [19]:

```
data_gamma = data_gamma[["ecut"] + [k for k in data_gamma if k.startswith("mode")]]
data_gamma
```

Out[19]:

ecut | mode0 | mode1 | mode2 | mode3 | mode4 | mode5 | |
---|---|---|---|---|---|---|---|

flow_alas_ecut_conv/w1/outdata/out_DDB | 4.0 | 0.0 | 0.0 | 0.0 | 0.043641 | 0.043641 | 0.043641 |

flow_alas_ecut_conv/w3/outdata/out_DDB | 6.0 | 0.0 | 0.0 | 0.0 | 0.044485 | 0.044485 | 0.044485 |

flow_alas_ecut_conv/w5/outdata/out_DDB | 8.0 | 0.0 | 0.0 | 0.0 | 0.044625 | 0.044625 | 0.044625 |

and we get some statistics about our data with:

In [20]:

```
data_gamma.describe()
```

Out[20]:

ecut | mode0 | mode1 | mode2 | mode3 | mode4 | mode5 | |
---|---|---|---|---|---|---|---|

count | 3.0 | 3.0 | 3.0 | 3.0 | 3.000000 | 3.000000 | 3.000000 |

mean | 6.0 | 0.0 | 0.0 | 0.0 | 0.044250 | 0.044250 | 0.044250 |

std | 2.0 | 0.0 | 0.0 | 0.0 | 0.000533 | 0.000533 | 0.000533 |

min | 4.0 | 0.0 | 0.0 | 0.0 | 0.043641 | 0.043641 | 0.043641 |

25% | 5.0 | 0.0 | 0.0 | 0.0 | 0.044063 | 0.044063 | 0.044063 |

50% | 6.0 | 0.0 | 0.0 | 0.0 | 0.044485 | 0.044485 | 0.044485 |

75% | 7.0 | 0.0 | 0.0 | 0.0 | 0.044555 | 0.044555 | 0.044555 |

max | 8.0 | 0.0 | 0.0 | 0.0 | 0.044625 | 0.044625 | 0.044625 |

`describe`

method already gives some useful info
about the convergence of the phonon modes.
Sometimes, however, we would like to visualize the data to have a better understanding of what's happening:

In [21]:

```
data_gamma.plot(x="ecut", y="mode3", style="-o");
```

Let's plot all the modes in different subplots with:

In [22]:

```
data_gamma.plot(x="ecut", y=[k for k in data_gamma if k.startswith("mode")], subplots=True, style="-o");
```

`ecut`

>= 6 Ha to get reasonably converged phonon frequencies at $\Gamma$.
In what follows, we assume that also the modes at the other $q$-points present a similar
convergence behaviour and we use `ecut`

= 6 Ha to keep the computational cost low.

For a quick introduction to Pandas, see:

Now we are finally ready for the calculation of the vibrational spectrum of $AlAs$.
We already managed to run DFPT calculations at $\Gamma$ with different values of `ecut`

and the
steps required to get a full band structure are not that different, provided that
the following differences are taken into account:

we need the dynamical matrix $D(q)$ on a homogeneous mesh so that it is possible to calculate $D(R)$ in anaddb via Fourier transform and then phonon frequencies for arbitrary q-points via Fourier interpolation

$AlAs$ is a polar semiconductor so we need to include the LO-TO splitting for $q \rightarrow 0$ that, in turns, requires the DFPT computation of the Born effective charges and of the dielectric constant.

In AbiPy, these concepts are translated in an easy-to-use API in which you pass an initial `AbinitInput`

object,
you specify the q-mesh for phonons in terms of `ph_nqpt`

and activate the computation of the
Born effective charges with the boolean flag `with_becs`

.

Let's have a look at the code (as usual there are more comments than lines of code):

In [23]:

```
from lesson_dfpt import build_flow_alas_phonons
abilab.print_source(build_flow_alas_phonons)
```

Out[23]:

```
def build_flow_alas_phonons(options):
"""
Build and return a Flow to compute the dynamical matrix on a (2, 2, 2) qmesh
as well as DDK and Born effective charges.
The final DDB with all perturbations will be merged automatically and placed
in the Flow `outdir` directory.
"""
scf_input = make_scf_input(ecut=6, ngkpt=(4, 4, 4))
return flowtk.PhononFlow.from_scf_input("flow_alas_phonons", scf_input,
ph_ngqpt=(2, 2, 2), with_becs=True)
```

We can finally construct the flow with:

In [24]:

```
flow_phbands = build_flow_alas_phonons(options=None)
```

and visualize the connections with:

In [25]:

```
flow_phbands.get_graphviz()
#flow_phbands.plot_networkx();
```

Out[25]:

Note that there are a lot of things happening under the hood here.

First of all, AbiPy generates `PhononTasks`

only for the $q$-points in the
irreducible wedge of the Brillouin zone corresponding to `ph_ngqpt`

.
Moreover, for a given $q$-point, only the irreducible atomic perturbations are explicitly computed
since the other atomic perturbations can be reconstructed by symmetry.
Fortunately you do not have to care about all these technical details as AbiPy and Abinit
will automate the whole procedure.

Remember that the $q$-point mesh cannot be chosen arbitrarily since all $q$ wavevectors should connect two $k$ points of the grid used for the electrons.

It is also worth stressing that the computational cost of the DFPT run depends on the q-point since only those symmetries that preserve the q-point as well as the direction of the perturbation can be employed (calculations at $\Gamma$ are therefore much faster than other q-points).

In [26]:

```
flow_phbands.get_vars_dataframe("rfphon", "rfatpol", "rfdir", "qpt", "kptopt")
```

Out[26]:

rfphon | rfatpol | rfdir | qpt | kptopt | class | |
---|---|---|---|---|---|---|

w0_t0 | None | None | None | None | None | ScfTask |

w1_t0 | None | None | (1, 0, 0) | (0, 0, 0) | 2 | DdkTask |

w1_t1 | None | None | (0, 1, 0) | (0, 0, 0) | 2 | DdkTask |

w1_t2 | None | None | (0, 0, 1) | (0, 0, 0) | 2 | DdkTask |

w1_t3 | 1 | [1, 1] | [1, 0, 0] | (0, 0, 0) | 2 | BecTask |

w1_t4 | 1 | [2, 2] | [1, 0, 0] | (0, 0, 0) | 2 | BecTask |

w2_t0 | 1 | [1, 1] | [1, 0, 0] | [0.5, 0.0, 0.0] | 3 | PhononTask |

w2_t1 | 1 | [1, 1] | [0, 1, 0] | [0.5, 0.0, 0.0] | 3 | PhononTask |

w2_t2 | 1 | [2, 2] | [1, 0, 0] | [0.5, 0.0, 0.0] | 3 | PhononTask |

w2_t3 | 1 | [2, 2] | [0, 1, 0] | [0.5, 0.0, 0.0] | 3 | PhononTask |

w3_t0 | 1 | [1, 1] | [1, 0, 0] | [0.5, 0.5, 0.0] | 3 | PhononTask |

w3_t1 | 1 | [2, 2] | [1, 0, 0] | [0.5, 0.5, 0.0] | 3 | PhononTask |

Now we can generate the directories and the input files of the `Flow`

.
Change lesson_dfpt.py so that the build_flow_alas_phonons function is called in main
instead of build_flow_alas_ecut_conv.
Run the script to generate the directory with the flow.
Finally, use

```
abirun.py flow_alas_phonons scheduler
```

to launch the entire calculation.

Our flow is completed and we have the final DDB file with all the $q$-points and all the independent atomic perturbations. Let's open this DDB file with:

In [27]:

```
ddb = abilab.abiopen("flow_alas_phonons/outdata/out_DDB")
print(ddb)
```

The `DdbFile`

object provides an easy-to-use interface that invokes `anaddb`

to post-process
the data stored in the DDB file.

`anacompare_phdos`

, for example, computes the phonon DOS with different $q$-meshes.
Each mesh is defined by a single integer, `nqsmall`

, that gives the number of
divisions used to sample the smallest reciprocal lattice vector.
The number of divisions along the other directions are chosen so that proportions are preserved:

In [28]:

```
c = ddb.anacompare_phdos(nqsmalls=[8, 10, 12, 14, 16])
```

In [29]:

```
c.plotter.combiplot();
```

A 16x16x16 $q$-mesh with the tethraedron method gives a well converged phonon DOS.

`anaget_phbst_and_phdos_files`

allows one to compute the phonon band structure on an automatically defined $q$-path as well as the the phonon DOS:

In [30]:

```
phbst_file, phdos_file = ddb.anaget_phbst_and_phdos_files(ndivsm=10, nqsmall=16, lo_to_splitting=False)
# Extract the phonon bands and the phonon DOS from phbst_file and phdos_file
phbands = phbst_file.phbands
phdos = phdos_file.phdos
```

Let's plot the bands with matplotlib:

In [31]:

```
phbands.plot();
```

and the high-symmetry q-path:

In [32]:

```
phbands.qpoints.plot();
```

In [33]:

```
phbst_file, phdos_file = ddb.anaget_phbst_and_phdos_files(ndivsm=10, nqsmall=16, lo_to_splitting=True)
# Extract the phonon bands and the phonon DOS from phbst_file and phdos_file
phbands = phbst_file.phbands
phdos = phdos_file.phdos
```

In [34]:

```
#phbst_file.phbands.qpoints.plot();
```

The band structure plot now correctly shows the non-analytical behaviour around $\Gamma$:

In [35]:

```
phbands.plot();
```

To plot bands and DOS on the same figure:

In [36]:

```
phbands.plot_with_phdos(phdos);
```

`PhdosFile`

contains the phonon frequencies, the displacement vectors
as well as the decomposition of the total DOS in terms of the contributions due to
the different types of atom in the unit cell.
This means that one can plot the type-projected phonon DOS with:

In [37]:

```
phdos_file.plot_pjdos_type(title="AlAs type-projected phonon DOS");
```

and it is even possible to plot fatbands and type-projected DOSes on the same figure with:

In [38]:

```
phbands.plot_fatbands(phdos_file=phdos_file);
```