select modes

README.md

@@ -32,18 +32,35 @@ phonopy -d --dim="1 1 1" -c unitcell.vasp
 rm POSCAR-*
 mv SPOSCAR POSCAR
 
-# 3. Run VASP using `IBRION=8` or `IBRION=6` and appropriate `POTIM`
+# 3. Run VASP using `IBRION=8` or `IBRION=6` with appropriate `POTIM`
 
 # 4. Extract force constants
 phonopy --fc --hdf5 vasprun.xml
 ```
 
 For infrared:
 ```bash
-# 5. Run a calculation with `LEPSILON = .TRUE.` on the unit cell
+# 1. Run a calculation with `LEPSILON = .TRUE.` **on the unit cell**
 
-# 6. Extract Born effective charges from calculations
+# 2. Extract Born effective charges from calculations
 phonopy-vasp-born vasprun.xml > BORN
+
+# 3. Get IR spectrum
+...
+```
+
+For Raman:
+```bash
+# 1. Get displaced geometries
+phonopy-vs-create-displaced
+
+# 2. Run a calculation with `LEPSILON = .TRUE.` on each geometry
+
+# 3. Collect dielectric constants
+...
+
+# 4. Get Raman spectrum
+...
 ```
 
 ## Who?

phonopy_vibspec/phonopy_phonons_analyzer.py

@@ -7,19 +7,22 @@
 from numpy.typing import NDArray
 from typing import Optional, List
 
-
 THZ_TO_INV_CM = 33.35641
 
+# [(value, coef), ...]
+# see https://en.wikipedia.org/wiki/Finite_difference_coefficient
+TWO_POINTS_STENCIL = [[-1, -.5], [1, .5]]  # two-points, centered
+
 
 class PhonopyPhononsAnalyzer:
     def __init__(self, phonon: phonopy.Phonopy):
-        self.phonon = phonon
+        self.phonotopy = phonon
         self.supercell = phonon.supercell
 
         # get eigenvalues and eigenvectors at gamma point
         # See https://github.com/phonopy/phonopy/issues/308#issuecomment-1769736200
-        self.phonon.symmetrize_force_constants()
-        self.phonon.run_mesh([1, 1, 1], with_eigenvectors=True)
+        self.phonotopy.symmetrize_force_constants()
+        self.phonotopy.run_mesh([1, 1, 1], with_eigenvectors=True)
 
         mesh_dict = phonon.get_mesh_dict()
 
@@ -32,7 +35,7 @@ def __init__(self, phonon: phonopy.Phonopy):
         self.eigendisps = (self.eigenvectors / sqrt_masses[numpy.newaxis, :]).reshape(-1, self.N, 3)  # in [Å]
 
         # get irreps
-        self.phonon.set_irreps([0, 0, 0])
+        self.phonotopy.set_irreps([0, 0, 0])
         self.irreps = phonon.get_irreps()
         self.irrep_labels = [None] * (self.N * 3)
 
@@ -63,7 +66,7 @@ def infrared_intensities(self, selected_modes: Optional[List[int]] = None) -> ND
         Intensities are given in [e²/AMU].
         """
 
-        born_tensor = self.phonon.nac_params['born']
+        born_tensor = self.phonotopy.nac_params['born']
 
         # select modes if any
         disps = self.eigendisps
@@ -79,12 +82,14 @@ def create_displaced_geometries(
         disp: float = 1e-2,
         modes: Optional[List[int]] = None,
         ref: Optional[str] = None,
-        stencil: List[float] = [-.5, .5]
+        stencil: List[List[float]] = TWO_POINTS_STENCIL
     ):
         """
         Create a set of displaced geometries of the unitcell in `directory`.
-        The goal is to use a two point stencil, so two geometries are generated per mode.
-        Tries to be clever by only considered one mode if they are generated (E, T).
+        The number of geometries that are generated depends on the stencil.
+
+        The `modes` is a 0-based list of mode.
+        If `mode` is `None`,  all mode are selected, except the acoustic ones, and only one version of degenerated ones.
         """
 
         # select modes
@@ -97,27 +102,27 @@ def create_displaced_geometries(
                 modes.append(dgset[0])
 
         # create geometries
-        base_geometry = self.phonon.unitcell
+        base_geometry = self.phonotopy.unitcell
         raman_disps_info = []
 
         for mode in modes:
             if mode < 0 or mode >= 3 * self.N:
                 raise IndexError(mode)
 
-            mode_displacement = disp
+            step = disp
             if ref == 'norm':
-                mode_displacement = disp / float(numpy.linalg.norm(self.eigendisps[mode]))
+                step = disp / float(numpy.linalg.norm(self.eigendisps[mode]))
 
-            raman_disps_info.append({'mode': mode, 'step': mode_displacement})
+            raman_disps_info.append({'mode': mode, 'step': step})
 
-            for i, value in enumerate(stencil):
+            for i, (value, _) in enumerate(stencil):
                 displaced_geometry = base_geometry.copy()
                 displaced_geometry.set_positions(
-                    base_geometry.positions + value * mode_displacement * self.eigendisps[mode]
+                    base_geometry.positions + value * step * self.eigendisps[mode]
                 )
 
                 calculator.write_crystal_structure(
-                    directory / 'unitcell_{:04d}_{:02d}.vasp'.format(mode, i),
+                    directory / 'unitcell_{:04d}_{:02d}.vasp'.format(mode + 1, i + 1),  # 1-based output
                     displaced_geometry,
                     interface_mode='vasp'
                 )

phonopy_vibspec/scripts/displaced_geometries.py

@@ -5,6 +5,8 @@
 import argparse
 import pathlib
 
+from typing import List
+
 from phonopy_vibspec.phonopy_phonons_analyzer import PhonopyPhononsAnalyzer
 from phonopy_vibspec.scripts import add_common_args
 
@@ -17,22 +19,35 @@ def valid_dir(inp: str) -> pathlib.Path:
     return path
 
 
+def list_of_modes(inp: str) -> List[int]:
+    try:
+        return [int(x) - 1 for x in inp.split()]
+    except ValueError:
+        raise argparse.ArgumentTypeError('invalid (space-separated) list of integers `{}` for modes'.format(inp))
+
+
 def main():
     parser = argparse.ArgumentParser(description=__doc__)
     add_common_args(parser)
 
     parser.add_argument('-d', '--displacement', type=float, help='Step size', default=1e-2)
     parser.add_argument('-r', '--ref', choices=['norm', 'none'], help='Reference for the step', default='none')
+    parser.add_argument('-m', '--modes', type=list_of_modes, help='List of modes (1-based)', default='')
     parser.add_argument('-o', '--output', type=valid_dir, help='Output directory', default='.')
 
     args = parser.parse_args()
 
-    results = PhonopyPhononsAnalyzer.from_phonopy(
+    phonons = PhonopyPhononsAnalyzer.from_phonopy(
         phonopy_yaml=args.phonopy,
         force_constants_filename=args.fc,
     )
 
-    results.create_displaced_geometries(pathlib.Path('.'), disp=args.displacement, ref=args.ref)
+    phonons.create_displaced_geometries(
+        pathlib.Path('.'),
+        disp=args.displacement,
+        ref=args.ref,
+        modes=args.modes if len(args.modes) > 0 else None
+    )
 
 
 if __name__ == '__main__':

tests/test_vibspec.py

@@ -2,48 +2,70 @@
 import pytest
 import yaml
 
-from phonopy_vibspec.phonopy_phonons_analyzer import PhonopyPhononsAnalyzer
+from phonopy_vibspec.phonopy_phonons_analyzer import PhonopyPhononsAnalyzer, TWO_POINTS_STENCIL
 
 
 def test_infrared_SiO2(context_SiO2):
-    results = PhonopyPhononsAnalyzer.from_phonopy(
+    phonons = PhonopyPhononsAnalyzer.from_phonopy(
         phonopy_yaml='phonopy_disp.yaml',
         force_constants_filename='force_constants.hdf5',
         born_filename='BORN'
     )
 
     # compute intensities
-    freqs = results.frequencies
-    ir_intensities = results.infrared_intensities()
+    freqs = phonons.frequencies
+    ir_intensities = phonons.infrared_intensities()
 
     # check that degenerate modes share the same frequencies and intensities
-    for dgset in (results.irreps._degenerate_sets):
+    for dgset in (phonons.irreps._degenerate_sets):
         if len(dgset) > 1:
             assert freqs[dgset[0]] == pytest.approx(freqs[dgset[1]], abs=1e-3)
             assert ir_intensities[dgset[0]] == pytest.approx(ir_intensities[dgset[1]], abs=1e-3)
 
     # SiO2 is D3, so A2 and E mode should be active, A1 should not!
-    for i, label in enumerate(results.irrep_labels):
+    for i, label in enumerate(phonons.irrep_labels):
         if label in ['A2', 'E']:
             assert ir_intensities[i] != pytest.approx(.0, abs=1e-8)
         else:
             assert ir_intensities[i] == pytest.approx(.0, abs=1e-8)
 
 
-def test_displaced_geometries_SiO2(context_SiO2, tmp_path):
-    results = PhonopyPhononsAnalyzer.from_phonopy(
+def test_create_displaced_geometries_SiO2(context_SiO2, tmp_path):
+    phonons = PhonopyPhononsAnalyzer.from_phonopy(
         phonopy_yaml='phonopy_disp.yaml',
         force_constants_filename='force_constants.hdf5',
         born_filename='BORN'
     )
 
-    results.create_displaced_geometries(tmp_path, ref='norm')
+    phonons.create_displaced_geometries(tmp_path, ref='norm')
 
     # check info:
     with (tmp_path / 'raman_disps.yml').open() as f:
         raman_disps = yaml.load(f, Loader=yaml.Loader)
 
-    assert raman_disps['stencil'] == [-.5, .5]
+    assert raman_disps['stencil'] == TWO_POINTS_STENCIL
 
     for mode_disp in raman_disps['modes']:
-        assert mode_disp['step'] == pytest.approx(0.01 / numpy.linalg.norm(results.eigendisps[mode_disp['mode']]))
+        assert mode_disp['step'] == pytest.approx(0.01 / numpy.linalg.norm(phonons.eigendisps[mode_disp['mode']]))
+
+
+def test_create_displaced_geometries_select_modes_SiO2(context_SiO2, tmp_path):
+    requested_modes = [3, 5, 6]
+
+    phonons = PhonopyPhononsAnalyzer.from_phonopy(
+        phonopy_yaml='phonopy_disp.yaml',
+        force_constants_filename='force_constants.hdf5',
+        born_filename='BORN'
+    )
+
+    phonons.create_displaced_geometries(tmp_path, modes=requested_modes)
+
+    with (tmp_path / 'raman_disps.yml').open() as f:
+        raman_disps = yaml.load(f, Loader=yaml.Loader)
+
+    # check only selected modes has been created
+    assert len(raman_disps['modes']) == len(requested_modes)
+
+    for mode, mode_disp in zip(requested_modes, raman_disps['modes']):
+        assert mode == mode_disp['mode']
+        assert mode_disp['step'] == pytest.approx(0.01)