QpointPhononModes

The QpointPhononModes object contains precalculated phonon frequencies and eigenvectors at certain q-points.

Reading From CASTEP
Reading From Phonopy
From Force Constants
Reordering frequencies
Plotting Dispersion
Calculating The Coherent Neutron Structure Factor
Calculating The Debye-Waller Exponent
Calculating Partial and Neutron-Weighted Density of States
Calculating Density of States
Docstring

Reading From CASTEP 

Phonon frequencies and eigenvectors can be read from a .phonon file using QpointPhononModes.from_castep.

from euphonic import QpointPhononModes

filename = 'quartz.phonon'
phonons = QpointPhononModes.from_castep(filename)

Reading From Phonopy 

Phonopy should be run with the --eigvecs flag, or EIGENVECTORS = .TRUE. to enable creation of a Euphonic QpointPhononModes object.

Using QpointPhononModes.from_phonopy Euphonic can read frequencies and eigenvectors from Phonopy files with the following default names:

mesh.yaml/mesh.hdf5
qpoints.yaml/qpoints.hdf5
bands.yaml/bands.hdf5

The file to be read can be specified with the phonon_name argument. Some of these files do not include the crystal information, so it must be read from a phonopy.yaml file, which can be specified with the summary_name argument. A path can also be specified.

from euphonic import QpointPhononModes

phonons = QpointPhononModes.from_phonopy(path='NaCl', phonon_name='mesh.hdf5')

From Force Constants 

See Force Constants

Reordering frequencies 

The stored frequencies can be reordered by comparing eigenvectors using QpointPhononModes.reorder_frequencies. This means that the same mode will have the same index across different q-points, so will be plotted as the same colour in a dispersion plot, allowing modes to be followed across q-space (see Plotting)

from euphonic import QpointPhononModes

phonons = QpointPhononModes.from_castep('quartz.phonon')
phonons.reorder_frequencies()

Plotting Dispersion 

See Plotting Dispersion

Calculating The Coherent Neutron Structure Factor 

The neutron structure factor can be calculated for each branch and q-point using the QpointPhononModes.calculate_structure_factor method, which returns a StructureFactor object. (See the docstring for algorithm details.)

Scattering lengths 

The coherent scattering length is a physical property of the nucleus. To provide this data explicitly, the scattering_lengths argument can be set as a dictionary mapping each atom identity to a pint.Quantity (see Units for details).

Alternatively, this argument may be a string referring to a data file with the coherent_scattering_length property. (See Reference Data for details.) By default, QpointPhononModes.calculate_structure_factor will use the "Sears1992" data set included in Euphonic. If you have a custom data file, this can be used instead.

Debye-Waller factor 

Inclusion of the Debye-Waller factor is optional, and can be provided in the dw keyword argument, see Calculating The Debye-Waller Exponent.

Example 

The following example shows a full calculation from the force constants to the structure factor with Debye-Waller:

import seekpath
import numpy as np
from euphonic import ureg, QpointPhononModes, ForceConstants
from euphonic.util import mp_grid

# Read the force constants
fc = ForceConstants.from_castep('quartz.castep_bin')

# Generate a recommended q-point path to calculate the structure factor on
# using seekpath
cell = fc.crystal.to_spglib_cell()
qpts = seekpath.get_explicit_k_path(cell)["explicit_kpoints_rel"]
# Calculate frequencies/eigenvectors for the q-point path
phonons = fc.calculate_qpoint_phonon_modes(qpts, asr='reciprocal')

# For the Debye-Waller calculation, generate and calculate
# frequencies/eigenvectors on a grid (generate a Monkhorst-Pack grid of
# q-points using the mp-grid helper function)
q_grid = mp_grid([5,5,5])
phonons_grid = fc.calculate_qpoint_phonon_modes(q_grid, asr='reciprocal')
# Now calculate the Debye-Waller exponent
temperature = 5*ureg('K')
dw = phonons_grid.calculate_debye_waller(temperature)

# Calculate the structure factor for each q-point in phonons. A
# StructureFactor object is returned
fm = ureg('fm')
scattering_lengths = {'Si': 4.1491*fm, 'O': 5.803*fm}
sf = phonons.calculate_structure_factor(scattering_lengths, dw=dw)

Calculating The Debye-Waller Exponent 

The Debye-Waller factor is an optional part of the structure factor calculation. The exponent part of the Debye-Waller factor is independent of Q and should be precalculated using QpointPhononModes.calculate_debye_waller (see the docstring for algorithm details). This requires a QpointPhononModes object calculated on a grid of q-points and a temperature, and returns a DebyeWaller object. The Debye-Waller exponent can be calculated by:

from euphonic import ureg, QpointPhononModes

phonons = QpointPhononModes.from_castep('quartz-grid.phonon')
temperature = 5*ureg('K')
dw = phonons.calculate_debye_waller(temperature)

Calculating Partial and Neutron-Weighted Density of States 

Partial and neutron-weighted density of states can be calculated using QpointPhononModes.calculate_pdos. Like DOS, this requires an array of energy bin edges in energy units. This returns a Spectrum1DCollection object, containing the atom-resolved partial density of states. If weighting is supplied, or the cross_sections are specified (this can be done in the same way as the scattering lengths) the neutron-weighted partial density of states is returned. This is also resolved per-atom, and units are area/energy per atom of sample.

The Spectrum1DCollection.metadata attribute labels each PDOS spectrum by species and index, and Spectrum1DCollection methods can be used to obtain per-species or total PDOS. For example, to produce both total and species-resolved coherent neutron-weighted PDOS:

import numpy as np
from euphonic import ureg, QpointPhononModes

phonons = QpointPhononModes.from_castep('quartz-grid.phonon')

ebins = np.arange(0, 165, 0.1)*ureg('meV')
pdos = phonons.calculate_pdos(ebins, weighting='coherent')  # atom resolved pdos

species_pdos = pdos.group_by('species')  # species resolved pdos
total_dos = pdos.sum()  # total dos

PDOS can also be adaptively broadened in the same way as DOS.

Calculating Density of States 

See Calculating DOS

Docstring 

class QpointPhononModes(crystal, qpts, frequencies, eigenvectors, weights=None)

A class to read and store vibrational data from model (e.g. CASTEP) output files

Variables

crystal – Lattice and atom information
n_qpts – Number of q-points in the object
qpts – Shape (n_qpts, 3) float ndarray. Q-point coordinates, in fractional coordinates of the reciprocal lattice
weights – Shape (n_qpts,) float ndarray. The weight for each q-point, for Brillouin Zone integration over symmetry-reduced q-points
frequencies – Shape (n_qpts, 3*crystal.n_atoms) float Quantity. Phonon frequencies per q-point and mode
eigenvectors – Shape (n_qpts, 3*crystal.n_atoms, crystal.n_atoms, 3) complex ndarray. The dynamical matrix eigenvectors in Cartesian coordinates, using the same Cartesian basis as the cell_vectors in the crystal object

__init__(crystal, qpts, frequencies, eigenvectors, weights=None)

Parameters

crystal (Crystal) – Lattice and atom information
qpts (ndarray) – Shape (n_qpts, 3) float ndarray. Q-point coordinates, in fractional coordinates of the reciprocal lattice
frequencies (Quantity) – Shape (n_qpts, 3*crystal.n_atoms) float Quantity. Phonon frequencies per q-point and mode
eigenvectors (ndarray) – Shape (n_qpts, 3*crystal.n_atoms, crystal.n_atoms, 3) complax ndarray. The dynamical matrix eigenvectors in Cartesian coordinates, using the same Cartesian basis as the cell_vectors in the crystal object
weights (Optional[ndarray]) – Shape (n_qpts,) float ndarray. The weight for each q-point, for Brillouin Zone integration over symmetry-reduced q-points. If None, equal weights are assumed

reorder_frequencies(reorder_gamma=True)

By doing a dot product of eigenvectors at adjacent q-points, determines which modes are most similar and reorders the frequencies at each q-point. This means that the same mode will have the same index across different q-points, so will be plotted as the same colour in a dispersion plot, and can be followed across q-space.

Parameters: reorder_gamma (bool) – Whether to reorder frequencies at gamma-equivalent points. If an analytical correction has been applied at the gamma points (i.e LO-TO splitting) mode assignments can be incorrect at adjacent q-points where the correction hasn’t been applied. So you might not want to reorder at gamma for some materials
Return type: None

calculate_structure_factor(scattering_lengths='Sears1992', dw=None)

Calculate the one phonon inelastic scattering for neutrons at each q-point

Parameters

scattering_lengths (Union[str, Dict[str, Quantity]]) –
Dataset of coherent scattering length for each element in the structure. This may be provided in 3 ways:
- A string naming an appropriate data collection packaged with Euphonic (including the default value ‘Sears1992’). This will be passed to the collection argument of euphonic.util.get_reference_data().
- A string filename for a user’s customised data file in the same format as those packaged with Euphonic.
- An explicit dictionary of float Quantity, giving spin- and isotope-averaged coherent scattering length for each element in the structure, e.g.:
```
{'O': 5.803*ureg('fm'), 'Zn': 5.680*ureg('fm')}
```
dw (Optional[DebyeWaller]) – Data for thermal motion effects. Typically this is computed over a converged Monkhort-Pack grid, which need not correspond to the q-points of this QpointPhononModes object.

Return type

StructureFactor

Returns

sf – An object containing the structure factor for each q-point and phonon mode

Notes

This method calculates the mode-resolved (not binned in energy) coherent one-phonon neutron scattering function \(S(Q, \omega_{\mathbf{q}\nu})\) per atom, as defined in 1:

\[S(Q, \omega_{\mathbf{q}\nu}) = \frac{1}{2N_{atom}} \ \left\lvert {\sum_{\kappa}{\frac{b_\kappa}{M_{\kappa}^{1/2}{\omega_{\mathbf{q}\nu}^{1/2}}} (\mathbf{Q}\cdot \mathbf{e}_{\mathbf{q}\nu\kappa})\exp\left(i\mathbf{Q}{\cdot}\mathbf{r}_{\kappa}\right) \exp\left(-W_{\kappa}\right)}} \right\rvert^2\]

Where \(\nu\) runs over phonon modes, \(\kappa\) runs over atoms, \(b_\kappa\) is the coherent neutron scattering length, \(M_{\kappa}\) is the atom mass, \(r_{\kappa}\) is the vector from the origin to atom \(\kappa\) in the unit cell, \(\mathbf{e}_{\mathbf{q}\nu\kappa}\) are the normalised Cartesian eigenvectors, \(\omega_{\mathbf{q}\nu}\) are the frequencies and \(\exp\left(-W_\kappa\right)\) is the Debye-Waller factor. \(N_{atom}\) is the number of atoms in the unit cell, so the returned structure factor is per atom of sample.

1: M.T. Dove, Structure and Dynamics, Oxford University Press, Oxford, 2003, 225-226

calculate_debye_waller(temperature, frequency_min=<Quantity(0.01, 'millielectron_volt')>, symmetrise=True)

Calculate the 3 x 3 Debye-Waller exponent for each atom over the q-points contained in this object

Parameters

temperature (Quantity) – Scalar float Quantity. The temperature to use in the Debye-Waller calculation
frequency_min (Quantity) – Scalar float Quantity in energy units. Excludes frequencies below this limit from the calculation, as the calculation contains a 1/frequency factor which would result in infinite values. This also allows negative frequencies to be excluded
symmetrise (bool) – Whether to symmetrise the Debye-Waller factor based on the crystal symmetry operations. Note that if the Debye-Waller exponent is not symmetrised, the results may not be the same for unfolded and symmetry-reduced q-point grids

Return type

DebyeWaller

Returns

dw – An object containing the 3x3 Debye-Waller exponent matrix for each atom

Notes

As part of the structure factor calculation, \(\exp\left(W_\kappa\right)\) is the anisotropic Debye-Waller factor for atom \(\kappa\), and the exponent can be written as:

\[W_\kappa = \sum\limits_{\alpha\beta}W^{\alpha\beta}_{\kappa} Q_{\alpha}Q_{\beta}\]

This function calculates the \(W^{\alpha\beta}_{\kappa}\) part of the exponent, which is a \(N_\mathrm{atom}\times3\times3\) matrix that is independent of Q, so for efficiency can be precalculated to be used in the structure factor calculation. The Debye-Waller exponent matrix is calculated by 2:

\[W^{\alpha\beta}_{\kappa} = \frac{\hbar}{4M_{\kappa}\sum\limits_{\mathbf{q}}{\mathrm{weight}_\mathbf{q}}} \sum\limits_{\mathbf{q}\nu \in{BZ}}\mathrm{weight}_\mathbf{q}\frac{e_{\mathbf{q}\nu\kappa\alpha}e^{*}_{\mathbf{q}\nu\kappa\beta}} {\omega_{\mathbf{q}\nu}} \mathrm{coth}\left(\frac{\hbar\omega_{\mathbf{q}\nu}}{2k_BT}\right)\]

Where the sum is over q-points and modes \(\nu\) in the first Brillouin Zone (BZ), \(\kappa\) runs over atoms, \(\alpha,\beta\) run over the Cartesian directions, \(M_{\kappa}\) is the atom mass, \(e_{q\nu\kappa\alpha}\) are the scalar components of the normalised Cartesian eigenvectors, \(\omega_{\mathbf{q}\nu}\) are the frequencies, and \(\mathrm{weight}_\mathbf{q}\) is the per q-point symmetry weight (if the q-points are not symmetry-reduced, all weights will be equal).

2: G.L. Squires, Introduction to the Theory of Thermal Neutron Scattering, Dover Publications, New York, 1996, 34-37

calculate_pdos(dos_bins, mode_widths=None, mode_widths_min=<Quantity(0.01, 'millielectron_volt')>, adaptive_method='reference', adaptive_error=0.01, adaptive_error_fit='cubic', weighting=None, cross_sections='BlueBook')

Calculates partial density of states for each atom in the unit cell

Parameters

dos_bins (Quantity) – Shape (n_e_bins + 1,) float Quantity. The energy bin edges to use for calculating the DOS.
mode_widths (Optional[Quantity]) – Shape (n_qpts, n_branches) float Quantity in energy units. The broadening width for each mode at each q-point, for adaptive broadening.
mode_widths_min (Quantity) – Scalar float Quantity in energy units. Sets a lower limit on the mode widths, as mode widths of zero will result in infinitely sharp peaks.
adaptive_method (Literal[‘reference’, ‘fast’]) – String. Specifies whether to use slow, reference adaptive method or faster, approximate method. Allowed options are ‘reference’ or ‘fast’, default is ‘reference’.
adaptive_error (Optional[float]) – Scalar float. Acceptable error for gaussian approximations when using the fast adaptive method, defined as the absolute difference between the areas of the true and approximate gaussians
adaptive_error_fit (Literal[‘cheby-log’, ‘cubic’]) – Select parametrisation of kernel width spacing to adaptive_error when using fast approximate method. ‘cheby-log’ is recommended: for backward-compatibilty, ‘cubic’ is the default.
weighting (Optional[str]) – One of {‘coherent’, ‘incoherent’, ‘coherent-plus-incoherent’}. If provided, produces a neutron-weighted DOS, weighted by either the coherent, incoherent, or sum of coherent and incoherent neutron scattering cross-sections.
cross_sections (Union[str, Dict[str, Quantity]]) –
A dataset of cross-sections for each element in the structure, it can be a string specifying a dataset, or a dictionary explicitly giving the cross-sections for each element.

If cross_sections is a string, it is passed to the collection argument of euphonic.util.get_reference_data(). This collection must contain the ‘coherent_cross_section’ or ‘incoherent_cross_section’ physical property, depending on the weighting argument. If weighting is None, this string argument is not used.

If cross sections is a dictionary, the weighting argument is ignored, and these cross-sections are used directly to calculate the neutron-weighted DOS. It must contain a key for each element in the structure, and each value must be a Quantity in the appropriate units, e.g:
```
{'La': 8.0*ureg('barn'), 'Zr': 9.5*ureg('barn')}
```

Return type

Spectrum1DCollection

Returns

dos – A collection of spectra, with the energy bins on the x-axis and PDOS for each atom in the unit cell on the y-axis, metadata describes the index and species of each PDOS. If the PDOS is not neutron-weighted, it is in units of modes per energy unit per atom, such that the integrated area of the total DOS is equal to 3. If the PDOS is weighted by the neutron cross-section, it is in units of area per unit energy per atom, see Notes for details.

Notes

The PDOS is calculated as:

\[PDOS_{\kappa}(\omega) = \frac{1}{N_\mathrm{atom}\sum\limits_{\mathbf{q}}{\mathrm{weight}_{\mathbf{q}}}} \sum_{\mathbf{q}\nu\alpha\in{BZ}} \mathrm{weight}_{\mathbf{q}} {e_{\mathbf{q}\nu\kappa\alpha}e^{*}_{\mathbf{q}\nu\kappa\alpha}} \delta(\omega - \omega_{\mathbf{q}\nu})\]

Where \(N_\mathrm{atom}\) is the number of atoms in the unit cell, \(\mathrm{weight}_\mathbf{q}\) is the per q-point weight, \(e_{\mathbf{q}\nu\kappa\alpha}\) is the normalised eigenvector at \(\mathbf{q}\) for mode \(\nu\), atom \(\kappa\) in Cartesian direction \(\alpha\) and \(\omega_{\mathbf{q}\nu}\) is the phonon frequency at \(\mathbf{q}\) for mode \(\nu\), and the sum \(BZ\) is over the 1st Brillouin Zone.

The neutron-weighted PDOS is calculated as:

\[PDOS^\mathrm{neutron}_{\kappa}(\omega) = \frac{\sigma_{\kappa}M_\mathrm{avg}}{M_\kappa} PDOS_{\kappa}(\omega)\]

Where \(\sigma_\kappa\) is the neutron scattering cross section, \(M_\mathrm{avg}\) is the average mass of atoms in the unit cell and \(M_\kappa\) is the mass of atom \(\kappa\), so the neutron-weighted PDOS is returned in units of area per unit energy.

to_dict()

Convert to a dictionary. See QpointPhononModes.from_dict for details on keys/values

Return type: Dict[str, Any]

to_qpoint_frequencies()

Create a QpointFrequencies object

Return type: QpointFrequencies

classmethod from_dict(d)

Convert a dictionary to a QpointPhononModes object

Parameters

d (dict) –

A dictionary with the following keys/values:

’crystal’: dict, see Crystal.from_dict
’qpts’: (n_qpts, 3) float ndarray
’frequencies’: (n_qpts, 3*crystal.n_atoms) float ndarray
’frequencies_unit’: str
’eigenvectors’: (n_qpts, 3*crystal.n_atoms, crystal.n_atoms, 3) complex ndarray

There are also the following optional keys:

’weights’: (n_qpts,) float ndarray

Return type

TypeVar(T, bound= QpointPhononModes)

classmethod from_json_file(filename)

Read from a JSON file. See QpointPhononModes.from_dict for required fields

Parameters: filename (str) – The file to read from
Return type: TypeVar(T, bound= QpointPhononModes)

classmethod from_castep(filename, average_repeat_points=True, prefer_non_loto=False)

Reads precalculated phonon mode data from a CASTEP .phonon file

Parameters

filename (str) – The path and name of the .phonon file to read
average_repeat_points (bool) – If multiple frequency/eigenvectors blocks are included with the same q-point index (i.e. for Gamma-point with LO-TO splitting), scale the weights such that these sum to the given weight
prefer_non_loto (bool) – If multiple frequency/eigenvector blocks are provided with the same q-point index and all-but-one include a direction vector, use the data from the point without a direction vector. (i.e. use the “exact” Gamma data without non-analytic correction.) This option takes priority over average_repeat_points.

Return type

TypeVar(T, bound= QpointPhononModes)

classmethod from_phonopy(path='.', phonon_name='band.yaml', phonon_format=None, summary_name='phonopy.yaml')

Reads precalculated phonon mode data from a Phonopy mesh/band/qpoints.yaml/hdf5 file. May also read from phonopy.yaml for structure information.

Parameters

path (str) – Path to directory containing the file(s)
phonon_name (str) – Name of Phonopy file including the frequencies
phonon_format (Optional[str]) – Format of the phonon_name file if it isn’t obvious from the phonon_name extension, one of {‘yaml’, ‘hdf5’}
summary_name (str) – Name of Phonopy summary file to read the crystal information from. Crystal information in the phonon_name file takes priority, but if it isn’t present, crystal information is read from summary_name instead

Return type

TypeVar(T, bound= QpointPhononModes)

calculate_dos(dos_bins, mode_widths=None, mode_widths_min=<Quantity(0.01, 'millielectron_volt')>, adaptive_method='reference', adaptive_error=0.01, adaptive_error_fit='cubic')

Calculates a density of states, in units of modes per atom per energy unit, such that the integrated area is equal to 3.

Parameters

dos_bins (Quantity) – Shape (n_e_bins + 1,) float Quantity in energy units. The energy bin edges to use for calculating the DOS
mode_widths (Optional[Quantity]) – Shape (n_qpts, n_branches) float Quantity in energy units. The broadening width for each mode at each q-point, for adaptive broadening
mode_widths_min (Quantity) – Scalar float Quantity in energy units. Sets a lower limit on the mode widths, as mode widths of zero will result in infinitely sharp peaks
adaptive_method (Literal[‘reference’, ‘fast’]) – String. Specifies whether to use slow, reference adaptive method or faster, approximate method.
adaptive_error (float) – Scalar float. Acceptable error for gaussian approximations when using the fast adaptive method, defined as the absolute difference between the areas of the true and approximate gaussians
adaptive_error_fit (Literal[‘cheby-log’, ‘cubic’]) – Select parametrisation of kernel width spacing to adaptive_error. ‘cheby-log’ is recommended: for backward-compatibilty, ‘cubic’ is the default.

Return type

Spectrum1D

Returns

dos – A spectrum containing the energy bins on the x-axis and DOS on the y-axis. The DOS is in units of modes per energy unit per atom, such that the integrated area is equal to 3.

Notes

The fast adaptive broadening method reduces computation time by reducing the number of Gaussian functions that have to be evaluated. Broadening kernels are only evaulated for a subset of mode width values with intermediate values approximated using interpolation.

The adaptive_error keyword argument is used to determine how many broadening kernels are computed exactly. The more exact kernels are used, the more accurate the Gaussian approximations will be, but computation time will also be increased.

calculate_dos_map(dos_bins, mode_widths=None, mode_widths_min=<Quantity(0.01, 'millielectron_volt')>)

Produces a bandstructure-like plot, using the DOS at each q-point

Parameters

dos_bins (Quantity) – Shape (n_e_bins + 1,) float Quantity in energy units. The energy bin edges to use for calculating the DOS
mode_widths (Optional[Quantity]) – Shape (n_qpts, n_branches) float Quantity in energy units. The broadening width for each mode at each q-point, for adaptive broadening
mode_widths_min (Quantity) – Scalar float Quantity in energy units. Sets a lower limit on the mode widths, as mode widths of zero will result in infinitely sharp peaks

Return type

Spectrum2D

Returns

dos_map – A 2D spectrum containing the q-point bins on the x-axis, energy bins on the y-axis and DOS on the z-axis

get_dispersion()

Creates a set of 1-D bands from mode data

Bands follow the same q-point order as in the qpts array, with x-axis spacing corresponding to the absolute distances between q-points. Discontinuities will appear as large jumps on the x-axis.

Return type: Spectrum1DCollection
Returns: dispersion – A sequence of mode bands with a common x-axis

to_json_file(filename)

Write to a JSON file. JSON fields are equivalent to from_dict keys

Parameters: filename (str) – Name of the JSON file to write to
Return type: None