Phonons & Phonopy: Pro Tips (2014)

An overview of the Phonopy (and Phono3py) lattice-dynamics codes, covering features, examples, applications and troubleshooting.

1. Phonons & Phonopy:
Pro Tips
J. M. Skelton
WMD Research Day
10
th
October 2014

2. Phonons and Lattice Dynamics
WMD Research Day: 10
th
Oct 2014 | Slide 2
Crystallography is generally concerned with the static properties of crystals,
describing features such as the average positions of atoms and the symmetry of a
crystal. Solid state physics takes a similar line as far as elementary electronic
properties are concerned.
We know, however, that atoms actually move around inside the crystal structure,
since it is these motions that give the concept of temperature […].
The static lattice model, which is only concerned with the average positions of atoms
and neglects their motions, can explain a large number of material features […].
There are, however, a number of properties that cannot be explained by a static
model…
Martin Dove, “Introduction to Lattice Dynamics”

3. Overview
WMD Research Day: 10
th
Oct 2014 | Slide 3
•Theory
The quantum harmonic oscillator; the 3D harmonic crystal
Ab initiothermodynamics
•Harmonic phonopy
Workflow
Calculating forces: options and things to watch out for (!)
Post processing
Input/output files; “hacking” phonopyfor other calculators
•Anharmonicity1: the quasi-harmonic approximation
Theory
phonopy-qha: workflow, output and example applications
•Anharmonicity2: phonon-phonon coupling
•Theory
•phono3py: workflow, setup and post processing
•Summary

4. The Quantum Harmonic Oscillator
WMD Research Day: 10
th
Oct 2014 | Slide 4
�
�=�+
1
2
ℏ��=
�
�
= Potential energy
�
�
�
�
�
= Spring constant
= Frequency
= Reduced mass
Where:
�=−�(�−�
0) �=
1
2
�(�−�
0)

5. Imaginary Frequencies and Phase Transitions
WMD Research Day: 10
th
Oct 2014 | Slide 5U
r-r
0 U
r-r
0
&#3627408472;>0 &#3627408472;<0
•If the system is on a potential-energy maximum, there is no restoring force along certain
modes -> these will have a negativeforce constant associated with them
•Since &#3627409172;=&#3627408472;/&#3627409159;, the mode must have an imaginary frequency(usually represented as a
negative frequency in phonon DOS/band structure curves
•The mechanism for some phase transitions is for one or more modes in the stable
structure to become imaginary at the transition temperature

6. The 3D Harmonic Crystal
WMD Research Day: 10
th
Oct 2014 | Slide 6
Φ
&#3627409148;&#3627409149;&#3627408470;&#3627408473;,&#3627408471;&#3627408473;′=
&#3627409173;
2
&#3627408440;
&#3627409173;&#3627408479;
&#3627409148;(&#3627408473;)&#3627409173;&#3627408479;
&#3627409149;(&#3627408473;
′
)
=−
&#3627409173;&#3627408441;
&#3627409148;(&#3627408470;&#3627408473;)
&#3627409173;&#3627408479;
&#3627409149;(&#3627408471;&#3627408473;′)
Φ
&#3627409148;&#3627409149;&#3627408470;&#3627408473;,&#3627408471;&#3627408473;′≈−
&#3627408441;
&#3627409148;(&#3627408470;&#3627408473;)
∆&#3627408479;
&#3627409149;(&#3627408471;&#3627408473;′)
&#3627408439;
&#3627409148;&#3627409149;&#3627408470;,&#3627408471;,??????=
1
&#3627408474;
&#3627408470;&#3627408474;
&#3627408471;

&#3627408473;′
Φ
&#3627409148;&#3627409149;&#3627408470;0,&#3627408471;&#3627408473;
′
exp[&#3627408470;??????.(&#3627408531;&#3627408471;&#3627408473;
′
−&#3627408531;(&#3627408470;0))]
•The force constant matrix Φ
&#3627409148;&#3627409149;(&#3627408470;&#3627408473;,&#3627408471;&#3627408473;′)can be obtained either from finite-displacement
calculations, or using DFPT
•The number of displacements which need to be evaluated to construct the dynamical
matrix can be reduced by symmetry
Force constant matrix:
From finite differences:
Dynamical matrix:
Sum over atom &#3627408471;in adjacent unit
cells &#3627408473;′-> supercellexpansion to
improve accuracy
&#3627408466;??????.Ω??????=&#3627408439;??????.&#3627408466;(??????)After diagonalisation:

7. The 3D Harmonic Crystal
WMD Research Day: 10
th
Oct 2014 | Slide 7
&#3627408439;??????=
&#3627408439;
&#3627408485;&#3627408485;(1,1,??????)&#3627408439;
&#3627408485;&#3627408486;(1,1,??????)&#3627408439;
&#3627408485;&#3627408487;(1,1,??????)…&#3627408439;
&#3627408485;&#3627408487;(1,??????,??????)
&#3627408439;
&#3627408486;&#3627408485;(1,1,??????)&#3627408439;
&#3627408486;&#3627408486;(1,1,??????)&#3627408439;
&#3627408486;&#3627408487;(1,1,??????)…&#3627408439;
&#3627408486;&#3627408487;(1,??????,??????)
&#3627408439;
&#3627408487;&#3627408485;(1,1,??????)&#3627408439;
&#3627408487;&#3627408486;(1,1,??????)&#3627408439;
&#3627408487;&#3627408487;(1,1,??????)…&#3627408439;
&#3627408487;&#3627408485;(1,??????,??????)
⋮ ⋮ ⋮ ⋱ ⋮
&#3627408439;
&#3627408487;&#3627408485;(??????,1,&#3627408530;)&#3627408439;
&#3627408487;&#3627408486;(??????,1,&#3627408478;)&#3627408439;
&#3627408487;&#3627408487;(??????,1,??????)…&#3627408439;
&#3627408487;&#3627408487;(??????,??????,??????)
????????????,λ=
&#3627408474;
1&#3627408479;
&#3627408485;(1,??????,λ)
&#3627408474;
1&#3627408479;
&#3627408486;(1,??????,λ)
&#3627408474;
1&#3627408479;
&#3627408487;(1,??????,λ)
⋮
&#3627408474;
??????&#3627408479;
&#3627408487;(??????,??????,λ)
Ω??????=
ω(λ
1,??????) . . . .
. ω(λ
2,??????) . . .
. . ω(λ
3,??????). .
. . . ⋱ ⋮
. . . .ω(λ
3??????,??????)

8. LO/TO Splitting
WMD Research Day: 10
th
Oct 2014 | Slide 8
•Two optic modes involve rows of atoms sliding past each other, while the third involves
separation of the ions
•The latter has an extra restoring force associated with it -> LO/TO splitting
•This can be modelled by a non-analytical correction to the phonon frequencies, using the
Born effective charge and macroscopic dielectric tensors
<111>

9. Ab InitioThermodynamics
WMD Research Day: 10
th
Oct 2014 | Slide 9
Helmholtz free energy:
For a solid:
&#3627408436;(&#3627408455;)=&#3627408456;(&#3627408455;)−&#3627408455;&#3627408454;(&#3627408455;)
&#3627408436;(&#3627408455;)=&#3627408456;
??????+&#3627408456;
??????&#3627408455;−&#3627408455;&#3627408454;
??????(&#3627408455;)
Where:&#3627408436;(&#3627408455;)
&#3627408456;
??????
&#3627408456;
??????(&#3627408455;)
&#3627408454;
??????(&#3627408455;)
= Helmholtz energy
= Lattice internal energy
= Vibrational internal energy
= Vibrational entropy
Equilibrium DFT
Phonons (!)
Thermodynamics requires phonons!
&#3627408456;
??????0K=????????????&#3627408440;

10. Partition function:&#3627408436;&#3627408455;=−&#3627408472;
??????&#3627408455;ln??????(&#3627408455;)
Ab InitioThermodynamics
WMD Research Day: 10
th
Oct 2014 | Slide 10
Helmholtz energy:&#3627408436;(&#3627408455;)=&#3627408456;
??????+&#3627408456;
??????&#3627408455;−&#3627408455;&#3627408454;
??????(&#3627408455;)
??????&#3627408455;=exp[−&#3627408456;
??????/&#3627408472;
??????&#3627408455;]
??????,λ
exp[ℏ&#3627409172;(??????,λ)/2&#3627408472;
??????&#3627408455;]
1−exp[ℏ&#3627409172;(??????,λ)/&#3627408472;
??????&#3627408455;]
&#3627408456;
??????&#3627408455;=
??????,λ
ℏ&#3627409172;(??????,λ)
1
2
+
1
expℏ&#3627409172;??????,λ/&#3627408472;
??????&#3627408455;−1
Vibrational energy:
&#3627408438;
??????=
&#3627409173;&#3627408456;
&#3627409173;&#3627408455;
??????
&#3627408454;=
&#3627409173;&#3627408436;
&#3627409173;&#3627408455;
Derivatives:
Phonon occupation
number

11. phonopy: Workflow
WMD Research Day: 10
th
Oct 2014 | Slide 11
Input Structure
Create
Displacements
Calculate Forces
Extract Forces
Post-Process
phonopy-d [--dim=“1 1 1”]
phonopy-f vasprun-{001..XXX}.xml
phonopy--fc vasprun.xml
phonopy[-t] [-p] [-s] Settings.conf

12. phonopy: Setup
WMD Research Day: 10
th
Oct 2014 | Slide 12
phonopy-d
[--dim=“1 1 1”]
[--tolerance=1e-5]
[-c POSCAR]
Set up the calculations
on an XxYxZsupercell
Position tolerance for
symmetry detection
Path to POSCAR-format structure
(“cell”), if not “POSCAR”
disp.yaml
POSCAR-001
POSCAR-002
POSCAR-003
...
SPOSCAR
POSCAR files containing single atomic displacements
“Unperturbed” supercellfor DFPT phonon
calculation
Information about the structure, supercell
and displaced atoms

13. phonopy: Calculating Forces
WMD Research Day: 10
th
Oct 2014 | Slide 13
ADDGRID = .TRUE.
EDIFF = 1E-8
ENCUT = 500-800 eV
LREAL = .FALSE.
PREC = High | Accurate
ADDGRID = .TRUE.
EDIFF = 1E-8
ENCUT = 500-800 eV
IBRION = 5|6|7|8
LREAL = .FALSE.
NSW = 1
PREC = High | Accurate
•Accurate forces are essential -> crank the standard settings right up
•LREAL = .FALSE.is essential, unless you manually adjust ROPT
•ADDGRID = .TRUE.doesn’t seem to be essential, but doesn’t cost much either
•For finite-difference/DFPT phonon calculations in VASP, set NSW = 1
Sample finite-differences INCAR: Sample VASP force-constants INCAR:

14. phonopy: Calculating Forces
WMD Research Day: 10
th
Oct 2014 | Slide 14
ADDGRID = .FALSE.
LREAL = Auto
ADDGRID = .FALSE.
LREAL = .FALSE.
ADDGRID = .TRUE.
LREAL = Auto
ADDGRID = .TRUE.
LREAL = .FALSE.

15. phonopy: Calculating Forces
WMD Research Day: 10
th
Oct 2014 | Slide 15
XC = PBEXC = LDA
XC = TPSS

16. phonopy: Calculating Forces
WMD Research Day: 10
th
Oct 2014 | Slide 16
4x4x4 Primitive Cell 4x4x4 Conventional Cell
•There are options in phonopyto project a calculation on the conventional cell back to
the primitive cell during post processing

17. phonopy: Some “Pro Tips”
WMD Research Day: 10
th
Oct 2014 | Slide 17
If using VASP FD/DFPT, set NWRITE = 3in the INCAR file, and you can run this
bash script on the OUTCAR to obtain a simulated IR spectrum “for free”:
http://homepage.univie.ac.at/david.karhanek/downloads.html#Entry02
If using DFPT with an LDA/GGA functional, set LEPSILON = .TRUE.in the
INCAR file to obtain the static dielectric constant, in particular the ionic-
relaxation part, for a small added cost
When using FD/DFPT, VASP tries to change the k-point set internally, which
requires NPAR = #Coresto be set in the INCAR file; setting ISYM = -1
avoids this, and although the number of displacements which need to be
evaluated may increase, the performance gained by using band parallelism can
quite easily offset this for low-symmetry systems (!)

18. phonopy: Post Processing
WMD Research Day: 10
th
Oct 2014 | Slide 18
phonopy-p -s Settings.conf
“Plot” “Save”
Settings file
DIM = 4 4 4
MP = 48 48 48
GAMMA_CENTER = .TRUE.
Sample phonon DOS settings file:
PbTe

19. phonopy: Post Processing
WMD Research Day: 10
th
Oct 2014 | Slide 19
phonopy-p -s Settings.conf
“Plot” “Save”
Settings file
DIM = 4 4 4
MP = 48 48 48
GAMMA_CENTER = .TRUE.
EIGENVECTORS = .TRUE .
PDOS = 1, 2
Sample phonon DOS settings file:
PbTe

20. phonopy: Post Processing
WMD Research Day: 10
th
Oct 2014 | Slide 20
“[Calculate] thermal
properties”
DIM = 4 4 4
MP = 48 48 48
GAMMA_CENTER = .TRUE.
Sample phonon DOS settings file:
phonopy-p -s -tSettings.conf
“Plot” “Save”
Settings file
PbTe

21. phonopy: Post Processing
WMD Research Day: 10
th
Oct 2014 | Slide 21
phonopy-p -s Settings.conf
“Plot” “Save”
Settings file
DIM = 4 4 4
BAND = 0.0 0.0 0.0 0.5 0.25 0.75
0.5 0.0 0.5 0.0 0.0 0.0
0.5 0.5 0.5
BAND_POINTS = 101
BAND_LABELS = \Gamma W X \Gamma L
[EIGENVECTORS = .TRUE.]
Sample phonon band structure settings file:
PbTe

22. phonopy: Post Processing
WMD Research Day: 10
th
Oct 2014 | Slide 22
phonopy-p -s Settings.conf
“Plot” “Save”
Settings file
DIM = 4 4 4
BAND = 0.0 0.0 0.0 0.5 0.25 0.75
0.5 0.0 0.5 0.0 0.0 0.0
0.5 0.5 0.5
BAND_POINTS = 101
BAND_LABELS = \Gamma W X \Gamma L
BAND_CONNECTION = .TRUE.
Sample phonon band structure settings file:
PbTe

23. phonopy: Non-Analytical Corrections
WMD Research Day: 10
th
Oct 2014 | Slide 23
EDIFF = 1E-8
ENCUT = 500-800 eV
LEPSILON = .TRUE.
LREAL = .FALSE.
NSW = 0
PREC = High | Accurate
INCAR for Born charges using DFPT:
EDIFF = 1E-8
ENCUT = 500-800 eV
LCALCEPS = .TRUE.
LREAL = .FALSE. ! Required?
NSW = 0
PREC = High | Accurate
[EFIELD_PEAD = Ex EyEz]
INCAR for Born charges using LCALCEPS:
•To apply a non-analytical correction (LO/TO splitting) to the phonon frequencies,
Phonopyneeds the Born effective charges and electronic-polarisation contribution to the
macroscopic dielectric constant
•For LDA/GGA functionals, these can be computed using DFPT; for others, they need to be
computed from the response to an electric field

24. phonopy: Non-Analytical Corrections
WMD Research Day: 10
th
Oct 2014 | Slide 24
outcar-born > BORN
<Conversion Factor>
ε
xxε
xyε
xzε
yxε
yyε
yzε
zxε
zyε
zz
Z
xxZ
xyZ
xzZ
yxZ
yyZ
yzZ
zxZ
zyZ
zz
Z
xxZ
xyZ
xzZ
yxZ
yyZ
yzZ
zxZ
zyZ
zz
Dielectric tensor
Born charge tensors
for unique atoms
Sample BORN file:
•Corrections are enabled by setting NAC = .TRUE.in the configutationfile, or passing
--nacas a command-line argument
•When this option is used, phonopyexpects to find a BORN file in the working directory

25. phonopy: Non-Analytical Corrections
WMD Research Day: 10
th
Oct 2014 | Slide 25
PbTe
NAC = .FALSE. NAC = .TRUE.

26. phonopy: Force-Constant Symmetrisation
WMD Research Day: 10
th
Oct 2014 | Slide 26
FC_SYMMETRY = 0 FC_SYMMETRY = 1
•Force-constant symmetrisation is enabled by setting FC_SYMMETRY = > 0in the
configuration file
•Note that the symmetrisation is done by default in most other codes (e.g. VASP) (!)
ZnS

27. phonopy: Output Files
WMD Research Day: 10
th
Oct 2014 | Slide 27
mesh: [ m
x, m
y, m
z]
nqpoint: 32000
natom: 8
phonon:
-q-position: [ q
x, q
y, q
z]
weight: w
1
band:
-# 1
frequency: ω
1
...
Sample mesh.yamlfile:
nqpoint: 808
npath: 8
natom: 8
phonon:
-q-position: [ q
x, q
y, q
z]
distance: d
1
band:
-# 1
frequency: ω
1
...
Sample band.yamlfile:
•If EIGENVECTORS = .TRUE. is set in the configuration file, the mode eigenvectors
will also appear in these files
•With BAND_CONNECTION = .TRUE. , the frequencies for each band in band.yaml
are ordered so that they connect across the band structure

28. phonopy: Output Files
WMD Research Day: 10
th
Oct 2014 | Slide 28
# Sigma = 0.053821
-0.5372... 0.0000...
-0.5103... 0.0000...
-0.4834... 0.0000...
-0.4564... 0.0000...
-0.4295... 0.0000...
-0.4026... 0.0000...
-0.3757... 0.0000...
-0.3488... 0.0000...
-0.3219... 0.0000...
...
Sample total_dos.dat file:
unit:
temperature: K
...
natom: 8
zero_point_energy: 18.9108676
high_T_entropy: 847.3220815
thermal_properties:
-temperature: 0.0000000
-free_energy: 18.9108676
-entropy: 0.0000000
-heat_capacity: 0.0000000
-energy: 18.9108676
...
Sample thermal_properties.yamlfile:
•The “partial_dos.dat” file generated
with EIGENVECTORS = .TRUE.
contains one column for each atom in
the primitive cell

29. “Hacking” phonopy
WMD Research Day: 10
th
Oct 2014 | Slide 29
128
2
1
d
1xd
1yd
1z
F
1xF
1yF
1z
F
2xF
2yF
2z
...
2
d
2xd
2yd
2z
F
1xF
1yF
1z
F
2xF
2yF
2z
...
Sample FORCE_SETS file:
128
11
Φ
xxΦ
xyΦ
xz
Φ
yxΦ
yyΦ
yz
Φ
zxΦ
zyΦ
zz
12
Φ
xxΦ
xyΦ
xz
Φ
yxΦ
yyΦ
yz
Φ
zxΦ
zyΦ
zz
...
Sample FORCE_CONSTANTS file:

30. “Hacking” phonopy: Phonopy-QE
WMD Research Day: 10
th
Oct 2014 | Slide 300.0
5.0
10.0
15.0
20.0
25.0
-5000 5001000150020002500300035004000
Phonon DOS / AU
v/ cm
-1
PBEsol PBE+D2 PBE
SulphamerazineForm-I

31. “Hacking” phonopy: Phonopy-Tinker
WMD Research Day: 10
th
Oct 2014 | Slide 31
IRMOF-10

32. Anharmonicity1: The QHA
WMD Research Day: 10
th
Oct 2014 | Slide 32
•In the harmonic approximation, &#3627408531;
0is
temperature independent -> cannot predict
thermal expansion
•At finite temperature, the system will
minimise its free energy (&#3627408436;or &#3627408442;), as opposed
to its lattice internal energy (&#3627408456;
??????)
•This can be modelled by computing &#3627408436;(&#3627408455;)as a
function of volume, from a sequence of
harmonic phonon calculations, and
performing an EoSfit to &#3627408436;at each
temperature
•Not really “anharmonic”, but not “purely
harmonic” either -> “quasi harmonic”
•Valid up to approx. 2&#3627408455;
&#3627408474;/3

33. Anharmonicity1: The QHA
WMD Research Day: 10
th
Oct 2014 | Slide 33
For a solid:&#3627408442;&#3627408455;=&#3627408456;
??????&#3627408457;+&#3627408456;
??????&#3627408455;,&#3627408457;+&#3627408477;&#3627408457;−&#3627408455;&#3627408454;
??????(&#3627408455;,&#3627408457;)
Gibbs energy:&#3627408442;&#3627408455;=&#3627408443;−&#3627408455;&#3627408454;=&#3627408456;+&#3627408477;&#3627408457;−&#3627408455;&#3627408454;
Within the QHA:&#3627408442;&#3627408455;,&#3627408477;=min
??????
[&#3627408436;&#3627408457;,&#3627408455;+&#3627408477;&#3627408457;]
Derived properties:&#3627408457;&#3627408455;,&#3627408437;(&#3627408455;)
α
??????(&#3627408455;)
-> From EoSfits
->
1
??????
????????????
????????????
??????
&#3627408454;(&#3627408455;) -> −
??????&#3627408442;
????????????
??????
&#3627408438;
??????(&#3627408455;) -> −
??????&#3627408443;
????????????
??????
Once we have &#3627408457;(&#3627408455;), we can compute any property for which the temperature
dependence is captured (to first approximation) by volume changes
In principle, QHA can
also model &#3627408477;dependence

34. phonopy-qha: Workflow
WMD Research Day: 10
th
Oct 2014 | Slide 34
EoSCurve
Harmonic
Phonopy
Thermal
Properties
Post-Process
phonopy-d --dim=“...”
phonopy-f vasprun-{001..XXX}.xml
phonopy-t -p -s Settings.conf
E/Vcurve -> “e-v.dat”
phonopy-qha
e-v.dat
thermal_properties-{001..XXX}.yaml
--tmax=980
[--tstep=10]
[--pressure=<p/GPa>]
[-p] [-s]
Needs to be two points smaller
than the maximum temperature
in the YAML files

35. phonopy-qha: Output
WMD Research Day: 10
th
Oct 2014 | Slide 35

36. phonopy-qha: Output
WMD Research Day: 10
th
Oct 2014 | Slide 36
bulk_modulus-temperature.dat
Cp-temperature.dat
Cp-temperature_polyfit.dat
Cv-volume.dat
dsdv-temperature.dat
entropy-volume.dat
gibbs-temperature.dat
gruneisen-temperature.dat
helmholtz-volume.dat
thermal_expansion.dat
volume_expansion.dat
volume-temperature.dat
<-Bis temperature dependent(!)
<-C
Vat each volume, at each temperature
<-S
Vat each volume, at each temperature
<-Average Gruneisenparameter (?)
<-Aat each volume, at each temperature
<-&#3627409148;
??????&#3627408455;(Volumetric)
<-&#3627409148;
??????&#3627408455;(Linear; ∆??????/??????
0, with ??????=
3
&#3627408457;)

37. phonopy-QHA: Examples
WMD Research Day: 10
th
Oct 2014 | Slide 37
PbTe
PbTe: 0 K -600 K
PbTe: 0 K -600 K

38. phonopy-QHA: Examples
WMD Research Day: 10
th
Oct 2014 | Slide 38
Functionalα
L/ 10
-6
K
-1
α
V/ 10
-6
K
-1
LDA 19.66 58.99
PW91 34.99 104.98
PBE 39.51 118.56
PBEsol 22.37 67.11
TPSS 29.04 87.13
revTPSS 25.85 77.56
Exp: α
L= 19.8/20.4 10
-6
K
-1
PbTe

39. phonopy-QHA: Examples
WMD Research Day: 10
th
Oct 2014 | Slide 39
SnS: 10 K -350 K
PbTe: 0 K -600 K

40. Anharmonicity2: Phonon-Phonon Coupling
WMD Research Day: 10
th
Oct 2014 | Slide 40
•The harmonic approximation treats phonons as being independent oscillators; this
means the phonons have infinite lifetimes
•In a real system, phonon-phonon coupling leads to scattering and a finite phonon
lifetime; this anharmoniceffect is required to model heat transport
&#3627409173;&#3627408467;??????,&#3627409158;
&#3627409173;??????
&#3627408479;=
&#3627409173;&#3627408467;??????,&#3627409158;
&#3627409173;??????
&#3627408479;
&#3627408465;&#3627408470;&#3627408467;&#3627408467;+
&#3627409173;&#3627408467;??????,&#3627409158;
&#3627409173;??????
&#3627408479;
&#3627408466;&#3627408485;&#3627408481;+
&#3627409173;&#3627408467;??????,&#3627409158;
&#3627409173;??????
&#3627408479;
&#3627408480;&#3627408464;??????&#3627408481;&#3627408481;=0
At equilibrium, the change in
occupation probability due to…
…diffusion, … …the external
heat current, ……and scattering ……must balance

41. Anharmonicity2: Phonon-Phonon Coupling
WMD Research Day: 10
th
Oct 2014 | Slide 41
&#3627409173;&#3627408467;??????,&#3627409158;
&#3627409173;??????
&#3627408531;=
&#3627409173;&#3627408467;??????,&#3627409158;
&#3627409173;??????
&#3627408531;
&#3627408465;&#3627408470;&#3627408467;&#3627408467;+
&#3627409173;&#3627408467;??????,&#3627409158;
&#3627409173;??????
&#3627408531;
&#3627408466;&#3627408485;&#3627408481;+
&#3627409173;&#3627408467;??????,&#3627409158;
&#3627409173;??????
&#3627408531;
&#3627408480;&#3627408464;??????&#3627408481;&#3627408481;=0BTE:
RTA:
&#3627409173;&#3627408467;??????,&#3627409158;
&#3627409173;??????
&#3627408531;
&#3627408480;&#3627408464;??????&#3627408481;&#3627408481;=
&#3627408467;??????,&#3627409158;−&#3627408467;
0
??????,&#3627409158;
τ(??????,&#3627409158;)
&#3627408467;??????,&#3627409158;−&#3627408467;
0
??????,&#3627409158;=−??????(??????,&#3627409158;)
&#3627409173;&#3627408467;
0
??????,&#3627409158;
&#3627409173;&#3627408455;
&#3627409147;&#3627408455;τ(??????,&#3627409158;)
Where: &#3627408467;??????,&#3627409158;
&#3627408467;
0
??????,&#3627409158;
????????????,&#3627409158;
??????(??????,&#3627409158;)
Mode occupation probability
Initial mode occupation probability
Mode group velocity
Mode relaxation time
Scattering assumed to be
related to phonon lifetimes

42. Anharmonicity2: Phonon-Phonon Coupling
WMD Research Day: 10
th
Oct 2014 | Slide 42
Phonon linewidths: ????????????,&#3627409158;=
1
Γ(??????,&#3627409158;)
Phonon-phonon interactions: Φ
&#3627409148;&#3627409149;&#3627409150;&#3627408470;,&#3627408471;,&#3627408472;=
&#3627409173;
3
&#3627408440;
&#3627409173;&#3627408479;
&#3627408470;,&#3627409148;&#3627409173;&#3627408479;
&#3627408471;,&#3627409149;&#3627409173;&#3627408479;
&#3627408472;,&#3627409150;
≈
−&#3627408441;
&#3627408470;,&#3627409148;
∆&#3627408479;
&#3627408471;,&#3627409149;∆&#3627408479;
&#3627408472;,&#3627409150;
…and finally: &#3627409157;
??????=
??????,??????
&#3627409172;??????,&#3627409158;
&#3627409173;&#3627408467;
0
??????,&#3627409158;
&#3627409173;&#3627408455;
????????????,&#3627409158;⨂????????????,&#3627409158;????????????,&#3627409158;
&#3627409157;
??????=
??????,??????
&#3627408438;
????????????,&#3627409158;????????????,&#3627409158;⨂????????????,&#3627409158;????????????,&#3627409158;
•In practice, the calculation involves a large set of single-point calculations to determine
Φ
&#3627409148;&#3627409149;&#3627409150;&#3627408470;,&#3627408471;,&#3627408472;, followed by a computationally-heavy post-processing step to get &#3627409157;
??????
•Worth noting that the BTE-RTA method may benefit from a cancellation of errors, and as
such works well despite ignoring various effects
Tensor product

43. phono3py: Workflow
WMD Research Day: 10
th
Oct 2014 | Slide 43
Input Structure
Create
Displacements
Calculate Forces
Extract Forces
Post-Process
phono3py -d --dim=“2 2 2”
phono3py --cf3 vasprun-{001..XXX}.xml
phono3py --tmax=1000 --tstep=10 --fc2 --fc3
--dim=“2 2 2” -v --mesh=“24 24 24”
--br--sigma=“0.1”
&#3627409157;
??????is driven by heat transport through acoustic
phonons; these have zero frequency at Γ, so
off-Γq-points are a must
phono3py --dim=“2 2 2”

44. phono3py: Setup
WMD Research Day: 10
th
Oct 2014 | Slide 44
phono3py -d --dim=“...”
[--cutpair=4]
[--dim2=“4 4 4”]/--fc2_extra=“4 4 4”]
phono3py --cf3 vasprun-{001..XXX}.xml
phono3py --dim=“...”
[--dim2...]
Set a pair-cutoffdistance for Φ
&#3627409148;&#3627409149;&#3627409150;&#3627408470;,&#3627408471;,&#3627408472;;
reduces number of displaced structures
Separate supercellsizes for calculating
Φ
&#3627409148;&#3627409149;&#3627408470;,&#3627408471;and Φ
&#3627409148;&#3627409149;&#3627409150;&#3627408470;,&#3627408471;,&#3627408472;
Extract forces from VASP output
Precalculate Φ
&#3627409148;&#3627409149;&#3627408470;,&#3627408471;and Φ
&#3627409148;&#3627409149;&#3627409150;&#3627408470;,&#3627408471;,&#3627408472;
to save time during post-processing
The --cutpairtag uses the same numbering for the displaced POSCAR files
as the full calculation; this means the cutoffcan be increased, and the extra
displacements added, systematically, to converge w.r.t. the interaction range

45. phono3py: Post Processing
WMD Research Day: 10
th
Oct 2014 | Slide 45
phono3py --tmax=1000 --tstep=10 --fc2 --fc3
--dim=“...” -v --mesh=“24 24 24”
--br--sigma=“0.1” [--nac]
[--dim2=“...”/--fc2_extra=“...”]
•The post-processing (mainly the phonon-lifetime calculations) takes a verylong time for
large supercells/large or low-symmetry structures
•It is possible to run the calculation on (ranges of) q-points separately, and then combine
them afterwards
•Various post-processing tags can be applied, e.g. to incorporate isotope effects
Read in pre-calculated
force constants
# of interactions per q-point becomes larger with mesh
size; cannot easily “max out” as for DOS calculations,
but needs to be converged

46. phono3py: Examples
WMD Research Day: 10
th
Oct 2014 | Slide 460
3
6
9
12
15
25125225325425525
κ
/ w m
-
1
K
-
1
T/ K 0
3
6
9
12
15
25125225325425525
κ
/ w m
-
1
K
-
1
T/ K 0
3
6
9
12
15
25125225325425525
κ
/ w m
-
1
K
-
1
T/ K
PbS PbSe
PbTe
Increasing V

47. Summary: So, What Can PhonopyDo?
WMD Research Day: 10
th
Oct 2014 | Slide 48
“Anharmonicity” Properties
phonopy Phonon (P)DOS,phonon dispersion, &#3627408436;(&#3627408455;), &#3627408456;
??????(&#3627408455;),
&#3627408454;
??????(&#3627408455;), &#3627408438;
??????&#3627408455;, phononproperties such as mode
eigenvectors, ??????(??????,&#3627409158;), IRs, thermal MSDs,etc.
phonopy-gruneisen ??????(??????,&#3627409158;)(mode Gruneisenparameters)
phonopy-qha &#3627408437;(&#3627408455;[,&#3627408477;]), &#3627408457;(&#3627408455;[,&#3627408477;]), &#3627408442;(&#3627408455;[,&#3627408477;]), &#3627408438;
??????(&#3627408455;[;&#3627408477;]), ??????(&#3627408455;[,&#3627408477;]), any
property where the &#3627408455;[,&#3627408477;]dependencecan be modelled
by changes in the lattice parameters
phono3py κ
??????(&#3627408455;); alsopossible toextract related quantities such as
Γ(??????,&#3627409158;),τ(??????,&#3627409158;), &#3627408438;
??????(??????,&#3627409158;), ??????(??????,&#3627409158;)

48. A Few Closing Remarks
WMD Research Day: 10
th
Oct 2014 | Slide 49
•In my experience, a well-chosen GGA functional (e.g. PBEsolfor bulk materials) gives
accurate forces
Provided tight convergence criterion are used, phonon frequencies and
thermodynamic properties show good agreement with experiment
•Using the QHA is a bit more expensive, but in return yields a lot of properties
Model the temperature dependence of properties without e.g. resorting to MD
averaging (although this certainly does have its merits)
[Cynical] Should end up with enough data for a decent PRB…
•phono3py produces very good values for &#3627409157;
??????, although it can be very expensive (“the GW
of lattice dynamics”?)
New code; still need to test various aspects of its functionality
Not many people doing this type of calculation at the moment, either with
phono3py or ShengBTE