ALAMODE

Users Guide

About

What is ALAMODE?

ALAMODE is an open source software designed for analyzing lattice anharmonicity and lattice thermal conductivity of solids. By using an external DFT package such as VASP and Quantum ESPRESSO, you can extract harmonic and anharmonic force constants straightforwardly with ALAMODE. Using the calculated anharmonic force constants, you can also estimate lattice thermal conductivity, phonon linewidth, and other anharmonic phonon properties from first principles.

Features

General
  • Extraction of harmonic and anharmonic force constants based on the supercell approach
  • Applicable to any crystal structures and low-dimensional systems
  • Accurate treatment of translational and rotational invariance
  • Interface to VASP, Quantum-ESPRESSO, and xTAPP codes
  • Mainly written in C++, parallelized with MPI+OpenMP
Harmonic properties
  • Phonon dispersion
  • Phonon DOS, atom-projected phonon DOS
  • Two-phonon DOS
  • Vibrational thermodynamic functions (heat capacity, entropy, free energy)
  • Mean-square displacement
  • Animation and visualization of phonon modes (requires VMD or XCrysDen)
  • 3-phonon scattering phase space
  • Phonon-isotope scattering rate
  • Participation ratio for analyzing localization of phonon modes
Anharmonic properties
  • Grüneisen parameter via cubic force constants
  • Lattice thermal conductivity by BTE-RTA
  • Cumulative thermal conductivity
  • Phonon linewidth due to 3-phonon interactions
  • Phonon frequency shift due to 3- and 4-phonon interactions
  • Temperature-dependent effective potential method

License

Copyright © 2014, 2015, 2016 Terumasa Tadano

This software is distributed under the MIT license. See the LICENSE.txt file for license rights and limitations.

How to Cite ALAMODE

Please cite the following article when you use ALAMODE:

T. Tadano, Y. Gohda, and S. Tsuneyuki, J. Phys.: Condens. Matter 26, 225402 (2014) [Link].

Acknowledgement

This project was supported by a Grant-in-Aid for Scientific Research on Innovative Areas ‘Materials Design through Computics: Complex Correlation and Non-Equilibrium Dynamics’. (http://computics-material.jp)

Author & Contact

Department of Applied Physics, The University of Tokyo, Japan

If you have any questions, suggestions, and problems regarding ALAMODE, please feel free to contact the author.

Download

You can download the latest and previous versions of ALAMODE at http://sourceforge.net/projects/alamode .

You can also download the package from the git repository as

If you choose the GitHub version, please use the ‘master’ branch.

Installation

Requirement

Mandatory requirements
  • C++ compiler (Intel compiler is recommended.)
  • LAPACK library
  • MPI library (Either OpenMPI, MPICH2, or IntelMPI)
  • Boost C++ library

In addition to the above requirements, users have to get and install a first-principles package (such as VASP, Wien2k, QUANTUM-ESPRESSO, or xTAPP) or another force field package (such as LAMMPS) by themselves in order to compute harmonic and anharmonic force constants.

Optional requirements
  • Python (> 2.6), Numpy, and Matplotlib
  • XcrySDen or VMD

We provide some small scripts written in Python (Python 2) for visualizing phonon dispersion relations, phonon DOSs, etc. To use these scripts, one need to install the above Python packages. Additionally, XcrySDen is necessary to visualize the normal mode directions and animate the normal mode. VMD may be more useful to make an animation, but it may be replaced by any other visualization software which support the XYZ format.

How to install

  1. Install the LAPACK, MPI, and Boost C++ libraries.

    To install the Boost C++ library, please download a source file from the webpage and unpack the file. Then, copy the ‘boost’ subdirectory to the include folder in the home directory (or anywhere you like). This can be done as follows:

    $ cd
    $ mkdir etc; cd etc
    (Download a source file and mv it to ~/etc)
    $ tar xvf boost_x_yy_z.tar.bz2
    $ cd ../
    $ mkdir include; cd include
    $ ln -s ../etc/boost_x_yy_z/boost .
    

In this example, we made a symbolic link to the ‘boost’ subdirectory in $HOME/include.

Instead of install from source, you can install the Boost library with Homebrew on Mac OSX.

  1. Download the package of ALAMODE from the download page or clone from the git repository.

  2. Change directory to the location of the file and untar the file as follows:

    $ tar -xvzf alamode-x.y.z.tar.gz
    

This will create a directory alamode-x.y.z containing the following subdirectories:

  • alm/ : Source files for alm (force constant calculation)
  • anphon/ : Source files for anphon (phonon calculation)
  • external/ : Third-party include files
  • include/ : Commonly-used include files
  • tools/ : Small auxiliary programs and scripts
  • docs/ : Source files for making documents
  • example/ : Example files
  1. Edit the Makefiles

In directories alm/ and anphon/, we provide sample Makefiles for Linux (with Intel compiler) and Mac OSX (with gcc). Please copy one of them as Makefile and modify it appropriately for your environment.

Here’s a typical setting for Linux with Intel compiler:

CXX = icpc
CXXFLAGS = -O2 -xHOST -openmp
INCLUDE = -I../include -I$(HOME)/include

CXXL = ${CXX}
LDFLAGS = -mkl

LAPACK =
LIBS = ${LAPACK}

To enable OpenMP parallelization, please add the -openmp (Intel) or -fopenmp (gcc) option in CXXFLAGS. In addition, the directory containing the boost/ subdirectory must be given in INCLUDE.

  1. Make executables by make command.

    If the compilation is successful, the binary file named alm (anphon) is created in the alm/ (anphon/) directory. To use some auxiliary scripts for post-processing and data conversion, please compile the programs in the tools directory as well. See README.md in the tools directory for details about the auxiliary programs.

Running ALAMODE

Program alm

Program alm estimates harmonic and anharmonic interatomic force constants (IFCs) based on the supercell approach.

  1. Perform usual SCF calculations for a primitive celll
Before performing phonon calculations, one needs to perform usual self-consistent field calculations and check the convergence with respect to the cutoff energy and the \(k\) point density. After that, please optimize the internal coordinate so that the atomic forces are negligibly small. Optimization of cell parameters may also be necessary, but please note that phonon properties are relatively sensitive to the cell parameters in polar materials such as perovskites.
  1. Decide the size of supercell
Next, please decide the size of a supercell. Here, one may use a conventional cell. When the primitive cell is fairly large (\(a \sim 10\) Å), one may proceed using the primitive cell.
  1. Prepare an input file for alm

Please make an input file for alm, say alm.in, which should contain &general, &interaction, &cutoff, and &position entries. For details of available input variables, please refer to here. Once the input file is properly prepared with MODE = suggest, necessary displacement patterns can be generated by executing alm as follows:

$ alm alm.in > alm.log

This produces the following files containing the pattern of atomic displacements.

  • PREFIX.HARMONIC_pattern
  • PREFIX.ANHARM?_pattern (If NORDER > 1)

In pattern files, all necessary displacement patterns are given in Cartesian coordinates.

Note

Pattern files just indicate the direction of displacements. The magnitude of displacements should be specified by each user. (\(\Delta u \sim 0.01\) Å is usual for calculating harmonic force constants.)

  1. Perform SCF calculations to generate displacement-force data set

Then, please prepare necessary input files for a DFT engine (or a classical force field engine) and calculate atomic forces for each displaced configuration. Once the atomic forces are calculated for all configurations, please collect the atomic displacements and atomic forces to separate files, say disp_all.dat and force_all.dat, in Rydberg atomic units. The detail of the file format is described in this page.

Note

We provide some auxiliary Python scripts to expedite the above procedure for VASP, Quantum ESPRESSO, and xTAPP users. The script files can be found in the tools/ directory. We are willing to support other software if necessary.

  1. Estimate IFCs by a least-squares fitting

In order to perform a fitting, please change the variable MODE of the input file alm.in to MODE = fitting. In addition please add the &fitting entry with appropriate NDATA, DFILE, and FFILE. (DFILE should be like DFILE = disp_all.dat.) Then, IFCs can be estimated by executing

$ alm alm.in > alm.log

which makes the following two files in the working directory.

  • PREFIX.fcs : The list of force constants
  • PREFIX.xml : XML file containing necessary information for subsequent phonon calculations

Program anphon

  1. Prepare an input file for anphon

To perform phonon calculations and thermal conductivity calculations, one needs to prepare another input file, say anphon.in, for the program anphon.

If one wants to perform (harmonic) phonon calculations, one should write MODE = phonons in the &general entry of anphon.in. Please make sure that FCSXML variable being set to the XML file generated by alm.

If one wants to conduct thermal conductivity calculations instead of usual phonon calculations, please switch to MODE = RTA with appropriate FCSXML containing cubic IFCs.

For details of input variables of anphon, please refer to the list of input variables for anphon.

  1. Execute anphon

Phonon properties and lattice thermal conductivity can be calculated via executing

$ anphon anphon.in > anphon.log

or

$ mpirun -np NPROCS anphon anphon.in > anphon.log

Here, NPROCS is the number of MPI threads. If the code is compiled with the OpenMP option, the OpenMP parallelization can also be used by setting the OMP_NUM_THREADS variable as

$ export OMP_NUM_THREADS=16

The number 16 should be modified appropriately for your environment.

Note

MPI+OpenMP hybrid parallelization is available when calculating thermal conductivity with MODE = RTA, in which anharmonic self-energies of all \(N_{\boldsymbol{q},irred}\times N_{j}\) phonon modes need to be calculated. Here \(N_{\boldsymbol{q},irred}\) and \(N_{j}\) are the number of irreducible \(q\) points and the number of phonon branches, respectively. These phonon modes are distributed across NPROCS MPI threads, and phonon self-energies are calculated in parallel. OpenMP is used for the double loop over the \(N_{j}\) branches inside the calculation of each phonon self-energy. Therefore, a good performance is expected when OMP_NUM_THREADS is a divisor of \(N_{j}^{2}\).

When the calculation finishes normally, various files are generated in the working directory.

  • PREFIX.bands : Phonon dispersion along designated Brillouin zone paths
  • PREFIX.dos : (Atom projected) phonon DOS
  • PREFIX.thermo : Thermodynamic functions
  • PREFIX.msd : Mean-square displacement of atoms
  • ...

The complete list of output files can be found here.

  1. Analyze the result

One can plot the phonon dispersion relation or phonon DOS using gnuplot. Alternatively, one can use a small script in the tools/ directory for visualizing these results. For example,

$ plotband.py target.bands

shows the phonon dispersion relation. Available command line options can be displayed by

$ plotband.py -h

We also provide a similar script for phonon DOS. Another script analyze_phonons.py may be useful to analyze the result of thermal conductivity calculations. For example, phonon lifetimes and mean-free-path at 300 K can be extracted by

$ analyze_phonons.py --calc tau --temp 300 target.result

It can also estimate a cumulative thermal conductivity by

$ analyze_phonons.py --calc cumulative --temp 300 --direction 1 target.result

For details, see the tutorial.

Tutorial

Input files prepared for this tutorial are located in the example/ directory of the ALAMODE package.

Silicon

_images/si222.png

Silicon. 2x2x2 conventional supercell

In the following, (anharmonic) phonon properties of bulk silicon (Si) are calculated by a 2x2x2 conventional cell containing 64 atoms.

  1. Get displacement pattern by alm
  2. Calculate atomic forces for the displaced configurations
  3. Estimate force constants by fitting
  4. Calculate phonon dispersion and phonon DOS
  5. Estimate anharmonic IFCs for thermal conductivity
  6. Calculate thermal conductivity
  7. Analyze results
1. Get displacement patterns by alm

Change directory to example/Si and open file si_alm.in. This file is an input for the code alm which estimate interatomic force constants (IFC) by least square fitting. In the file, the crystal structure of a 2x2x2 conventional supercell of Si is specified in the &cell and the &position fields as the following:

&general
  PREFIX = si222
  MODE = suggest
  NAT = 64; NKD = 1
  KD = Si
/

&interaction
  NORDER = 1  # 1: harmonic, 2: cubic, ..
/

&cell
  20.406 # factor in Bohr unit
  1.0 0.0 0.0 # a1
  0.0 1.0 0.0 # a2
  0.0 0.0 1.0 # a3
/

&cutoff 
  Si-Si None
/

&position
  1 0.0000000000000000 0.0000000000000000 0.0000000000000000   
  1 0.0000000000000000 0.0000000000000000 0.5000000000000000
  1 0.0000000000000000 0.2500000000000000 0.2500000000000000
  1 0.0000000000000000 0.2500000000000000 0.7500000000000000
  1 0.0000000000000000 0.5000000000000000 0.0000000000000000
  1 0.0000000000000000 0.5000000000000000 0.5000000000000000
  1 0.0000000000000000 0.7500000000000000 0.2500000000000000

Replace the lattice constant of the supercell (20.406 Bohr) by your own value.

Then, execute alm

$ alm si_alm.in > si_alm.log1

which creates a file si222.pattern_HARMONIC in the working directory. In the pattern file, suggested displacement patterns are defined in Cartesian coordinates. As you can see in the file, there is only one displacement pattern for harmonic IFCs of bulk Si.

2. Calculate atomic forces for the displaced configurations

Next, calculate atomic forces for all the displaced configurations defined in si222.pattern_HARMONIC. To do so, you first need to decide the magnitude of displacements \(\Delta u\), which should be small so that anharmonic contributions are negligible. In ordinary case, \(\Delta u \sim 0.01\) Å is a reasonable choice.

Then, prepare input files necessary to run an external DFT code for each configuration. Since this procedure is a little tiresome, we provide a subsidiary Python script for VASP, Quantum-ESPRESSO (QE), and xTAPP. Using the script displace.py in the tools/ directory, you can generate the necessary input files as follows:

QE

$ python displace.py --QE=si222.pw.in --mag=0.02 si222.pattern_HARMONIC

VASP

$ python displace.py --VASP=POSCAR.orig --mag=0.02 si222.pattern_HARMONIC

xTAPP

$ python displace.py --xTAPP=si222.cg --mag=0.02 si222.pattern_HARMONIC

The --mag option specifies the displacement length in units of Angstrom. You need to specify an input file with equilibrium atomic positions either by the --QE, --VASP, or --xTAPP.

Then, calculate atomic forces for all the configurations. This can be done with a simple shell script as follows:

!#/bin/bash

# Assume we have 20 displaced configurations for QE [disp01.pw.in,..., disp20.pw.in].

for ((i=1;i<=20;i++))
do
    num=`echo $i | awk '{printf("%02d",$1)}'`
    mkdir ${num}
    cd ${num}
    cp ../disp${num}.pw.in .
    pw.x < disp${num}.pw.in > disp${num}.pw.out
    cd ../
done

Important

In QE, please set tprnfor=.true. to get atomic forces.

The next step is to collect the displacement data and force data by the Python script extract.py (also in the tools/ directory). This script can extract atomic displacements, atomic forces, and total energies from multiple output files as follows:

QE

$ python extract.py --QE=si222.pw.in --get=disp *.pw.out > disp.dat
$ python extract.py --QE=si222.pw.in --get=force *.pw.out > force.dat

VASP

$ python extract.py --VASP=POSCAR.orig --get=disp vasprun*.xml > disp.dat
$ python extract.py --VASP=POSCAR.orig --get=force vasprun*.xml > force.dat

xTAPP

$ python extract.py --xTAPP=si222.cg --get=disp *.str > disp.dat
$ python extract.py --xTAPP=si222.cg --get=force *.str > force.dat

In the above examples, atomic displacements of all the configurations are merged as disp.dat, and the corresponding atomic forces are saved in the file force.dat. These files will be used in the following fitting procedure as DFILE and FFILE. (See Format of DFILE and FFILE).

Note

For your convenience, we provide the disp.dat and force.dat files in the reference/ subdirectory. These files were generated by the Quantum-ESPRESSO package with --mag=0.02. You can proceed to the next step by copying these files to the working directory.

3. Estimate force constants by fitting

Edit the file si_alm.in to perform least-square fitting. Change the MODE = suggest to MODE = fitting as follows:

&general
  PREFIX = si222
  MODE = fitting   # <-- here
  NAT = 64; NKD = 1
  KD = Si
/

Also, add the &fitting field as:

&fitting
  NDATA = 1
  DFILE = disp.dat
  FFILE = force.dat
/

Then, execute alm again

$ alm si_alm.in > si_alm.log2

This time alm extract harmonic IFCs from the given displacement-force data set (disp.dat and force.dat above).

You can find files si222.fcs and si222.xml in the working directory. The file si222.fcs contains all IFCs in Rydberg atomic units. You can find symmetrically irreducible sets of IFCs in the first part as:

 *********************** Force Constants (FCs) ***********************
 *        Force constants are printed in Rydberg atomic units.       *
 *        FC2: Ry/a0^2     FC3: Ry/a0^3     FC4: Ry/a0^4   etc.      *
 *        FC?: Ry/a0^?     a0 = Bohr radius                          *
 *                                                                   *
 *        The value shown in the last column is the distance         *
 *        between the most distant atomic pairs.                     *
 *********************************************************************

 ----------------------------------------------------------------------
      Index              FCs         P        Pairs     Distance [Bohr]
 (Global, Local)              (Multiplicity)                           
 ----------------------------------------------------------------------

  *FC2
       1       1     2.7617453e-01   1     1x     1x       0.000
       2       2    -1.2195057e-03   2     1x     2x      10.203
       3       3    -6.2496693e-04   2     1x    33x      10.203
       4       4     8.5935598e-03   1     1x     3x       7.215
       5       5     1.9655898e-03   1     1x     3y       7.215
       6       6    -4.2490872e-03   1     1x    17x       7.215
       7       7    -3.9172885e-03   1     1x    17z       7.215
       8       8     7.9671376e-03   4     1x     6x      14.429
       9       9    -2.3601630e-04   4     1x    34x      14.429
      10      10    -6.7104831e-02   1     1x     9x       4.418
      11      11    -4.7380801e-02   1     1x     9y       4.418
      12      12     1.3282532e-03   1     1x    10x       8.460
      13      13    -8.2116209e-04   1     1x    10y       8.460
      14      14    -6.9561306e-04   1     1x    10z       8.460
      15      15     4.5271736e-04   1     1x    26x       8.460
      16      16    -4.0571255e-03   1     1x    11x      11.118
      17      17     5.0517278e-04   1     1x    11y      11.118
      18      18    -2.0691115e-04   1     1x    25x      11.118
      19      19    -4.7362015e-04   1     1x    25z      11.118
      20      20    -1.0042599e-03   2     1x    20x      12.496
      21      21     1.2863153e-03   2     1x    20y      12.496
      22      22    -1.6874372e-04   2     1x    20z      12.496
      23      23    -2.0122242e-04   2     1x    35x      12.496
      24      24     6.0327036e-04   1     1x    28x      13.254
      25      25    -5.1714558e-04   1     1x    28y      13.254

Harmonic IFCs \(\Phi_{\mu\nu}(i,j)\) in the supercell are given in the third column, and the multiplicity \(P\) is the number of times each interaction \((i, j)\) occurs within the given cutoff radius. For example, \(P = 2\) for the pair \((1x, 2x)\) because the distance \(r_{1,2}\) is exactly the same as the distance \(r_{1,2'}\) where the atom 2’ is a neighboring image of atom 2 under the periodic boundary condition. If you compare the magnitude of IFCs, the values in the third column should be divided by \(P\).

In the log file si_alm2.log, the fitting error is printed. Try

$ grep "Fitting error" si_alm2.log
Fitting error (%) : 1.47766

The other file si222.xml contains crystal structure, symmetry, IFCs, and all other information necessary for subsequent phonon calculations.

4. Calculate phonon dispersion and phonon DOS

Open the file si_phband.in and edit it for your system.

&general
  PREFIX = si222
  MODE = phonons
  FCSXML =si222.xml

  NKD = 1; KD = Si
  MASS = 28.0855
/

&cell
  10.203
  0.0 0.5 0.5
  0.5 0.0 0.5
  0.5 0.5 0.0
/

&kpoint
  1  # KPMODE = 1: line mode
  G 0.0 0.0 0.0 X 0.5 0.5 0.0 51
  X 0.5 0.5 1.0 G 0.0 0.0 0.0 51
  G 0.0 0.0 0.0 L 0.5 0.5 0.5 51
/

Please specify the XML file you obtained in step 3 by the FCSXML-tag as above. In the &cell-field, you need to define the lattice vector of a primitive cell. In phonon dispersion calculations, the first entry in the &kpoint-field should be 1 (KPMODE = 1).

Then, execute anphon

$ anphon si_phband.in > si_phband.log

which creates a file si222.bands in the working directory. In this file, phonon frequencies along the given reciprocal path are printed in units of cm-1 as:

#  G X G L
# 0.000000 0.615817 1.486715 2.020028
# k-axis, Eigenvalues [cm^-1]
0.000000   1.097373e-10   1.097373e-10   1.097373e-10   5.053714e+02   5.053714e+02   5.053714e+02
0.012316   6.132848e+00   6.132848e+00   1.033887e+01   5.052979e+02   5.052979e+02   5.053619e+02
0.024633   1.225469e+01   1.225469e+01   2.066997e+01   5.050779e+02   5.050779e+02   5.053331e+02
0.036949   1.835448e+01   1.835448e+01   3.098558e+01   5.047130e+02   5.047130e+02   5.052842e+02
0.049265   2.442116e+01   2.442116e+01   4.127802e+01   5.042056e+02   5.042056e+02   5.052137e+02
0.061582   3.044357e+01   3.044357e+01   5.153973e+01   5.035592e+02   5.035592e+02   5.051196e+02
0.073898   3.641049e+01   3.641049e+01   6.176326e+01   5.027783e+02   5.027783e+02   5.049995e+02

You can plot the phonon dispersion relation with gnuplot or any other plot software.

For visualizing phonon dispersion relations, we provide a Python script plotband.py in the tools/ directory (Matplotlib is required.). Try

$ python plotband.py si222.bands

Then, the phonon dispersion is shown by a pop-up window as follows:

_images/si_phband.png

You can save the figure as png, eps, or other formats from this window. You can also change the energy unit of phonon frequency from cm-1 to THz or meV by the --unit option. For more detail of the usage of plotband.py, type

$ python plotband.py -h

Next, let us calculate the phonon DOS. Copy si_phband.in to si_phdos.in and edit the &kpoint field as follows:

&kpoint
  2  # KPMODE = 2: uniform mesh mode
  20 20 20
/

Then, execute anphon

$ anphon si_phdos.in > si_phdos.log

This time anphon creates files si222.dos and si222.thermo in the working directory, which contain phonon DOS and thermodynamic functions, respectively. For visualizing phonon DOS and projected DOSs, we provide a Python script plotdos.py in the tools/ directory (Matplotlib is required.). The command

$ python plotdos.py --emax 550 --nokey si222.dos

will show the phonon DOS of Si by a pop-up window:

_images/si_phdos.png

To obtain more sharp DOS, try again with a denser \(k\) grid and a smaller DELTA_E value.

5. Estimate cubic IFCs for thermal conductivity

Copy file si_alm.in to si_alm2.in. Edit the &general, &interaction, and &cutoff fields of si_alm2.in as the following:

&general
  PREFIX = si222_cubic
  MODE = suggest
  NAT = 64; NKD = 1
  KD = Si
/

Change the PREFIX from si222 to si222_cubic and set MODE to suggest.

&interaction
  NORDER = 2
/

Change the NORDER tag from NORDER = 1 to NORDER = 2 to include cubic IFCs. Here, we consider cubic interaction pairs up to second nearest neighbors by specifying the cutoff radii as:

&cutoff
  Si-Si None 7.3
/

7.3 Bohr is slightly larger than the distance of second nearest neighbors (7.21461 Bohr). Change the cutoff value appropriately for your own case. (Atomic distance can be found in the file si_alm.log.)

Then, execute alm

$ alm si_alm2.in > si_alm2.log

which creates files si222_cubic.pattern_HARMONIC and si222_cubic.pattern_ANHARM3.

Then, calculate atomic forces of displaced configurations given in the file si222_cubic.pattern_ANHARM3, and collect the displacement (force) data to a file disp3.dat (force3.dat) as you did for harmonic IFCs in Steps 3 and 4.

Note

Since making disp3.dat and force3.dat requires moderate computational resources, you can skip this procedure by copying files reference/disp3.dat and reference/force3.dat to the working directory. The files we provide were generated by the Quantum-ESPRESSO package with --mag=0.04.

In si_alm2.in, change MODE = suggest to MODE = fitting and add the following:

&fitting
  NDATA = 20
  DFILE = disp3.dat
  FFILE = force3.dat
  FC2XML = si222.xml # Fix harmonic IFCs
/

By the FC2XML tag, harmonic IFCs are fixed to the values in si222.xml. Then, execute alm again

$ alm si_alm2.in > si_alm2.log2

which creates files si222_cubic.fcs and si222_cubic.xml. This time cubic IFCs are also included in these files.

Note

In the above example, we obtained cubic IFCs by least square fitting with harmonic IFCs being fixed to the value of the previous harmonic calculation. You can also estimate both harmonic and cubic IFCs simultaneously instead. To do this, merge disp.dat and disp3.dat ( and force files) as

$ cat disp.dat disp3.dat > disp_merged.dat
$ cat force.dat force3.dat > force_merged.dat

and change the &fitting field as the following:

&fitting
  NDATA = 21
  DFILE = disp_merged.dat
  FFILE = force_merged.dat
/
6. Calculate thermal conductivity

Copy file si_phdos.in to si_RTA.in and edit the MODE and FCSXML tags as follows:

&general
  PREFIX = si222
  MODE = RTA
  FCSXML = si222_cubic.xml

  NKD = 1; KD = Si
  MASS = 28.0855
/

In addition, change the \(k\) grid density as:

&kpoint
  2
  10 10 10
/

Then, execute anphon as a background job

$ anphon si_RTA.in > si_RTA.log &

Please be patient. This can take a while. When the job finishes, you can find a file si222.kl in which the lattice thermal conductivity is saved. You can plot this file using gnuplot (or any other plotting softwares) as follows:

$ gnuplot
gnuplot> set logscale xy
gnuplot> set xlabel "Temperature (K)"
gnuplot> set ylabel "Lattice thermal conductivity (W/mK)"
gnuplot> plot "si222.kl" usi 1:2 w lp
_images/si_kappa.png

Calculated lattice thermal conductivity of Si (click to enlarge)

As you can see, the thermal conductivity diverges in \(T\rightarrow 0\) limit. This occurs because we only considered intrinsic phonon-phonon scatterings in the present calculation and neglected phonon-boundary scatterings which are dominant in the low temperature range. The effect of the boundary scattering can be included using the python script analyze_phonons.py in the tools directory:

$ analyze_phonons.py --calc kappa_boundary --size 1.0e+6 si222.result > si222_boundary_1mm.kl

In this script, the phonon lifetimes are altered using the Matthiessen’s rule

\[\frac{1}{\tau_{q}^{total}} = \frac{1}{\tau_{q}^{p-p}} + \frac{2|\boldsymbol{v}_{q}|}{L}.\]

Here, the first term on the right-hand side of the equation is the scattering rate due to the phonon-phonon scattering and the second term is the scattering rate due to a grain boundary of size \(L\). The size \(L\) must be specified using the --size option in units of nm. The result is also shown in the above figure and the divergence is cured with the boundary effect.

Note

When a calculation is performed with a smearing method (ISMEAR=0 or 1) instead of the tetrahedron method (ISMEAR=-1), the thermal conductivity may have a peak structure in the very low temperature region even without the boundary effect. This peak occurs because of the finite smearing width \(\epsilon\) used in the smearing methods. As we descrease the \(\epsilon\) value, the peak value of \(\kappa\) should disapper. In addition, a very dense \(q\) grid is necessary for describing phonon excitations and thermal transport in the low temperature region (regardless of the ISMEAR value).

7. Analyze results

There are many ways to analyze the result for better understandings of nanoscale thermal transport. Some selected features are introduced below:

Phonon lifetime

The file si222.result contains phonon linewidths at irreducible \(k\) points. You can extract phonon lifetime from this file as follows:

$ analyze_phonons.py --calc tau --temp 300 si222.result > tau300K.dat
$ gnuplot
gnuplot> set xrange [1:]
gnuplot> set logscale y
gnuplot> set xlabel "Phonon frequency (cm^{-1})"
gnuplot> set ylabel "Phonon lifetime (ps)"
gnuplot> plot "tau300K.dat" using 3:4 w p
_images/si_tau.png

Phonon lifetime of Si at 300 K (click to enlarge)

In the above figure, phonon lifetimes calculated with \(20\times 20\times 20\ q\) points are also shown by open circles.

Cumulative thermal conductivity

Following the procedure below, you can obtain the cumulative thermal conductivity:

$ analyze_phonons.py --calc cumulative --temp 300 --length 10000:5 si222.result > cumulative_300K.dat
$ gnuplot
gnuplot> set logscale x
gnuplot> set xlabel "L (nm)"
gnuplot> set ylabel "Cumulative kappa (W/mK)"
gnuplot> plot "cumulative_300K.dat" using 1:2 w lp
_images/si_cumulative.png

Cumulative thermal conductivity of Si at 300 K (click to enlarge)

To draw a smooth line, you will have to use a denser \(q\) grid as shown in the figure by the orange line, which are obtained with \(20\times 20\times 20\ q\) points.

Thermal conductivity spectrum

To calculate the spectrum of thermal conductivity, modify the si_RTA.in as follows:

&general
  PREFIX = si222
  MODE = RTA
  FCSXML = si222_cubic.xml

  NKD = 1; KD = Si
  MASS = 28.0855

  EMIN = 0; EMAX = 550; DELTA_E = 1.0 # <-- frequency range
/

&cell
10.203
 0.0 0.5 0.5
 0.5 0.0 0.5
 0.5 0.5 0.0
/

&kpoint
 2
 10 10 10
/

&analysis
 KAPPA_SPEC = 1 # compute spectrum of kappa
/

The frequency range is specified with the EMIN, EMAX, and DELTA_E tags, and the KAPPA_SPEC = 1 is set in the &analysis field. Then, execute anphon again

$ anphon si_RTA.in > si_RTA2.log

After the calculation finishes, you can find the file si222.kl_spec which contains the spectra of thermal conductivity at various temperatures. You can plot the data at room temperature as follows:

$ awk '{if ($1 == 300.0) print $0}' si222.kl_spec > si222_300K.kl_spec
$ gnuplot
gnuplot> set xlabel "Frequency (cm^{-1})"
gnuplot> set ylabel "Spectrum of kappa (W/mK/cm^{-1})"
gnuplot> plot "si222_300K.kl_spec" using 2:3 w l lt 2 lw 2
_images/si_kappa_spec.png

Spectrum of thermal conductivity of Si at 300 K (click to enlarge)

In the above figure, a computational result with \(20\times 20\times 20\ q\) points is also shown by dashed line. From the figure, we can see that low energy phonons below 200 cm\(^{-1}\) account for more than 80% of the total thermal conductivity at 300 K.

Making input files for alm

Format of input files

Each input file should consist of entry fields. Available entry fields are

&general, &interaction, &cutoff, &cell, &position, and &fitting.

Each entry field starts from the key label &field and ends at the terminate character “/”. (This is equivalent to Fortran namelist.)

For example, &general entry field of program alm should be given like

&general
  # Comment line
  PREFIX = prefix
  MODE = fitting
/

Multiple entries can be put in a single line. Also, characters put on the right of sharp (“#”) will be neglected. Therefore, the above example is equivalent to the following:

&general
  PREFIX = prefix; MODE = fitting  # Comment line
/

Each variable should be written inside the appropriate entry field.

List of input variables

“&general”-field
  • PREFIX-tag : Job prefix to be used for names of output files
Default:None
Type:String

  • MODE-tag = fitting | suggest
fitting
Perform fittings to estimate harmonic and anharmonic IFCs.
This mode requires appropriate DFILE and FFILE.
suggest
This mode suggests the displacement patterns necessary
to estimate harmonic and anharmonic IFCS.
Default:None
Type:String

  • NAT-tag : Number of atoms in the supercell
Default:None
Type:Integer

  • NKD-tag : Number of atomic species
Default:None
Type:Integer

  • KD-tag = Name[1], ... , Name[NKD]
Default:None
Type:Array of strings
Example:In the case of GaAs with NKD = 2, it should be KD = Ga As.

  • NSYM-tag = 0 | 1 | nsym
0
The program automatically generates the crystal symmetry
operations (rotational and translational parts).
When PRINTSYM = 1, symmetry operations will be saved
in the file “SYMM_INFO”.
1
Only the identity operation will be considered.
nsym
“nsym” symmetry operations will be read from “SYMM_INFO”
file.
Default:0
Type:Integer

  • TOLERANCE-tag : Tolerance for finding symmetry operations
Default:1.0e-6
Type:Double

  • PRINTSYM-tag = 0 | 1
0 Symmetry operations won’t be saved in “SYMM_INFO”
1 Symmetry operations will be saved in “SYMM_INFO”
Default:0
type:Integer

  • PERIODIC-tag = PERIODIC[1], PERIODIC[2], PERIODIC[3]
0
Do not consider periodic boundary conditions when
searching for interacting atoms.
1
Consider periodic boundary conditions when
searching for interacting atoms.
Default:1 1 1
type:Array of integers
Description:This tag is useful for generating interacting atoms in low dimensional systems. When PERIODIC[i] is zero, periodic boundary condition is turned off along the direction of the lattice vector \(\boldsymbol{a}_{i}\).

  • HESSIAN-tag = 0 | 1
0 Do not save the Hessian matrix
1 Save the entire Hessian matrix of the supercell as PREFIX.hessian.
Default:0
type:Integer

“&interaction”-field
  • NORDER-tag : The order of force constants to be calculated. Anharmonic terms up to \((m+1)\)th order will be considered with NORDER = \(m\).
Default:None
Type:Integer
Example:NORDER should be 1 for harmonic calculations, and 2 to include cubic terms.

  • NBODY-tag : Entry for excluding multiple-body interactions from anharmonic force constants
Default:NBODY = [2, 3, 4, ..., NORDER + 1]
Type:Array of integers
Example:If one wants to exclude three-body interactions from cubic force constants, one should explicitly give NBODY = 2 2.

“&cutoff”-field

In this entry field, one needs to specify cutoff radii of interaction for each order in units of Bohr. In the current implementation, cutoff radii should be defined for every possible pair of atomic elements. For example, the cutoff entry for a harmonic calculation (NORDER = 1) of Si (NKD = 1) should be like

&cutoff
 Si-Si 10.0
/

This means that the cutoff radii of 10 \(a_{0}\) will be used for harmonic Si-Si terms. Please note that the first column should be two character strings, which are contained in the KD-tag, connected by a hyphen (’-’).

When one wants to consider cubic terms (NORDER = 2), please specify the cutoff radius for cubic terms in the third column as the following:

&cutoff
 Si-Si 10.0 5.6 # Pair r_{2} r_{3}
/

Instead of giving specific cutoff radii, one can write “None” as follows:

&cutoff
 Si-Si None 5.6
/

which means that all possible harmonic terms between Si-Si atoms will be included.

Caution

Setting ‘None’ for anharmonic terms can greatly increase the number of parameters and thereby increase the computational cost.

When there are more than two atomic elements, please specify the cutoff radii between every possible pair of atomic elements. In the case of MgO (NKD = 2), the cutoff entry should be like

&cutoff
 Mg-Mg 8.0
 O-O 8.0
 Mg-O 10.0
/

which can equivalently be written by using the wild card (’*’) as

&cutoff
 *-* 8.0
 Mg-O 10.0 # Overwrite the cutoff radius for Mg-O harmonic interactions
/

Important

Cutoff radii specified by an earlier entry will be overwritten by a new entry that comes later.

Once the cutoff radii are properly given, harmonic force constants \(\Phi_{i,j}^{\mu,\nu}\) satisfying \(r_{ij} \le r_{c}^{\mathrm{KD}[i]-\mathrm{KD}[j]}\) will be searched.

In the case of cubic terms, force constants \(\Phi_{ijk}^{\mu\nu\lambda}\) satisfying \(r_{ij} \le r_{c}^{\mathrm{KD}[i]-\mathrm{KD}[j]}\), \(r_{ik} \le r_{c}^{\mathrm{KD}[i]-\mathrm{KD}[k]}\), and \(r_{jk} \le r_{c}^{\mathrm{KD}[j]-\mathrm{KD}[k]}\) will be searched and determined by fitting.


“&cell”-field

Please give the cell parameters in this entry in units of Bohr as the following:

&cell
 a
 a11 a12 a13
 a21 a22 a23
 a31 a32 a33
/

The cell parameters are then given by \(\vec{a}_{1} = a \times (a_{11}, a_{12}, a_{13})\), \(\vec{a}_{2} = a \times (a_{21}, a_{22}, a_{23})\), and \(\vec{a}_{3} = a \times (a_{31}, a_{32}, a_{33})\).


“&position”-field

In this field, one needs to specify the atomic element and fractional coordinate of atoms in the supercell. Each line should be

ikd xf[1] xf[2] xf[3]

where ikd is an integer specifying the atomic element (ikd = 1, ..., NKD) and xf[i] is the fractional coordinate of an atom. There should be NAT such lines in the &position entry field.


“&fitting”-field

This field is necessary when MODE = fitting.

  • DFILE-tag : File name containing atomic displacements in Cartesian coordinate
Default:None
Type:String
Description:The format of DFILE can be found here

  • FFILE-tag : File name containing atomic forces in Cartesian coordinate
Default:None
Type:String
Description:The format of FFILE can be found here

  • NDATA-tag : Number of displacement-force data sets
Default:None
Type:Integer
Description:DFILE and FFILE should contain at least NDATA\(\times\) NAT lines.

  • NSTART, NEND-tags : Specifies the range of data to be used for fitting
Default:NSTART = 1, NEND = NDATA
Type:Integer
Example:To use the data in the range of [20:30] out of 50 entries, the tags should be NSTART = 20 and NEND = 30.

  • ICONST-tag = 0 | 1 | 2 | 3
0 No constraints
1 Constraints for translational invariance will be imposed between IFCs.
2
In addition to ICONST = 1, constraints for rotational invariance
will be imposed up to (NORDER + 1)th order.
3
In addition to ICONST = 2, constraints for rotational invariance
between (NORDER + 1)th order and (NORDER + 2)th order, which
are zero, will be considered.
11
Same as ICONST = 1 but the constraint is imposed algebraically
rather than numerically.
Default:1
Type:Integer

  • ROTAXIS-tag : Rotation axis used to estimate constraints for rotational invariance. This entry is necessary when ICONST = 2, 3.
Default:None
Type:String
Example:When one wants to consider the rotational invariance around the \(x\)-axis, one should give ROTAXIS = x. If one needs additional constraints for the rotation around the \(y\)-axis, ROTAXIS should be ROTAXIS = xy.

  • FC2XML-tag : XML file to which the harmonic terms will be fixed upon fitting
Default:None
Type:String
Description:When FC2XML-tag is given, harmonic force constants will be fixed to the values stored in the FC2XML file. This may be useful for optimizing cubic and higher-order terms without changing the harmonic terms. Please make sure that the number of harmonic terms in the new computational condition is the same as that in the FC2XML file.

  • FC3XML-tag : XML file to which the cubic terms will be fixed upon fitting
Default:None
Type:String
Description:Same as the FC2XML-tag, but FC3XML is to fix cubic force constants.

Format of DFILE and FFILE

The displacement-force data sets obtained by first-principles (or classical force-field) calculations have to be saved to DFILE and FFILE to estimate IFCs with MODE = fitting. In DFILE, please explicitly specify the atomic displacements \(u_{\alpha}(\ell\kappa)\) in units of Bohr as follows:

\[\begin{eqnarray*} u_{x}(1) & u_{y}(1) & u_{z}(1) \\ u_{x}(2) & u_{y}(2) & u_{z}(2) \\ & \vdots & \\ u_{x}(\mathrm{NAT}) & u_{y}(\mathrm{NAT}) & u_{z}(\mathrm{NAT}) \end{eqnarray*}\]

When there are NAT atoms in the supercell and NDATA data sets, there should be NAT \(\times\) NDATA lines in the DFILE without blank lines. In FFILE, please specify the corresponding atomic forces in units of Ryd/Bohr.

Making input files for anphon

Format of input files

Each input file should consist of entry fields. Available entry fields are

&general, &cell, &analysis, and &kpoint.

The format of the input file is the same as that of alm which can be found here.

List of input variables

“&general”-field
  • PREFIX-tag : Job prefix to be used for names of output files
Default:None
Type:String

  • MODE-tag = phonons | RTA
phonons
Calculate phonon dispersion relations, phonon DOS,
Grüneisen parameters etc.
RTA
Calculate phonon lifetimes and lattice thermal conductivity
based on the Boltzmann transport equation (BTE)
with the relaxation time approximation (RTA).
Default:None
Type:String

  • NKD-tag : Number of atomic species
Default:None
Type:Integer

  • KD-tag = Name[1], ... , Name[NKD]
Default:None
Type:Array of strings
Example:In the case of GaAs with NKD = 2, it should be KD = Ga As.

  • MASS-tag = mass[1], ... , mass[NKD]
Default:None
Type:Array of double
Example:In the case of Bi2Te3 with NKD = 2, MASS should be MASS = 208.98 127.60.

  • FCSXML-tag : XML file containing force constants generated by the program alm
Default:None
Type:String

  • FC2XML-tag : XML file containing harmonic force constants for different size of supercell
Default:None
Type:String
Description:When FC2XML is given, the harmonic force constants in this file will be used for calculating dynamical matrices. It is possible to use different size of supercell for harmonic and anharmonic terms, which are specified by FC2XML and FCSXML respectively.

  • NSYM-tag = 0 | 1 | nsym
0
The program automatically generates the crystal symmetry
operations (rotational and translational parts).
When PRINTSYM = 1, symmetry operations will be saved
in the file “SYMM_INFO_PRIM”.
1
Only the identity operation will be considered.
nsym
“nsym” symmetry operations will be read from “SYMM_INFO_PRIM” file.
Default:0
Type:Integer

  • TOLERANCE-tag : Tolerance for finding symmetry operations
Default:1.0e-6
Type:Double

  • PRINTSYM-tag = 0 | 1
0 Symmetry operations won’t be saved in “SYMM_INFO_PRIM”
1 Symmetry operations will be saved in “SYMM_INFO_PRIM”
Default:0
type:Integer

  • NONALAYTIC-tag = 0 | 1 | 2
0
Non-analytic correction is not considered.
1
Include the non-analytic correction by the damping method proposed by Parlinski.
2
Include the non-analytic correction by the mixed-space approach
Default:0
Type:Integer
Description:When NONANALYTIC > 0, appropriate NA_SIGMA and BORNINFO have to be given.

  • NA_SIGMA-tag : Damping factor for the non-analytic term
Default:None
Type:Double
Description:This entry is necessary when NONANALYTIC = 1. The definition of NA_SIGMA is described in the formalism section.

  • BORNINFO-tag : File containing the macroscopic dielectric tensor and Born effective charges for the non-analytic correction
Default:None
Type:String
Description:The details of the file format can be found here.

  • TMIN, TMAX, DT-tags : Temperature range and its stride in units of Kelvin
Default:TMIN = 0, TMAX = 1000, DT = 10
Type:Double

  • EMIN, EMAX, DELTA_E-tags : Energy range and its stride in units of kayser (cm-1)
Default:EMIN = 0, EMAX = 1000, DELTA_E = 10
Type:Double

  • ISMEAR-tag = -1 | 0 | 1
-1 Tetrahedron method
0 Lorentzian smearing with width of EPSILON
1 Gaussian smearing with width of EPSILON
Default:-1
Type:Integer
Description:ISMEAR specifies the method for Brillouin zone integration

  • EPSILON-tag : Smearing width in units of Kayser (cm-1)
Default:10.0
Type:Double
Description:This variable is neglected when ISMEAR = -1

  • TRISYM-tag : Flag to use symmetry operations to reduce the number of triples of \(k\) points for self-energy calculations
0 Symmetry will not be used
1 Use symmetry to reduce triples of \(k\) points
Default:1
Type:Integer
Description:This variable is used only when MODE = RTA.

Note

TRISYM = 1 can reduce the computational cost, but phonon linewidth stored to the file PREFIX.result needs to be averaged at points of degeneracy. For that purpose, a subsidiary program analyze_phonons.py* should be used.


  • RESTART-tag : Flag to restart the calculation when MODE = RTA
0 Calculate from scratch
1 Restart from the existing file
Default:1 if there is a file named PREFIX.result; 0 otherwise
Type:Integer

“&cell”-field

Please specify the cell parameters of the primitive cell as:

&cell
 a
 a11 a12 a13
 a21 a22 a23
 a31 a32 a33
/

The cell parameters are then given by \(\vec{a}_{1} = a \times (a_{11}, a_{12}, a_{13})\), \(\vec{a}_{2} = a \times (a_{21}, a_{22}, a_{23})\), and \(\vec{a}_{3} = a \times (a_{31}, a_{32}, a_{33})\).

Note

The lattice constant \(a\) must be consistent with the value used for the program alm. For example, if one used \(a = 20.4 a_{0}\) for a 2x2x2 supercell of Si, one should use \(a = 10.2 a_{0}\) here for the primitive cell.


“&kpoint”-field

This entry field is used to specify the list of \(k\) points to be calculated. The first entry KPMODE specifies the types of calculation which is followed by detailed entries.

  • KPMODE = 0 : Calculate phonon frequencies at given \(k\) points

For example, if one wants to calculate phonon frequencies at Gamma (0, 0, 0) and X (0, 1/2, 1/2) of an FCC crystal, the &kpoint entry should be written as

&kpoint
 0
 0.000 0.000 0.000
 0.000 0.500 0.500
/
  • KPMODE = 1 : Band dispersion calculation

For example, if one wants to calculate phonon dispersion relations along G-K-X-G-L of a FCC crystal, the &kpoint entry should be written as follows:

&kpoint
 1
 G 0.000 0.000 0.000  K 0.375 0.375 0.750 51
 K 0.375 0.375 0.750  X 0.500 0.500 1.000 51
 X 0.000 0.500 0.500  G 0.000 0.000 0.000 51
 G 0.000 0.000 0.000  L 0.500 0.500 0.500 51
/

The 1st and 5th columns specify the character of Brillouin zone edges, which are followed by fractional coordinates of each point. The last column indicates the number of sampling points.

  • KPMODE = 2 : Uniform \(k\) grid for phonon DOS and thermal conductivity

In order to perform a calculation with 20x20x20 \(k\) grid, the entry should be

&kpoint
 2
 20 20 20
/

“&analysis”-field
  • GRUNEISEN-tag = 0 | 1
0 Grüneisen parameters will not be calculated
1 Grüneisen parameters will be stored
Default:0
Type:Integer
Description:When MODE = phonons and GRUNEISEN = 1, Grüneisen parameters will be stored in PREFIX.gru (KPMODE = 1) or PREFIX.gru_all (KPMODE = 2).

Note

To compute Grüneisen parameters, cubic force constants must be contained in the FCSXML file.


  • PRINTEVEC-tag = 0 | 1
0 Do not print phonon eigenvectors
1 Print phonon eigenvectors in the PREFIX.evec file
Default:0
Type:Integer

  • PRINTXSF-tag = 0 | 1
0 Do not save an AXSF file
1 Create an AXSF file PREFIX.axsf
Default:0
Type:Integer
Description:This is to visualize the direction of vibrational modes at gamma (0, 0, 0) by XCrySDen. This option is valid only when MODE = phonons.

  • PRINTVEL-tag = 0 | 1
0 Do not print group velocity
1 Store phonon velocities to a file
Default:0
Type:Integer
Description:When MODE = phonons and PRINTVEL = 1, group velocities of phonons will be stored in PREFIX.phvel (KPMODE = 1) or PREFIX.phvel_all (KPMODE = 2).

  • PRINTMSD-tag = 0 | 1
0 Do not print mean-square-displacement (MSD) of atoms
1 Save MSD of atoms to the file PREFIX.mds
Default:0
Type:Integer
Description:This flag is available only when MODE = phonons and KPMODE = 2.

  • PDOS-tag = 0 | 1
0 Only the total DOS will be printed in PREFIX.dos
1 Atom-projected phonon DOS will be stored in PREFIX.dos
Default:0
Type:Integer
Description:This flag is available only when MODE = phonons and KPMODE = 2.

  • TDOS-tag = 0 | 1
0 Do not compute two-phonon DOS
1 Two-phonon DOSs will be stored in PREFIX.tdos
Default:0
Type:Integer
Description:This flag is available only when MODE = phonons and KPMODE = 2.

Note

Calculation of two-phonon DOS is computationally expensive.


  • SPS-tag = 0 | 1 | 2
0 Do not compute scattering phase space
1 Total and mode-decomposed scattering phase space involving three-phonon processes will be stored in PREFIX.sps
2 Three-phonon scattering phase space with the Bose factor will be stored in PREFIX.sps_Bose
Default:0
Type:Integer
Description:This flag is available only when MODE = phonons and KPMODE = 2.

  • PRINTPR-tag = 0 | 1
0 Do not compute the (atomic) participation ratio
1 Compute participation ratio and atomic participation ratio, which will be stored in PREFIX.pr and PREFIX.apr respectively.
Default:0
Type:Integer
Description:This flag is available when MODE = phonons.

  • KAPPA_SPEC-tag = 0 | 1
0 Do not compute the thermal conductivity spectra
1 Compute the thermal conductivity spectra, which will be stored in PREFIX.kappa_spec .
Default:0
Type:Integer
Description:This flag is available when MODE = RTA.

  • ISOTOPE-tag = 0 | 1
0 Do not consider phonon-isotope scatterings
1 Consider phonon-isotope scatterings
2 Consider phonon-isotope scatterings as in ISOTOPE = 1 and the corresponding selfenergy will be stored in PREFIX.gamma_isotope
Default:0
Type:Integer
Description:When MODE = RTA and ISOTOPE = 1 or 2, phonon scatterings due to isotopes will be considered perturbatively. ISOFACT should be properly given.

  • ISOFACT-tag = isofact[1], ... , isofact[NKD]
Default:0
Type:Array of doubles
Description:Isotope factor is a dimensionless value defined by \(\sum_{i} f_{i} (1 - m_{i}/\bar{m})^{2}\). Here, \(f_{i}\) is the fraction of the \(i\)th isotope of an element having mass \(m_{i}\), and \(\bar{m}=\sum_{i}f_{i}m_{i}\) is the average mass, respectively. This quantity is equivalent to \(g_{2}\) appearing in the original paper by S. Tamura [Phys. Rev. B, 27, 858.].

  • ANIME-tag = k1, k2, k3
Default:None
Type:Array of doubles
Description:This tag is to animate vibrational mode. k1, k2, and k3 specify the momentum of phonon modes to animate, which should be given in units of the reciprocal lattice vector. For example, ANIME = 0.0 0.0 0.5 will animate phonon modes at (0, 0, 1/2). When ANIME is given, ANIME_CELLSIZE is also necessary. You can choose the format of animation files, either AXSF or XYZ, by ANIME_FORMAT tag.

  • ANIME_CELLSIZE-tag = L1, L2, L3
Default:None
Type:Array of integers
Description:This tag specifies the cell size for animation. L1, L2, and L3 should be large enough to be commensurate with the reciprocal point given by the ANIME tag.

  • ANIME_FORMAT = xsf | xyz
Default:xyz
Type:String
Description:When ANIME_FORMAT = xsf, PREFIX.anime???.axsf files are created for XcrySDen. When ANIME_FORMAT = xyz, PREFIX.anime???.xyz files are created for VMD (and any other supporting software such as Jmol).

Format of BORNINFO

When one wants to consider the LO-TO splitting near the \(\Gamma\) point, it is necessary to set NONANALYTIC = 1 and provide BORNINFO file containing dielectric tensor \(\epsilon^{\infty}\) and Born effective charge \(Z^{*}\). In BORNINFO file, the dielectric tensor should be written in first 3 lines which are followed by Born effective charge tensors for each atom as the following.

\[\begin{eqnarray*} \epsilon_{xx}^{\infty} & \epsilon_{xy}^{\infty} & \epsilon_{xz}^{\infty} \\ \epsilon_{yx}^{\infty} & \epsilon_{yy}^{\infty} & \epsilon_{yz}^{\infty} \\ \epsilon_{zx}^{\infty} & \epsilon_{zy}^{\infty} & \epsilon_{zz}^{\infty} \\ Z_{1,xx}^{*} & Z_{1,xy}^{*} & Z_{1,xz}^{*} \\ Z_{1,yx}^{*} & Z_{1,yy}^{*} & Z_{1,zz}^{*} \\ Z_{1,zx}^{*} & Z_{1,zy}^{*} & Z_{1,zz}^{*} \\ & \vdots & \\ Z_{N_p,xx}^{*} & Z_{N_p,xy}^{*} & Z_{N_p,xz}^{*} \\ Z_{N_p,yx}^{*} & Z_{N_p,yy}^{*} & Z_{N_p,zz}^{*} \\ Z_{N_p,zx}^{*} & Z_{N_p,zy}^{*} & Z_{N_p,zz}^{*} \\ \end{eqnarray*}\]

Here, \(N_p\) is the number of atoms contained in the primitive cell.

Attention

Please pay attention to the order of Born effective charges.

List of output files

Output files of alm

  • PREFIX.pattern_HARMONIC, PREFIX.pattern_ANHARM?
In these files, displacement patterns are printed in units of \(\boldsymbol{e}_{x,y,z}\). These files are created when MODE = suggest. Patterns for anharmonic force constants are printed only when NORDER > 1.
  • PREFIX.fcs
Harmonic and anharmonic force constants in Rydberg atomic units. In the first section, only symmetry-reduced force constants are printed. All symmetry-related force constants are shown in the following section with the symmetry prefactor (\(\pm 1\)). Created when MODE = fitting.
  • PREFIX.xml
A XML file containing necessary information for performing phonon calculations. The files can be read by anphon using the FCSXML-tag. Created when MODE = fitting.

Output files of anphon

  • PREFIX.bands
Phonon dispersion along given \(k\) paths in units of cm -1. Created when MODE = phonons with KPMODE = 1.
  • PREFIX.dos
Phonon density of states (DOS). Atom projected phonon DOSs are also printed when PDOS = 1. Created when MODE = phonons with KPMODE = 2.
  • PREFIX.tdos
Two-phonon density of states for all irreducible :math:’k’ points. Created when MODE = phonons with KPMODE = 2 and TDOS = 1.
  • PREFIX.thermo
Constant volume heat capacity, vibrational entropy, internal energy, and vibrational free energy. Created when MODE = phonons with KPMODE = 2.
  • PREFIX.msd
Mean-square-displacements of atoms. Created when MODE = phonons with KPMODE = 2 and PRINTMSD = 1.
  • PREFIX.sps
Total and mode-decomposed scattering phase space. Created when MODE = phonons with KPMODE = 2 and SPS = 1.
  • PREFIX.pr
Participation ratio of every phonon modes. Created when MODE = phonons and PRINTPR = 1.
  • PREFIX.apr
Atomic participation ratio of every phonon modes. Created when MODE = phonons and PRINTPR = 1.
  • PREFIX.phvel
Phonon group velocity along given \(k\) paths. Created when MODE = phonons with KPMODE = 1 and PRINTVEL = 1.
  • PREFIX.phvel_all
Magnitude of group velocity \(|\boldsymbol{v}|\) of all phonon modes at the uniform \(k\) grid. Created when MODE = phonons with KPMODE = 2 and PRINTVEL = 1.
  • PREFIX.evec
Eigenvalues and eigenvectors of dynamical matrices. Eigenvalues are printed in Rydberg atomic units. Created when MODE = phonons with PRINTEVEC = 1.
  • PREFIX.gru
Grüneisen parameters along given \(k\) paths. Created when MODE = phonons with KPMODE = 1 and GRUNEISEN = 1.
  • PREFIX.gru_all
Grüneisen parameters of all phonon modes at the uniform \(k\) grid. Created when MODE = phonons with KPMODE = 2 and GRUNEISEN = 1.
  • PREFIX.axsf
Zone-center phonon modes with directions indicated by arrows. This file can be visualized by XcrySDen. Created when MODE = phonons with PRINTXSF = 1.
  • PREFIX.anime???.axsf and PREFIX.anime???.xyz
Files for animating phonon modes. ??? is the mode number. Created when MODE = phonons with a proper ANIME-tag. If ANIME_FORMAT = xsf, axsf files will be created which can be displayed by XcrySDen. If ANIME_FORMAT = xyz, xyz files will be created which can be visualized by VMD, Jmol, etc.

  • PREFIX.result
In this file, phonon frequency, group velocity, and anharmonic phonon linewidths are printed. This file is updated during thermal conductivity calculations (MODE = RTA). In addition, this file is read when the restart mode is turned on (RESTART = 1).
  • PREFIX.kl

Lattice thermal conductivity tensor. Created when MODE = RTA.

  • PREFIX.kl_spec

Spectra of lattice thermal conductivity. Only diagonal components will be saved. Created when MODE = RTA and KAPPA_SPEC = 1.

  • PREFIX.gamma_isotope
Phonon selfenergy due to isotope scatterings calculated by the Tamura’s formula. Created when MODE = RTA and ISOTOPE = 2.

Formalism for program alm

Interatomic force constants (IFCs)

The starting point of the computational methodology is to approximate the potential energy of interacting atoms by a Taylor expansion with respect to atomic displacements by

(1)\[\begin{split}&U - U_{0} = \sum_{n=1}^{N} U_{n} = U_{1} + U_{2} + U_{3} + \cdots, \\ &U_{n} = \frac{1}{n!} \sum_{\substack{\ell_{1}\kappa_{1}, \dots, \ell_{n}\kappa_{n} \\ \mu_{1},\dots, \mu_{n}}} \Phi_{\mu_{1}\dots\mu_{n}}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n}) \; u_{\mu_{1}}(\ell_{1}\kappa_{1})\cdots u_{\mu_{n}}(\ell_{n}\kappa_{n}).\end{split}\]

Here, \(u_{\mu}(\ell\kappa)\) is the atomic displacement of \(\kappa\)th atom in the \(\ell\)th unit cell along \(\mu\)th direction, and \(\Phi_{\mu_{1}\dots\mu_{n}}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n})\) is the \(n\)th-order interatomic force constant (IFC).

Symmetry relationship between IFCs

The are several relationships between IFCs which may be used to reduce the number of independence IFCs.

  • Permutation

Firstly, IFCs should be invariant under the exchange of triplet \((\ell,\kappa,\mu)\), e.g.

\[\Phi_{\mu_{1}\mu_{2}\mu_{3}}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2};\ell_{3}\kappa_{3})=\Phi_{\mu_{1}\mu_{3}\mu_{2}}(\ell_{1}\kappa_{1};\ell_{3}\kappa_{3};\ell_{2}\kappa_{2})=\dots.\]
  • Periodicity

Secondly, since IFCs should depend on interatomic distances, they are invariant under a translation in units of lattice vector, namely

\[\Phi_{\mu_{1}\mu_{2}\dots\mu_{n}}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2};\dots;\ell_{n}\kappa_{n})=\Phi_{\mu_{1}\mu_{2}\dots\mu_{n}}(0\kappa_{1};\ell_{2}-\ell_{1}\kappa_{2};\dots;\ell_{n}-\ell_{1}\kappa_{n}).\]
  • Crystal symmetry

A crystal symmetry operation maps an atom \(\vec{r}(\ell\kappa)\) to another equivalent atom \(\vec{r}(LK)\) by rotation and translation. Since the potential energy is invariant under any crystal symmetry operations, IFCs should transform under a symmetry operation as follows:

(2)\[\sum_{\nu_{1},\dots,\nu_{n}}\Phi_{\nu_{1}\dots\nu_{n}}(L_{1}K_{1};\dots;L_{n}K_{n}) O_{\mu_{1}\nu_{1}}\cdots O_{\mu_{n}\nu_{n}} = \Phi_{\mu_{1}\dots\mu_{n}}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n}),\]

where \(O\) is the rotational matrix of the symmetry operation. Let \(N_{s}\) be the number of symmetry operations, there are \(N_{s}\) relationships between IFCs which may be used to find independent IFCs.

Note

In the current implementation of alm, independent IFCs are searched in Cartesian coordinate where the matrix element \(O_{\mu\nu}\) is 0 or \(\pm1\) in all symmetry operations except for those of hexagonal (trigonal) lattice. Also, except for hexagonal (trigonal) systems, the product \(O_{\mu_{1}\nu_{1}}\cdots O_{\mu_{n}\nu_{n}}\) in the left hand side of equation (2) becomes non-zero only for a specific pair of \(\{\nu\}\) (and becomes 0 for all other \(\{\nu\}\)s). Therefore, let \(\{\nu^{\prime}\}\) be such a pair of \(\{\nu\}\), the equation (2) can be reduced to

(3)\[\Phi_{\nu_{1}^{\prime}\dots\nu_{n}^{\prime}}(L_{1}K_{1};\dots;L_{n}K_{n}) = s \Phi_{\mu_{1}\dots\mu_{n}}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n}),\]

where \(s=\pm1\). The code employs equation (3) instead of equation (2) to reduce the number of IFCs. If IFCs of the left-hand side and the right-hand side of equation (3) are equivalent and the coupling coefficient is \(s=-1\), the IFC is removed since it becomes zero. For hexagonal (trigonal) systems, there can be symmetry operations where multiple terms in the left-hand side of equation (2) become non-zero. For such cases, equation (2) is not used to reduce the number of IFCs. Alternatively, the corresponding symmetry relationships are imposed as constraints between IFCs in solving fitting problems.

Constraints between IFCs

Since the potential energy is invariant under rigid translation and rotation, it may be necessary for IFCs to satisfy corresponding constraints.

The constraints for translational invariance are given by

(4)\[\sum_{\ell_{1}\kappa_{1}}\Phi_{\mu_{1}\mu_{2}\dots\mu_{n}}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2};\dots;\ell_{n}\kappa_{n}) = 0,\]

which should be satisfied for arbitrary pairs of \(\ell_{2}\kappa_{2},\dots,\ell_{n}\kappa_{n}\) and \(\mu_{1},\dots,\mu_{n}\). The code alm imposes equation (4) by default (ICONST = 1).

The constraints for rotational invariance are

\[\begin{split}&\sum_{\ell^{\prime}\kappa^{\prime}} (\Phi_{\mu_{1}\dots\mu_{n}\nu}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n};\ell^{\prime}\kappa^{\prime}) r_{\mu}(\ell^{\prime}\kappa^{\prime}) - \Phi_{\mu_{1}\dots\mu_{n}\mu}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n};\ell^{\prime}\kappa^{\prime}) r_{\nu}(\ell^{\prime}\kappa^{\prime})) \\ & \hspace{10mm} + \sum_{\lambda = 1}^{n}\sum_{\mu_{\lambda}^{\prime}} \Phi_{\mu_{1}\dots\mu_{\lambda}^{\prime}\dots\mu_{n}}(\ell_{1}\kappa_{1};\dots;\ell_{\lambda}\kappa_{\lambda};\dots;\ell_{n}\kappa_{n}) (\delta_{\mu,\mu_{\lambda}}\delta_{\nu,\mu_{\lambda}^{\prime}} - \delta_{\nu,\mu_{\lambda}}\delta_{\mu,\mu_{\lambda}^{\prime}}) = 0,\end{split}\]

which must be satisfied for arbitrary pairs of \((\ell_{1}\kappa_{1},\dots,\ell_{n}\kappa_{n};\mu_{1},\dots,\mu_{n};\mu,\nu)\). This is complicated since \((n+1)\)th-order IFCs (first line) are related to \(n\)th-order IFCs (second line).

For example, the constraints for rotational invariance related to harmonic terms can be found as

(5)\[\begin{split}&\sum_{\ell_{2}\kappa_{2}} (\Phi_{\mu_{1}\nu}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2})r_{\mu}(\ell_{2}\kappa_{2})-\Phi_{\mu_{1}\mu}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2})r_{\nu}(\ell_{2}\kappa_{2})) \notag \\ & \hspace{10mm} + \Phi_{\nu}(\ell_{1}\kappa_{1})\delta_{\mu,\mu_{1}} - \Phi_{\mu}(\ell_{1}\kappa_{1})\delta_{\nu,\mu_{1}} = 0,\end{split}\]

and

(6)\[\begin{split}& \sum_{\ell_{3}\kappa_{3}} (\Phi_{\mu_{1}\mu_{2}\nu}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2};\ell_{3}\kappa_{3}) r_{\mu}(\ell_{3}\kappa_{3}) \notag - \Phi_{\mu_{1}\mu_{2}\mu}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2};\ell_{3}\kappa_{3}) r_{\nu}(\ell_{3}\kappa_{3})) \\ & \hspace{10mm} + \Phi_{\nu\mu_{2}}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2})\delta_{\mu,\mu_{1}} - \Phi_{\mu\mu_{2}}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2})\delta_{\nu,\mu_{1}} \notag \\ & \hspace{10mm} + \Phi_{\mu_{1}\nu}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2})\delta_{\mu,\mu_{2}} - \Phi_{\mu_{1}\mu}(\ell_{1}\kappa_{1};\ell_{2}\kappa_{2})\delta_{\nu,\mu_{2}} = 0.\end{split}\]

When NORDER = 1, equation (5) will be considered if ICONST = 2, whereas equation (6) will be neglected. To further consider equation (6), please use ICONST = 3, though it may enforce a number of harmonic IFCs to be zero since cubic terms don’t exist in harmonic calculations (NORDER = 1).

Estimate IFCs by least-square fitting

The code alm extracts harmonic and anharmonic IFCs from a displacement-force data set by solving the following linear least-square problem:

(7)\[\text{minimize} \ \ \chi^{2} = \sum_{t}^{m} \sum_{i} \|F_{i,t}^{\mathrm{DFT}} - F_{i,t}^{\mathrm{ALM}} \|^{2}.\]

Here, \(m\) is the number of atomic configurations and the index \(i = (\ell,\kappa,\mu)\) is the triplet of coordinates. The model force \(F_{i,t}^{\mathrm{ALM}}\) is a linear function of IFCs \(\{\Phi\}\) which can be obtained by differentiating \(U\) (Eq. (1)) by \(u_{i}\). The parameters (IFCs) are determined so as to best mimic the atomic forces obtained by DFT calculations, \(F_{i,t}^{\mathrm{DFT}}\).

To evaluate goodness of fit, alm reports the relative error \(\sigma\) defined by

(8)\[\sigma = \sqrt{\frac{\chi^{2}}{\sum_{t}^{m}\sum_{i} (F_{i,t}^{\mathrm{DFT}})^{2}}},\]

where the numerator is the residual of fit and the denominator is the square sum of DFT forces.

Formalism for program anphon

Dynamical matrix

The dynamical matrix is given by

(9)\[D_{\mu\nu}(\kappa\kappa^{\prime};\boldsymbol{q}) = \frac{1}{\sqrt{M_{\kappa}M_{\kappa^{\prime}}}} \sum_{\ell^{\prime}}\Phi_{\mu\nu}(\ell\kappa;\ell^{\prime}\kappa^{\prime})\exp{\left[i\boldsymbol{q}\cdot(\boldsymbol{r}(\ell^{\prime})-\boldsymbol{r}(\ell))\right]},\]

where \(M_{\kappa}\) is the atomic mass of atom \(\kappa\). By diagonalizing the dynamical matrix, one can obtain \(m\) (\(=3N_{\kappa}\)) eigenvalues \(\omega_{\boldsymbol{q}j}^{2}\) (\(j = 1, 2, \dots, m\)) and corresponding eigenvectors \(\boldsymbol{e}_{\boldsymbol{q}j}\) for each \(\boldsymbol{q}\) point. Here, \(\boldsymbol{e}_{\boldsymbol{q}j}\) is a column vector consisting of atomic polarization \(e_{\mu}(\kappa;\boldsymbol{q}j)\). Let \(D(\boldsymbol{q})\) denote a matrix form of equation (9), the eigenvalues may be written as

(10)\[\omega_{\boldsymbol{q}j}^{2} = (\boldsymbol{e}_{\boldsymbol{q}j}^{*})^{\mathrm{T}} D(\boldsymbol{q})\boldsymbol{e}_{\boldsymbol{q}j}.\]

Next, we introduce \(m\times m\) matrices \(\Lambda\) and \(W\) which are defined as \(\Lambda(\boldsymbol{q}) = \mathrm{diag} (\omega_{\boldsymbol{q}1}^{2},\dots,\omega_{\boldsymbol{q}m}^{2})\) and \(W(\boldsymbol{q}) = (\boldsymbol{e}_{\boldsymbol{q}1},\dots,\boldsymbol{e}_{\boldsymbol{q}m})\), respectively. Then, equation (10) can be denoted as

\[\Lambda(\boldsymbol{q}) = W^{\dagger}(\boldsymbol{q})D(\boldsymbol{q})W(\boldsymbol{q}).\]

When one needs to capture the LO-TO splitting near the zone-center by the supercell approach, it is necessary to add the non-analytic part of the dynamical matrix defined by

\[D_{\mu\nu}^{\mathrm{NA}}(\kappa\kappa^{\prime};\boldsymbol{q}) = \frac{1}{\sqrt{M_{\kappa}M_{\kappa^{\prime}}}} \frac{4\pi e^{2}}{\Omega} \frac{(Z_{\kappa}^{*}\boldsymbol{q})_{\mu}(Z_{\kappa^{\prime}}^{*}\boldsymbol{q})_{\nu}}{\boldsymbol{q}\cdot\epsilon^{\infty}\boldsymbol{q}},\]

where \(\Omega\) is the volume of the primitive cell, \(Z_{\kappa}^{*}\) is the Born effective charge tensor of atom \(\kappa\), and \(\epsilon^{\infty}\) is the dielectric constant tensor, respectively. In program anphon, either the Parlinski’s way [1] or the mixed-space approach [2] can be used. In the Parlinski’s approach (NONANALYTIC = 1), the total dynamical matrix is given by

\[D(\boldsymbol{q}) + D^{\textrm{NA}}(\boldsymbol{q})\exp{(-q^{2}/\sigma^{2})},\]

where \(\sigma\) is a damping factor. \(\sigma\) must be chosen carefully so that the non-analytic contribution becomes negligible at Brillouin zone boundaries. In the mixed-space approach (NONANALYTIC = 2), the total dynamical matrix is given by

\[D(\boldsymbol{q}) + D^{\textrm{NA}}(\boldsymbol{q})\frac{1}{N}\sum_{\ell^{\prime}}\exp{\left[i\boldsymbol{q}\cdot(\boldsymbol{r}(\ell^{\prime})-\boldsymbol{r}(\ell))\right]}.\]

The second term vanishes at commensurate \(\boldsymbol{q}\) points other than \(\Gamma\) point (\(\boldsymbol{q} = 0\)).

To include the non-analytic term, one needs to set NONANALYTIC > 0 and give appropriate BORNINFO and NA_SIGMA tags.

Group velocity

The group velocity of phonon mode \(\boldsymbol{q}j\) is given by

\[\boldsymbol{v}_{\boldsymbol{q}j} = \frac{\partial \omega_{\boldsymbol{q}j}}{\partial \boldsymbol{q}}.\]

To evaluate the group velocity numerically, we employ a central difference where \(\boldsymbol{v}\) may approximately be given by

\[\boldsymbol{v}_{\boldsymbol{q}j} \approx \frac{\omega_{\boldsymbol{q}+\Delta\boldsymbol{q}j} - \omega_{\boldsymbol{q}-\Delta\boldsymbol{q}j}}{2\Delta\boldsymbol{q}}.\]

If one needs to save the group velocities, please turn on the PRINTVEL-tag.

Thermodynamics functions

The specific heat at constant volume \(C_{\mathrm{v}}\), the internal energy \(U\), the vibrational entropy \(S\), and the Helmholtz free energy \(F\) of individual harmonic oscillator are given as follows:

\[\begin{align*} U &= \frac{1}{N_{q}}\sum_{\boldsymbol{q},j} \hbar\omega_{\boldsymbol{q}j} \left[\frac{1}{e^{\hbar\omega_{\boldsymbol{q}j}/kT} - 1} + \frac{1}{2}\right], \\ C_{\mathrm{v}} &= \frac{k}{N_{q}}\sum_{\boldsymbol{q},j} \left(\frac{\hbar\omega_{\boldsymbol{q}j}}{2kT}\right)^{2} \mathrm{cosech}^{2}\left(\frac{\hbar\omega_{\boldsymbol{q}j}}{2kT}\right),\\ S &= \frac{k}{N_{q}}\sum_{\boldsymbol{q},j} \left[\frac{\hbar\omega_{\boldsymbol{q}j}}{kT} \frac{1}{e^{\hbar\omega_{\boldsymbol{q}j}/kT} - 1} - \log{\left( 1 - e^{-\hbar\omega_{\boldsymbol{q}j}/kT}\right)}\right], \\ F &= \frac{1}{N_{q}}\sum_{\boldsymbol{q},j}\left[ \frac{\hbar\omega_{\boldsymbol{q}j}}{2} + kT\log{\left( 1 - e^{-\hbar\omega_{\boldsymbol{q}j}/kT}\right)} \right]. \end{align*}\]

Here, \(k\) is the Boltzmann constant. These quantities will be saved in the PREFIX.thermo file.

Mean square displacement

The mean square displacement tensor of atom \(\kappa\) is given by

\[\begin{align} \left< u_{\mu}(\kappa)u_{\nu}(\kappa) \right> & = \frac{\hbar}{2M_{\kappa}N_{q}}\sum_{\boldsymbol{q},j} \frac{1}{2\omega_{\boldsymbol{q}j}}\left(e_{\mu}(\kappa;\boldsymbol{q}j)e_{\nu}^{*}(\kappa;\boldsymbol{q}j)+ e_{\mu}^{*}(\kappa;\boldsymbol{q}j)e_{\nu}(\kappa;\boldsymbol{q}j)\right) \notag \\ & \hspace{25mm} \times \coth{\left(\frac{\hbar\omega_{\boldsymbol{q}j}}{2kT}\right)}. \end{align}\]

When PRINTMSD is turned on, the code print the diagonal part of the mean square displacement tensor

\[\begin{split}\left< u_{\mu}^{2}(\kappa)\right> = \frac{\hbar}{M_{\kappa}N_{q}}\sum_{\boldsymbol{q},j}\frac{1}{\omega_{\boldsymbol{q}j}} |e_{\mu}(\kappa;\boldsymbol{q}j)|^{2} \left(n_{\boldsymbol{q}j}+\frac{1}{2}\right),\end{split}\]

where \(n_{\boldsymbol{q}j} = 1/(e^{\hbar\omega_{\boldsymbol{q}j}/kT}-1)\) is the Bose-Einstein distribution function.

Phonon DOS

When KPMODE = 2, the program anphon saves the (one) phonon density of states (DOS) to the file PREFIX.dos. The one-phonon DOS is given by

\[\mathrm{DOS}(\omega) = \frac{1}{N_{q}}\sum_{\boldsymbol{q},j}\delta(\omega - \omega_{\boldsymbol{q}j}).\]

If PDOS = 1 is given, the program also prints the atom-projected phonon DOS which is given by

\[\mathrm{PDOS}(\kappa;\omega) = \frac{1}{N_{q}}\sum_{\boldsymbol{q},j}|\boldsymbol{e}(\kappa;\boldsymbol{q}j)|^{2}\delta(\omega - \omega_{\boldsymbol{q}j}).\]

In addition, TDOS-tag is available to compute the two-phonon DOS defined by

\[\mathrm{DOS2}(\omega;\boldsymbol{q};\pm) = \frac{1}{N_{q}}\sum_{\boldsymbol{q}_{1},\boldsymbol{q}_{2}, j_{1}, j_{2}} \delta(\omega\pm\omega_{\boldsymbol{q}_{1}j_{1}}-\omega_{\boldsymbol{q}_{2}j_{2}})\delta_{\boldsymbol{q}\pm\boldsymbol{q}_{1},\boldsymbol{q}_{2}+\boldsymbol{G}},\]

where \(\boldsymbol{G}\) is a reciprocal lattice vector. The sign \(\pm\) correspond to absorption and emission processes, respectively. Please note that the computation of the two-phonon DOS can be expensive especially when \(N_{q}\) or \(N_{\kappa}\) is large.

(Atomic) participation ratio

Participation ratio (PR) and atomic participation ratio (APR) defined in the following may be useful to analyze the localized nature of the phonon mode \(\boldsymbol{q}j\).

  • Participation ratio (PR)
\[PR_{\boldsymbol{q}j} = \left(\sum_{\kappa}^{N_{\kappa}} \frac{|\boldsymbol{e}(\kappa;\boldsymbol{q}j)|^{2}}{M_{\kappa}}\right)^{2} \Bigg/ N_{\kappa} \sum_{\kappa}^{N_{\kappa}} \frac{|\boldsymbol{e}(\kappa;\boldsymbol{q}j)|^{4}}{M_{\kappa}^{2}}\]
  • Atomic participation ratio (APR)
\[APR_{\boldsymbol{q}j,\kappa} = \frac{|\boldsymbol{e}(\kappa;\boldsymbol{q}j)|^{2}}{M_{\kappa}} \Bigg/ \left( N_{\kappa} \sum_{\kappa}^{N_{\kappa}} \frac{|\boldsymbol{e}(\kappa;\boldsymbol{q}j)|^{4}}{M_{\kappa}^{2}} \right)^{1/2}\]

For an extended eigenmode, the PR value is of order 1, whereas for a localized eigenmodes PR is of order \(1/N_{\kappa}\) [3]. APR is an atomic decomposition of PR that satisfies \(PR_{\boldsymbol{q}j} = \sum_{\kappa} (APR_{\boldsymbol{q}j,\kappa})^{2}\). To print the PR and APR, please set MODE = phonons and PRINTPR = 1 in the &analysis entry field.

Scattering phase space

When KPMODE = 2 and SPS = 1, the three-phonon scattering phase space \(P_{3}\) is calculated and saved to the file PREFIX.sps. \(P_{3}\) is defined as

\[P_{3}(\boldsymbol{q}j) = \frac{1}{3m^{3}} (2P_{3}^{(+)}(\boldsymbol{q}j) + P_{3}^{(-)}(\boldsymbol{q}j)),\]

where \(m\) is the number of phonon branches and

\[P_{3}^{(\pm)}(\boldsymbol{q}j) = \frac{1}{N_{q}}\sum_{\boldsymbol{q}_{1},\boldsymbol{q}_{2}, j_{1}, j_{2}}\delta(\omega_{\boldsymbol{q}j}\pm\omega_{\boldsymbol{q}_{1}j_{1}}-\omega_{\boldsymbol{q}_{2}j_{2}})\delta_{\boldsymbol{q}\pm\boldsymbol{q}_{1},\boldsymbol{q}_{2}+\boldsymbol{G}}.\]

anphon also print the total scattering phase space

\[P_{3} = \frac{1}{N_{q}}\sum_{\boldsymbol{q}} P_{3}(\boldsymbol{q}j).\]

When SPS = 2, the three-phonon scattering phase space with the occupation factor \(W_{3}^{(\pm)}\) will be calculated and saved to the file PREFIX.sps_Bose. \(W_{3}^{(\pm)}\) is defined as

\[\begin{split}W_{3}^{(\pm)}(\boldsymbol{q}j) = \frac{1}{N_{q}}{\sum_{\boldsymbol{q}_{1},\boldsymbol{q}_{2}, j_{1}, j_{2}}} \left\{ \begin{array}{c} n_{2} - n_{1} \\ n_{1} + n_{2} + 1 \end{array} \right\} \delta(\omega_{\boldsymbol{q}j}\pm\omega_{\boldsymbol{q}_{1}j_{1}}-\omega_{\boldsymbol{q}_{2}j_{2}})\delta_{\boldsymbol{q}\pm\boldsymbol{q}_{1},\boldsymbol{q}_{2}+\boldsymbol{G}}.\end{split}\]

Here, \(n_{1}=n(\omega_{\boldsymbol{q}_{1}j_{1}})\) and \(n_{2}=n(\omega_{\boldsymbol{q}_{2}j_{2}})\) where \(n(\omega) = \frac{1}{e^{\hbar\omega/k_B T}-1}\) is the Bose-Einstein distribution function. Since \(n(\omega)\) is temperature dependent, \(W_{3}^{(\pm)}\) is also temperature dependent. The file PREFIX.sps_Bose contains \(W_{3}^{(\pm)}\) for all phonon modes at various temperatures specified with TMIN, TMAX, and DT tags.

Grüneisen parameter

The mode Grüneisen parameter, defined as \(\gamma_{\boldsymbol{q}j} = - \frac{\partial \log{\omega_{\boldsymbol{q}j}}}{\partial \log{V}}\), is calculated by

\[\gamma_{\boldsymbol{q}j}= -\frac{(\boldsymbol{e}_{\boldsymbol{q}j}^{*})^{\mathrm{T}} \delta D(\boldsymbol{q})\boldsymbol{e}_{\boldsymbol{q}j}}{6\omega_{\boldsymbol{q}j}},\]

where \(\delta D(\boldsymbol{q})\) is a change in the dynamical matrix due to a volume change \(\delta V\), which is given by

\[\begin{align} \delta D_{\mu\nu}(\kappa\kappa^{\prime};\boldsymbol{q}) &= \frac{1}{\sqrt{M_{\kappa}M_{\kappa^{\prime}}}} \sum_{\ell^{\prime}}\delta\Phi_{\mu\nu}(\ell\kappa;\ell^{\prime}\kappa^{\prime})\exp{\left[i\boldsymbol{q}\cdot(\boldsymbol{r}(\ell^{\prime})-\boldsymbol{r}(\ell))\right]},\\ \delta\Phi_{\mu\nu}(\ell\kappa;\ell^{\prime}\kappa^{\prime}) &= \sum_{\ell^{\prime\prime},\kappa^{\prime\prime},\lambda}\Phi_{\mu\nu\lambda}(\ell\kappa;\ell^{\prime}\kappa^{\prime};\ell^{\prime\prime}\kappa^{\prime\prime})r_{\lambda}(\ell^{\prime\prime}\kappa^{\prime\prime}). \end{align}\]

Please set GRUNEISEN = 1 and give an appropriate FCSXML file containing cubic IFCs to print Grüneisen parameters.

Anharmonic self-energy

The anharmonic self-energy due to cubic anharmonicity to the lowest order is given by

(11)\[\begin{split}\Sigma_{\boldsymbol{q}j}(i\omega_m) &= \frac{1}{2\hbar^{2}}\sum_{\boldsymbol{q}_{1},\boldsymbol{q}_{2}}\sum_{j_{1},j_{2}} |V^{(3)}_{-\boldsymbol{q}j,\boldsymbol{q}_{1}j_{1},\boldsymbol{q}_{2}j_{2}}|^{2} \notag \\ & \times \left[ \frac{n_{1}+n_{2} + 1}{i\omega_{m} + \omega_{1} + \omega_{2}} - \frac{n_{1}+n_{2} + 1}{i\omega_{m} - \omega_{1} - \omega_{2}} + \frac{n_{1}-n_{2}}{i\omega_{m} - \omega_{1} + \omega_{2}} - \frac{n_{1}-n_{2}}{i\omega_{m} + \omega_{1} - \omega_{2}} \right],\end{split}\]

where \(i\omega_{m}\) is the Matsubara frequency. In equation (11), we simply denoted \(\omega_{\boldsymbol{q}_{i}j_{i}}\) as \(\omega_{i}\). The matrix element \(V^{(3)}\) is given by

\[\begin{split}V^{(3)}_{\boldsymbol{q}j,\boldsymbol{q}^{\prime}j^{\prime},\boldsymbol{q}^{\prime\prime}j^{\prime\prime}} & = \left( \frac{\hbar}{2N_{q}}\right)^{\frac{3}{2}} \frac{1}{\sqrt{\omega_{\boldsymbol{q}n}\omega_{\boldsymbol{q}^{\prime}j^{\prime}}\omega_{\boldsymbol{q}^{\prime\prime}j^{\prime\prime}}}} \sum_{\ell,\ell^{\prime},\ell^{\prime\prime}} \exp{\left[\mathrm{i}(\boldsymbol{q}\cdot\boldsymbol{r}(\ell)+\boldsymbol{q}^{\prime}\cdot\boldsymbol{r}(\ell^{\prime})+\boldsymbol{q}^{\prime\prime}\cdot\boldsymbol{r}(\ell^{\prime\prime}))\right]} \notag \\ & \times \sum_{\kappa,\kappa^{\prime},\kappa^{\prime\prime}} \frac{1}{\sqrt{M_{\kappa}M_{\kappa^{\prime}}M_{\kappa^{\prime\prime}}}} \sum_{\mu,\nu,\lambda} \Phi_{\mu\nu\lambda}(\ell\kappa;\ell^{\prime}\kappa^{\prime};\ell^{\prime\prime}\kappa^{\prime\prime}) e_{\mu}(\kappa;\boldsymbol{q}j)e_{\nu}(\kappa^{\prime};\boldsymbol{q}^{\prime}j^{\prime})e_{\lambda}(\kappa^{\prime\prime};\boldsymbol{q}^{\prime\prime}j^{\prime\prime}) \; ,\end{split}\]

which becomes zero unless \(\boldsymbol{q}+\boldsymbol{q}^{\prime}+\boldsymbol{q}^{\prime\prime}\) is an integral multiple of \(\boldsymbol{G}=n_{1}\boldsymbol{b}_{1}+n_{2}\boldsymbol{b}_{2}+n_{3}\boldsymbol{b}_{3}\). Phonon linewidth \(\Gamma_{\boldsymbol{q}j}\), which is the imaginary part of the phonon self-energy, can be obtained by the analytic continuation to the real axis (\(i\omega_{m}\to \omega + i0^{+}\)) as

(12)\[\begin{split} \Gamma^{\mathrm{anh}}_{\boldsymbol{q}j}(\omega) &= \frac{\pi}{2\hbar^{2}}\sum_{\boldsymbol{q}_{1},\boldsymbol{q}_{2}}\sum_{j_{1},j_{2}} |V^{(3)}_{-\boldsymbol{q}j,\boldsymbol{q}_{1}j_{1},\boldsymbol{q}_{2}j_{2}}|^{2} \notag \\ & \times \left[ -(n_{1}+n_{2} + 1)\delta{(\omega + \omega_{1} + \omega_{2})} + (n_{1}+n_{2} + 1) \delta{(\omega - \omega_{1} - \omega_{2})} \right. \notag \\ & \left. \hspace{12mm} - (n_{1}-n_{2})\delta{(\omega - \omega_{1} + \omega_{2})} + (n_{1}-n_{2})\delta{(\omega + \omega_{1} - \omega_{2})} \right].\end{split}\]

The computation of equation (12) is the most expensive part of the thermal conductivity calculations. Therefore, we employ the crystal symmetry to reduce the number of triplet pairs \((\boldsymbol{q}j,\boldsymbol{q}^{\prime}j^{\prime},\boldsymbol{q}^{\prime\prime}j^{\prime\prime})\) of \(V^{(3)}\) to calculate. To disable the reduction, please set TRISYM = 0.

Isotope scattering

The effect of isotope scatterings can be considered by the mass perturbation approach proposed by S. Tamura [4] by the ISOTOPE-tag. The corresponding phonon linewidth is given by

\[\Gamma_{\boldsymbol{q}j}^{\mathrm{iso}}(\omega)= \frac{\pi}{4N_{q}} \omega_{\boldsymbol{q}j}^{2}\sum_{\boldsymbol{q}_{1},j_{1}}\delta(\omega-\omega_{\boldsymbol{q}_{1}j_{1}}) \sum_{\kappa}g_{2}(\kappa)|\boldsymbol{e}^{*}(\kappa;\boldsymbol{q}_{1}\boldsymbol{j}_{1})\cdot\boldsymbol{e}(\kappa;\boldsymbol{q}\boldsymbol{j})|^{2},\]

where \(g_{2}\) is a dimensionless factor given by

\[g_{2}(\kappa)=\sum_{i}f_{i}(\kappa)\left(1 - \frac{m_{i}(\kappa)}{M_{\kappa}}\right)^{2}.\]

Here, \(f_{i}\) is the fraction of the \(i\)th isotope of an element having mass \(m_i\), and \(M_{\kappa}=\sum_{i}f_{i}m_{i}(\kappa)\) is the average mass, respectively. The \(g_{2}\) values should be provided by the ISOFACT-tag. The average mass \(M_{\kappa}\) is substituted by the value specified in the MASS-tag.

Lattice thermal conductivity

The lattice thermal conductivity tensor \(\kappa_{\mathrm{ph}}^{\mu\nu}(T)\) is estimated within the relaxation-time approximation as

\[\kappa_{\mathrm{ph}}^{\mu\nu}(T) = \frac{1}{\Omega N_{q}} \sum_{\boldsymbol{q},j}c_{\boldsymbol{q}j}(T)v_{\boldsymbol{q}j}^{\mu}v_{\boldsymbol{q}j}^{\nu}\tau_{\boldsymbol{q}j}(T),\]

where \(c_{\boldsymbol{q}j} = \hbar\omega_{\boldsymbol{q}j}\partial n_{\boldsymbol{q}j}/\partial T\) and \(\tau_{\boldsymbol{q}j}(T)\) is the phonon lifetime. The phonon lifetime is estimated using the Matthiessen’s rule as

\[\tau_{\boldsymbol{q}j}^{-1}(T) = 2 (\Gamma_{\boldsymbol{q}j}^{\mathrm{anh}}(T) + \Gamma_{\boldsymbol{q}j}^{\mathrm{iso}}).\]

The lattice thermal conductivity will be written to the file PREFIX.kl.

The spectra of the lattice thermal conductivity \(\kappa_{\mathrm{ph}}^{\mu\mu}(\omega)\) can also be calculated by setting KAPPA_SPEC = 1 in the &analysis field. \(\kappa_{\mathrm{ph}}^{\mu\mu}(\omega)\) is defined as

\[\kappa_{\mathrm{ph}}^{\mu\mu}(\omega) = \frac{1}{\Omega N_{q}}\sum_{\boldsymbol{q},j}c_{\boldsymbol{q}j}v_{\boldsymbol{q}j}^{\mu}v_{\boldsymbol{q}j}^{\mu}\tau_{\boldsymbol{q}j} \delta(\omega-\omega_{\boldsymbol{q}j}).\]

If we integrate this quantity over \(\omega\), we then obtain the bulk thermal conductivity, namely \(\kappa_{\mathrm{ph}}^{\mu\mu} = \int_{0}^{\infty} \kappa_{\mathrm{ph}}^{\mu\mu}(\omega) \; \mathrm{d}\omega\).

Cumulative thermal conductivity

The accumulative lattice thermal conductivity \(\kappa_{\mathrm{ph,acc}}^{\mu\nu}(L)\) is defined as

\[\kappa_{\mathrm{ph,acc}}^{\mu\mu}(L) = \frac{1}{\Omega N_{q}} \sum_{\boldsymbol{q},j}c_{\boldsymbol{q}j}v_{\boldsymbol{q}j}^{\mu}v_{\boldsymbol{q}j}^{\mu}\tau_{\boldsymbol{q}j}\Theta (L-|\boldsymbol{v}_{\boldsymbol{q}j}|\tau_{\boldsymbol{q}j}),\]

where \(\Theta(x)\) is the step function. This quantity can be calculated by using the script analyze_phonons.py with --calc cumulative flag. One can also use another definition for the accumulative thermal conductivity:

\[\kappa_{\mathrm{ph,acc}}^{\mu\nu}(L) = \frac{1}{\Omega N_{q}} \sum_{\boldsymbol{q},j}c_{\boldsymbol{q}j}v_{\boldsymbol{q}j}^{\mu}v_{\boldsymbol{q}j}^{\nu}\tau_{\boldsymbol{q}j}\Theta (L-|v_{\boldsymbol{q}j}^{\mu}|\tau_{\boldsymbol{q}j}).\]

In this case, the contribution to the total thermal conductivity is limited only from phonon modes whose mean-free-path along the \(\mu\)-direction is smaller than \(L\). To calculate this, please use the --calc cumulative2 flag and specify the direction \(\mu\) by the --direction option.

Delta function

In order to compute the phonon DOSs and the imaginary part of phonon self-energies, it is necessary to evaluate the Brillouin-zone integration containing Dirac’s delta function. For that purpose, we provide 3 options through the ISMEAR-tag.

When ISMEAR = 0, the delta function is replaced by the Lorentzian function as

\[\delta(\omega) \approx \frac{1}{\pi}\frac{\epsilon^{2}}{\omega^{2}+\epsilon^{2}}.\]

When ISMEAR = 1, the delta function is replaced by the Gaussian function as

\[\delta(\omega) \approx \frac{1}{\sqrt{\pi}\epsilon}\exp{(-\omega^{2}/\epsilon^{2})},\]

which decays faster than the Lorentzian function. For both cases, \(\epsilon\) should be given by the EPSILON-tag, which must be chosen carefully to avoid any unscientific results. \(\epsilon\) should be small enough to capture detailed phonon structures such as phonon DOS or energy conservation surface related to three-phonon processes, but it should be large enough to avoid unscientific oscillations. Choosing an appropriate value for \(\epsilon\) is not a trivial task since it may depend on the phonon structure and the density of \(\boldsymbol{q}\) points.

To avoid such issues, the program anphon employs the tetrahedron method [5] by default (ISMEAR = -1) for numerical evaluations of Brillouin zone integration containing \(\delta(\omega)\). When the tetrahedron method is used, the EPSILON-tag is neglected. We recommend using the tetrahedron method whenever possible.


[1]K. Parlinski, Z. Q. Li, and Y. Kawazoe, Phys. Rev. Lett. 81, 3298 (1998).
[2]Y. Wang et al., J. Phys.: Condens. Matter 22, 202201 (2010).
[3]J. Hafner and M. Krajci, J. Phys.: Condens. Matter 5, 2489 (1993).
[4]S. -I. Tamura, Phys. Rev. B 27, 858 (1983).
[5]P. E. Blöchl, O. Jepsen, and O. K. Andersen, Phys. Rev. B 49, 1450555 (1994).

Indices and tables