Phonon Properties of Diamond-Structure Crystals (DFPT)

by Sebastian Tillack for exciting fluorine

Purpose: In this tutorial, you will learn how to perform a phonon calculation using density-functional perturbation theory (DFPT). You will work through the example of diamond.

Note: We strongly recommend to do the super-cell phonon tutorial Phonon properties of diamond-structure crystals (Super-Cell) first to familiarize yourself with the standard phonon input parameters. All phonon-derived properties from the super-cell tutorial can be applied in the same way to DFPT phonons.

---

0. Prerequisites

This tutorial assumes that you have already set the relevant environment variables. Otherwise, please have a look at How to set environment variables for tutorials scripts. We further assume that you have already completed Phonon properties of diamond-structure crystals (Super-Cell).

In order to start, create a new working directory, e.g., /home/tutorials/diamond-phonons-dfpt.

$ cd /home/tutorials
$ mkdir diamond-phonons-dfpt
$ cd diamond-phonons-dfpt

---

1. Theoretical Background

---

☛ Theory of linear response for phonons

Within density-functional perturbation theory (DFPT) we calculate the linear response of the Kohn-Sham (KS) system upon an external perturbation $\delta$. We express the quantities of the perturbed system as an expansion in the perturbation operator and write down the KS equations for the perturbed system

(1)

\begin{align} (\hat{\bf h}\phantom{}^0 + \lambda \,\delta \hat{\bf h} + \dots) (\psi_{n{\bf k}}^0({\bf r}) + \lambda\, \delta \psi_{n{\bf k}}({\bf r}) + \dots) = (\epsilon_{n{\bf k}}^0 + \lambda \,\delta \epsilon_{n{\bf k}} + \dots) (\psi_{n{\bf k}}^0({\bf r}) + \lambda \,\delta \psi_{n{\bf k}}({\bf r}) + \dots) \;. \end{align}

Collecting terms linear in $\lambda$ we find the so called Sternheimer equation

(2)

\begin{align} [\hat{\bf h}\phantom{}^0 - \epsilon_{n{\bf k}}^0] \,\delta \psi_{n{\bf k}}({\bf r}) = -[\delta \hat{\bf h} - \delta \epsilon_{n{\bf k}}] \,\psi_{n{\bf k}}^0({\bf r}) \;, \end{align}

from which the linear response of the wavefunctions $\delta \psi_{n{\bf k}}$ can be obtained. Equivalently to the KS equations, the Sternheimer equation has to be solved self consistently, since the Hamiltonian response implicitly depends on the wave function response via the density and potential response:

(3)

\begin{align} \delta \psi_{n{\bf k}} \rightarrow \delta n({\bf r}) \rightarrow \delta V({\bf r}) \rightarrow \delta \hat{\bf h} \rightarrow \delta \psi_{n{\bf k}} \;. \end{align}

In the case of DFPT for phonons, the external perturbation is a collective phonon-like displacement of the nuclei from their equilibrium positions $\pmb{\tau}_{\kappa{\bf R}}$

(4)

\begin{align} \delta = \delta^{\bf q}_{\kappa \alpha} \quad:\quad \pmb{\tau}_{\kappa {\bf R}} \longrightarrow \pmb{\tau}_{\kappa {\bf R}} + {\rm e}^{{\rm i}{\bf q}\cdot{\bf R}} \delta \pmb{\tau}_{\kappa} \;, \end{align}

where $\kappa$ and ${\bf R}$ label atoms and unit-cells, respectively, ${\bf q}$ is the phonon wave vector and $\delta \pmb{\tau}_{\kappa}$ is an infinitesimal displacement vector.

In order to exploit symmetries, we do not displace every atom individually in the three Cartesian directions (canonical basis) but we use special symmetry adapted displacement patterns (irreducible representations, irreps):

(5)

\begin{align} \delta^{\bf q}_{Id} = \sum_{\kappa,\alpha} p_{\kappa \alpha, Id}({\bf q})\, \delta^{\bf q}_{\kappa \alpha} \;, \end{align}

where $I$ labels an irrep (for a given ${\bf q}$) and $d$ labels a so-called member of the irrep.

Once we have converged the Sternheimer equation for the density and potential response, we can use the results to compute in a single step the linear response of the atomic forces ${\bf F}_{\kappa}$ which directly correspond to the entries of the dynamical matrix

(6)

\begin{align} \delta^{\bf q}_{\kappa \alpha} F_{\kappa' \beta} = \sqrt{M_\kappa\, M_{\kappa'}}\, D_{\kappa \alpha, \kappa' \beta}({\bf q}) \;. \end{align}

---

2. Running a DFPT phonon calculation

In comparison to the input file from Phonon properties of diamond-structure crystals (Super-Cell) two essential changes in the input file are necessary.

First, in the groundstate element, we replace the attribute xctype by the new element libxc. By that, we specify to take the exchange-correlation potential from the libxc library instead of the exciting internal implementation. In the case of the PBE xc-functional, the corresponding line in the input file reads

<libxc correlation="XC_GGA_C_PBE" exchange="XC_GGA_X_PBE"/>

The reason for this change is that a DFPT calculation requires the xc-kernel (the second derivatives of the xc-energy) which is not natively implemented in exciting and is always obtained from the libxc library. For full consistency, we recommend to also use the xc-potential (first derivative of the xc-energy) from libxc. However, in principle you can also use the exciting internal implementations of the xc-potential using the attribute xctype. In this case, only the xc-kernel of the corresponding potential is obtained from libxc.

The second change is to switch the method used to compute phonons to "dfpt" to trigger a DFPT calculation.

...
   <phonons 
      do="fromscratch" 
      method="dfpt"
      ngridq="2 2 2">
      ...
   </phonons>
...

The complete input file for this calculation is given in the following.

<input>
 
   <title>Diamond: DFPT Phonons</title>
 
   <structure speciespath="$EXCITINGROOT/species">
 
      <crystal scale="6.7468">
         <basevect>0.5 0.5 0.0</basevect>
         <basevect>0.5 0.0 0.5</basevect>
         <basevect>0.0 0.5 0.5</basevect>
      </crystal>
 
      <species speciesfile="C.xml" rmt="1.25">
         <atom coord="0.00 0.00 0.00" />
         <atom coord="0.25 0.25 0.25" />
      </species>
 
   </structure>
 
   <groundstate 
      do="fromscratch" 
      ngridk="2 2 2" 
      swidth="0.0001"
      rgkmax="4.0"
      gmaxvr="14"
      maxscl="25"
      >
      <libxc correlation="XC_GGA_C_PBE" exchange="XC_GGA_X_PBE"/>
   </groundstate>
 
   <phonons 
      do="fromscratch" 
      method="dfpt"
      ngridq="2 2 2">
      <qpointset>
         <qpoint> 0.0 0.0 0.0 </qpoint>
         <qpoint> 0.5 0.5 0.0 </qpoint>
         <qpoint> 0.5 0.5 0.5 </qpoint>
      </qpointset>
   </phonons>
 
</input>

Run the calculation as usual with

$ SETUP-excitingroot.sh
$ time exciting_smp &

---

3. Output of a DFPT phonon calculation

3.1 Common output files

The most important output files such as PHONON.OUT and the dynamical matrices in DYN_Q####_####_####_S##_A###_P#.OUT are the same as in Phonon properties of diamond-structure crystals (Super-Cell). Hence, you still can see the resulting phonon frequencies using

$ CONVERT-au2invcm.py

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

--------------------------------------------
Total number of atoms    :     2
Total number of q-points :     3
--------------------------------------------

q-point   1  is  0.0 0.0 0.0
   mode   1      frequency:     0.0000  cm-1
   mode   2      frequency:     0.0000  cm-1
   mode   3      frequency:     0.0000  cm-1
   mode   4      frequency:  1537.2055  cm-1
   mode   5      frequency:  1537.2055  cm-1
   mode   6      frequency:  1537.2055  cm-1

q-point   2  is  0.5 0.5 0.0
   mode   1      frequency:   677.4866  cm-1
   mode   2      frequency:   677.4866  cm-1
   mode   3      frequency:  1005.5120  cm-1
   mode   4      frequency:  1005.5120  cm-1
   mode   5      frequency:  1202.3322  cm-1
   mode   6      frequency:  1202.3322  cm-1

q-point   3  is  0.5 0.5 0.5
   mode   1      frequency:   452.2903  cm-1
   mode   2      frequency:   452.2903  cm-1
   mode   3      frequency:   894.6191  cm-1
   mode   4      frequency:  1290.8184  cm-1
   mode   5      frequency:  1290.8184  cm-1
   mode   6      frequency:  1377.7608  cm-1

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

Note the differences in the resulting frequencies compared to the super-cell approach. These are due to the different methods and become much smaller for more converged numerical settings.

In order to get the phonon density of states and the full phonon dispersion curves, one can proceed as in tutorial Phonon properties of diamond-structure crystals (Super-Cell). The resulting phonon density-of-states will look like:

The obtained phonon dispersion curves are shown below.

---

☛ show converged plot

Better converged results (using an 8$\times$8$\times$8 k-grid, rgkmax = "7.0", and 4 q-points in each direction) can be found below.

---

3.2 DFPT specific output

A DFPT phonon calculation produces different subdirectories than a super-cell calculation. In a DFPT calculation we use special symmetry adapted displacement patterns, the irreps, instead of the canonical displacements of each atom in each of the three Cartesian directions as done in the super-cell calculation. Hence, the names of the subdirectories in a DFPT calculation are Q####_####_####_I###, where the first part specifies the q-point in the grid and the second part labels the irrep at that specific q-point. Inside of each subdirectory, you will find the following files:

filename description PHONON_INFO_fileext.OUT General information about the run. Similar to INFO.OUT for the groundstate. DYN_fileext_#.OUT Force response (or dynamical matrix row) for the #-th member of the irrep. PHONON_DRHO_fileext_#.OUT Electron density response for the #-th member of the irrep. PHONON_DVEFF_fileext_#.OUT Effective potential response for the #-th member of the irrep.

where fileext is the same as the subdirectory name.

In the following, we take a closer look to the file PHONON_INFO_Q0000_0000_0000_I001.OUT.

The first part of the file contains information about the use of symmetries, parallelization settings and the whole phonon run.

********************************************************************************
* symmetries and parallelization settings                                      *
********************************************************************************

Irreducible representations are used.
k-points are reduced by symmetries.
Response quantities are symmetrized.

q-points and symmetries:

  iq  q-vector (lattice coordinates)  N_sym  N_k  N_irrep (size)
--------------------------------------------------------------------------------
   1 ( 0.000000, 0.000000, 0.000000)     48    3        2 (3*,3)
   2 ( 0.500000, 0.000000, 0.000000)     12    4        4 (2,1,1,2)
   3 ( 0.500000, 0.500000, 0.000000)     16    4        3 (2,2,2)

current task / total (remaining) number of tasks:         1/       9 (9)
load of current task / total (remaining) load:            9/      66 (66)
procs for current task / total number of procs:           1/       1
average / maximum load per proc:                         66/      66
number of idle procs:                                     0

The first lines tell that the irreps have been used as displacement patterns and that the computational cost has been reduced by exploiting symmetries.

The following table displays all (irreducible) q-points that are computed in that phonon calculation. N_sym is the number of crystal symmetries in the small group of q. N_k indicates the number of k-points in the irreducible Brillouin zone obtained with the symmetries in the small group. N_irrep denotes the number of irreps for the given q-point and the values in parenthesis give their respective size (or equivalently the number of members per irrep). The q-point and irrep this file corresponds to is marked with an asterisk.

From the last lines we can extract that the whole calculation consists of 9 independent tasks from which all 9 are remaining and that this particular file corresponds to task 1. The meaning of the other lines will be discussed below in the Section 4.2.

In the next part, the displacement patterns of the irreps are given.

********************************************************************************
* irreducible representations / atomic displacement patterns                   *
********************************************************************************

irrep member 1
  C     1    +0.7071067812 +0.0000000000i
             +0.0000000000 +0.0000000000i
             +0.0000000000 +0.0000000000i
  C     2    -0.7071067812 +0.0000000000i
             +0.0000000000 +0.0000000000i
             +0.0000000000 +0.0000000000i

irrep member 2
  C     1    +0.0000000000 +0.0000000000i
             +0.7071067812 +0.0000000000i
             +0.0000000000 +0.0000000000i
  C     2    +0.0000000000 +0.0000000000i
             -0.7071067812 +0.0000000000i
             +0.0000000000 +0.0000000000i

irrep member 3
  C     1    +0.0000000000 +0.0000000000i
             +0.0000000000 +0.0000000000i
             +0.7071067812 +0.0000000000i
  C     2    +0.0000000000 +0.0000000000i
             +0.0000000000 +0.0000000000i
             -0.7071067812 +0.0000000000i

In this example, the irrep has three members and for each member, the three components of the (potentially complex) displacement vector for each atom are given.

The remainder of the file contains information about the self-consistency cycle of the Sternheimer equation and shows the change in the potential response for each iteration:

********************************************************************************
* scf loop for Sternheimer equation started                                    *
********************************************************************************
Initial density response read from file.

Sternheimer SCF loop iteration   1
maximum change in potential response:                0.101255    
Fermi energy response:                                0.00000    
Time spent for iteration (seconds):                          0.77

Sternheimer SCF loop iteration   2
maximum change in potential response:                0.697163E-01
Fermi energy response:                                0.00000    
Time spent for iteration (seconds):                          0.72

...

---

4. Useful tips for real-life calculations

4.1 Resume from a previous run

Similar to a ground-state calculation, also a DFPT phonon calculation can be resumed from a previous run. This might prove helpful in cases where an expensive calculation did not finish, e.g., due to given run-time limits or following an unintended crash. In order to restart a calculation from a previous run, set the attribute do to "fromfile" in the phonons element. This will read the density and potential response from the files PHONON_DRHO_fileext_#.OUT and PHONON_DVEFF_fileext_#.OUT and restart the iterative solution of the Sternheimer equation. For that to work, however, it is important that the respective files have been written in the first place. To achieve this, we recommend to make use of the attribute nwrite from the groundstate element in the case of more elaborate calculations. This attribute specifies after how many iterations the current density and potential response are written to file.

Note that subdirectories in which the DYN_fileext_#.OUT files are present will be skipped in following runs.

4.2 Parallel execution using multiple processes and MPI

---

☛ Parallelization explained

In general, DFPT phonon calculations can be highly parallelized. There are two levels of parallelization.

1. The outer level parallelizes over pairs (q,I) of q-points and irreps. Each such pair corresponds to one of the subdirectories and constitutes an independent part of the full calculation. The (q,I) pairs do not depend on one another and can be calculated in parallel.
2. The inner level parallelizes over k-points and irrep members d within each of the (q,I) pairs. MPI parallelization over this inner level is also implemented. The product of the number of k-points and irrep members determines the computational load of each of the (q,I) pairs. The computational load is directly proportional to the time needed for a single SCF iteration. It also is the upper limit for the number of processes that can be utilized.

Important information about the computational cost of the full calculation and the parallelization settings that are in use can be found in the file PHONON_RUN_INFO.OUT. Let's discuss the content of the file for the example above.

We already know the first part from the PHONON_INFO_fileext.OUT file.

================================================================================
= DFPT PHONON CALCULATION INFORMATION                                          =
================================================================================

Irreducible representations will be used.
k-points will be reduced by symmetries.
Response quantities will be symmetrized.

********************************************************************************
* q-points and displacements                                                   *
********************************************************************************

number of q-points:               3
number of phonon modes:           6

  iq  q-vector (lattice coordinates)  N_sym  N_k  N_irrep (size)
--------------------------------------------------------------------------------
   1 ( 0.000000, 0.000000, 0.000000)     48    3        2 (3,3)
   2 ( 0.500000, 0.000000, 0.000000)     12    4        4 (2,1,1,2)
   3 ( 0.500000, 0.500000, 0.000000)     16    4        3 (2,2,2)

It shows some information about the usage of symmetries and prints all the q-points and irreps with their size (number of irrep members) and the number of (irreducible) k-points for each q-point.

The next part displays information about the computational load of the calculation.

********************************************************************************
* independent parts and computational load                                     *
********************************************************************************

total (remaining) number of independent parts:      9 (9)
total (remaining) computational load:              66 (66)

remaining parts
  ip   iq irrep dim load procs
------------------------------
   1    1     1   3    9     1
   2    1     2   3    9     1
   3    2     1   2    8     1
   4    2     2   1    4     1
   5    2     3   1    4     1
   6    2     4   2    8     1
   7    3     1   2    8     1
   8    3     2   2    8     1
   9    3     3   2    8     1

done parts
  iq irrep dim load
-------------------

In this example, there are 9 independent parts (= number of (q,I) pairs = number of subdirectories). This number equals the sum of the column N_irrep in the table above. The total computational load of the full calculation is 66. I.e., 66 is the largest number of processes that can be effectively utilized in this example.

The two tables below list all the remaining and already completed (none in this example) parts. For each part ip, it gives the index of the q-point iq and the index of the irreducible representation irrep and its dimension dim (= number of the members of irrep). The last two columns give the computational load (= product of number of k-points and irrep members) and the number of processes that will work on this part.

The last part of the file gives further information about the parallelization settings.

********************************************************************************
* parallelization settings                                                     *
********************************************************************************

number of processes:                        1
optimal average load per process:        66.0
minimum load per process:                  66
maximum load per process:                  66
number of idle processes:                   0

proc load (imbalance)
----------------------
   0   66 (+0.0)

execution schedule
 proc|        5|       10|       15|       20|       25|       30|       35|       40|       45|       50|       55|       60|       65|
    0|###P4###:::P5:::#######P3#######:::::::P6:::::::#######P7#######:::::::P8:::::::#######P9#######::::::::P1::::::::########P2########

In this example, there will be only one process. Both the optimal, minimum and maximum computational load per process is 66 in this example and there are no idle processes. The table below shows the load assigned to each process. In parenthesis, the imbalance, i.e., the difference to the optimal average load is provided.

The last part is a "graphical" visualization of the execution schedule. It shows, which process will work on which part and in which order (from left to right). The scale above the schedule is in units of computational load.

---

Note that one unit of computational load is different for each calculation and it is mainly determined by the number of (L)APW+LO basis functions (= matrix size). Keep reading for a more meaningful example.

---

☛ Choosing parallelization settings

In order to find out how the calculation will be parallelized for a given number of processes you have the option of performing a dry run. Let's assume, we plan to execute the example above with 8 processes. Before we start the actual calculation, we want to find out how the calculation will be parallelized. To this extend, we change the attribute do to "dry" in the phonons element and add the new attribute drynumprocs = "8".

...
   <phonons 
      do="dry" 
      method="dfpt"
      ngridq="2 2 2"
      drynumprocs="8">
      ...
   </phonons>
...

Running exciting, this case, will only produce the file PHONON_RUN_INFO.OUT. (You might have to delete the DYN_fileext_#.OUT files if you have already completed the phonon calculation above in the same directory.) This time, the file looks like the following.

================================================================================
= DFPT PHONON CALCULATION INFORMATION                                          =
================================================================================

Irreducible representations will be used.
k-points will be reduced by symmetries.
Response quantities will be symmetrized.

********************************************************************************
* q-points and displacements                                                   *
********************************************************************************

number of q-points:               3
number of phonon modes:           6

  iq  q-vector (lattice coordinates)  N_sym  N_k  N_irrep (size)
--------------------------------------------------------------------------------
   1 ( 0.000000, 0.000000, 0.000000)     48    3        2 (3,3)
   2 ( 0.500000, 0.000000, 0.000000)     12    4        4 (2,1,1,2)
   3 ( 0.500000, 0.500000, 0.000000)     16    4        3 (2,2,2)

********************************************************************************
* independent parts and computational load                                     *
********************************************************************************

total (remaining) number of independent parts:      9 (9)
total (remaining) computational load:              66 (66)

remaining parts
  ip   iq irrep dim load procs
------------------------------
   1    1     1   3    9     5
   2    1     2   3    9     3
   3    2     1   2    8     8
   4    2     2   1    4     4
   5    2     3   1    4     4
   6    2     4   2    8     8
   7    3     1   2    8     8
   8    3     2   2    8     8
   9    3     3   2    8     8

done parts
  iq irrep dim load
-------------------

********************************************************************************
* parallelization settings                                                     *
********************************************************************************

number of processes:                        8
optimal average load per process:         8.2
minimum load per process:                   6
maximum load per process:                   9
number of idle processes:                   0

proc load (imbalance)
----------------------
   0    9 (+0.8)
   1    9 (+0.8)
   2    9 (+0.8)
   3    6 (-2.2)
   4    6 (-2.2)
   5    9 (+0.8)
   6    9 (+0.8)
   7    9 (+0.8)

execution schedule
 proc|        5|
    0|##::##::##::####  
    1|##::##::##P4####  
    2|##::##::##::#P1#  
    3|P3P6P7P8P9::####  
    4|##::##::##%%####  
    5|##::##::##P5::::::
    6|##::##::##%%::P2::
    7|##::##::##%%::::::

The table in the second section of the file now shows how many processes will be working on which part (e.g., 3 processes will work on part 2). From the third section we can tell that 8 is a reasonable number of processes for this example, since the optimal average load is 8.2 and the maximum load per process is 9 which is the smallest possible integer greater than the average load and meaning only very little "waste" of resources.

From the "graphical" visualization of the execution schedule, we read that all processes will start working together on part 3 (P3), followed by a joint execution of P6, P7, P8, and P9. After that, the processes split work and processes 0 to 3 will work on P4 while processes 4 to 7 will work on P5. After both P4 and P5 are finished, procs 0 to 4 will work on P1 while procs 5 to 7 will work on P2.

Note the prioritization of the inner level parallelization. All processes start to jointly work on a single part instead of individually starting many parts in parallel.

---