primitive cell should be used for irreps

phonopy_vibspec/phonons_analyzer.py

@@ -29,7 +29,7 @@ class PhononsAnalyzer:
 
     def __init__(self, phonon: phonopy.Phonopy):
         self.phonotopy = phonon
-        self.unitcell = phonon.unitcell
+        self.structure = phonon.primitive
 
         # get eigenvalues and eigenvectors at gamma point
         # See https://github.com/phonopy/phonopy/issues/308#issuecomment-1769736200
@@ -41,7 +41,7 @@ def __init__(self, phonon: phonopy.Phonopy):
 
         mesh_dict = phonon.get_mesh_dict()
 
-        self.N = self.unitcell.get_number_of_atoms()
+        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⁻¹]
 
@@ -52,18 +52,22 @@ def __init__(self, phonon: phonopy.Phonopy):
         self.eigenvectors = mesh_dict['eigenvectors'][0].real.T  # in  [Å sqrt(AMU)]
 
         # compute displacements with Eq. 4 of 10.1039/C7CP01680H
-        sqrt_masses = numpy.repeat(numpy.sqrt(self.unitcell.masses), 3)
+        sqrt_masses = numpy.repeat(numpy.sqrt(self.structure.masses), 3)
         self.eigendisps = (self.eigenvectors / sqrt_masses[numpy.newaxis, :]).reshape(-1, self.N, 3)  # in [Å]
 
         # get irreps
-        self.phonotopy.set_irreps([0, 0, 0])
-        self.irreps = phonon.get_irreps()
-        self.irrep_labels = [''] * (self.N * 3)
+        self.irrep_labels = ['A'] * (self.N * 3)
 
-        # TODO: that's internal API, so subject to change!
-        for label, dgset in zip(self.irreps._get_ir_labels(), self.irreps._degenerate_sets):
-            for j in dgset:
-                self.irrep_labels[j] = label
+        try:
+            self.phonotopy.set_irreps([0, 0, 0])
+            self.irreps = phonon.get_irreps()
+
+            # TODO: that's internal API, so subject to change!
+            for label, dgset in zip(self.irreps._get_ir_labels(), self.irreps._degenerate_sets):
+                for j in dgset:
+                    self.irrep_labels[j] = label
+        except RuntimeError as e:
+            l_logger.warn('Error while computing irreps ({}). Incorrect labels will be assigned.'.format(e))
 
     @classmethod
     def from_phonopy(
@@ -150,7 +154,7 @@ def prepare_raman(
         mode_calcs = []
         steps = []
 
-        base_geometry = self.unitcell
+        base_geometry = self.structure
 
         for mode in modes:
             if mode < 0 or mode >= 3 * self.N:
@@ -211,6 +215,8 @@ def make_vesta_for_modes(
     ):
         """Make a VESTA file for each `modes` (or all except acoustic if `mode` is None) containing a vector for
         each atom, corresponding to the eigenvector.
+
+        Note: this use the primitive cell, which might not be equal to the unit cell, so this might be confusing.
         """
 
         l_logger.info('Make VESTA files for each mode')
@@ -219,7 +225,7 @@ def make_vesta_for_modes(
         if modes is None:
             modes = list(range(3, 3 * self.N))
 
-        cell = self.unitcell.cell
+        cell = self.structure.cell
         norms = numpy.linalg.norm(cell, axis=1)
         cart_to_cell = numpy.linalg.inv(cell)
 
@@ -243,7 +249,7 @@ def make_vesta_for_modes(
             with (directory / self.VESTA_MODE_TEMPLATE.format(mode + 1)).open('w') as f:
                 make_vesta_file(
                     f,
-                    self.unitcell,
+                    self.structure,
                     vectors,
                     title='Mode {} ({:.3f} cm⁻¹, {})'.format(mode + 1, self.frequencies[mode], self.irrep_labels[mode])
                 )

phonopy_vibspec/vesta.py

@@ -256,9 +256,12 @@ def make_vesta_file(f: TextIO, structure: PhonopyAtoms, vectors: Optional[List[V
     # add bond search parameters the best we can (e.g., sum of covalent radii for each pair)
     f.write('SBOND\n')
     symbols_set = list(set(structure.symbols))
-    idx = 0
+    idx = -1
     for i, si in enumerate(symbols_set):
         for sj in symbols_set[i:]:
+            if si == sj and si != 'C':  # avoid bonds between the same element, as it seems to confuse VESTA
+                continue
+
             idx += 1
 
             f.write('{:>4} {:>4s} {:>4} {:10.6f} {:10.6f} {:2} {:2} {:2} {:2} {:2} {:10.6f}  {:10.6f}'

tests/conftest.py

@@ -12,3 +12,9 @@ def context_SiO2(monkeypatch):
 def context_SiO2_supercell(monkeypatch):
     dir_ = pathlib.Path(__file__).parent / 'tests_files/SiO2_supercell'
     monkeypatch.chdir(dir_)
+
+
+@pytest.fixture
+def context_CaO(monkeypatch):
+    dir_ = pathlib.Path(__file__).parent / 'tests_files/CaO'
+    monkeypatch.chdir(dir_)

tests/test_vibspec.py

@@ -179,3 +179,20 @@ def test_SiO2_vesta(context_SiO2, tmp_path):
 
     for mode in modes:
         assert (tmp_path / phonons.VESTA_MODE_TEMPLATE.format(mode + 1)).exists()
+
+
+def test_CaO_vesta(context_CaO, tmp_path):
+    # Irreps analysis requires primitive cell, which in this case is not equal to the unit cell!!
+
+    phonons = PhononsAnalyzer.from_phonopy(
+        phonopy_yaml='phonopy_disp.yaml',
+        force_constants_filename='force_constants.hdf5',
+    )
+
+    # just check that everything goes without error
+    modes = [0, 3]
+    assert phonons.irrep_labels != ['A'] * (3 * phonons.N)
+    phonons.make_vesta_for_modes(tmp_path, modes)
+
+    for mode in modes:
+        assert (tmp_path / phonons.VESTA_MODE_TEMPLATE.format(mode + 1)).exists()

tests/tests_files/CaO/BORN

@@ -0,0 +1,4 @@
+# epsilon and Z* of atoms 1 5
+  21.21789450    0.00000000    0.00000000    0.00000000   21.21789450    0.00000000    0.00000000    0.00000000   21.21789450 
+   2.59124958    0.00000000    0.00000000    0.00000000    2.59124958    0.00000000    0.00000000    0.00000000    2.59124958 
+  -2.59124958    0.00000000    0.00000000    0.00000000   -2.59124958    0.00000000    0.00000000    0.00000000   -2.59124958 

tests/tests_files/CaO/phonopy_disp.yaml

@@ -0,0 +1,132 @@
+phonopy:
+  version: "2.20.0"
+  frequency_unit_conversion_factor: 15.633302
+  symmetry_tolerance: 1.00000e-05
+  configuration:
+    cell_filename: "../unitcell.vasp"
+    create_displacements: ".true."
+    primitive_axes: "auto"
+    dim: "1 1 1"
+
+physical_unit:
+  atomic_mass: "AMU"
+  length: "angstrom"
+  force_constants: "eV/angstrom^2"
+
+space_group:
+  type: "Fm-3m"
+  number: 225
+  Hall_symbol: "-F 4 2 3"
+
+primitive_matrix:
+- [  0.000000000000000,  0.500000000000000,  0.500000000000000 ]
+- [  0.500000000000000,  0.000000000000000,  0.500000000000000 ]
+- [  0.500000000000000,  0.500000000000000,  0.000000000000000 ]
+
+supercell_matrix:
+- [   1,   0,   0 ]
+- [   0,   1,   0 ]
+- [   0,   0,   1 ]
+
+primitive_cell:
+  lattice:
+  - [     0.000000000000000,     2.419633000000000,     2.419633000000000 ] # a
+  - [     2.419633000000000,     0.000000000000000,     2.419633000000000 ] # b
+  - [     2.419633000000000,     2.419633000000000,     0.000000000000000 ] # c
+  points:
+  - symbol: Ca # 1
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 40.078000
+  - symbol: O  # 2
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.500000000000000 ]
+    mass: 15.999400
+  reciprocal_lattice: # without 2pi
+  - [    -0.206642908242696,     0.206642908242696,     0.206642908242696 ] # a*
+  - [     0.206642908242696,    -0.206642908242696,     0.206642908242696 ] # b*
+  - [     0.206642908242696,     0.206642908242696,    -0.206642908242696 ] # c*
+
+unit_cell:
+  lattice:
+  - [     4.839266000000000,     0.000000000000000,     0.000000000000000 ] # a
+  - [     0.000000000000000,     4.839266000000000,     0.000000000000000 ] # b
+  - [     0.000000000000000,     0.000000000000000,     4.839266000000000 ] # c
+  points:
+  - symbol: Ca # 1
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 40.078000
+    reduced_to: 1
+  - symbol: Ca # 2
+    coordinates: [  0.000000000000000,  0.500000000000000,  0.500000000000000 ]
+    mass: 40.078000
+    reduced_to: 1
+  - symbol: Ca # 3
+    coordinates: [  0.500000000000000,  0.000000000000000,  0.500000000000000 ]
+    mass: 40.078000
+    reduced_to: 1
+  - symbol: Ca # 4
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.000000000000000 ]
+    mass: 40.078000
+    reduced_to: 1
+  - symbol: O  # 5
+    coordinates: [  0.500000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 15.999400
+    reduced_to: 5
+  - symbol: O  # 6
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.500000000000000 ]
+    mass: 15.999400
+    reduced_to: 5
+  - symbol: O  # 7
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.500000000000000 ]
+    mass: 15.999400
+    reduced_to: 5
+  - symbol: O  # 8
+    coordinates: [  0.000000000000000,  0.500000000000000,  0.000000000000000 ]
+    mass: 15.999400
+    reduced_to: 5
+
+supercell:
+  lattice:
+  - [     4.839266000000000,     0.000000000000000,     0.000000000000000 ] # a
+  - [     0.000000000000000,     4.839266000000000,     0.000000000000000 ] # b
+  - [     0.000000000000000,     0.000000000000000,     4.839266000000000 ] # c
+  points:
+  - symbol: Ca # 1
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 40.078000
+    reduced_to: 1
+  - symbol: Ca # 2
+    coordinates: [  0.000000000000000,  0.500000000000000,  0.500000000000000 ]
+    mass: 40.078000
+    reduced_to: 1
+  - symbol: Ca # 3
+    coordinates: [  0.500000000000000,  0.000000000000000,  0.500000000000000 ]
+    mass: 40.078000
+    reduced_to: 1
+  - symbol: Ca # 4
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.000000000000000 ]
+    mass: 40.078000
+    reduced_to: 1
+  - symbol: O  # 5
+    coordinates: [  0.500000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 15.999400
+    reduced_to: 5
+  - symbol: O  # 6
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.500000000000000 ]
+    mass: 15.999400
+    reduced_to: 5
+  - symbol: O  # 7
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.500000000000000 ]
+    mass: 15.999400
+    reduced_to: 5
+  - symbol: O  # 8
+    coordinates: [  0.000000000000000,  0.500000000000000,  0.000000000000000 ]
+    mass: 15.999400
+    reduced_to: 5
+
+displacements:
+- atom:    1
+  displacement:
+    [   0.0100000000000000,  0.0000000000000000,  0.0000000000000000 ]
+- atom:    5
+  displacement:
+    [   0.0100000000000000,  0.0000000000000000,  0.0000000000000000 ]