material

MaterialParameters(N=None, S=None, L=None, L_dot_S=None, L_tens_S=None, lambda_S=None, lambda_L=None, m_psi=None)

Class for DM-relevant material properties, such as the number of fermions, spin, orbital angular momentum, etc.

Attributes:

Name	Type	Description
N	dict	Fermion numbers.
S	dict	Spin vectors.
L	dict	Orbital angular momentum vectors.
L_dot_S	dict	\(L \cdot S\)
L_tens_S	dict	Spin orbit coupling tensor \(L \otimes S\)
lambda_S	ArrayLike	spin-coefficient for magnons
lambda_L	ArrayLike	orbital angular mom.-coefficient for magnons
m_psi	dict	Dictionary of masses for different particles.

Methods:

Name	Description
validate_for_phonons	Validates that the material properties are suitable for phonon calculations.
validate_for_magnons	Validates that the material properties are suitable for magnon calculations.

Parameters:

Name	Type	Description	Default
N	dict	Fermion numbers.	None
S	dict	Spin vectors.	None
L	dict	Orbital angular momentum vectors.	None
L_dot_S	dict	\(L \cdot S\)	None
L_tens_S	dict	Spin orbit coupling tensor \(L \otimes S\)	None
lambda_S	ArrayLike	spin-coefficient for magnons	None
lambda_L	ArrayLike	orbital angular mom.-coefficient for magnons	None
m_psi	dict	Masses of the fermions. Defaults to NIST values.	None

Source code in darkmagic/material.py

def __init__(
    self,
    N: dict = None,
    S: dict = None,
    L: dict = None,
    L_dot_S: dict = None,
    L_tens_S: dict = None,
    lambda_S: ArrayLike = None,
    lambda_L: ArrayLike = None,
    m_psi: dict = None,
):
    r"""
    Material properties constructor. All dicts have keys "n", "p", "e" for neutron, proton and electron. Any missing values are instantiated to 0.

    Args:
        N (dict, optional): Fermion numbers.
        S (dict, optional): Spin vectors.
        L (dict, optional): Orbital angular momentum vectors.
        L_dot_S (dict, optional): $L \cdot S$
        L_tens_S (dict, optional): Spin orbit coupling tensor $L \otimes S$
        lambda_S (ArrayLike, optional): spin-coefficient for magnons
        lambda_L (ArrayLike, optional): orbital angular mom.-coefficient for magnons
        m_psi (dict, optional): Masses of the fermions. Defaults to NIST values.
    """
    # Phonons
    self.N = N
    self.S = S
    self.L = L
    self.L_dot_S = L_dot_S
    self.L_tens_S = L_tens_S
    # Magnons
    self.lambda_S = lambda_S
    self.lambda_L = lambda_L
    # Mass of the particles
    self.m_psi = m_psi

validate_for_phonons(n_atoms)

Validates that the material properties are suitable for phonons. Namely, at least one of N, S, L, L_dot_S or L_tens_S must be defined.

Parameters:

Name	Type	Description	Default
n_atoms	int	Number of atoms in the material.	required

Raises:

Type	Description
AssertionError	If any of the required material properties for phonons are missing or have incorrect dimensions.

Source code in darkmagic/material.py

def validate_for_phonons(self, n_atoms: int) -> None:
    """
    Validates that the material properties are suitable for phonons.
    Namely, at least one of N, S, L, L_dot_S or L_tens_S must be defined.

    Args:
        n_atoms (int): Number of atoms in the material.

    Raises:
        AssertionError: If any of the required material properties for phonons are missing or have incorrect dimensions.

    """
    assert any([self.N, self.S, self.L, self.L_dot_S, self.L_tens_S])
    for d in [self.N, self.S, self.L, self.L_dot_S, self.L_tens_S]:
        if d:
            assert any(np.any(v) for v in d.values())

    self._validate_input(n_atoms)

validate_for_magnons(n_atoms)

Validates that the material properties are suitable for magnons. Namely, at least one of lambda_S and lambda_L must be defined.

Parameters:

Name	Type	Description	Default
n_atoms	int	Number of atoms in the material.	required

Raises:

Type	Description
AssertionError	If any of the required material properties for magnons are missing or have incorrect dimensions.

Source code in darkmagic/material.py

def validate_for_magnons(self, n_atoms: int) -> None:
    """
    Validates that the material properties are suitable for magnons. Namely, at least one of lambda_S and lambda_L must be defined.

    Args:
        n_atoms (int): Number of atoms in the material.

    Raises:
        AssertionError: If any of the required material properties for magnons are missing or have incorrect dimensions.

    """
    assert any([np.any(self.lambda_S), np.any(self.lambda_L)])
    # TODO: not nice to have so many return values
    self._validate_input(n_atoms)

Material(name, properties, structure, m_atoms)

Represents a generic material with its structural and atomic properties.

Attributes:

Name	Type	Description
name	str	The name of the material.
properties	MaterialProperties	The properties of the material.
real_frac_to_cart	ndarray	The transformation matrix from fractional to Cartesian coordinates (units 1/eV), in real space.
real_cart_to_frac	ndarray	The transformation matrix from Cartesian (units 1/eV) to fractional coordinates, in real space.
recip_frac_to_cart	ndarray	The transformation matrix from fractional to Cartesian coordinates (units eV), in k-space.
recip_cart_to_frac	ndarray	The transformation matrix from Cartesian (units eV) to fractional coordinates, in k-space.
m_atoms	ArrayLike	an array of atomic masses, in eV.
m_cell	ndarray	The total mass of the atoms in the material, in eV.
xj	ndarray	The Cartesian coordinates (units 1/eV) of the atoms in the material.
structure	Structure	the crystal structure pymatgen Structure object.
n_atoms	int	The number of atoms in the material.

Parameters:

Name	Type	Description	Default
name	str	The name of the material.	required
properties	MaterialProperties	The properties of the material.	required
structure	Structure	The structure of the material.	required
m_atoms	ArrayLike	atomic masses in eV.	required

Source code in darkmagic/material.py

def __init__(
    self,
    name: str,
    properties: MaterialParameters,
    structure: Structure,
    m_atoms: ArrayLike,
):
    """
    Constructor for a generic Material object

    Args:
        name (str): The name of the material.
        properties (MaterialProperties): The properties of the material.
        structure (Structure): The structure of the material.
        m_atoms (ArrayLike): atomic masses in eV.
    """
    # Material properties
    self.name = name
    self.properties = properties

    # Define transformation matrices
    self.real_frac_to_cart = structure.lattice.matrix.T
    self.real_cart_to_frac = LA.inv(self.real_frac_to_cart)
    self.recip_frac_to_cart = structure.lattice.reciprocal_lattice.matrix.T
    self.recip_cart_to_frac = LA.inv(self.recip_frac_to_cart)

    # Atomic and structural properties
    self.m_atoms = m_atoms
    self.m_cell = np.sum(m_atoms)
    self.xj = structure.cart_coords
    self.structure = structure
    self.n_atoms = len(structure.species)

    # Internal variables
    self._max_dE = None
    self._q_cut = None

PhononMaterial(name, properties, phonopy_yaml_path)

Bases: Material

A class for materials with phonons.

Attributes:

Name	Type	Description
phonopy_file	Phonopy	The Phonopy object for the material's phonons
n_modes	int	The number of phonon modes in the material.
born	ndarray	The born effective charges
epsilon	ndarray	The dielectric tensor

Parameters:

Name	Type	Description	Default
name	str	The name of the material.	required
properties	MaterialProperties	The properties of the material.	required
phonopy_yaml_path	str	The path to the Phonopy YAML file.	required

Source code in darkmagic/material.py

def __init__(
    self, name: str, properties: MaterialParameters, phonopy_yaml_path: str
):
    """
    Constructor for PhononMaterial objects.

    Args:
        name (str): The name of the material.
        properties (MaterialProperties): The properties of the material.
        phonopy_yaml_path (str): The path to the Phonopy YAML file.

    """
    # TODO: Need a check for when phonopy_yaml does not have NAC
    phonopy_file = phonopy.load(phonopy_yaml=phonopy_yaml_path, is_nac=True)
    # TODO: should be a dict that has the correct factor for all codes
    length_factor = const.bohr_to_Ang if phonopy_file.calculator == "qe" else 1.0
    self.phonopy_file = phonopy_file
    n_atoms = phonopy_file.primitive.get_number_of_atoms()
    self.n_modes = 3 * n_atoms

    properties.validate_for_phonons(n_atoms)

    m_atoms = phonopy_file.primitive.masses * const.amu_to_eV

    # NAC parameters (born effective charges and dielectric tensor)
    self.born = np.array(
        phonopy_file.nac_params.get("born", np.zeros((n_atoms, 3, 3)))
    )
    self.epsilon = np.array(
        phonopy_file.nac_params.get("dielectric", np.identity(3))
    )

    # Create a Structure object
    # At some point should make careful assessment of primitive vs unit_cell
    # PhonoDark uses primitive, but what about when it's different from unit_cell?
    positions = phonopy_file.primitive.scaled_positions
    lattice = (
        np.array(phonopy_file.primitive.cell) * const.Ang_to_inveV * length_factor
    )
    species = phonopy_file.primitive.symbols

    structure = Structure(lattice, species, positions)

    super().__init__(name, properties, structure, m_atoms)

max_dE: float property

Returns omega_ph_max = max(omega_ph) if there are optical modes, otherwise returns the average over the entire Brillouin zone. The quantities are obviously not the same but should be the same order. See theoretical framework paper, paragraph in middle of page 24 (of published version).

TODO: clarify this

Returns:

Name	Type	Description
float	float	the maximum energy deposition

q_cut: float property

The Debye-Waller factor suppresses the rate at larger q beyond q ~ np.sqrt(m_atom * omega_ph). This method calculates an estimate for the cutoff value of q.

Returns:

Name	Type	Description
float	float	The cutoff value of q.

get_eig(k_points, with_eigenvectors=True)

Calculates the phonon frequencies and eigenvectors for the given k-points.

Parameters:

Name	Type	Description	Default
k_points	ArrayLike	Numpy array of k-points in fractional coordinates.	required
with_eigenvectors	bool	Flag indicating whether to calculate eigenvectors.	True

Returns:

Type	Description
Tuple[ndarray, ndarray]	A tuple containing the phonon frequencies and eigenvectors. The phonon frequencies are represented as a numpy array of shape (n_k,n_modes) The eigenvectors are represented as a numpy array of shape (n_k, n_atoms, n_modes, 3) where n_k is the number of k-points, n_modes is the number of modes, n_atoms is the number of atoms, and the last index is for the x, y, z components of the eigenvectors.

Source code in darkmagic/material.py

def get_eig(
    self, k_points: ArrayLike, with_eigenvectors: bool = True
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Calculates the phonon frequencies and eigenvectors for the given k-points.

    Args:
        k_points (ArrayLike): Numpy array of k-points in fractional coordinates.
        with_eigenvectors (bool, optional): Flag indicating whether to calculate eigenvectors.

    Returns:
        A tuple containing the phonon frequencies and eigenvectors.

            * The phonon frequencies are represented as a numpy array of shape (n_k,n_modes)

            * The eigenvectors are represented as a numpy array of shape (n_k, n_atoms, n_modes, 3)

            where n_k is the number of k-points, n_modes is the number of modes,
            n_atoms is the number of atoms, and the last index is
            for the x, y, z components of the eigenvectors.
    """
    # run phonopy in mesh mode
    self.phonopy_file.run_qpoints(k_points, with_eigenvectors=with_eigenvectors)

    mesh_dict = self.phonopy_file.get_qpoints_dict()

    eigenvectors_pre = mesh_dict.get("eigenvectors", None)
    # print(eigenvectors_pre)
    # convert frequencies to correct units
    omega = const.THz_to_eV * mesh_dict["frequencies"]

    eigenvectors = np.zeros(
        (len(k_points), self.n_modes, self.n_atoms, 3), dtype=complex
    )
    # Need to reshape the eigenvectors from (n_k, n_modes, n_modes)
    # to (n_k, n_modes, n_atoms, 3) # TODO: is this correct?
    if with_eigenvectors:
        # TODO: Should rewrite this with a reshape...
        for q in range(len(k_points)):
            for nu in range(self.n_modes):
                eigenvectors[q, nu] = np.array_split(
                    eigenvectors_pre[q].T[nu], self.n_atoms
                )

    return omega, eigenvectors

get_W_tensor(grid)

Computes the W tensor for the given Monkhorst-Pack grid. The W tensor for atom \(j\) is given by: $$ \mathbf{W}j = \frac{\Omega}{4 m_j} \sum\nu \int_\text{1BZ} \frac{d^3k}{(2\pi)^3} \frac{\epsilon_{\nu j \bm{k}} \otimes \epsilon_{\nu j \bm{k}}^*}{\omega_{\nu \bm{k}}} $$

The Debye-Waller factor can be computed from the W tensor as: $$ W_j(\bm{q}) = \bm{q} \cdot (\mathbf{W}_j \bm{q}) $$

Parameters:

Name	Type	Description	Default
grid	MonkhorstPackGrid	The Monkhorst-Pack grid.	required

Returns:

Type	Description
ndarray	np.ndarray: The W tensor.

Source code in darkmagic/material.py

def get_W_tensor(self, grid: MonkhorstPackGrid) -> np.ndarray:
    r"""
    Computes the W tensor for the given Monkhorst-Pack grid. The W tensor for atom $j$ is given by:
    $$
    \mathbf{W}_j = \frac{\Omega}{4 m_j} \sum_\nu \int_\text{1BZ} \frac{d^3k}{(2\pi)^3} \frac{\epsilon_{\nu j \bm{k}} \otimes \epsilon_{\nu j \bm{k}}^*}{\omega_{\nu \bm{k}}}
    $$

    The Debye-Waller factor can be computed from the W tensor as:
    $$
    W_j(\bm{q}) = \bm{q} \cdot (\mathbf{W}_j \bm{q})
    $$

    Args:
        grid (MonkhorstPackGrid): The Monkhorst-Pack grid.

    Returns:
        np.ndarray: The W tensor.

    """
    omega, epsilon = self.get_eig(grid.k_frac)
    # epsilon is (n_k, n_modes, n_atoms, 3)
    eps_tensor = np.einsum("...i,...j->...ij", epsilon, np.conj(epsilon))

    # Sum over all modes and divide by the frequency
    W = (
        1
        / (4 * self.m_atoms[None, :, None, None])
        * np.sum(eps_tensor / omega[..., None, None, None], axis=1)
    )
    # Integrate over the BZ
    return np.sum(W * grid.weights[:, None, None, None], axis=0) / np.sum(
        grid.weights
    )

MagnonMaterial(name, properties, hamiltonian, m_cell, nodmi=False, noaniso=False)

Bases: Material

A class for materials with magnons

Attributes:

Name	Type	Description
hamiltonian	SpinHamiltonian	The spin Hamiltonian of the material.
n_modes	int	The number of magnon modes.
dispersion	MagnonDispersion	The magnon dispersion.

In the current implementation, the hamiltonian only contains the magnetic atoms and their interactions. So m_cell needs to be specified separately

Parameters:

Name	Type	Description	Default
name	str	The name of the material.	required
properties	MaterialParameters	The properties of the material.	required
hamiltonian	SpinHamiltonian	The spin Hamiltonian	required
m_cell	float	the total mass of all ions in the cell	required
nodmi	bool	Whether to include DM interactions.	False
noaniso	bool	Whether to include anisotropic exchange.	False

Source code in darkmagic/material.py

def __init__(
    self,
    name: str,
    properties: MaterialParameters,
    hamiltonian: SpinHamiltonian,
    m_cell: float,
    nodmi: bool = False,
    noaniso: bool = False,
):
    """
    Constructor for a magnon material

    In the current implementation, the hamiltonian only
    contains the magnetic atoms and their interactions.
    So m_cell needs to be specified separately

    Args:
        name: The name of the material.
        properties: The properties of the material.
        hamiltonian: The spin Hamiltonian
        m_cell: the total mass of all ions in the cell
        nodmi: Whether to include DM interactions.
        noaniso: Whether to include anisotropic exchange.
    """
    # Ensure the hamiltonian is in the correct units
    # hamiltonian.cell *= const.Ang_to_inveV
    # In the future we should have a check that it comes in in units of A
    # And convert it here
    self.hamiltonian = hamiltonian
    n_atoms = len(hamiltonian.magnetic_atoms)
    self.n_modes = n_atoms
    self.dispersion = MagnonDispersion(
        hamiltonian, phase_convention="tanner", nodmi=nodmi, noaniso=noaniso
    )

    n_atoms = len(hamiltonian.magnetic_atoms)  # Number of magnetic atoms
    properties.validate_for_magnons(n_atoms)
    # Atom positions in cartesian coordinates (units of 1/eV)
    self.xj = np.array(
        [
            hamiltonian.get_atom_coordinates(atom, relative=False)
            for atom in hamiltonian.magnetic_atoms
        ]
    )
    # sqrt(Sj/2)
    # TODO: make this an internal variable
    self.sqrt_spins_2 = np.sqrt(
        np.array([atom.spin for atom in hamiltonian.magnetic_atoms]) / 2
    )
    # The vectors for rotating to local coordiante system
    # TODO: make this an internal variable
    self.rj = self.dispersion.u

    positions = np.array([a.position for a in hamiltonian.magnetic_atoms])
    lattice = hamiltonian.cell
    species = [a.type for a in hamiltonian.magnetic_atoms]
    structure = Structure(lattice, species, positions)

    m_atoms = [m_cell / n_atoms] * n_atoms
    super().__init__(name, properties, structure, m_atoms)

max_dE: float property

TODO: this needs improvement

Returns the maximum dE possible for the material. For magnons, we estimate this as roughly 3 times the highest magnon frequency at the Brillouin zone (BZ) boundary. If there are no gapped modes at the gamma point, the maximum dE will be 0.

Returns:

Name	Type	Description
float	float	The maximum dE value.

Notes

This calculation should be an average over the Brillouin zone (BZ).

q_cut: float property

For magnons there is no q_cut, so we just set this to a very large number.

Returns:

Name	Type	Description
q_cut	float	a very large number.

get_eig(k, G)

Calculate the eigenvalues and magnon polarization vectors for a given k-point and G-vector.

Parameters:

Name	Type	Description	Default
k	ArrayLike	Single k-point, cartesian coordinates (units of eV).	required
G	ArrayLike	Single G-vector, cartesian coordinates (units of eV).	required

Returns:

Type	Description
Tuple[ndarray, ndarray]	Tuple[np.ndarray, np.ndarray]: A tuple containing the eigenvalues (omega_nu_k) and eigenvectors (epsilon_nu_k_G). - omega_nu_k: (N,) an array of complex numbers representing the eigenvalues in eV. - epsilon_nu_k_G: (N,3) array of complex numbers representing the eigenvectors (magnon polarization vectors) in eV/?? N is the number of magnon modes.

Source code in darkmagic/material.py

def get_eig(self, k: ArrayLike, G: ArrayLike) -> Tuple[np.ndarray, np.ndarray]:
    """
    Calculate the eigenvalues and magnon polarization vectors for a given k-point and G-vector.

    Args:
        k (ArrayLike): Single k-point, cartesian coordinates (units of eV).
        G (ArrayLike): Single G-vector, cartesian coordinates (units of eV).

    Returns:
        Tuple[np.ndarray, np.ndarray]: A tuple containing the eigenvalues (omega_nu_k) and eigenvectors (epsilon_nu_k_G).
            - omega_nu_k: (N,) an array of complex numbers representing the eigenvalues in eV.
            - epsilon_nu_k_G: (N,3) array of complex numbers representing the eigenvectors (magnon polarization vectors) in eV/??
            N is the number of magnon modes.

    """
    # Calculate the prefactor
    prefactor = self.sqrt_spins_2 * np.exp(1j * np.dot(self.xj, G))

    # See Tanner's Disseration and RadTools doc for explanation of the 1/2
    N = self.n_atoms
    omega_nu_k, Tk = self._get_omega_T(self.dispersion.h(k) / 2)
    Uk_conj = np.conjugate(Tk[:N, :N])  # This is U_{j,nu,k})*
    V_minusk = np.conjugate(Tk[N:, :N])  # This is ((V_{j,nu,-k})*)*

    epsilon_nu_k_G = (prefactor[:, None] * V_minusk).T @ np.conjugate(self.rj) + (
        prefactor[:, None] * Uk_conj
    ).T @ self.rj

    return omega_nu_k, epsilon_nu_k_G  # (n,) and (n,3) array of complex numbers