Phonon calculations¶
Module for calculating vibrational normal modes for periodic systems using the socalled small displacement method (see e.g. [Alfe]). So far, spacegroup symmetries are not exploited to reduce the number of atomic displacements that must be calculated and subsequent symmetrization of the force constants.
For polar materials the dynamical matrix at the zone center acquires a nonanalytical contribution that accounts for the LOTO splitting. This contribution requires additional functionality to evaluate and is not included in the present implementation. Its implementation in conjunction with the small displacement method is described in [Wang].
Example¶
Simple example showing how to calculate the phonon dispersion for bulk aluminum using a 7x7x7 supercell within effective medium theory:
from ase.build import bulk
from ase.calculators.emt import EMT
from ase.dft.kpoints import ibz_points, bandpath
from ase.phonons import Phonons
# Setup crystal and EMT calculator
atoms = bulk('Al', 'fcc', a=4.05)
calc = EMT()
# Phonon calculator
N = 7
ph = Phonons(atoms, calc, supercell=(N, N, N), delta=0.05)
ph.run()
# Read forces and assemble the dynamical matrix
ph.read(acoustic=True)
# Highsymmetry points in the Brillouin zone
points = ibz_points['fcc']
G = points['Gamma']
X = points['X']
W = points['W']
K = points['K']
L = points['L']
U = points['U']
point_names = ['$\Gamma$', 'X', 'U', 'L', '$\Gamma$', 'K']
path = [G, X, U, L, G, K]
# Band structure in meV
path_kc, q, Q = bandpath(path, atoms.cell, 100)
omega_kn = 1000 * ph.band_structure(path_kc)
# Calculate phonon DOS
omega_e, dos_e = ph.dos(kpts=(50, 50, 50), npts=5000, delta=5e4)
omega_e *= 1000
# Plot the band structure and DOS
import matplotlib.pyplot as plt
plt.figure(1, (8, 6))
plt.axes([.1, .07, .67, .85])
for n in range(len(omega_kn[0])):
omega_n = omega_kn[:, n]
plt.plot(q, omega_n, 'k', lw=2)
plt.xticks(Q, point_names, fontsize=18)
plt.yticks(fontsize=18)
plt.xlim(q[0], q[1])
plt.ylabel("Frequency ($\mathrm{meV}$)", fontsize=22)
plt.grid('on')
plt.axes([.8, .07, .17, .85])
plt.fill_between(dos_e, omega_e, y2=0, color='lightgrey', edgecolor='k', lw=1)
plt.ylim(0, 35)
plt.xticks([], [])
plt.yticks([], [])
plt.xlabel("DOS", fontsize=18)
plt.show()
Mode inspection:
# Write modes for specific qvector to trajectory files
ph.write_modes([l/2 for l in L], branches=[2], repeat=(8, 8, 8), kT=3e4)
[Alfe]  D. Alfe, PHON: A program to calculate phonons using the small displacement method, Comput. Phys. Commun. 180, 2622 (2009) 
[Wang]  Y. Wang et al., A mixedspace approach to firstprinciples calculations of phonon frequencies for polar materials, J. Phys.: Cond. Matter 22, 202201 (2010) 
List of all Methods¶

class
ase.phonons.
Phonons
(*args, **kwargs)[source]¶ Class for calculating phonon modes using the finite displacement method.
The matrix of force constants is calculated from the finite difference approximation to the firstorder derivative of the atomic forces as:
2 nbj nbj nbj d E F  F+ C =  ~  , mai dR dR 2 * delta mai nbj
where F+/F denotes the force in direction j on atom nb when atom ma is displaced in direction +i/i. The force constants are related by various symmetry relations. From the definition of the force constants it must be symmetric in the three indices mai:
nbj mai bj ai C = C > C (R ) = C (R ) . mai nbj ai n bj n
As the force constants can only depend on the difference between the m and n indices, this symmetry is more conveniently expressed as shown on the right handside.
The acoustic sumrule:
_ _ aj \ bj C (R ) =  ) C (R ) ai 0 /__ ai m (m, b) != (0, a)
Ordering of the unit cells illustrated here for a 1dimensional system (in case
refcell=None
in constructor!):m = 0 m = 1 m = 2 m = 1        * b  *  *  *        * a  *  *  *       
Example:
>>> from ase.build import bulk >>> from ase.phonons import Phonons >>> from gpaw import GPAW, FermiDirac >>> atoms = bulk('Si', 'diamond', a=5.4) >>> calc = GPAW(kpts=(5, 5, 5), h=0.2, occupations=FermiDirac(0.)) >>> ph = Phonons(atoms, calc, supercell=(5, 5, 5)) >>> ph.run() >>> ph.read(method='frederiksen', acoustic=True)
Initialize with base class args and kwargs.

apply_cutoff
(D_N, r_c)[source]¶ Zero elements for interatomic distances larger than the cutoff.
Parameters:
 D_N: ndarray
 Dynamical/force constant matrix.
 r_c: float
 Cutoff in Angstrom.

band_structure
(path_kc, modes=False, born=False, verbose=True)[source]¶ Calculate phonon dispersion along a path in the Brillouin zone.
The dynamical matrix at arbitrary qvectors is obtained by Fourier transforming the realspace force constants. In case of negative eigenvalues (squared frequency), the corresponding negative frequency is returned.
Frequencies and modes are in units of eV and Ang/sqrt(amu), respectively.
Parameters:
 path_kc: ndarray
 List of kpoint coordinates (in units of the reciprocal lattice vectors) specifying the path in the Brillouin zone for which the dynamical matrix will be calculated.
 modes: bool
 Returns both frequencies and modes when True.
 born: bool
 Include nonanalytic part given by the Born effective charges and the static part of the highfrequency dielectric tensor. This contribution to the force constant accounts for the splitting between the LO and TO branches for q > 0.
 verbose: bool
 Print warnings when imaginary frequncies are detected.

dos
(kpts=(10, 10, 10), npts=1000, delta=0.001, indices=None, verbose=True)[source]¶ Calculate phonon dos as a function of energy.
Parameters:
 qpts: tuple
 Shape of MonkhorstPack grid for sampling the Brillouin zone.
 npts: int
 Number of energy points.
 delta: float
 Broadening of Lorentzian lineshape in eV.
 indices: list
 If indices is not None, the atomicpartial dos for the specified atoms will be calculated.
 verbose: bool
 Print warnings when imaginary frequncies are detected.

read
(method=’Frederiksen’, symmetrize=3, acoustic=True, cutoff=None, born=False, **kwargs)[source]¶ Read forces from pickle files and calculate force constants.
Extra keyword arguments will be passed to
read_born_charges
.Parameters:
 method: str
 Specify method for evaluating the atomic forces.
 symmetrize: int
 Symmetrize force constants (see doc string at top) when
symmetrize != 0
(default: 3). Since restoring the acoustic sum rule breaks the symmetry, the symmetrization must be repeated a few times until the changes a insignificant. The integer gives the number of iterations that will be carried out.  acoustic: bool
 Restore the acoustic sum rule on the force constants.
 cutoff: None or float
 Zero elements in the dynamical matrix between atoms with an interatomic distance larger than the cutoff.
 born: bool
 Read in Born effective charge tensor and highfrequency static dielelctric tensor from file.

read_born_charges
(name=None, neutrality=True)[source]¶ Read Born charges and dieletric tensor from pickle file.
The charge neutrality sumrule:
_ _ \ a ) Z = 0 /__ ij a
Parameters:
 neutrality: bool
 Restore charge neutrality condition on calculated Born effective charges.

write_modes
(q_c, branches=0, kT=0.02585199101165164, born=False, repeat=(1, 1, 1), nimages=30, center=False)[source]¶ Write modes to trajectory file.
Parameters:
 q_c: ndarray
 qvector of the modes.
 branches: int or list
 Branch index of modes.
 kT: float
 Temperature in units of eV. Determines the amplitude of the atomic displacements in the modes.
 born: bool
 Include nonanalytic contribution to the force constants at q > 0.
 repeat: tuple
 Repeat atoms (l, m, n) times in the directions of the lattice vectors. Displacements of atoms in repeated cells carry a Bloch phase factor given by the qvector and the cell lattice vector R_m.
 nimages: int
 Number of images in an oscillation.
 center: bool
 Center atoms in unit cell if True (default: False).
