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 reordering can be seen the plotting dispersion (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 = 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:

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
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

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

Parameters

crystal (Crystal) – Lattice and atom information
qpts (ndarray) – Shape (n_qpts, 3) float ndarray. The Q-point coordinates
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) complex ndarray. The dynamical matrix eigenvectors
weights (Optional[ndarray]) – Shape (n_qpts,) float ndarray. The weight for each q-point. 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

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) one-phonon neutron scattering function \(S(Q, \omega_{q\nu})\) per atom, as defined in 1. Note that internally Euphonic uses atomic units so \(\hbar\) has been omitted from the formulation:

\[S(Q, \omega_{q\nu}) = \frac{1}{2N_{atom}} \ \left\lvert \ \sum\limits_\kappa\frac{b_\kappa}{M_{\kappa}^{1/2}\omega_{q\nu}^{1/2}} \ [Q\cdot\epsilon_{q\nu\kappa\alpha}]e^{iQ{\cdot}r_\kappa}e^{-W} \ \right\rvert^2\]

Where \(\nu\) runs over phonon modes, \(\kappa\) runs over atoms, \(\alpha\) runs over the Cartesian directions, \(b_\kappa\) is the coherent neutron scattering length, \(M_{\kappa}\) is the atom mass, \(r_{\kappa}\) is the vector to atom \(\kappa\) in the unit cell, \(\epsilon_{q\nu\kappa\alpha}\) are the eigevectors, \(\omega_{q\nu}\) are the frequencies and \(e^{-W}\) 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 for each atom

Notes

As part of the structure factor calculation, the anisotropic Debye-Waller factor is defined as:

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

The Debye-Waller exponent is defined as \(W^{\kappa}_{\alpha\beta}\) and is independent of Q, so for efficiency can be precalculated to be used in the structure factor calculation. The Debye-Waller exponent is calculated by 2. Note that internally Euphonic uses atomic units so \(\hbar\) has been omitted from the formulation:

\[W^{\kappa}_{\alpha\beta} = \frac{1}{4M_{\kappa}\sum\limits_{q}{weight_q}} \sum\limits_{q\nu}weight_q\frac{\epsilon_{q\nu\kappa\alpha}\epsilon^{*}_{q\nu\kappa\beta}} {\omega_{q\nu}} coth(\frac{\omega_{q\nu}}{2k_BT})\]

Where \(\nu\) runs over phonon modes, \(\kappa\) runs over atoms, \(\alpha,\beta\) run over the Cartesian directions, \(M_{\kappa}\) is the atom mass, \(\epsilon_{q\nu\kappa\alpha}\) are the eigenvectors, \(\omega_{q\nu}\) are the frequencies, and \(weight_q\) is the per q-point weight. The q-points should be distributed over the 1st Brillouin Zone.

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')>, 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.
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. If weighting is None, the y-axis is in 1/energy units. If weighting is specified or cross_sections are supplied, the y-axis is in area/energy units per average atom.

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

~T

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: ~T

classmethod from_castep(filename)¶

Reads precalculated phonon mode data from a CASTEP .phonon file

Parameters: filename (str) – The path and name of the .phonon file to read
Return type: ~T

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

~T

calculate_dos(dos_bins, mode_widths=None, mode_widths_min=<Quantity(0.01, 'millielectron_volt')>)¶

Calculates a density of states

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

Spectrum1D

Returns

dos – A spectrum containing the energy bins on the x-axis and dos on the y-axis

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