# coding: utf-8
# Distributed under the terms of the MIT license.
""" This file implements classes to store and manipulate electronic and
vibrational bandstructures, with or without projection data.
"""
import numpy as np
from matador.orm.spectral.spectral import Spectral
from matador.utils.chem_utils import INVERSE_CM_TO_EV, KELVIN_TO_EV
EPS = 1e-4
[docs]class Dispersion(Spectral):
"""Parent class for continuous spectra in reciprocal space, i.e.
electronic and vibrational bandstructures.
Note:
This class speaks of "k-points" as general reciprocal space points
used to display the dispersion curves; these correspond to CASTEP's
phonon_kpoints or spectral_kpoints, and not the k-points
used to generate the underlying wavefunction or dynamical matrix.
"""
[docs] def find_full_kpt_branch(self):
"""Find all branch indices from branch start indices."""
branch_inds = []
for ind, start_ind in enumerate(self.kpoint_branch_start[:-1]):
branch_inds.append(
list(range(start_ind, self.kpoint_branch_start[ind + 1]))
)
branch_inds.append(list(range(self.kpoint_branch_start[-1], self.num_kpoints)))
if not sum([len(branch) for branch in branch_inds]) == self.num_kpoints:
raise RuntimeError(
"Error parsing kpoints: number of kpoints does "
"not match number in branches"
)
return branch_inds
[docs] def set_branches_and_spacing(self):
"""Set the relevant kpoint spacing and branch attributes."""
branch_start, spacing = self.find_kpoint_branches()
self._data["kpoint_path_spacing"] = spacing
self._data["kpoint_branch_start"] = branch_start
[docs] def find_kpoint_branches(self):
"""Separate a kpoint path into discontinuous branches,
Returns:
list[list[int]]: list of lists containing the indices of
the discontinuous kpoint branches.
float: estimated k-point spacing from median of their
separations.
"""
kpt_diffs = np.linalg.norm(np.diff(self.kpoint_path_cartesian, axis=0), axis=1)
spacing = np.median(kpt_diffs)
# add 0 as its the start of the first path, then add all indices
# have to add 1 to the where to get the start rather than end of the branch
branch_start = [0] + (np.where(kpt_diffs > 3 * spacing)[0] + 1).tolist()
return branch_start, spacing
[docs] def linearise_path(self, preserve_kspace_distance=False):
"""For a given k-point path, normalise the spacing between points, mapping
it onto [0, 1].
Keyword arguments:
preserve_kspace_distance (bool): if True, point separation
will be determined by their actual separation in
reciprocal space, otherwise they will be evenly spaced.
Returns:
np.ndarray: 3xN array containing k-points mapped onto [0, 1].
"""
path = [0]
for branch in self.kpoint_branches:
for ind, kpt in enumerate(self.kpoint_path[branch]):
if ind != len(branch) - 1:
if preserve_kspace_distance:
diff = np.sqrt(
np.sum((kpt - self.kpoint_path[branch[ind + 1]]) ** 2)
)
else:
diff = 1.0
path.append(path[-1] + diff)
path = np.asarray(path)
path /= np.max(path)
if len(path) != self.num_kpoints - len(self.kpoint_branches) + 1:
raise RuntimeError("Linearised kpoint path has wrong number of kpoints!")
return path
[docs] def reorder_bands(self):
"""Reorder the bands of this Dispersion object directly."""
self._data["eigs_s_k"] = self.get_band_reordering(
self.eigs, self.kpoint_branches
)
[docs] @staticmethod
def get_band_reordering(eigs, kpoint_branches):
"""Recursively reorder eigenvalues such that bands join up correctly,
based on local gradients.
Parameters:
dispersion (numpy.ndarray): array containing eigenvalues as a
function of q/k
branches (:obj:`list` of :obj:`int`): list containing branches of
k/q-point path
Returns:
numpy.ndarray: reordered branches.
"""
sorted_eigs = np.array(eigs, copy=True)
num_bands = np.shape(sorted_eigs)[1]
for channel_ind, channel in enumerate(eigs):
eigs = channel
for branch_ind, branch in enumerate(kpoint_branches):
eigs_branch = channel[:, branch]
converged = False
counter = 0
i_cached = 0
while not converged and counter < len(branch):
counter += 1
for i in range(i_cached + 1, len(branch) - 1):
guess = 2 * eigs_branch[:, i] - eigs_branch[:, i - 1]
argsort_guess = np.argsort(guess)
if np.any(
np.argsort(guess) != np.argsort(eigs_branch[:, i + 1])
):
tmp_copy = np.array(channel, copy=True)
for ind, mode in enumerate(
np.argsort(eigs_branch[:, i]).tolist()
):
eigs_branch[mode, i + 1 :] = tmp_copy[:, branch][
argsort_guess[ind], i + 1 :
]
for other_branch in kpoint_branches[branch_ind:]:
eigs_other_branch = channel[:, other_branch]
for ind, mode in enumerate(
np.argsort(channel[:, i]).tolist()
):
eigs_other_branch[mode] = tmp_copy[:, other_branch][
argsort_guess[ind]
]
channel[:, other_branch] = eigs_other_branch
channel[:, branch] = eigs_branch
i_cached = i
break
else:
converged = True
sorted_eigs[channel_ind] = channel.reshape(1, num_bands, len(channel[0]))
return sorted_eigs
[docs] def plot_dispersion(self, **kwargs):
"""Make a plot of the band structure, with projections, if found."""
from matador.plotting.spectral_plotting import plot_spectral
_kwargs = {
"plot_dos": False,
"plot_bandstructure": True,
"plot_pdis": "projector_weights" in self,
"phonons": "Vibrational" in self.__class__.__name__,
}
_kwargs.update(kwargs)
plot_spectral(self, **_kwargs)
def _reshaped_eigs(self, eigs, shape):
"""Attempts to reshape the eigenvalues into the desired shape.
Parameters:
eigs (np.ndarray): the eigs to reshape.
shape (tuple): the desired shape.
Returns:
np.ndarray: the reshaped eigs.
"""
raise NotImplementedError(
"Wrong eigenvalue shape passed, and reshape function is not yet implemented. "
"Eigs should have shape {}, not {}".format(shape, np.shape(eigs))
)
[docs]class ElectronicDispersion(Dispersion):
"""Class that stores electronic dispersion data. Attributes are
all implemented as properties based on underlying raw data.
Attributes:
source (list): list of source files.
num_kpoints (int): number of kpoints.
num_spins (int): number of spin channels.
num_bands (int): number of bands.
num_electrons (int): number of bands.
eigs_s_k (numpy.ndarray): eigenvalue array of shape
(num_spins, num_bands, num_kpoints).
kpoint_path (numpy.ndarray): array of shape (num_kpoints, 3)
containing the k-point path in fractional coordinates.
kpoint_path_cartesian (numpy.ndarray): as above, in Cartesian
coordinates.
fermi_energy (float): Fermi energy in eV (takes average of spin
channels if more than one is present).
spin_fermi_energy (list[float]): Fermi energy for each spin channel
in eV.
band_gap (float): smallest band gap across spin channels/momenta.
spin_band_gap (list[float]): band gap per spin channel.
projectors (list[tuple]): list of projector labels in format
(`element`, `l-channel`).
projector_weights (numpy.ndarray): array of projector_weights with shape
(num_kpoints, num_bands, num_projectors)
Note:
projector_weights will only work with one spin channel.
"""
required_keys = [
"num_kpoints",
"num_spins",
"num_bands",
"num_electrons",
"eigs_s_k",
"kpoint_path",
"lattice_cart",
"fermi_energy",
]
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
if kwargs.get("projection_data") is not None:
projection_data = kwargs.get("projection_data")
self._data["projectors"] = projection_data["projectors"]
self._data["projector_weights"] = projection_data["projector_weights"]
else:
self._data["projectors"] = self._data.get("projectors")
self._data["projector_weights"] = self._data.get("projector_weights")
if self._data.get("projectors") is not None:
if self.num_spins != 1:
raise NotImplementedError(
"Projected dispersion not implemented" " for multiple spin channels"
)
# only want to take projectors and projector_weights from this data
proj_shape = np.shape(self._data["projector_weights"])
expected_shape = (self.num_kpoints, self.num_bands, self.num_projectors)
if proj_shape != expected_shape:
raise RuntimeError(
f"Incompatible shape of projector weights: {proj_shape}, was expecting {expected_shape}"
)
shape = (self.num_spins, self.num_bands, self.num_kpoints)
if np.shape(self._data["eigs_s_k"]) != shape:
self._data["eigs_s_k"] = self._reshaped_eigs(self._data["eigs_s_k"], shape)
@property
def num_spins(self):
"""Number of spin channels in spectrum."""
return self._data["num_spins"]
@property
def num_electrons(self):
"""Number of electrons."""
return self._data["num_electrons"]
@property
def eigs_s_k(self):
"""Array of electronic eigenvalues with shape
(num_spins, num_bands, num_kpoints).
"""
return self._data["eigs_s_k"]
@property
def eigs(self):
"""Alias for `self.eigs_s_k`."""
return self.eigs_s_k
@property
def fermi_energy(self):
"""Return the Fermi energy as described in the raw data."""
return self._data["fermi_energy"] or np.mean(
self._data.get("spin_fermi_energy")
)
@property
def spin_fermi_energy(self):
"""Return the Fermi energy as described in the raw data."""
return self._data.get("spin_fermi_energy")
@property
def band_gap(self):
"""Return the band gap of the system."""
if not self._data.get("band_gap"):
self.set_gap_data()
return self._data["band_gap"]
@property
def band_gap_path_inds(self):
"""Return the indices of the k-points that comprise the smallest
band gap.
"""
if not self._data.get("band_gap_path_inds"):
self.set_gap_data()
return self._data["band_gap_path_inds"]
@property
def spin_band_gap(self):
"""Return the band gap for each spin channel."""
if not self._data.get("spin_band_gap"):
self.set_gap_data()
return self._data["spin_band_gap"]
@property
def spin_band_gap_path_inds(self):
"""Return the indices of the k-points that comprise the smallest
band gap for each spin channel.
"""
if not self._data.get("spin_band_gap_path_inds"):
self.set_gap_data()
return self._data["spin_band_gap_path_inds"]
[docs] def set_gap_data(self):
"""Loop over bands to set the band gap, VBM, CBM, their
positions and the smallest direct gap inside self._data,
for each spin channel. Sets self.band_gap to be the smallest of
the band gaps across all spin channels.
"""
spin_keys = [
"spin_band_gap",
"spin_band_gap_path",
"spin_band_gap_path_inds",
"spin_valence_band_min",
"spin_conduction_band_max",
"spin_gap_momentum",
"spin_direct_gap",
"spin_direct_gap_path",
"spin_direct_gap_path_inds",
"spin_direct_valence_band_min",
"spin_direct_conduction_band_max",
]
for key in spin_keys:
self._data[key] = self.num_spins * [None]
for ispin in range(self.num_spins):
vbm = -1e10
cbm = 1e10
cbm_pos = []
vbm_pos = []
# calculate indirect gap
for _, branch in enumerate(self.kpoint_branches):
for nb in range(self.num_bands):
band = (
self.eigs_s_k[ispin][nb][branch] - self.spin_fermi_energy[ispin]
)
argmin = np.argmin(band)
argmax = np.argmax(band)
if vbm + EPS < band[argmax] < 0:
vbm = band[argmax]
vbm_pos = [branch[argmax]]
elif vbm - EPS <= band[argmax] < 0:
vbm = band[argmax]
vbm_pos.extend([branch[argmax]])
if cbm - EPS > band[argmin] > 0:
cbm = band[argmin]
cbm_pos = [branch[argmin]]
elif cbm + EPS >= band[argmin] > 0:
cbm = band[argmin]
cbm_pos.extend([branch[argmin]])
if band[argmin] + EPS / 2 < 0 < band[argmax] - EPS / 2:
vbm = 0
cbm = 0
vbm_pos = [0]
cbm_pos = [0]
break
smallest_diff = 1e10
for _cbm_pos in cbm_pos:
for _vbm_pos in vbm_pos:
if abs(_vbm_pos - _cbm_pos) < smallest_diff:
tmp_cbm_pos = _cbm_pos
tmp_vbm_pos = _vbm_pos
smallest_diff = abs(_vbm_pos - _cbm_pos)
cbm_pos = tmp_cbm_pos
vbm_pos = tmp_vbm_pos
self._data["spin_valence_band_min"][ispin] = (
vbm + self.spin_fermi_energy[ispin]
)
self._data["spin_conduction_band_max"][ispin] = (
cbm + self.spin_fermi_energy[ispin]
)
self._data["spin_band_gap"][ispin] = cbm - vbm
self._data["spin_band_gap_path"][ispin] = [
self.kpoint_path[cbm_pos],
self.kpoint_path[vbm_pos],
]
self._data["spin_band_gap_path_inds"][ispin] = [cbm_pos, vbm_pos]
self._data["spin_gap_momentum"][ispin] = np.linalg.norm(
self.kpoint_path_cartesian[cbm_pos]
- self.kpoint_path_cartesian[vbm_pos]
)
# calculate direct gap
direct_gaps = np.zeros(self.num_kpoints)
direct_cbms = np.zeros(self.num_kpoints)
direct_vbms = np.zeros(self.num_kpoints)
for ind, _ in enumerate(self.kpoint_path):
direct_cbm = 1e10
direct_vbm = -1e10
for nb in range(self.num_bands):
band_eig = (
self.eigs_s_k[ispin][nb][ind] - self.spin_fermi_energy[ispin]
)
if direct_vbm <= band_eig < EPS:
direct_vbm = band_eig
if direct_cbm >= band_eig > EPS:
direct_cbm = band_eig
direct_gaps[ind] = direct_cbm - direct_vbm
direct_cbms[ind] = direct_cbm
direct_vbms[ind] = direct_vbm
self._data["spin_direct_gap"][ispin] = np.min(direct_gaps)
self._data["spin_direct_conduction_band_max"][ispin] = (
direct_cbms[np.argmin(direct_gaps)] + self.spin_fermi_energy[ispin]
)
self._data["spin_direct_valence_band_min"][ispin] = (
direct_vbms[np.argmin(direct_gaps)] + self.spin_fermi_energy[ispin]
)
self._data["spin_direct_gap"][ispin] = np.min(direct_gaps)
self._data["spin_direct_gap_path"][ispin] = 2 * [
self.kpoint_path[np.argmin(direct_gaps)]
]
self._data["spin_direct_gap_path_inds"][ispin] = 2 * [
np.argmin(direct_gaps)
]
if (
np.abs(
self._data["spin_direct_gap"][ispin]
- self._data["spin_band_gap"][ispin]
)
< EPS
):
self._data["spin_valence_band_min"][ispin] = (
direct_vbm + self.spin_fermi_energy[ispin]
)
self._data["spin_conduction_band_max"][ispin] = (
direct_cbm + self.spin_fermi_energy[ispin]
)
self._data["spin_band_gap_path_inds"][ispin] = self._data[
"spin_direct_gap_path_inds"
][ispin]
cbm_pos = self._data["spin_direct_gap_path_inds"][ispin][0]
vbm_pos = self._data["spin_direct_gap_path_inds"][ispin][1]
self._data["spin_band_gap_path"][ispin] = self._data[
"spin_direct_gap_path"
][ispin]
self._data["spin_gap_momentum"][ispin] = np.linalg.norm(
self.kpoint_path_cartesian[cbm_pos]
- self.kpoint_path_cartesian[vbm_pos]
)
spin_gap_index = np.argmin(self._data["spin_band_gap"])
# use smallest spin channel gap data for standard non-spin data access
for key in spin_keys:
self._data[key.replace("spin_", "")] = self._data[key][spin_gap_index]
[docs] def new_from_trimmed_path(self, k_start_ind=0, k_end_ind=None):
"""Returns a new ElectronicDispersion object with the kpoint
path trimmed by the provided indices.
"""
_new_data_dict = {}
if k_end_ind is None:
k_end_ind = len(self.kpoint_path)
_new_data_dict["kpoint_path"] = self.kpoint_path[k_start_ind:k_end_ind]
_new_data_dict["eigs_s_k"] = self.eigs[:, :, k_start_ind:k_end_ind]
_new_data_dict["num_bands"] = self.num_bands
_new_data_dict["num_electrons"] = self.num_electrons
_new_data_dict["num_spins"] = self.num_spins
_new_data_dict["num_kpoints"] = len(_new_data_dict["kpoint_path"])
_new_data_dict["lattice_cart"] = self.lattice_cart
_new_data_dict["fermi_energy"] = self.fermi_energy
_new_data_dict["spin_fermi_energy"] = self.spin_fermi_energy
return ElectronicDispersion(**_new_data_dict)
[docs]class VibrationalDispersion(Dispersion):
"""Class that stores vibrational dispersion data. Attributes are
all implemented as properties based on underlying raw data.
Attributes:
source (list): list of source files.
num_kpoints (int): number of kpoints.
num_atoms (int): number of atoms.
num_modes (int): number of phonon modes.
eigs (numpy.ndarray): eigenvalue array of shape
(1, num_modes, num_kpoints), in frequency units below, with
first index denoting the single "spin channel" for phonons.
freq_unit (str): human-readable frequency unit used for eig array.
kpoint_path (numpy.ndarray): array of shape (num_kpoints, 3)
containing the k-point path in fractional coordinates.
kpoint_path_cartesian (numpy.ndarray): as above, in Cartesian
coordinates.
"""
required_keys = ["num_kpoints", "num_modes", "eigs_q", "kpoint_path"]
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
shape = (1, self.num_modes, self.num_qpoints)
if np.shape(self._data["eigs_q"]) != shape:
self._data["eigs_q"] = self._reshaped_eigs(self._data["eigs_q"], shape)
@property
def num_atoms(self):
"""Number of atoms in cell."""
return self._data["num_atoms"]
@property
def num_modes(self):
"""Number of phonon modes."""
return self._data["num_modes"]
@property
def num_bands(self):
"""Alias for number of modes."""
return self._data["num_modes"]
@property
def eigs_q(self):
"""Eigenvalues in frequency units `self.freq_unit`, with shape
(1, num_modes, num_kpoints).
"""
return np.asarray(self._data["eigs_q"])
@property
def softest_mode_freq(self):
"""The frequency of the softest mode in the calculation.
Negative modes correspond to imaginary frequencies.
"""
return np.min(self.eigs)
@property
def debye_temperature(self):
"""Returns the Debye temperature in K."""
return self.debye_freq * INVERSE_CM_TO_EV / KELVIN_TO_EV
@property
def debye_freq(self):
"""Returns the Debye frequency in cm^-1."""
return np.max(self.eigs)