ForceConstants¶
The ForceConstants object contains the force constants, supercell, and
crystal structure information required to calculate phonon frequencies and
eigenvectors at any arbitrary q via Fourier interpolation.
Reading From CASTEP¶
The force constants matrix and other required information can be read from a
.castep_bin or .check file with
ForceConstants.from_castep:
from euphonic import ForceConstants
filename = 'quartz/quartz.castep_bin'
fc = ForceConstants.from_castep(filename)
By default CASTEP may not write the force constants, if you receive an error
saying the force constants could not be read, in the .param file ensure a
PHONON_FINE_METHOD has been chosen e.g. PHONON_FINE_METHOD: interpolate,
and set PHONON_WRITE_FORCE_CONSTANTS: true, then rerun CASTEP to trigger the
force constants to be written.
Reading From Phonopy¶
When using Phonopy with Euphonic, it is recommended that all the required data
(force constants, crystal structure, born charges if applicable) be collected
in a single phonopy.yaml file. This can be done by running Phonopy with the
--include-all flag or with INCLUDE_ALL = .TRUE.
(phonopy >= 2.5.0 only).
Required information is read from Phonopy output files using
ForceConstants.from_phonopy.
A path keyword argument can be supplied (if the files are in another
directory), and by default phonopy.yaml is read, but the filename can be
changed with the summary_name keyword argument:
from euphonic import ForceConstants
fc = ForceConstants.from_phonopy(path='NaCl',
summary_name='phonopy_fc.yaml')
If you are using an older version of Phonopy, the force constants and born
charges can also be read from Phonopy plaintext or hdf5 files by specifying the
fc_name and born_name keyword arguments:
from euphonic import ForceConstants
fc = ForceConstants.from_phonopy(path='NaCl',
fc_name='force_constants.hdf5',
born_name='BORN')
Calculating phonon frequencies and eigenvectors¶
Phonon frequencies and eigenvectors are calculated using
ForceConstants.calculate_qpoint_phonon_modes
(see the docstring for algorithm details). A Numpy array of q-points of shape
(n_qpts, 3) must be provided, and a
QpointPhononModes object is returned. A
recommended q-point path for plotting bandstructures can be generated using
seekpath:
import seekpath
import numpy as np
from euphonic import ForceConstants
# Read quartz data from quartz.castep_bin
filename = 'quartz/quartz.castep_bin'
fc = ForceConstants.from_castep(filename)
# Generate a recommended q-point path using seekpath
cell = fc.crystal.to_spglib_cell()
qpts = seekpath.get_explicit_k_path(cell)["explicit_kpoints_rel"]
# Calculate frequencies/eigenvectors
phonons = fc.calculate_qpoint_phonon_modes(qpts, asr='reciprocal')
Calculating phonon frequencies only¶
This uses the same algorithm as for calculating both the frequencies and
eigenvectors, only with lower memory requirements as the eigenvectors
are not stored. This is done using
ForceConstants.calculate_qpoint_frequencies
which returns a QpointFrequencies object.
Docstring¶
-
class
ForceConstants(crystal, force_constants, sc_matrix, cell_origins, born=None, dielectric=None)¶ A class to read and store the data required for a phonon interpolation calculation from model (e.g. CASTEP) output, and calculate phonon frequencies/eigenvectors at arbitrary q-points via Fourier interpolation
- Variables
crystal – Lattice and atom information
force_constants – Shape (n_cells_in_sc, 3*n_atoms, 3*n_atoms) float Quantity in energy/length**2 units. The force constants matrix
sc_matrix – Shape (3, 3) int ndarray. The supercell matrix
n_cells_in_sc – Number of cells in the supercell
cell_origins – Shape (n_cells_in_sc, 3) int ndarray. The locations of the unit cells within the supercell
born – Shape (n_atoms, 3, 3) float Quantity in charge units or None. The Born charges for each atom
dielectric – Shape (3, 3) float Quantity in charge**2/(length*energy) units or None. The dielectric permittivity tensor
-
__init__(crystal, force_constants, sc_matrix, cell_origins, born=None, dielectric=None)¶ - Parameters
crystal (
Crystal) – Lattice and atom informationforce_constants (
Quantity) – Shape (n_cells_in_sc, 3*n_atoms, 3*n_atoms) float Quantity in energy/length**2 units. The force constants matrixsc_matrix (
ndarray) – Shape (3, 3) int ndarray. The supercell matrixcell_origins (
ndarray) – Shape (n_cells_in_sc, 3) int ndarray. The locations of the unit cells within the supercellborn (
Optional[Quantity]) – Shape (n_atoms, 3, 3) float Quantity in charge units. The Born charges for each atomdielectric (
Optional[Quantity]) – Shape (3, 3) float Quantity in charge**2/(length*energy) units. The dielectric permittivity tensor
-
calculate_qpoint_phonon_modes(qpts, weights=None, asr=None, dipole=True, dipole_parameter=1.0, eta_scale=1.0, splitting=True, insert_gamma=False, reduce_qpts=True, use_c=None, n_threads=None, return_mode_gradients=False, return_mode_widths=False)¶ Calculate phonon frequencies and eigenvectors at specified q-points from a force constants matrix via Fourier interpolation
- Parameters
qpts (
ndarray) – Shape (n_qpts, 3) float ndarray. The q-points to interpolate onto in reciprocal cell vector unitsweights (
Optional[ndarray]) – Shape (n_qpts,) float ndarray. The weight for each q-point. If not given, equal weights are appliedasr (
Optional[str]) – One of {‘realspace’, ‘reciprocal’}. Which acoustic sum rule correction to apply. ‘realspace’ applies the correction to the force constant matrix in real space. ‘reciprocal’ applies the correction to the dynamical matrix at every q-pointdipole (
bool) – Whether to calculate the dipole tail correction to the dynamical matrix at each q-point using the Ewald sum, if the Born charges and dielectric permitivitty tensor are present.dipole_parameter (
float) – Changes the cutoff in real/reciprocal space for the dipole Ewald sum. A higher value uses more reciprocal terms. If tuned correctly this can result in performance improvements. See euphonic-optimise-dipole-parameter program for help on choosing a good dipole_parameter.eta_scale (
float) –Deprecated since version 0.6.0.
Please use dipole_parameter instead
splitting (
bool) – Whether to calculate the LO-TO splitting at the gamma points. Only applied if dipole is True and the Born charges and dielectric permitivitty tensor are present.insert_gamma (
bool) – If splitting is True, this will insert gamma points into qpts to store the extra split frequencies. Note this means that the length of qpts in the output QpointPhononModes object will not necessarily be the same as the input qpts. If qpts already contains double gamma points where you want split frequencies, leave this as False.reduce_qpts (
bool) – Whether to use periodicity to reduce all q-points and only calculate for unique q-points within the 1st BZ. This won’t change the output but could increase performance.use_c (
Optional[bool]) – Whether to use C instead of Python to calculate and diagonalise the dynamical matrix. By default this is None and will use the C extension if it is installed, and fall back to Python if not. If use_c=True, this will force use of the C extension and an error will be raised if it is not installed. If use_c=False, this will force use of Python, even if the C extension is installedn_threads (
Optional[int]) – The number of OpenMP threads to use when looping over q-points in C. By default this is None, in which case the environment variable EUPHONIC_NUM_THREADS will be used to determine number of threads, if this is not set then the value returned from multiprocessing.cpu_count() will be usedreturn_mode_gradients (
bool) – Whether to also return the vector mode gradients (in Cartesian coordinates). These can be converted to mode widths and used in adaptive broadening for DOS. For details on how these are calculated see the Notes sectionreturn_mode_widths (
bool) –Deprecated since version 0.5.2.
The mode widths as calculated were only applicable for adaptive broadening of DOS, this argument will be removed in favour of the more flexible return_mode_gradients, which will allow the calculation of direction-specific mode widths in the future, for example. The mode widths can still be obtained from the mode gradients using euphonic.util.mode_gradients_to_widths
- Return type
Union[QpointPhononModes,Tuple[QpointPhononModes,Quantity]]- Returns
qpoint_phonon_modes – A QpointPhononModes object containing the interpolated frequencies and eigenvectors at each q-point. Note that if there is LO-TO splitting, and insert_gamma=True, the number of input q-points may not be the same as in the output object
mode_gradients – Optional shape (n_qpts, n_branches, 3) float Quantity, the vector mode gradients dw/dQ in Cartesian coordinates. Is only returned if return_mode_gradients is true
- Raises
ImportCError – If use_c=True but the C extension cannot be imported
Notes
Phonon Frequencies/Eigenvectors Calculation
Phonon frequencies/eigenvectors are calculated at any q-point by Fourier interpolation of a force constants matrix. The force constants matrix is defined as 1:
\[\phi_{\alpha, {\alpha}'}^{\kappa, {\kappa}'} = \frac{\delta^{2}E}{{\delta}u_{\kappa,\alpha}{\delta}u_{{\kappa}',{\alpha}'}}\]Which gives the Dynamical matrix at q:
\[D_{\alpha, {\alpha}'}^{\kappa, {\kappa}'}(q) = \frac{1}{\sqrt{M_\kappa M_{\kappa '}}} \sum_{a}\phi_{\alpha, \alpha '}^{\kappa, \kappa '}e^{-iq\cdot r_a}\]The eigenvalue equation for the dynamical matrix is then:
\[D_{\alpha, {\alpha}'}^{\kappa, {\kappa}'}(q) \epsilon_{q\nu\kappa\alpha} = \omega_{q\nu}^{2} \epsilon_{q\nu\kappa\alpha}\]Where \(\nu\) runs over phonon modes, \(\kappa\) runs over atoms, \(\alpha\) runs over the Cartesian directions, \(a\) runs over unit cells in the supercell, \(u_{\kappa, \alpha}\) is the displacement of atom \(\kappa\) in direction \(\alpha\), \(M_{\kappa}\) is the mass of atom \(\kappa\), \(r_{a}\) is the vector to the origin of cell \(a\) in the supercell, \(\epsilon_{q\nu\kappa\alpha}\) are the eigevectors, and \(\omega_{q\nu}^{2}\) are the frequencies squared.
In polar materials, there is an additional long-ranged correction to the force constants matrix (applied if dipole=True) and a non-analytical correction at the gamma point 2 (applied if splitting=True).
- 1
Dove, Introduction to Lattice Dynamics, Cambridge University Press, Cambridge, 1993, 83-87
- 2
Gonze, K. C. Charlier, D. C. Allan, M. P. Teter, Phys. Rev. B, 1994, 50, 13035-13038
Mode Gradients Calculation
The mode gradients can be used to calculate mode widths to be used in the adaptive broadening scheme 3 when computing a DOS - this broadens each mode contribution individually according to its mode width. The mode widths are proportional to the mode gradients and can be estimated using
euphonic.util.mode_gradients_to_widthsThe mode gradients \(\frac{d\omega_{q\nu}}{dQ}\) at each Q are calculated as the same time as the phonon frequencies and eigenvectors as follows.
Firstly, the eigenvalue equation above can be written in matrix form as:
\[E(q)\Omega(q) = D(q)E(q)\]\[\Omega(q) = E^{-1}(q)D(q)E(q)\]Where \(\Omega(q)\) is the diagonal matrix of phonon frequencies squared \(\omega_{q\nu}^{2}\) and \(E(q)\) is the matrix containing eigenvectors for all modes. \(\frac{d\omega_{q\nu}}{dQ}\), can then be obtained by differentiating the above equation with respect to Q using the product rule:
\[\frac{d\Omega}{dQ} = 2\omega_{q\nu}\frac{d\omega_{q\nu}}{dQ}\delta_{\nu, \nu^{\prime}}\]\[\frac{d\omega_{q\nu}}{dQ} = \frac{1}{2\omega_{q\nu}}{( \frac{d{E^{-1}}}{dQ}DE + E^{-1}\frac{dD}{dQ}E + E^{-1}D\frac{dE}{dQ})}\]Given that eigenvectors are normalised and orthogonal, an identity can be employed:
\[\frac{d}{d\lambda}E^{\prime}E = 0\]So terms involving the first derivative of the eigenvector matrix can be set to zero, resulting in:
\[\frac{d\omega_{q\nu}}{dQ} = \frac{1}{2\omega_{q\nu}}{(E^{-1}\frac{dD}{dQ}E)}\]\(\frac{dD}{dQ}\) can be obtained by differentiating the Fourier equation above:
\[\frac{dD}{dQ} = \frac{-i r_a}{\sqrt{M_\kappa M_{\kappa '}}} \sum_{a}\phi_{\alpha, \alpha '}^{\kappa, \kappa '}e^{-iq\cdot r_a}\]- 3
Yates, X. Wang, D. Vanderbilt and I. Souza, Phys. Rev. B, 2007, 75, 195121
-
calculate_qpoint_frequencies(qpts, weights=None, asr=None, dipole=True, dipole_parameter=1.0, eta_scale=1.0, splitting=True, insert_gamma=False, reduce_qpts=True, use_c=None, n_threads=None, return_mode_gradients=False, return_mode_widths=False)¶ Calculate phonon frequencies (without eigenvectors) at specified q-points. See ForceConstants.calculate_qpoint_phonon_modes for argument and algorithm details
- Return type
Union[QpointFrequencies,Tuple[QpointFrequencies,Quantity]]
-
to_dict()¶ Convert to a dictionary. See ForceConstants.from_dict for details on keys/values
- Return type
Dict[str,Any]- Returns
dict
-
to_json_file(filename)¶ Write to a JSON file. JSON fields are equivalent to ForceConstants.from_dict keys
- Parameters
filename (
str) – Name of the JSON file to write to- Return type
None
-
classmethod
from_dict(d)¶ Convert a dictionary to a ForceConstants object
- Parameters
d (
Dict[str,Any]) –A dictionary with the following keys/values:
’crystal’: dict, see Crystal.from_dict
’force_constants’: (n_cells_in_sc, 3*crystal.n_atoms, 3*crystal.n_atoms) float ndarray
’force_constants_unit’: str
’sc_matrix’: (3,3) int ndarray
’cell_origins’: (n_cells_in_sc, 3) int ndarray
There are also the following optional keys:
’born’: (3*crystal.n_atoms, 3, 3) float ndarray
’born_unit’: str
’dielectric’: (3, 3) float ndarray
’dielectric_unit’: str
- Return type
~T
- Returns
forceconstants
-
classmethod
from_json_file(filename)¶ Read from a JSON file. See ForceConstants.from_dict for required fields
- Parameters
filename (
str) – The file to read from- Return type
~T
- Returns
forceconstants
-
classmethod
from_castep(filename)¶ Reads from a .castep_bin or .check file
- Parameters
filename (
str) – The path and name of the file to read- Return type
~T
- Returns
forceconstants
-
classmethod
from_phonopy(path='.', summary_name='phonopy.yaml', born_name=None, fc_name='FORCE_CONSTANTS', fc_format=None)¶ Reads data from the phonopy summary file (default phonopy.yaml) and optionally born and force constants files. Only attempts to read from born or force constants files if these can’t be found in the summary file.
- Parameters
path (
str) – Path to directory containing the file(s)summary_name (
str) – Filename of phonopy summary file, default phonopy.yaml. By default any information (e.g. force constants) read from this file takes priorityborn_name (
Optional[str]) – Name of the Phonopy file containing born charges and dielectric tensor (by convention in Phonopy this would be called BORN). Is only read if Born charges can’t be found in the summary_name filefc_name (
str) – Name of file containing force constants. Is only read if force constants can’t be found in summary_namefc_format (
Optional[str]) – One of {‘phonopy’, ‘hdf5’}. Format of file containing force constants data, if it can’t be determined from the extension. FORCE_CONSTANTS is type ‘phonopy’
- Return type
~T
- Returns
forceconstants