remove degenerate modes

phonopy_vibspec/phonons_analyzer.py

@@ -106,7 +106,10 @@ def prepare_raman(
         stencil: Optional[list] = None
     ) -> RamanSpectrum:
         """
-        Prepare a Raman spectrum by creating a set of displaced geometries of the unitcell in `directory`.
+        Prepare a Raman spectrum by preparing numerical differentiation, and thus creating a set of displaced
+        geometries of the unitcell in `directory`.
+        Degenerate modes are removed from the numerical differentiation (i.e., only the first mode of
+        each degenerate set will be computed).
         The number of geometries that are generated per mode depends on the stencil.
 
         The `modes` is a 0-based list of mode to include.
@@ -123,13 +126,27 @@ def prepare_raman(
         irrep_labels = [self.irrep_labels[m] for m in modes]
 
         # create geometries
+        dgsets = {}
+        for dgset in self.irreps._degenerate_sets:
+            for i in dgset:
+                dgsets[i] = dgset[0]
+
+        mode_equiv = []
+        mode_calcs = []
         steps = []
+
         base_geometry = self.phonotopy.unitcell
 
         for mode in modes:
             if mode < 0 or mode >= 3 * self.N:
                 raise IndexError(mode)
 
+            mode_equiv.append(dgsets[mode])
+            if dgsets[mode] != mode:
+                continue
+            else:
+                mode_calcs.append(mode)
+
             step = disp
             if ref == 'norm':
                 step = disp / float(numpy.linalg.norm(self.eigendisps[mode]))
@@ -149,11 +166,15 @@ def prepare_raman(
                 )
 
         calculator = RamanSpectrum(
+            # input
             modes=modes,
             frequencies=frequencies,
             irrep_labels=irrep_labels,
-            steps=steps,
-            stencil=stencil
+            # nd:
+            nd_stencil=stencil,
+            nd_mode_equiv=mode_equiv,
+            nd_modes=mode_calcs,
+            nd_steps=steps,
         )
 
         return calculator

phonopy_vibspec/spectra.py

@@ -14,26 +14,31 @@ def __init__(
         frequencies: List[float],
         irrep_labels: Optional[List[str]] = None,
         dalpha_dq: Optional[NDArray[float]] = None,
-        # calculations
-        stencil: Optional[list] = None,
-        steps: List[float] = None,
-        dielectrics: Optional[NDArray] = None,
+        # nd
+        nd_stencil: Optional[list] = None,
+        nd_mode_equiv: Optional[List[int]] = None,
+        nd_modes: Optional[List[int]] = None,
+        nd_steps: List[float] = None,
+        nd_dielectrics: Optional[NDArray] = None,
     ):
         assert irrep_labels is None or len(irrep_labels) == len(modes)
         assert len(frequencies) == len(modes)
         assert dalpha_dq is None or dalpha_dq.shape[0] == len(modes)
 
-        assert steps is None or len(steps) == len(modes)
-        assert dielectrics is None or dielectrics.shape[0] == len(modes)
+        assert nd_mode_equiv is None or len(nd_mode_equiv) == len(modes)
+        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.modes = modes
         self.frequencies = frequencies
         self.irrep_labels = irrep_labels
         self.dalpha_dq = dalpha_dq
 
-        self.steps = steps
-        self.dielectrics = dielectrics
-        self.stencil = stencil
+        self.nd_stencil = nd_stencil
+        self.nd_mode_equiv = nd_mode_equiv
+        self.nd_modes = nd_modes
+        self.nd_steps = nd_steps
+        self.dielectrics = nd_dielectrics
 
     def __len__(self):
         return len(self.modes)
@@ -53,20 +58,25 @@ def from_hdf5(cls, path: pathlib.Path):
             if 'input/dalpha_dq' in f:
                 dalpha_dq = f['input/dalpha_dq'][:]
 
-            stencil = None
-            steps = None
-            dielectrics = None
+            nd_stencil = None
+            nd_steps = None
+            nd_dielectrics = None
 
-            if 'numerical_differentiation' in f:
-                stencil = f['numerical_differentiation/stencil'][:]
-                steps = f['numerical_differentiation/steps'][:]
+            if 'nd' in f:
+                g_nd = f['nd']
 
-                if 'numerical_differentiation/dielectrics' in f:
-                    dielectrics = f['numerical_differentiation/dielectrics'][:]
+                nd_stencil = g_nd['stencil'][:]
+                nd_mode_equiv = g_nd['mode_equiv'][:]
+                nd_modes = g_nd['modes'][:]
+                nd_steps = g_nd['steps'][:]
+
+                if 'dielectrics' in g_nd:
+                    nd_dielectrics = g_nd['dielectrics'][:]
 
             return cls(
                 modes=modes, frequencies=frequencies, irrep_labels=irrep_labels, dalpha_dq=dalpha_dq,
-                stencil=stencil, steps=steps, dielectrics=dielectrics
+                nd_stencil=nd_stencil, nd_mode_equiv=nd_mode_equiv, nd_modes=nd_modes, nd_steps=nd_steps,
+                nd_dielectrics=nd_dielectrics
             )
 
     def to_hdf5(self, path: pathlib.Path):
@@ -84,14 +94,22 @@ def to_hdf5(self, path: pathlib.Path):
             if self.dalpha_dq:
                 g_input.create_dataset('dalpha_dq', data=self.dalpha_dq)
 
-            g_calcs = f.create_group('numerical_differentiation')
-            g_calcs.create_dataset('stencil', data=numpy.array(self.stencil))
+            g_nd = f.create_group('nd')
+
+            if self.nd_stencil:
+                g_nd.create_dataset('stencil', data=numpy.array(self.nd_stencil))
+
+            if self.nd_mode_equiv:
+                g_nd.create_dataset('mode_equiv', data=self.nd_mode_equiv)
+
+            if self.nd_modes:
+                g_nd.create_dataset('modes', data=self.nd_modes)
 
-            if self.steps:
-                g_calcs.create_dataset('steps', data=self.steps)
+            if self.nd_steps:
+                g_nd.create_dataset('steps', data=self.nd_steps)
 
             if self.dielectrics is not None:
-                g_calcs.create_dataset('dielectrics', data=self.dielectrics)
+                g_nd.create_dataset('dielectrics', data=self.dielectrics)
 
 
 class InfraredSpectrum:

tests/test_vibspec.py

@@ -13,11 +13,12 @@ def test_infrared_SiO2(context_SiO2):
     )
 
     # get spectrum
-    infrared_spectrum = phonons.infrared_spectrum()
-    ir_intensities = infrared_spectrum.compute_intensities()
+    spectrum = phonons.infrared_spectrum()
+    ir_intensities = spectrum.compute_intensities()
 
-    assert len(infrared_spectrum.modes) == 24  # skip acoustic
-    assert numpy.allclose(infrared_spectrum.frequencies, phonons.frequencies[3:])
+    assert len(spectrum.modes) == 24  # skip acoustic
+    assert numpy.allclose(spectrum.frequencies, phonons.frequencies[3:])
+    assert spectrum.irrep_labels == phonons.irrep_labels[3:]
 
     # check that degenerate modes share the same intensities
     for dgset in (phonons.irreps._degenerate_sets):
@@ -28,7 +29,7 @@ def test_infrared_SiO2(context_SiO2):
             assert ir_intensities[dgset[0] - 3] == pytest.approx(ir_intensities[dgset[1] - 3], abs=1e-3)
 
     # SiO2 is D3, so A2 and E mode should be active, A1 should not!
-    for i, label in enumerate(infrared_spectrum.irrep_labels):
+    for i, label in enumerate(spectrum.irrep_labels):
         if label in ['A2', 'E']:
             assert ir_intensities[i] != pytest.approx(.0, abs=1e-8)
         else:
@@ -63,14 +64,24 @@ def test_prepare_raman_SiO2(context_SiO2, tmp_path):
         force_constants_filename='force_constants.hdf5'
     )
 
-    spectrum = phonons.prepare_raman(tmp_path, ref='norm')
-    assert len(spectrum.modes) == 24
-    assert spectrum.stencil == TWO_POINTS_STENCIL
+    spectrum = phonons.prepare_raman(tmp_path)
 
-    for i in range(len(spectrum.modes)):
-        assert spectrum.steps[i] == pytest.approx(
-            0.01 / numpy.linalg.norm(phonons.eigendisps[spectrum.modes[i]])
-        )
+    assert len(spectrum.modes) == 24  # skip acoustic
+    assert numpy.allclose(spectrum.frequencies, phonons.frequencies[3:])
+    assert spectrum.irrep_labels == phonons.irrep_labels[3:]
+
+    assert len(spectrum.nd_modes) == 16  # remove degenerates
+    assert spectrum.nd_stencil == TWO_POINTS_STENCIL
+    assert spectrum.nd_steps == [1e-2] * len(spectrum.nd_modes)
+
+    # check if files have been created
+    for i in range(len(spectrum)):
+        mode = spectrum.modes[i]
+        f = tmp_path / 'unitcell_{:04d}_{:02d}.vasp'.format(mode + 1, 1)
+        if spectrum.nd_mode_equiv[i] != mode:
+            assert not f.exists()
+        else:
+            assert f.exists()
 
 
 def test_prepare_raman_select_modes_SiO2(context_SiO2, tmp_path):
@@ -81,10 +92,14 @@ def test_prepare_raman_select_modes_SiO2(context_SiO2, tmp_path):
         force_constants_filename='force_constants.hdf5'
     )
 
-    spectrum = phonons.prepare_raman(tmp_path, modes=requested_modes)
+    spectrum = phonons.prepare_raman(tmp_path, modes=requested_modes, ref='norm')
 
     assert spectrum.modes == requested_modes
-    assert spectrum.steps == [1e-2] * len(requested_modes)
+
+    for i in range(len(spectrum.nd_modes)):
+        assert spectrum.nd_steps[i] == pytest.approx(
+            0.01 / numpy.linalg.norm(phonons.eigendisps[requested_modes[i]])
+        )
 
 
 def test_raman_spectrum_save_SiO2(context_SiO2, tmp_path):
@@ -102,5 +117,7 @@ def test_raman_spectrum_save_SiO2(context_SiO2, tmp_path):
     assert numpy.allclose(spectrum.frequencies, spectrum_read.frequencies)
     assert numpy.allclose(spectrum.modes, spectrum_read.modes)
 
-    assert numpy.allclose(spectrum.stencil, spectrum_read.stencil)
-    assert numpy.allclose(spectrum.steps, spectrum_read.steps)
+    assert numpy.allclose(spectrum.nd_stencil, spectrum_read.nd_stencil)
+    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)