Merge 836caa5 into ec611b9

pymatgen/io/phonopy.py

@@ -5,17 +5,20 @@
 """
 Module for interfacing with phonopy, see https://atztogo.github.io/phonopy/
 """
+from typing import Dict, List, Any
 
 import numpy as np
 from monty.dev import requires
 from monty.serialization import loadfn
+from scipy.interpolate import InterpolatedUnivariateSpline
 
 from pymatgen.core import Lattice, Structure
 from pymatgen.phonon.bandstructure import (
     PhononBandStructure,
     PhononBandStructureSymmLine,
 )
 from pymatgen.phonon.dos import CompletePhononDos, PhononDos
+from pymatgen.phonon.gruneisen import GruneisenParameter, GruneisenPhononBandStructureSymmLine
 from pymatgen.symmetry.bandstructure import HighSymmKpath
 
 try:
@@ -24,9 +27,11 @@
     from phonopy.structure.atoms import PhonopyAtoms
 except ImportError:
     Phonopy = None
+    write_disp_yaml = None
+    PhonopyAtoms = None
 
 
-@requires(Phonopy, "phonopy not installed!")
+@requires(Phonopy, "phonopy not installed!")  # type: ignore
 def get_pmg_structure(phonopy_structure):
     """
     Convert a PhonopyAtoms object to pymatgen Structure object.
@@ -51,7 +56,7 @@ def get_pmg_structure(phonopy_structure):
     )
 
 
-@requires(Phonopy, "phonopy not installed!")
+@requires(Phonopy, "phonopy not installed!")  # type: ignore
 def get_phonopy_structure(pmg_structure):
     """
     Convert a pymatgen Structure object to a PhonopyAtoms object.
@@ -430,3 +435,246 @@ def get_phonon_band_structure_symm_line_from_fc(
     labels_dict = {a: k for a, k in zip(labels, kpoints) if a != ""}
 
     return PhononBandStructureSymmLine(kpoints, frequencies, structure.lattice, labels_dict=labels_dict)
+
+
+def get_gruneisenparameter(gruneisen_path, structure=None, structure_path=None) -> GruneisenParameter:
+    """
+    Get Gruneisen object from gruneisen.yaml file, as obtained from phonopy (Frequencies in THz!).
+    The order is structure > structure path > structure from gruneisen dict.
+    Newer versions of phonopy include the structure in the yaml file,
+    the structure/structure_path is kept for compatibility.
+
+    Args:
+        gruneisen_path: Path to gruneisen.yaml file (frequencies have to be in THz!)
+        structure: pymatgen Structure object
+        structure_path: path to structure in a file (e.g., POSCAR)
+
+    Returns: GruneisenParameter object
+
+    """
+
+    gruneisen_dict = loadfn(gruneisen_path)
+
+    if structure_path and structure is None:
+        structure = Structure.from_file(structure_path)
+    else:
+        try:
+            structure = get_structure_from_dict(gruneisen_dict)
+        except ValueError:
+            raise ValueError("\nPlease provide a structure.\n")
+
+    qpts, multiplicities, frequencies, gruneisen = ([] for _ in range(4))
+    phonopy_labels_dict = {}
+
+    for p in gruneisen_dict["phonon"]:
+        q = p["q-position"]
+        qpts.append(q)
+        if "multiplicity" in p:
+            m = p["multiplicity"]
+        else:
+            m = 1
+        multiplicities.append(m)
+        bands, gruneisenband = ([] for _ in range(2))
+        for b in p["band"]:
+            bands.append(b["frequency"])
+            if "gruneisen" in b:
+                gruneisenband.append(b["gruneisen"])
+        frequencies.append(bands)
+        gruneisen.append(gruneisenband)
+        if "label" in p:
+            phonopy_labels_dict[p["label"]] = p["q-position"]
+
+    qpts_np = np.array(qpts)
+    multiplicities_np = np.array(multiplicities)
+    # transpose to match the convention in PhononBandStructure
+    frequencies_np = np.transpose(frequencies)
+    gruneisen_np = np.transpose(gruneisen)
+
+    return GruneisenParameter(
+        gruneisen=gruneisen_np,
+        qpoints=qpts_np,
+        multiplicities=multiplicities_np,
+        frequencies=frequencies_np,
+        structure=structure,
+    )
+
+
+def get_gs_ph_bs_symm_line_from_dict(
+    gruneisen_dict, structure=None, structure_path=None, labels_dict=None, fit=False
+) -> GruneisenPhononBandStructureSymmLine:
+    r"""
+    Creates a pymatgen GruneisenPhononBandStructure object from the dictionary
+    extracted by the gruneisen.yaml file produced by phonopy. The labels
+    will be extracted from the dictionary, if present. If the 'eigenvector'
+    key is found the eigendisplacements will be calculated according to the
+    formula::
+
+        exp(2*pi*i*(frac_coords \\dot q) / sqrt(mass) * v
+
+    and added to the object. A fit algorithm can be used to replace diverging
+    Gruneisen values close to gamma.
+
+    Args:
+        gruneisen_dict (dict): the dictionary extracted from the gruneisen.yaml file
+        structure (Structure): pymatgen structure object
+        structure_path: path to structure file
+        labels_dict (dict): dict that links a qpoint in frac coords to a label.
+            Its value will replace the data contained in the band.yaml.
+        fit (bool): Substitute Grueneisen parameters close to the gamma point
+            with points obtained from a fit to a spline if the derivate from
+            a smooth curve (i.e. if the slope changes by more than 200% in the
+            range of 10% around the gamma point).
+            These derivations occur because of very small frequencies
+            (and therefore numerical inaccuracies) close to gamma.
+    """
+
+    if structure_path and structure is None:
+        structure = Structure.from_file(structure_path)
+    else:
+        try:
+            structure = get_structure_from_dict(gruneisen_dict)
+        except ValueError:
+            raise ValueError("\nPlease provide a structure.\n")
+
+    qpts, frequencies, gruneisenparameters = ([] for _ in range(3))
+    phonopy_labels_dict = {}  # type: Dict[Any,Any]
+
+    if fit:
+        for pa in gruneisen_dict["path"]:
+            phonon = pa["phonon"]  # This is a list
+            start = pa["phonon"][0]
+            end = pa["phonon"][-1]
+
+            if start["q-position"] == [0, 0, 0]:  # Gamma at start of band
+                qpts_temp, frequencies_temp, gruneisen_temp, distance = (
+                    [] for _ in range(4)
+                )  # type: List[Any],List[Any],List[Any],List[Any]
+                for i in range(pa["nqpoint"]):
+                    bands, gruneisenband = ([] for _ in range(2))  # type: List[Any], List[Any]
+                    for b in phonon[pa["nqpoint"] - i - 1]["band"]:
+                        bands.append(b["frequency"])
+                        # Fraction of leftover points in current band
+                        gruen = _extrapolate_grun(b, distance, gruneisen_temp, gruneisenband, i, pa)
+                        gruneisenband.append(gruen)
+                    q = phonon[pa["nqpoint"] - i - 1]["q-position"]
+                    qpts_temp.append(q)
+                    d = phonon[pa["nqpoint"] - i - 1]["distance"]
+                    distance.append(d)
+                    frequencies_temp.append(bands)
+                    gruneisen_temp.append(gruneisenband)
+                    if "label" in phonon[pa["nqpoint"] - i - 1]:
+                        phonopy_labels_dict[phonon[pa["nqpoint"] - i - 1]]["label"] = phonon[pa["nqpoint"] - i - 1][
+                            "q-position"
+                        ]
+
+                qpts.extend(list(reversed(qpts_temp)))
+                frequencies.extend(list(reversed(frequencies_temp)))
+                gruneisenparameters.extend(list(reversed(gruneisen_temp)))
+
+            elif end["q-position"] == [0, 0, 0]:  # Gamma at end of band
+                distance = []
+                for i in range(pa["nqpoint"]):
+                    bands, gruneisenband = ([] for _ in range(2))
+                    for b in phonon[i]["band"]:
+                        bands.append(b["frequency"])
+                        gruen = _extrapolate_grun(b, distance, gruneisenparameters, gruneisenband, i, pa)
+                        gruneisenband.append(gruen)
+                    q = phonon[i]["q-position"]
+                    qpts.append(q)
+                    d = phonon[i]["distance"]
+                    distance.append(d)
+                    frequencies.append(bands)
+                    gruneisenparameters.append(gruneisenband)
+                    if "label" in phonon[i]:
+                        phonopy_labels_dict[phonon[i]["label"]] = phonon[i]["q-position"]
+
+            else:  # No Gamma in band
+                distance = []
+                for i in range(pa["nqpoint"]):
+                    bands, gruneisenband = ([] for _ in range(2))
+                    for b in phonon[i]["band"]:
+                        bands.append(b["frequency"])
+                        gruneisenband.append(b["gruneisen"])
+                    q = phonon[i]["q-position"]
+                    qpts.append(q)
+                    d = phonon[i]["distance"]
+                    distance.append(d)
+                    frequencies.append(bands)
+                    gruneisenparameters.append(gruneisenband)
+                    if "label" in phonon[i]:
+                        phonopy_labels_dict[phonon[i]["label"]] = phonon[i]["q-position"]
+
+    else:
+        for pa in gruneisen_dict["path"]:
+            for p in pa["phonon"]:
+                q = p["q-position"]
+                qpts.append(q)
+                bands, gruneisen_bands = ([] for _ in range(2))
+                for b in p["band"]:
+                    bands.append(b["frequency"])
+                    gruneisen_bands.append(b["gruneisen"])
+                frequencies.append(bands)
+                gruneisenparameters.append(gruneisen_bands)
+                if "label" in p:
+                    phonopy_labels_dict[p["label"]] = p["q-position"]
+
+    qpts_np = np.array(qpts)
+    # transpose to match the convention in PhononBandStructure
+    frequencies_np = np.transpose(frequencies)
+    gruneisenparameters_np = np.transpose(gruneisenparameters)
+
+    rec_latt = structure.lattice.reciprocal_lattice
+    labels_dict = labels_dict or phonopy_labels_dict
+    return GruneisenPhononBandStructureSymmLine(
+        qpoints=qpts_np,
+        frequencies=frequencies_np,
+        gruneisenparameters=gruneisenparameters_np,
+        lattice=rec_latt,
+        labels_dict=labels_dict,
+        structure=structure,
+        eigendisplacements=None,
+    )
+
+
+def _extrapolate_grun(b, distance, gruneisenparameter, gruneisenband, i, pa):
+    leftover_fraction = (pa["nqpoint"] - i - 1) / pa["nqpoint"]
+    if leftover_fraction < 0.1:
+        diff = abs(b["gruneisen"] - gruneisenparameter[-1][len(gruneisenband)]) / abs(
+            gruneisenparameter[-2][len(gruneisenband)] - gruneisenparameter[-1][len(gruneisenband)]
+        )
+        if diff > 2:
+            x = list(range(len(distance)))
+            y = [i[len(gruneisenband)] for i in gruneisenparameter]
+            y = y[-len(x) :]  # Only elements of current band
+            extrapolator = InterpolatedUnivariateSpline(x, y, k=5)
+            g_extrapolated = extrapolator(len(distance))
+            gruen = float(g_extrapolated)
+        else:
+            gruen = b["gruneisen"]
+    else:
+        gruen = b["gruneisen"]
+    return gruen
+
+
+def get_gruneisen_ph_bs_symm_line(gruneisen_path, structure=None, structure_path=None, labels_dict=None, fit=False):
+    r"""
+    Creates a pymatgen GruneisenPhononBandStructure from a band.yaml file.
+    The labels will be extracted from the dictionary, if present.
+    If the 'eigenvector' key is found the eigendisplacements will be
+    calculated according to the formula:
+    \\exp(2*pi*i*(frac_coords \\dot q) / sqrt(mass) * v
+     and added to the object.
+
+    Args:
+        gruneisen_path: path to the band.yaml file
+        structure: pymaten Structure object
+        structure_path: path to a structure file (e.g., POSCAR)
+        labels_dict: dict that links a qpoint in frac coords to a label.
+        fit: Substitute Grueneisen parameters close to the gamma point
+            with points obtained from a fit to a spline if the derivate from
+            a smooth curve (i.e. if the slope changes by more than 200% in the
+            range of 10% around the gamma point).
+            These derivations occur because of very small frequencies
+            (and therefore numerical inaccuracies) close to gamma.
+    """
+    return get_gs_ph_bs_symm_line_from_dict(loadfn(gruneisen_path), structure, structure_path, labels_dict, fit)

pymatgen/io/tests/test_phonopy.py

@@ -1,5 +1,4 @@
 import os
-import sys
 import unittest
 from pathlib import Path
 
@@ -17,7 +16,6 @@
     print(ex)
     Phonopy = None
 
-
 test_dir = os.path.join(PymatgenTest.TEST_FILES_DIR, "phonopy")
 
 
@@ -181,5 +179,43 @@ def test_get_phonon_band_structure_symm_line_from_fc(self):
         self.assertAlmostEqual(bs.bands[2][10], 2.869229797603161)
 
 
+@unittest.skipIf(Phonopy is None, "Phonopy not present")
+class TestGruneisen(unittest.TestCase):
+    def test_ph_bs_symm_line(self):
+        self.bs_symm_line_1 = get_gruneisen_ph_bs_symm_line(
+            gruneisen_path=os.path.join(PymatgenTest.TEST_FILES_DIR, "gruneisen/gruneisen_band_Si.yaml"),
+            structure_path=os.path.join(PymatgenTest.TEST_FILES_DIR, "gruneisen/eq/POSCAR_Si"),
+            fit=True,
+        )
+        self.bs_symm_line_2 = get_gruneisen_ph_bs_symm_line(
+            gruneisen_path=os.path.join(PymatgenTest.TEST_FILES_DIR, "gruneisen/gruneisen_band_Si.yaml"),
+            structure_path=os.path.join(PymatgenTest.TEST_FILES_DIR, "gruneisen/eq/POSCAR_Si"),
+            fit=False,
+        )
+
+        # check if a bit of the gruneisen parameters happens
+
+        self.assertNotEqual(self.bs_symm_line_1.gruneisen[0][0], self.bs_symm_line_2.gruneisen[0][0])
+        with self.assertRaises(ValueError):
+            self.bs_symm_line_2 = get_gruneisen_ph_bs_symm_line(
+                gruneisen_path=os.path.join(PymatgenTest.TEST_FILES_DIR, "gruneisen/gruneisen_eq_plus_minus_InP.yaml")
+            )
+
+    def test_gruneisen_parameter(self):
+        self.gruneisenobject_Si = get_gruneisenparameter(
+            os.path.join(PymatgenTest.TEST_FILES_DIR, "gruneisen/gruneisen_mesh_Si.yaml"),
+            structure_path=os.path.join(PymatgenTest.TEST_FILES_DIR, "gruneisen/eq/POSCAR_Si"),
+        )
+
+        self.assertAlmostEqual(self.gruneisenobject_Si.frequencies[0][0], 0.2523831291)
+        self.assertAlmostEqual(self.gruneisenobject_Si.gruneisen[0][0], -0.1190736091)
+
+        # catch the exception when no structure is present
+        with self.assertRaises(ValueError):
+            get_gruneisenparameter(
+                os.path.join(PymatgenTest.TEST_FILES_DIR, "gruneisen/gruneisen_mesh_InP_without_struct.yaml")
+            )
+
+
 if __name__ == "__main__":
     unittest.main()