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

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.

      ngridq="2 2 2">

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

   <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>
      <species speciesfile="C.xml" rmt="1.25">
         <atom coord="0.00 0.00 0.00" />
         <atom coord="0.25 0.25 0.25" />
      ngridk="2 2 2" 
      <libxc correlation="XC_GGA_C_PBE" exchange="XC_GGA_X_PBE"/>
      ngridq="2 2 2">
         <qpoint> 0.0 0.0 0.0 </qpoint>
         <qpoint> 0.5 0.5 0.0 </qpoint>
         <qpoint> 0.5 0.5 0.5 </qpoint>

Run the calculation as usual with

$ 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



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.


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

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License