extract dielectrics

phonopy_vibspec/phonons_analyzer.py

@@ -17,6 +17,8 @@
 
 
 class PhononsAnalyzer:
+    """Use Phonopy to extract phonon frequencies and eigenmodes, as well as irreps
+    """
     def __init__(self, phonon: phonopy.Phonopy):
         self.phonotopy = phonon
         self.supercell = phonon.supercell
@@ -167,6 +169,7 @@ def prepare_raman(
 
         calculator = RamanSpectrum(
             # input
+            cell_volume=base_geometry.volume,
             modes=modes,
             frequencies=frequencies,
             irrep_labels=irrep_labels,

phonopy_vibspec/spectra.py

@@ -1,15 +1,72 @@
 import pathlib
-from typing import List, Optional
-
 import h5py
 import numpy
 from numpy.typing import NDArray
+from typing import List, Optional
+
+from phonopy.interface.vasp import VasprunxmlExpat
+
+
+class InfraredSpectrum:
+    def __init__(
+        self,
+        # input
+        modes: List[int],
+        frequencies: List[float],
+        dmu_dq: NDArray[float],
+        irrep_labels: Optional[List[str]] = None,
+    ):
+        assert irrep_labels is None or len(irrep_labels) == len(modes)
+        assert len(frequencies) == len(modes)
+        assert dmu_dq.shape[0] == len(modes)
+
+        self.modes = modes
+        self.frequencies = frequencies
+        self.irrep_labels = irrep_labels
+        self.dmu_dq = dmu_dq
+
+    def __len__(self):
+        return len(self.modes)
+
+    @classmethod
+    def from_hdf5(cls, path: pathlib.Path):
+        with h5py.File(path) as f:
+            assert f.attrs['spectrum'] == 'IR'
+
+            if f.attrs['version'] > 1:
+                raise Exception('unsupported version')
+
+            if 'input' not in f:
+                raise Exception('missing `input` group')
+
+            modes = f['input/modes'][:]
+            frequencies = f['input/frequencies'][:]
+            dmu_dq = f['input/dmu_dq'][:]
+            irrep_labels = None
+            if 'input/irrep_labels' in f:
+                irrep_labels = f['input/irrep_labels'][:]
+
+            return cls(modes=modes, frequencies=frequencies, irrep_labels=irrep_labels, dmu_dq=dmu_dq)
+
+    def to_hdf5(self, path: pathlib.Path):
+        with h5py.File(path, 'w') as f:
+            f.attrs['spectrum'] = 'IR'
+            f.attrs['version'] = 1
+
+            g_input = f.create_group('input')
+            g_input.create_dataset('modes', data=self.modes)
+            g_input.create_dataset('frequencies', data=self.frequencies)
+            g_input.create_dataset('dmu_dq', data=self.dmu_dq)
+
+    def compute_intensities(self) -> NDArray[float]:
+        return (self.dmu_dq ** 2).sum(axis=1)
 
 
 class RamanSpectrum:
     def __init__(
         self,
         # input
+        cell_volume: float,
         modes: List[int],
         frequencies: List[float],
         irrep_labels: Optional[List[str]] = None,
@@ -29,6 +86,7 @@ def __init__(
         assert nd_steps is None or len(nd_steps) == len(nd_modes)
         assert nd_dielectrics is None or nd_dielectrics.shape[0] == len(nd_modes)
 
+        self.cell_volume = cell_volume
         self.modes = modes
         self.frequencies = frequencies
         self.irrep_labels = irrep_labels
@@ -49,6 +107,7 @@ def from_hdf5(cls, path: pathlib.Path):
             if 'input' not in f:
                 raise Exception('missing `input` group')
 
+            cell_volume = f['input/cell_volume'][0]
             modes = f['input/modes'][:]
             frequencies = f['input/frequencies'][:]
             irrep_labels = None
@@ -74,6 +133,7 @@ def from_hdf5(cls, path: pathlib.Path):
                     nd_dielectrics = g_nd['dielectrics'][:]
 
             return cls(
+                cell_volume=cell_volume,
                 modes=modes, frequencies=frequencies, irrep_labels=irrep_labels, dalpha_dq=dalpha_dq,
                 nd_stencil=nd_stencil, nd_mode_equiv=nd_mode_equiv, nd_modes=nd_modes, nd_steps=nd_steps,
                 nd_dielectrics=nd_dielectrics
@@ -85,6 +145,7 @@ def to_hdf5(self, path: pathlib.Path):
             f.attrs['version'] = 1
 
             g_input = f.create_group('input')
+            g_input.create_dataset('cell_volume', data=[self.cell_volume])
             g_input.create_dataset('modes', data=self.modes)
             g_input.create_dataset('frequencies', data=self.frequencies)
 
@@ -111,57 +172,48 @@ def to_hdf5(self, path: pathlib.Path):
             if self.dielectrics is not None:
                 g_nd.create_dataset('dielectrics', data=self.dielectrics)
 
+    def extract_dielectrics(self, files: List[pathlib.Path]):
+        """Extract dielectric tensors from `files`.
+        Assumes that the files are in the EXACT same order as created (i.e., increasing mode and step).
+        """
 
-class InfraredSpectrum:
-    def __init__(
-        self,
-        # input
-        modes: List[int],
-        frequencies: List[float],
-        dmu_dq: NDArray[float],
-        irrep_labels: Optional[List[str]] = None,
-    ):
-        assert irrep_labels is None or len(irrep_labels) == len(modes)
-        assert len(frequencies) == len(modes)
-        assert dmu_dq.shape[0] == len(modes)
+        nsteps = len(self.nd_stencil)
+        assert len(files) == nsteps * len(self.nd_modes)
 
-        self.modes = modes
-        self.frequencies = frequencies
-        self.irrep_labels = irrep_labels
-        self.dmu_dq = dmu_dq
+        dielectrics = numpy.zeros((len(self.nd_modes), 2, 3, 3))
 
-    def __len__(self):
-        return len(self.modes)
+        for i, path in enumerate(files):
+            imode = i // nsteps
+            step = i % nsteps
 
-    @classmethod
-    def from_hdf5(cls, path: pathlib.Path):
-        with h5py.File(path) as f:
-            assert f.attrs['spectrum'] == 'IR'
+            with path.open('rb') as f:
+                vasprun = VasprunxmlExpat(f)
+                vasprun.parse()
 
-            if f.attrs['version'] > 1:
-                raise Exception('unsupported version')
+                dielectric_tensor = vasprun.epsilon
+                if dielectric_tensor is None:
+                    raise Exception('no dielectric tensor in `{}`'.format(path))
 
-            if 'input' not in f:
-                raise Exception('missing `input` group')
+                dielectrics[imode, step] = dielectric_tensor
 
-            modes = f['input/modes'][:]
-            frequencies = f['input/frequencies'][:]
-            dmu_dq = f['input/dmu_dq'][:]
-            irrep_labels = None
-            if 'input/irrep_labels' in f:
-                irrep_labels = f['input/irrep_labels'][:]
+        self.dielectrics = dielectrics
+        self._compute_dalpha_dq()
 
-            return cls(modes=modes, frequencies=frequencies, irrep_labels=irrep_labels, dmu_dq=dmu_dq)
+    def _compute_dalpha_dq(self):
+        assert self.dielectrics is not None
 
-    def to_hdf5(self, path: pathlib.Path):
-        with h5py.File(path, 'w') as f:
-            f.attrs['spectrum'] = 'IR'
-            f.attrs['version'] = 1
+        # compute:
+        dadqs = numpy.zeros((len(self.nd_modes), 3, 3))
+        for i, mode in enumerate(self.nd_modes):
+            dadqs[i] = numpy.sum(
+                [c * d for (v, c), d in zip(self.nd_stencil, self.dielectrics[i])], axis=0) / self.nd_steps[i]
 
-            g_input = f.create_group('input')
-            g_input.create_dataset('modes', data=self.modes)
-            g_input.create_dataset('frequencies', data=self.frequencies)
-            g_input.create_dataset('dmu_dq', data=self.dmu_dq)
+        print(dadqs)
 
-    def compute_intensities(self) -> NDArray[float]:
-        return (self.dmu_dq ** 2).sum(axis=1)
+        # distribute:
+        mode_to_index = dict((m, i) for i, m in enumerate(self.nd_modes))
+        dalpha_dq = numpy.zeros((len(self), 3, 3))
+        for i, mode in enumerate(self.modes):
+            dalpha_dq[i] = dadqs[mode_to_index[self.nd_mode_equiv[i]]]
+
+        self.dalpha_dq = dalpha_dq

tests/test_vibspec.py

@@ -1,3 +1,5 @@
+import pathlib
+
 import numpy
 import pytest
 
@@ -66,6 +68,7 @@ def test_prepare_raman_SiO2(context_SiO2, tmp_path):
 
     spectrum = phonons.prepare_raman(tmp_path)
 
+    assert spectrum.cell_volume == phonons.phonotopy.unitcell.volume
     assert len(spectrum.modes) == 24  # skip acoustic
     assert numpy.allclose(spectrum.frequencies, phonons.frequencies[3:])
     assert spectrum.irrep_labels == phonons.irrep_labels[3:]
@@ -113,6 +116,7 @@ def test_raman_spectrum_save_SiO2(context_SiO2, tmp_path):
 
     spectrum_read = RamanSpectrum.from_hdf5(tmp_path / 'raman.hdf5')
 
+    assert spectrum.cell_volume == spectrum_read.cell_volume
     assert numpy.allclose(spectrum.modes, spectrum_read.modes)
     assert numpy.allclose(spectrum.frequencies, spectrum_read.frequencies)
     assert numpy.allclose(spectrum.modes, spectrum_read.modes)
@@ -121,3 +125,19 @@ def test_raman_spectrum_save_SiO2(context_SiO2, tmp_path):
     assert numpy.allclose(spectrum.nd_mode_equiv, spectrum_read.nd_mode_equiv)
     assert numpy.allclose(spectrum.nd_modes, spectrum_read.nd_modes)
     assert numpy.allclose(spectrum.nd_steps, spectrum_read.nd_steps)
+
+
+def test_raman_spectrum_extract_dielectrics(context_SiO2):
+    calc_directory = pathlib.Path.cwd() / 'calc_dielectrics'
+    spectrum = RamanSpectrum.from_hdf5(pathlib.Path(calc_directory / 'raman.hdf5'))
+
+    nsteps = len(spectrum.nd_stencil)
+    files = []
+    for mode in spectrum.nd_modes:
+        for i in range(nsteps):
+            files.append(calc_directory / 'unitcell_{:04d}_{:02d}'.format(mode + 1, i + 1) / 'vasprun.xml')
+
+    spectrum.extract_dielectrics(files)
+
+    assert spectrum.dielectrics is not None
+    assert spectrum.dalpha_dq is not None