Phonon calculations

Phonon calculations#

Module for calculating vibrational normal modes for periodic systems using the so-called small displacement method (see e.g. [Alfe]). So far, space-group symmetries are not exploited to reduce the number of atomic displacements that must be calculated and subsequent symmetrization of the force constants.

A key step in the calculation of phonon modes and energies through small displacements is the construction of the \(\mathbf{q}\)-dependent dynamical matrix from harmonic force constants and phase factors. These phase factors are determined by the input \(\mathbf{q}\)-vector and the atom-to-atom vectors which go from an atom \(a_1\) in the reference unit cell to an atom \(a_2\) in cell \(N\) in the supercell; however, due to the periodicity of the supercell, the shortest vector from \(a_1\) to atom \((N, a_2)\) is not necessarily the obvious displacement vector between atoms in the supercell. Instead, the shortest vector can be found by checking all the nearest-neighbour periodic images of \((N, a_2)\) and choosing the shortest vector. Furthermore, there may be multiple images of \((N, a_2)\) that yield the shortest vector, and in such cases, it will be beneficial to use the averaged phase factor of all the shortest atom-to-atom image vectors. For further details on mean minimum image phases see equations (47) and (48) and surrounding text in [Togo]. In the ASE Phonons class, the usage of mean minimum image phases can be enabled by setting use_mean_minimum_images=True at class initialization. This is highly recommended and tends to have much better converged phonon energies for smaller supercells. The tolerance for detecting mean minimum images can be adjusted with the minimum_image_tol keyword.

For polar materials the dynamical matrix at the zone center acquires a non-analytical contribution that accounts for the LO-TO 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].

List of all methods#

class ase.phonons.Phonons(atoms: Atoms, calc=None, supercell: tuple[int, int, int] = (1, 1, 1), select_atoms: list[int] | list[str] | None = None, name: str = 'phonon', delta: float = 0.01, center_refcell: bool = False, use_mean_minimum_images: bool | None = None, minimum_image_tol: float = 1e-05, comm=None)[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 first-order 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 hand-side.

The acoustic sum-rule:

            _ _
 aj         \    bj
C  (R ) = -  )  C  (R )
 ai  0      /__  ai  m
           (m, b)
             !=
           (0, a)

Ordering of the unit cells illustrated here for a 1-dimensional system (in case center_refcell=False in constructor!):

    m = 0        m = 1        m = -2        m = -1
-----------------------------------------------------
|            |            |            |            |
|        * b |        *   |        *   |        *   |
|            |            |            |            |
|   * a      |   *        |   *        |   *        |
|            |            |            |            |
-----------------------------------------------------

Examples

>>> from ase.build import bulk
>>> from ase.phonons import Phonons
>>> from ase.calculators.emt import EMT

>>> atoms = bulk('Al', 'fcc', a=4.05)
>>> calc = EMT()
>>> ph = Phonons(
...     atoms,
...     calc,
...     supercell=(7, 7, 7),
...     delta=0.05,
...     use_mean_minimum_images=True
... )
>>> ph.run()
>>> ph.read()
>>> path = atoms.cell.bandpath('GXULGK', npoints=100)
>>> bs = ph.get_band_structure(path, verbose=False)
>>> dos = ph.get_dos(kpts=(20, 20, 20)).sample_grid(npts=100, width=1e-3)

Init with an instance of Atoms and a calculator.

Parameters:

  • atoms – The atoms to work on.

  • calc – Calculator for the supercell calculation.

  • supercell – Size of supercell given by the number of repetitions (n1, n2, n3) of the small unit cell in each direction.

  • select_atoms – Select which atoms to generate displacements for. By default, displacements are generated for all atoms.

  • name – Base name to use for files.

  • delta – Magnitude of displacement in Ang.

  • center_refcell – Reference cell in which the atoms will be displaced. If False, then corner cell in supercell is used. If True, then cell in the center of the supercell is used.

  • use_mean_minimum_images – Use averaged phase factor of all the shortest atom-to-atom image vectors. This will be the default in the future and is strongly recommended as it requires smaller supercells to converge the phonon energies.

  • minimum_image_tol – Tolerance (in Å) for finding equivalent atoms when using mean minimum image phase factors.

  • comm – MPI communicator for the phonon calculation. Default is to use world.

acoustic(C_avNav: ndarray) → None[source]#

Restore acoustic sumrule on force constants.

apply_cutoff(C_avNav: ndarray, r_c: float) → None[source]#

Zero elements for interatomic distances larger than the cutoff.

Parameters:

  • C_avNav (ndarray) – Matrix of force constants.

  • r_c (float) – Cutoff radius in Ångstrom.

band_structure(path_kc, modes=False, verbose=True)[source]#

Calculate phonon dispersion along a path in the Brillouin zone.

The dynamical matrix at arbitrary q-vectors is obtained by Fourier transforming the real-space force constants. In case of negative eigenvalues (squared frequency), the corresponding negative frequency is returned.

Frequencies and modes are in units of eV and 1/sqrt(amu), respectively.

Parameters:

  • path_kc (ndarray) – List of k-point 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.

  • verbose (bool) – Print warnings when imaginary frequncies are detected.

Returns:

  • omega_kl (np.ndarray) – Phonon band energies.

  • u_klav (np.ndarray) – Eigenvectors (only when modes is True).

calculate_dynamical_matrix(q_c: ndarray, D_avNav: ndarray) → ndarray[source]#

Computation of the dynamical matrix in momentum space D_XX(q). This is a Fourier transform from real-space force constants for a given momentum vector q.

If the phonons class was initialized with use_mean_minimum_images=True, the q-dependent phase factor for a given set of atoms a1 and (N, a2) is averaged over all the shortest atom-to-atom vectors from atom1 to a periodic image of atom (N, a2).

Parameters:

  • q_c – Phonon wave vector in scaled coordinates.

  • D_avNav – Force constant matrix scaled by atom masses.

Returns:

D_XX – Dynamical matrix, a complex-valued array D(q) with shape=(3 * natoms, 3 * natoms).

Return type:

np.ndarray

check_eq_forces()[source]#

Check maximum size of forces in the equilibrium structure.

compute_dynamical_matrix(q_scaled: ndarray, D_N: ndarray) → ndarray[source]#

Computation of the dynamical matrix in momentum space D_ab(q). This is a Fourier transform from real-space dynamical matrix D_N for a given momentum vector q.

Note

Deprecated. Please use calculate_dynamical_matrix().

Parameters:

  • q_scaled (np.ndarray) – q vector in scaled coordinates.

  • D_N (np.ndarray) – Dynamical matrix in real-space.

Returns:

D_q – 2D complex-valued array D(q) with shape=(3 * natoms, 3 * natoms).

Return type:

np.ndarray

get_band_structure(path, modes: bool = False, verbose: bool = True)[source]#

Calculate and return the phonon band structure.

This method computes the phonon band structure for a given path in reciprocal space. It is a wrapper around the internal band_structure() method of the Phonons class. The method can optionally calculate and return phonon modes.

Frequencies and modes are in units of eV and \(1/\sqrt{\mathrm{amu}}\), respectively.

Parameters:

  • path (BandPath object) – The BandPath object defining the path in the reciprocal space over which the phonon band structure is calculated.

  • modes (bool, optional) – If True, phonon modes will also be calculated and returned. Defaults to False.

  • verbose (bool, optional) – If True, enables verbose output during the calculation. Defaults to True.

Returns:

If modes is False, returns a BandStructure object containing the phonon band structure. If modes is True, returns a tuple, where the first element is the BandStructure object and the second element is an ndarray of phonon modes.

If modes are returned, the array is of shape (k-point, bands, atoms, 3) and the norm-squared of the mode is 1 / m_{eff}, where m_{eff} is the effective mass of the mode.

Return type:

BandStructure or tuple of (BandStructure, ndarray)

get_dos(kpts: tuple[int, int, int] = (10, 10, 10), indices: list | None = None, verbose: bool = True)[source]#

Return a phonon density of states.

Parameters:

  • kpts (tuple) – Shape of Monkhorst-Pack grid for sampling the Brillouin zone.

  • indices (list) – If indices is not None, the amplitude-weighted atomic-partial DOS for the specified atoms will be calculated.

  • verbose (bool) – Print warnings when imaginary frequncies are detected.

Returns:

Density of states.

Return type:

RawDOSData

get_force_constant() → ndarray[source]#

Deprecated version of get_force_constants, can safely be removed by summer 2027

Note

Please use get_force_constants().

get_force_constants() → ndarray[source]#

Return matrix of force constants.

read(method: str = 'Frederiksen', symmetrize: int = 3, acoustic: bool = True, cutoff: float | None = None)[source]#

Read forces from json files and calculate force constants.

Parameters:

  • method – Specify method for evaluating the atomic forces, method=’Frederiksen’ imposes momentum conservation.

  • symmetrize – 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 – Restore the acoustic sum rule on the force constants.

  • cutoff – Zero elements in the dynamical matrix between atoms with an interatomic distance larger than the cutoff.

read_born_charges(name='born', neutrality=True)[source]#

Read Born charges and dieletric tensor from JSON file.

The charge neutrality sum-rule:

_ _
\    a
 )  Z   = 0
/__  ij
 a

Deprecated since version 3.22.1: Current implementation of non-analytical correction is likely incorrect, see issue: #941

Parameters:

  • neutrality (bool) – Restore charge neutrality condition on calculated Born effective charges.

  • name (str) – Key used to identify the file with Born charges for the unit cell in the JSON cache.

symmetrize(C_avNav: ndarray) → ndarray[source]#

Symmetrize force constant matrix.

write_modes(q_c, branches=0, kT: float = 0.02585199101165164, repeat: tuple[int, int, int] = (1, 1, 1), nimages: int = 30, center: bool = False) → None[source]#

Write modes to trajectory file.

Note

To exaggerate the amplitudes for better visualization, multiply kT by the square of the desired factor.

Parameters:

  • q_c (ndarray of shape (3,)) – q-vector 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.

  • 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 q-vector 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).

References#

[Alfe]

D. Alfe, PHON: A program to calculate phonons using the small displacement method, Comput. Phys. Commun. 180, 2622 (2009)

[Togo]

A. Togo, et al. Implementation strategies in phonopy and phono3py, Journal of Physics: Condensed Matter 35.35 (2023): 353001.

[Wang]

Y. Wang et al., A mixed-space approach to first-principles calculations of phonon frequencies for polar materials, J. Phys.: Cond. Matter 22, 202201 (2010)