Merge pull request #763 from zerothi/phonon-re-order

changed phonon displacement shape

CHANGELOG.md

@@ -83,6 +83,7 @@ we hit release version 1.0.0.
 - removed `Selector` and `TimeSelector`, they were never used internally
 
 ### Changed
+- changed `Eigenmode.displacement` return shape, please read the documentation
 - bumped minimal Python version to 3.9, #640
 - documentation build system on RTD is updated, #745
 - `gauge` arguments now accept 'cell' and 'orbital' in replacements for 'R' and 'r', respectively

src/sisl/physics/phonon.py

@@ -295,7 +295,7 @@ def PDOS(self, E, distribution="gaussian"):
         """
         return PDOS(E, self.mode, self.hw, distribution)
 
-    def displacement(self):
+    def displacement(self, atol: float = 1e-9):
         r"""Calculate real-space displacements for a given mode (in units of the characteristic length)
 
         The displacements per mode may be written as:
@@ -311,26 +311,26 @@ def displacement(self):
 
         Parameters
         ----------
-        mode : array_like
-            vectors describing the phonon modes, 2nd dimension contains the modes. In case of degenerate
-            modes the vectors *may* be rotated upon return.
-        hw : array_like
-            frequencies of the modes, for any negative frequency the returned displacement will be 0.
-        mass : array_like
-            masses for the atoms (has to have length ``mode.shape[1] // 3``
+        atol :
+            absolute tolerance for whether a phonon is 0 or not.
+            Since the phonon energy is used in the calculation of the displacement vector
+            we have to remove phonon modes with 0 energy.
+            The displacements for phonon modes with an absolute energy below `atol` will
+            be 0.
 
         Returns
         -------
         numpy.ndarray
-            displacements per mode with final dimension ``(3, mode.shape[0])``, displacements are in Ang
+            displacements per mode with final dimension ``(len(self), self.parent.na, 3)``, displacements are in Ang
         """
-        idx = (self.c == 0).nonzero()[0]
+        # get indices for the zero modes
+        idx = (np.fabs(self.c) <= atol).nonzero()[0]
         mode = self.mode
         U = mode.copy()
         U[idx, :] = 0.0
 
         # Now create the remaining displacements
-        idx = delete(_a.arangei(U.shape[0]), idx)
+        idx = delete(_a.arange(U.shape[0]), idx)
 
         # Generate displacement factor
         factor = _displacement_const / fabs(self.c[idx]).reshape(-1, 1) ** 0.5
@@ -339,5 +339,4 @@ def displacement(self):
         U[idx] = (mode[idx, :] * factor).reshape(
             len(idx), -1, 3
         ) / self._geometry().mass.reshape(1, -1, 1) ** 0.5
-        U = np.swapaxes(U, 0, 2)
         return U

src/sisl/physics/tests/test_dynamical_matrix.py

@@ -125,6 +125,12 @@ def test_dos_pdos_velocity(self, setup):
         em = D.eigenmode()
         assert np.allclose(em.DOS(E), em.PDOS(E).sum(0))
 
+    def test_displacement(self, setup):
+        D = setup.D.copy()
+        D.construct(setup.func)
+        em = D.eigenmode()
+        assert em.displacement().shape == (len(em), D.geometry.na, 3)
+
     def test_pickle(self, setup):
         import pickle as p