In this example, we will calculate the phonon dispersion relations of Cu₂O at the DFT-PBE0/TZVP level of theory.

Table of contents

Prerequisites

To calculate the phonon dispersion relations of Cu₂O at the DFT-PBE0/TZVP level of theory, we need the optimized crystal structure of Cu₂O at the same level of theory. For details, see the page Geometry optimization of Cu₂O.

Phonopy reads in the optimized geometry from a CRYSTAL output file. There are several possible choices for the output file:

1. Output from IR and Raman spectra calculation (recommended). The output includes dielectric tensor and Born effective charges, which are required to calculate the non-analytical term correction to phonon dispersion relations. If you don't need the Raman spectrum, you can save time by leaving out keyword INTRAMAN.
2. Output from harmonic frequency calculation at the Γ-point (recommended only if you are sure that you will not need the non-analytical term correction). 
3. Output from geometry optimization (not recommended). It is always a good practice to check that the structure is a true local minimum based on the Γ-point frequencies before carrying out much more expensive phonon dispersion relation calculations.

Copy the CRYSTAL output file to a new directory and rename the file to crystal.o. This is the default filename that Phonopy expects for a CRYSTAL calculation.

Phonon supercell dimensions

In the directory with the file crystal.o, create a file called DIM that contains the dimensions of the phonon supercell. Here, we use a rather small 2x2x2 phonon supercell and the contents of the DIM file would be the following line:

2 2 2

TEMPLATE file for phonon supercell calculations

By default, Phonopy creates just the phonon supercell geometries and you need to manually complete the input files. If you provide a file called TEMPLATE in the directorywhere you are running the phonon calculations, Phonopy will prepare complete input files automatically. The file must contain the basis set block for CRYSTAL input and the computational parameters block without final END. Here is an example for Cu₂O:

8 7
0 0 6 2.0 1.0
  27032.382631      0.21726302465E-03
  4052.3871392      0.16838662199E-02
  922.32722710      0.87395616265E-02
  261.24070989      0.35239968808E-01
  85.354641351      0.11153519115
  31.035035245      0.25588953961
0 0 2 2.0 1.0
  12.271113873      0.39768730901
  4.9159842006      0.24627849430
0 0 1 0.0 1.0
 0.90086482370      1.0000000000
0 1 1 0.0 1.0
 0.25000000000      1.0 1.0
0 2 4 4.0 1.0
  75.300554155      0.60685103418E-02
  17.743733858      0.41912575824E-01
  5.5355828651      0.16153841088
  2.0685535103      0.35706951311
0 2 1 0.0 1.0
 0.78238772422      1.0000000000
0 3 1 0.0 1.0
 1.2000000000       1.0000000000
29 12
0 0 8 2.0 1.0
 377518.79923      0.22811766128E-03
 56589.984311      0.17688035931E-02
 12878.711706      0.91993460227E-02
 3645.3752143      0.37411016434E-01
 1187.0072945      0.12189873737
 426.46421902      0.28983900714
 165.70660164      0.41531872174
 65.598942707      0.21905799287
0 0 4 2.0 1.0
  414.41265811     -0.24682525053E-01
  128.32056039     -0.11716827406
  20.622089750      0.55301315941
  8.7821226045      0.52242718609
0 0 2 2.0 1.0
 13.741372006      -0.22736061821
 2.2431246833      0.71761210873
0 0 1 1.0 1.0
 0.89370549079     1.0000000000
0 0 1 0.0 1.0
  0.35              1.0000000000
0 1 1 0.0 1.0
  0.14              1.0 1.0
0 2 6 6.0 1.0
 2034.7596692       0.23524822298E-02
 481.90468106      0.19134070751E-01
 154.67482963      0.90171105278E-01
 57.740576969      0.26063284735
 23.099052811      0.42093485770
 9.3882478591      0.24344615121
0 2 3 6.0 1.0
 37.596171210      -0.28991094530E-01
 5.1240690810      0.54919083831
 2.0119996085      0.93793330488
0 2 1 0.0 1.0
 0.73860686002     1.0000000000
0 3 4 10.0 1.0
 74.129460637      0.14363216676E-01
 21.359842587      0.86628177096E-01
 7.4995240537      0.25631430541
 2.7601394169      0.40374062368
0 3 1 0.0 1.0
 0.95362061236     0.39427042447
0 3 1 0.0 1.0
 0.28420862520     0.23091146816
99 0
ENDBAS
DFT
PBE0
END
SHRINK
0 4
4 4 4
TOLINTEG
8 8 8 8 16
MAXCYCLE
120
FMIXING
80
EXCHSIZE
30000000
BIPOSIZE
30000000

Creating phonon supercells

Phonopy is operated from the command line with the command phonopy and appropriate flags. In the directory with the files crystal.o, DIM, and TEMPLATE, execute:

phonopy -v --crystal --dim="$(<DIM)" -d

Here we tell phonopy to invoke CRYSTAL mode (--crystal), tell it to show extended output (-v, stands for verbose), create the necessary displacements of atoms (-d, displacements) and do all this in a supercell that is of size 2x2x2 the primitive cell of Cu₂O. The choice of 2x2x2 is chosen to be computationally effective, see Useful tips 3.

The output printed in the terminal (see below) will show the lattice parameters and atom positions of the supercell as well as say "disp.yaml and supercells have been created." The folder now contains the input files for all the supercells that need to be calculated.

Figure 1. Output of Phonopy for Cu₂O when creating 2x2x2 sized supercells. Figure: Jarno Linnera.

Completing the phonon supercell input files

By default, the prepared .d12-files contain only the keywords "GRADCAL", "EXTERNAL" and "TOLDEE 10", telling CRYSTAL to perform force calculations on a geometry that is specified in an .ext-file and to use a strict SCF energy convergence threshold. Information about the basis set and DFT parameters should then be added to the CRYSTAL input files. This can be automated, as the interface also looks for a file called TEMPLATE in the folder. By having the basis set and DFT parameters in the file TEMPLATE, they are added in all the supercell-***.d12-files which are then ready to be calculated right off the bat. The TEMPLATE file used in this example can be found here.

The supercell calculation outputs at the DFT-PBE0/TZVP level are provided in the Optional files above if you do not wish to calculate them yourself. Please note that if you use the output files provided here, they may differ from the ones that would be produced by the TEMPLATE file given here. The example outputs were produced and calculated with older versions of the CRYSTAL-Phonopy interface and CRYSTAL.

Submitting the phonon supercell input files

The phonon supercells are typically fairly large and may require several gigabytes of memory per CPU. Here we use 3 GB per core:

jsub -np 12 -smp -mem 3G crycrystal supercell-001.d12

jsub -np 12 -smp -mem 3G crycrystal supercell-002.d12

Processing the phonon supercell outputs

When all the force calculations have succesfully finished, we feed them to Phonopy with the command

phonopy -v --crystal -f supercell-*o

Phonopy reads the CRYSTAL output files and creates the file FORCE_SETS.

[jarno@wihuri wiki_test]$ phonopy -v --crystal -f supercell-*o
        _
  _ __ | |__   ___  _ __   ___   _ __  _   _
 | '_ \| '_ \ / _ \| '_ \ / _ \ | '_ \| | | |
 | |_) | | | | (_) | | | | (_) || |_) | |_| |
 | .__/|_| |_|\___/|_| |_|\___(_) .__/ \__, |
 |_|                            |_|    |___/
                                      1.13.2

1. Drift force of "supercell-001.o" to be subtracted
  0.00000000   0.00000000   0.00000000
2. Drift force of "supercell-002.o" to be subtracted
  0.00000000   0.00000000  -0.00000000
FORCE_SETS has been created.
                 _
   ___ _ __   __| |
  / _ \ '_ \ / _` |
 |  __/ | | | (_| |
  \___|_| |_|\__,_|

Calculating force constants

You are now ready to calculate and plot the phonon dispersion. It is, however, a good idea to calculate and write the force constants to a file, so they can be read from the file whenever needed instead of computing again. This is accomplished with the command

 phonopy -v --crystal --dim="2 2 2" --writefc

The output should contain the structure and the following block in the end.

----------------------------------------------------------------------------
Computing force constants...
Max drift of force constants: 0.001669 (zz) -0.000000 (zx)
Force constants are written into FORCE_CONSTANTS.
----------------------------------------------------------------------------
 One of the following run modes may be specified for phonon calculations.
 - Mesh sampling (MP, --mesh)
 - Q-points (QPOINTS, --qpoints)
 - Band structure (BAND, --band)
 - Animation (ANIME, --anime)
 - Modulation (MODULATION, --modulation)
 - Characters of Irreps (IRREPS, --irreps)
 - Create displacements (-d)
----------------------------------------------------------------------------
                 _
   ___ _ __   __| |
  / _ \ '_ \ / _` |
 |  __/ | | | (_| |
  \___|_| |_|\__,_|

Phonopy also prints the maximum drift of force constants for given directions. It can be described as a measure of accuracy for the calculated force constants. If the drift is large, the results may not be reliable. The drift shown here for Cu₂O is very small and warrants no action. If the system shows large drift, that may be overcome with symmetrization of the force constants. For more information, see the Phonopy manual.

Calculating phonon dispersion relations

As Phonopy suggests, we may use one of the listed run modes. We will be using the flag --band for the phonon dispersion (phonon band structure). We also need to provide Phonopy with information on the reciprocal space paths along which the phonon dispersion is calculated and plotted. This can be given in the command when running the program, but it is much more convenient to write the information in a file, such as band.conf used here. An example for the contents of band.conf for Cu₂O is given below and a matching file is found here. The nomenclature of the band points is based on the paper of Setyawan and Curtarolo. 

DIM = 2 2 2
BAND = 0 0 0    0 1/2 0    1/2 1/2 0    0 0 0    1/2 1/2 1/2    0 1/2 0    1/2 1/2 0    1/2 1/2 1/2
BAND_LABELS = $\Gamma$ X M $\Gamma$ R X M R
BAND_POINTS = 101

Here we specify the supercell dimensions with DIM, the path in the reciprocal space with BAND, names of the points where band segments meet with BAND_LABELS and the number of points sampled between each band segment with BAND_POINTS. Now the command

phonopy -v --crystal --readfc -p -s --factor=521.47083 band.conf

will tell Phonopy to calculate and plot the data (-p), save it into a PDF (-s), use a factor of 521.47083 for converting energy units to cm^-1 (default is THz), read the force constants from the file rather than re-calculate them (--readfc) and finally to read all the other keywords from band.conf. Phonopy again prints out the structure and the path in the reciprocal space, as shown below.

Max drift of force constants: 0.001669 (xx) -0.000000 (zx)
Reciprocal space paths in reduced coordinates:
[ 0.00  0.00  0.00] --> [ 0.00  0.50  0.00]
[ 0.00  0.50  0.00] --> [ 0.50  0.50  0.00]
[ 0.50  0.50  0.00] --> [-0.00 -0.00  0.00]
[ 0.00  0.00  0.00] --> [ 0.50  0.50  0.50]
[ 0.50  0.50  0.50] --> [-0.00  0.50 -0.00]
[ 0.00  0.50  0.00] --> [ 0.50  0.50  0.00]
[ 0.50  0.50  0.00] --> [ 0.50  0.50  0.50]
                 _
   ___ _ __   __| |
  / _ \ '_ \ / _` |
 |  __/ | | | (_| |
  \___|_| |_|\__,_|

Plotting tools and useful tips

Now the folder should contain the files band.yaml, which contains the information about the phonon dispersion, and the plot in band.pdf. The .pdf should look like this.

Figure 2. Phonon dispersion of Cu₂O. Figure: Jarno Linnera.

If you wish to modify the plot, you can write your own scripts to read the data from the files or use the plotting script bandplot. For usage and functionalities, enter

 phonopy-bandplot -h

Useful tips:

1. SCF convergence issues with displaced supercells

Even though the displacements are small, sometimes the supercells can suffer from SCF convergence issues. This can be, in most cases, salvaged by giving CRYSTAL a better initial guess for the wavefunction. A good candidate for it is obtained by doing a single-point calculation for the regular supercell, where no atoms are displaced (remove the keyword GRADCAL, forces are not needed). Copying the wavefunction file supercell.w as supercell-***.w for the needed cell(s) and inserting the keyword GUESSP in the SCF-section of supercell-***.d12 will result in CRYSTAL using the pre-calculated wavefunction as the initial guess.

2. Reducing the CPU time of the force calculations

Following the previous advice, one can use the wavefunction from the undisplaced supercell as an initial guess for all of the displaced supercells. Whether or not this is useful, depends on the situation. If the number of cells is large, shaving some cycles from the SCF procedure can result in a decent reduction of CPU time.

3. Ensuring convergence with respect to interatomic distances

The supercell size of 2x2x2 has been tested to produce the correct dispersion for Cu₂O although the lattice constant is only roughly 8.6 Å. Generally, in order to ensure the force constants between a pair of atoms vanish within the supercell, distances of 10 - 15 Å are recommended. This is to make sure all interactions of relevant magnitude are taken into account. To make sure the results are converged, up to a 4x4x4 supercell would be preferable to check for Cu₂O.

Summary of the commands

Expect something to go wrong if you don't know what you're doing

1. phonopy -v --crystal --dim="2 2 2" -d

2. Supercell calculations

3. phonopy -v --crystal -f supercell-*o

4. phonopy -v --crystal --dim="2 2 2" -p -s band.conf