Merge pull request #1 from pierre-24/add_kp

Allow to compute the spectra and displacements at other points in the Brilouin zone

Author: pierre-24

phonopy_vibspec/phonons_analyzer.py

@@ -1,17 +1,18 @@
 import pathlib
 import numpy
 
+from numpy.typing import NDArray
+
 import phonopy
 from phonopy.interface import calculator as phonopy_calculator
 
-from typing import Optional, List, Tuple
+from typing import Optional, List, Tuple, Union
 
 from phonopy_vibspec import logger
+from phonopy.units import VaspToCm
 from phonopy_vibspec.spectra import RamanSpectrum, InfraredSpectrum
 from phonopy_vibspec.vesta import VestaVector, make_vesta_file
 
-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
@@ -21,35 +22,37 @@
 
 
 class PhononsAnalyzer:
-    """Use Phonopy to extract phonon frequencies and eigenmodes, as well as irreps
+    """Use Phonopy to extract phonon frequencies and eigenmodes, as well as irreps, at a given q-point
+    (default is Gamma).
     """
 
     DC_GEOMETRY_TEMPLATE = 'dielec_mode{:04d}_step{:02d}.vasp'
     VESTA_MODE_TEMPLATE = 'mode{:04d}.vesta'
 
-    def __init__(self, phonon: phonopy.Phonopy):
-        self.phonotopy = phonon
+    def __init__(self, phonon: phonopy.Phonopy, q: Union[NDArray, Tuple[float, float, float]] = (.0, .0, .0)):
+        self.phonopy = phonon
+        self.q = q
         self.structure = phonon.primitive
 
         # get eigenvalues and eigenvectors at gamma point
         # See https://github.com/phonopy/phonopy/issues/308#issuecomment-1769736200
         l_logger.info('Symmetrize force constant')
-        self.phonotopy.symmetrize_force_constants()
-
-        l_logger.info('Run mesh')
-        self.phonotopy.run_mesh([1, 1, 1], with_eigenvectors=True)
+        self.phonopy.symmetrize_force_constants()
 
-        mesh_dict = phonon.get_mesh_dict()
+        l_logger.info('Fetch dynamical matrix at q=({})'.format(', '.join('{:.3f}'.format(x) for x in q)))
+        self.phonopy.dynamical_matrix.run(self.q)
+        dm = self.phonopy.dynamical_matrix.dynamical_matrix
+        eigv, eigf = numpy.linalg.eigh(dm)
 
         self.N = self.structure.get_number_of_atoms()
         l_logger.info('Analyze {} modes (including acoustic)'.format(3 * self.N))
-        self.frequencies = mesh_dict['frequencies'][0] * THZ_TO_INV_CM  # in [cm⁻¹]
+        self.frequencies = numpy.sqrt(numpy.abs(eigv.real)) * numpy.sign(eigv.real) * VaspToCm  # in [cm⁻¹]
 
-        l_logger.info('The 5 first modes are {}'.format(
+        l_logger.info('The 5 first modes are at (in cm⁻¹) {}'.format(
             ', '.join('{:.3f}'.format(x) for x in self.frequencies[:5]))
         )
 
-        self.eigenvectors = mesh_dict['eigenvectors'][0].real.T  # in  [Å sqrt(AMU)]
+        self.eigenvectors = eigf.real.T  # in  [Å sqrt(AMU)]
 
         # compute displacements with Eq. 4 of 10.1039/C7CP01680H
         sqrt_masses = numpy.repeat(numpy.sqrt(self.structure.masses), 3)
@@ -59,7 +62,7 @@ def __init__(self, phonon: phonopy.Phonopy):
         self.irrep_labels = ['A'] * (self.N * 3)
 
         try:
-            self.phonotopy.set_irreps([0, 0, 0])
+            self.phonopy.set_irreps(q)
             self.irreps = phonon.get_irreps()
 
             # TODO: that's internal API, so subject to change!
@@ -74,7 +77,8 @@ def from_phonopy(
         cls,
         phonopy_yaml: str = 'phonopy_disp.yaml',
         force_constants_filename: str = 'force_constants.hdf5',
-        born_filename: Optional[str] = None
+        born_filename: Optional[str] = None,
+        q: Union[NDArray, Tuple[float, float, float]] = (.0, .0, .0)
     ) -> 'PhononsAnalyzer':
         """
         Use the Python interface of Phonopy, see https://phonopy.github.io/phonopy/phonopy-module.html.
@@ -86,7 +90,7 @@ def from_phonopy(
             phonopy_yaml=phonopy_yaml,
             force_constants_filename=force_constants_filename,
             born_filename=born_filename,
-        ))
+        ), q=q)
 
     def infrared_spectrum(self, modes: Optional[List[int]] = None) -> InfraredSpectrum:
         """
@@ -96,11 +100,11 @@ def infrared_spectrum(self, modes: Optional[List[int]] = None) -> InfraredSpectr
 
         l_logger.info('Create IR spectrum object')
 
-        born_tensor = self.phonotopy.nac_params['born']
+        born_tensor = self.phonopy.nac_params['born']
 
         # select modes if any
         if modes is None:
-            modes = list(range(3, 3 * self.N))
+            modes = list(range(3 if numpy.allclose(self.q, [.0, .0, .0]) else 0, 3 * self.N))
         else:
             for mode in modes:
                 if mode < 0 or mode >= 3 * self.N:
@@ -139,7 +143,11 @@ def prepare_raman(
 
         # select modes if any
         if modes is None:
-            modes = list(range(3, 3 * self.N))
+            modes = list(range(3 if numpy.allclose(self.q, [.0, .0, .0]) else 0, 3 * self.N))
+        else:
+            for mode in modes:
+                if mode < 0 or mode >= 3 * self.N:
+                    raise IndexError(mode)
 
         frequencies = [self.frequencies[m] for m in modes]
         irrep_labels = [self.irrep_labels[m] for m in modes]

phonopy_vibspec/scripts/__init__.py

@@ -9,6 +9,17 @@ def list_of_modes(inp: str) -> List[int]:
         raise argparse.ArgumentTypeError('invalid (space-separated) list of integers `{}` for modes'.format(inp))
 
 
+def q_point(inp: str) -> Tuple[float, float, float]:
+    chunks = inp.split()
+    if len(chunks) != 3:
+        raise argparse.ArgumentTypeError('invalid (space-separated) q-point `{}`, must contains 3 elements'.format(inp))
+
+    try:
+        return tuple(float(x) for x in chunks)
+    except ValueError:
+        raise argparse.ArgumentTypeError('invalid (space-separated) list of floats `{}` for q-point'.format(inp))
+
+
 def add_common_args(parser: argparse.ArgumentParser):
     parser.add_argument(
         '-c', '--phonopy', type=str, help='Phonopy YAML containing the cells', default='phonopy_disp.yaml')
@@ -17,6 +28,12 @@ def add_common_args(parser: argparse.ArgumentParser):
 
     parser.add_argument('-m', '--modes', type=list_of_modes, help='List of modes (1-based)', default='')
 
+    parser.add_argument(
+        '-q',
+        type=q_point,
+        help='q-point at which this should be computed (default is gamma)',
+        default='0 0 0')
+
 
 def interval(s_interval: str) -> Tuple[float, float]:
     """get interval

phonopy_vibspec/scripts/prepare_raman.py

@@ -23,6 +23,7 @@ def main():
     phonons = PhononsAnalyzer.from_phonopy(
         phonopy_yaml=args.phonopy,
         force_constants_filename=args.fc,
+        q=args.q
     )
 
     raman_spectrum = phonons.prepare_raman(

phonopy_vibspec/scripts/spectrum_ir.py

@@ -21,7 +21,8 @@ def main():
     phonons = PhononsAnalyzer.from_phonopy(
         phonopy_yaml=args.phonopy,
         force_constants_filename=args.fc,
-        born_filename=args.born
+        born_filename=args.born,
+        q=args.q
     )
 
     ir_spectrum = phonons.infrared_spectrum(modes=args.modes if len(args.modes) > 0 else None)

phonopy_vibspec/scripts/vesta_modes.py

@@ -47,6 +47,7 @@ def main():
     phonons = PhononsAnalyzer.from_phonopy(
         phonopy_yaml=args.phonopy,
         force_constants_filename=args.fc,
+        q=args.q
     )
 
     phonons.make_vesta_for_modes(

tests/test_vibspec.py

@@ -69,7 +69,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 spectrum.cell_volume == phonons.phonopy.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:]