From 003a88df97f9130ddb92d2f23edc76593a034322 Mon Sep 17 00:00:00 2001 From: Janosh Riebesell Date: Sat, 16 Dec 2023 10:32:47 -0800 Subject: [PATCH 1/7] fix typos and improve code readability --- .../electronic_structure/bandstructure.py | 2 +- pymatgen/phonon/bandstructure.py | 75 +++++++++---------- pymatgen/phonon/plotter.py | 1 + 3 files changed, 37 insertions(+), 41 deletions(-) diff --git a/pymatgen/electronic_structure/bandstructure.py b/pymatgen/electronic_structure/bandstructure.py index 95047e48e8e..5dc883f3c2c 100644 --- a/pymatgen/electronic_structure/bandstructure.py +++ b/pymatgen/electronic_structure/bandstructure.py @@ -779,7 +779,7 @@ def get_equivalent_kpoints(self, index): TODO: now it uses the label we might want to use coordinates instead (in case there was a mislabel) """ - # if the kpoint has no label it can"t have a repetition along the band + # if the kpoint has no label it can't have a repetition along the band # structure line object if self.kpoints[index].label is None: diff --git a/pymatgen/phonon/bandstructure.py b/pymatgen/phonon/bandstructure.py index c052e00ee8e..0ac8d671ed9 100644 --- a/pymatgen/phonon/bandstructure.py +++ b/pymatgen/phonon/bandstructure.py @@ -24,13 +24,13 @@ def get_reasonable_repetitions(n_atoms: int) -> tuple[int, int, int]: according to the number of atoms in the system. """ if n_atoms < 4: - return (3, 3, 3) + return 3, 3, 3 if 4 <= n_atoms < 15: - return (2, 2, 2) + return 2, 2, 2 if 15 <= n_atoms < 50: - return (2, 2, 1) + return 2, 2, 1 - return (1, 1, 1) + return 1, 1, 1 def eigenvectors_from_displacements(disp, masses) -> np.ndarray: @@ -137,8 +137,8 @@ def __init__( self.nb_qpoints = len(self.qpoints) # normalize directions for nac_frequencies and nac_eigendisplacements - self.nac_frequencies = [] - self.nac_eigendisplacements = [] + self.nac_frequencies: list[tuple[list[float], np.ndarray]] = [] + self.nac_eigendisplacements: list[tuple[list[float], np.ndarray]] = [] if nac_frequencies is not None: for freq in nac_frequencies: self.nac_frequencies.append(([idx / np.linalg.norm(freq[0]) for idx in freq[0]], freq[1])) @@ -177,10 +177,10 @@ def get_nac_frequencies_along_dir(self, direction: Sequence) -> np.ndarray | Non the frequencies as a numpy array o(3*len(structure), len(qpoints)). None if not found. """ - versor = [i / np.linalg.norm(direction) for i in direction] - for d, f in self.nac_frequencies: - if np.allclose(versor, d): - return f + versor = [idx / np.linalg.norm(direction) for idx in direction] + for dist, freq in self.nac_frequencies: + if np.allclose(versor, dist): + return freq return None @@ -195,10 +195,10 @@ def get_nac_eigendisplacements_along_dir(self, direction) -> np.ndarray | None: the eigendisplacements as a numpy array of complex numbers with shape (3*len(structure), len(structure), 3). None if not found. """ - versor = [i / np.linalg.norm(direction) for i in direction] - for d, e in self.nac_eigendisplacements: - if np.allclose(versor, d): - return e + versor = [idx / np.linalg.norm(direction) for idx in direction] + for dist, eigen_disp in self.nac_eigendisplacements: + if np.allclose(versor, dist): + return eigen_disp return None @@ -426,45 +426,40 @@ def get_equivalent_qpoints(self, index: int) -> list[int]: TODO: now it uses the label we might want to use coordinates instead (in case there was a mislabel) """ - # if the qpoint has no label it can"t have a repetition along the band + # if the qpoint has no label it can't have a repetition along the band # structure line object if self.qpoints[index].label is None: return [index] list_index_qpoints = [] - for i in range(self.nb_qpoints): - if self.qpoints[i].label == self.qpoints[index].label: - list_index_qpoints.append(i) + for idx in range(self.nb_qpoints): + if self.qpoints[idx].label == self.qpoints[index].label: + list_index_qpoints.append(idx) return list_index_qpoints - def get_branch(self, index: int) -> list[dict]: - r"""Returns in what branch(es) is the qpoint. There can be several - branches. + def get_branch(self, index: int) -> list[dict[str, str | int]]: + r"""Returns in what branch(es) is the qpoint. There can be several branches. Args: - index: the qpoint index + index (int): the qpoint index Returns: - A list of dictionaries [{"name","start_index","end_index","index"}] - indicating all branches in which the qpoint is. It takes into - account the fact that one qpoint (e.g., \\Gamma) can be in several - branches + list[dict[str, str | int]]: [{"name","start_index","end_index","index"}] + indicating all branches in which the qpoint is. It takes into + account the fact that one qpoint (e.g., \\Gamma) can be in several + branches """ - to_return = [] - for i in self.get_equivalent_qpoints(index): - for b in self.branches: - if b["start_index"] <= i <= b["end_index"]: - to_return.append( - { - "name": b["name"], - "start_index": b["start_index"], - "end_index": b["end_index"], - "index": i, - } + lst = [] + for pt_idx in self.get_equivalent_qpoints(index): + for branch in self.branches: + start_idx, end_idx = branch["start_index"], branch["end_index"] + if start_idx <= pt_idx <= end_idx: + lst.append( + {"name": branch["name"], "start_index": start_idx, "end_index": end_idx, "index": pt_idx} ) - return to_return + return lst def write_phononwebsite(self, filename: str | PathLike) -> None: """Write a json file for the phononwebsite: @@ -606,7 +601,7 @@ def from_dict(cls, dct: dict) -> PhononBandStructureSymmLine: eigendisplacements = ( np.array(dct["eigendisplacements"]["real"]) + np.array(dct["eigendisplacements"]["imag"]) * 1j ) - structure = Structure.from_dict(dct["structure"]) if "structure" in dct else None + struct = Structure.from_dict(dct["structure"]) if "structure" in dct else None return cls( dct["qpoints"], np.array(dct["bands"]), @@ -614,5 +609,5 @@ def from_dict(cls, dct: dict) -> PhononBandStructureSymmLine: dct["has_nac"], eigendisplacements, dct["labels_dict"], - structure=structure, + structure=struct, ) diff --git a/pymatgen/phonon/plotter.py b/pymatgen/phonon/plotter.py index 0cbc5e32eab..4f0f3c024b0 100644 --- a/pymatgen/phonon/plotter.py +++ b/pymatgen/phonon/plotter.py @@ -285,6 +285,7 @@ def _make_ticks(self, ax: Axes) -> Axes: ax.set_xticks(ticks_labels[0]) ax.set_xticklabels(ticks_labels[1]) + # plot vertical lines at each of the ticks for idx, label in enumerate(ticks["label"]): if label is not None: ax.axvline(ticks["distance"][idx], color="black") From 32fe1ecd874b9e0a1253caee9dd6089aaa2df84e Mon Sep 17 00:00:00 2001 From: Janosh Riebesell Date: Sat, 16 Dec 2023 10:37:59 -0800 Subject: [PATCH 2/7] clarify has_nac doc str --- pymatgen/phonon/bandstructure.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pymatgen/phonon/bandstructure.py b/pymatgen/phonon/bandstructure.py index 0ac8d671ed9..7df000c7f94 100644 --- a/pymatgen/phonon/bandstructure.py +++ b/pymatgen/phonon/bandstructure.py @@ -158,7 +158,9 @@ def has_imaginary_freq(self, tol: float = 1e-5) -> bool: @property def has_nac(self) -> bool: - """True if nac_frequencies are present.""" + """True if nac_frequencies are present (i.e. the band structure has been + calculated taking into account Born-charge-derived non-analytical corrections at Gamma). + """ return len(self.nac_frequencies) > 0 @property From ff94e58819112ad4d85959762f6d43f82aa1b18b Mon Sep 17 00:00:00 2001 From: Janosh Riebesell Date: Sat, 16 Dec 2023 10:38:29 -0800 Subject: [PATCH 3/7] add PhononBandStructure.has_imaginary_gamma_freq() to check for imaginary modes at the gamma point --- pymatgen/phonon/bandstructure.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/pymatgen/phonon/bandstructure.py b/pymatgen/phonon/bandstructure.py index 7df000c7f94..534b2a7ee2e 100644 --- a/pymatgen/phonon/bandstructure.py +++ b/pymatgen/phonon/bandstructure.py @@ -156,6 +156,18 @@ def has_imaginary_freq(self, tol: float = 1e-5) -> bool: """True if imaginary frequencies are present in the BS.""" return self.min_freq()[1] + tol < 0 + def has_imaginary_gamma_freq(self, tol: float = 1e-5) -> bool: + """Checks if there are imaginary modes specifically at the gamma point. + + Args: + tol: Tolerance for determining if a frequency is imaginary. + + Returns: + True if there are imaginary modes at the gamma point, False otherwise. + """ + gamma_freqs = self.bands[:, 0] # frequencies at the Gamma point + return any(freq < -tol for freq in gamma_freqs) + @property def has_nac(self) -> bool: """True if nac_frequencies are present (i.e. the band structure has been From fd408b85ff0950d603e63d91d6b57bbf054cdd45 Mon Sep 17 00:00:00 2001 From: Janosh Riebesell Date: Sat, 16 Dec 2023 10:39:03 -0800 Subject: [PATCH 4/7] improve tests for phBS.has_imaginary_freq and add new ones for has_imaginary_gamma_freq --- tests/phonon/test_bandstructure.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/tests/phonon/test_bandstructure.py b/tests/phonon/test_bandstructure.py index a82f115bbec..e26d9ce1801 100644 --- a/tests/phonon/test_bandstructure.py +++ b/tests/phonon/test_bandstructure.py @@ -38,8 +38,6 @@ def test_basic(self): assert list(self.bs.min_freq()[0].frac_coords) == [0, 0, 0] assert self.bs.min_freq()[1] == approx(-0.03700895020) - assert self.bs.has_imaginary_freq() - assert not self.bs.has_imaginary_freq(tol=0.5) assert_allclose(self.bs.asr_breaking(), [-0.0370089502, -0.0370089502, -0.0221388897]) assert self.bs.nb_bands == 6 @@ -47,6 +45,20 @@ def test_basic(self): assert_allclose(self.bs.qpoints[1].frac_coords, [0.01, 0, 0]) + def test_has_imaginary_freq(self): + for tol in (0, 0.01, 0.02, 0.03, 0.04, 0.05): + assert self.bs.has_imaginary_freq(tol=tol) == (tol < 0.04) + + for tol in (0, 0.01, 0.02, 0.03, 0.04, 0.05): + assert self.bs2.has_imaginary_freq(tol=tol) == (tol < 0.01) + + # test Gamma point imaginary frequency detection + for tol in (0, 0.01, 0.02, 0.03, 0.04, 0.05): + assert self.bs.has_imaginary_gamma_freq(tol=tol) == (tol < 0.04) + + for tol in (0, 0.01, 0.02, 0.03, 0.04, 0.05): + assert self.bs2.has_imaginary_gamma_freq(tol=tol) == (tol < 0.01) + def test_nac(self): assert self.bs.has_nac assert not self.bs2.has_nac From 02b557980d19b00c401409eb6259a168b0f687d4 Mon Sep 17 00:00:00 2001 From: Janosh Riebesell Date: Sat, 16 Dec 2023 10:46:59 -0800 Subject: [PATCH 5/7] add PhononDos.mae() method to compare two DOSs --- pymatgen/phonon/dos.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/pymatgen/phonon/dos.py b/pymatgen/phonon/dos.py index c08e7f8aca7..990719f6454 100644 --- a/pymatgen/phonon/dos.py +++ b/pymatgen/phonon/dos.py @@ -329,6 +329,28 @@ def zero_point_energy(self, structure: Structure | None = None) -> float: return zpe + def mae(self, other: PhononDos, two_sided: bool = True) -> float: + """Mean absolute error between two DOSs. + + Args: + other: Another DOS object. + two_sided: Whether to calculate the two-sided MAE meaning interpolate each DOS to the + other's frequencies and averaging the two MAEs. Defaults to True. + + Returns: + float: Mean absolute error. + """ + # Interpolate other.densities to align with self.frequencies + self_interpolated = np.interp(self.frequencies, other.frequencies, other.densities) + self_mae = np.abs(self.densities - self_interpolated).mean() + + if two_sided: + other_interpolated = np.interp(other.frequencies, self.frequencies, self.densities) + other_mae = np.abs(other.densities - other_interpolated).mean() + return (self_mae + other_mae) / 2 + + return self_mae + class CompletePhononDos(PhononDos): """This wrapper class defines a total dos, and also provides a list of PDos. From c96a163c572615c58919fbd8c3b0f0896fbe1af3 Mon Sep 17 00:00:00 2001 From: Janosh Riebesell Date: Sat, 16 Dec 2023 10:47:17 -0800 Subject: [PATCH 6/7] add TestPhononDos.test_mae() --- tests/phonon/test_dos.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/tests/phonon/test_dos.py b/tests/phonon/test_dos.py index 9c1d7886944..a94271ebaf2 100644 --- a/tests/phonon/test_dos.py +++ b/tests/phonon/test_dos.py @@ -3,6 +3,7 @@ import json import re +import pytest from pytest import approx from pymatgen.core import Element @@ -85,6 +86,19 @@ def test_eq(self): assert self.dos != 2 * self.dos assert 2 * self.dos == self.dos + self.dos + def test_mae(self): + assert self.dos.mae(self.dos) == 0 + assert self.dos.mae(self.dos + 1) == 1 + assert self.dos.mae(self.dos - 1) == 1 + assert self.dos.mae(2 * self.dos) == pytest.approx(0.786546967) + assert (2 * self.dos).mae(self.dos) == pytest.approx(0.786546967) + + # test two_sided=False after shifting DOS freqs so MAE requires interpolation + dos2 = PhononDos(self.dos.frequencies + 0.01, self.dos.densities) + assert self.dos.mae(dos2 + 1, two_sided=False) == pytest.approx(0.999999999) + assert self.dos.mae(dos2 - 1, two_sided=False) == pytest.approx(1.00000000000031) + assert self.dos.mae(2 * dos2, two_sided=False) == pytest.approx(0.786546967) + class TestCompletePhononDos(PymatgenTest): def setUp(self): From 9c59b59c660d48cd0542a69ee50b87e39d48f36f Mon Sep 17 00:00:00 2001 From: Janosh Riebesell Date: Sat, 16 Dec 2023 11:44:00 -0800 Subject: [PATCH 7/7] increase default has_imaginary_freq tol from 1e-5 to 1e-3 --- pymatgen/phonon/bandstructure.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/pymatgen/phonon/bandstructure.py b/pymatgen/phonon/bandstructure.py index 534b2a7ee2e..52d60cc3a08 100644 --- a/pymatgen/phonon/bandstructure.py +++ b/pymatgen/phonon/bandstructure.py @@ -152,18 +152,20 @@ def min_freq(self) -> tuple[Kpoint, float]: return self.qpoints[idx[1]], self.bands[idx] - def has_imaginary_freq(self, tol: float = 1e-5) -> bool: - """True if imaginary frequencies are present in the BS.""" + def has_imaginary_freq(self, tol: float = 1e-3) -> bool: + """True if imaginary frequencies are present anywhere in the band structure. Always True if + has_imaginary_gamma_freq is True. + + Args: + tol: Tolerance for determining if a frequency is imaginary. Defaults to 1e-3. + """ return self.min_freq()[1] + tol < 0 - def has_imaginary_gamma_freq(self, tol: float = 1e-5) -> bool: - """Checks if there are imaginary modes specifically at the gamma point. + def has_imaginary_gamma_freq(self, tol: float = 1e-3) -> bool: + """Checks if there are imaginary modes at the gamma point. Args: - tol: Tolerance for determining if a frequency is imaginary. - - Returns: - True if there are imaginary modes at the gamma point, False otherwise. + tol: Tolerance for determining if a frequency is imaginary. Defaults to 1e-3. """ gamma_freqs = self.bands[:, 0] # frequencies at the Gamma point return any(freq < -tol for freq in gamma_freqs)