diff --git a/UserGuide/gettingstarted.rst b/UserGuide/gettingstarted.rst
index 5cf1c49b..88beed99 100644
--- a/UserGuide/gettingstarted.rst
+++ b/UserGuide/gettingstarted.rst
@@ -612,7 +612,7 @@ The two most common file formats are the Phonopy and the ASE calculation.
 
 For phonopy, we expect you generate a FORCE_CONSTANTS file and a phonopy.yaml file containing the information about the structure and the atoms in the supercell.
 
-NOTE: Up to version 1.0 the Phonopy importer assumes the FORCE_CONSTANTS and the structure are written in Rydberg atomic units. This will probably change in the near future.
+The importer reads the physical units from the physical_unit block of the phonopy.yaml file and converts them to the internal conventions. Without that block the phonopy defaults (Angstrom and eV/Angstrom^2) are assumed.
 
 .. code::
 
diff --git a/cellconstructor/Phonons.py b/cellconstructor/Phonons.py
index 18ab7a8c..09813e66 100644
--- a/cellconstructor/Phonons.py
+++ b/cellconstructor/Phonons.py
@@ -1367,209 +1367,66 @@ def load_phonopy(self, yaml_filename = "phonopy.yaml", fc_filename = None):
         LOAD FROM PHONOPY FORCE CONSTANTS
         =================================
 
-        This subroutine load the dynamical matrix from the phonopy FORCE_CONSTANT file.
-        It needs two files: the file with the structure information,
-        and the file with the force constant matrix.
+        Load the dynamical matrix from phonopy output files.
+        It needs two files: the phonopy.yaml file with the structure
+        and the FORCE_CONSTANTS file with the real space force
+        constants (both the full and the compact format are supported).
 
-        TODO: Not working!!
+        Units are read from the physical_unit block of the yaml file
+        and converted to the internal conventions: Angstrom for the
+        positions and Ry/bohr^2 for the force constants. If the block
+        is missing, the phonopy defaults (Angstrom and eV/Angstrom^2)
+        are assumed.
 
         Parameters
         ----------
             yaml_filename : string
-                Path to the YAML file, this contains the info of the structure and the supercell.
+                Path to the phonopy.yaml file, contains the structure,
+                the supercell and the physical units.
             fc_filename: string
                 Path to the FORCE_CONSTANTS file. If None, a file called FORCE_CONSTANTS in the same directory
                 as phonopy.yaml will be looked for.
         """
-        raise NotImplementedError("Error, this subroutine contains bugs. Use CC.Methods.phonopy_fc2_to_tensor2 instead")
-        warnings.warn("This subroutine is not tested yet, use it with care.")
-
-        unit_cell = np.zeros((3,3), dtype = np.double)
-        supercell = np.zeros(3, dtype = np.intc)
-        coords = []
-        atoms = []
-        masses = {}
-
-        superstruct = None
-        unit_cell_itau = []
-
-        with open(yaml_filename, "r") as fp:
-
-            read_primitive_cell = False
-            read_coord = False
-            read_lattice = False
-            read_supercell = False
-            read_superstruct = False
-            counter = 0
-            for line in fp.readlines():
-
-                line = line.strip()
-                if not line:
-                    continue
-
-                data = line.replace(",","").split()
-
-                if line == "supercell_matrix:":
-                    read_supercell = True
-                    counter = 0
-                    continue
-
-                if read_supercell and len(data) == 6:
-                    supercell[counter] = int(data[2 + counter])
-                    counter += 1
-
-                    if counter == 3:
-                        counter = 0
-                        read_supercell = False
-
-                if line == "unit_cell:":
-                    read_primitive_cell = True
-                    continue
+        units, supercell_matrix, cells = _read_phonopy_yaml(yaml_filename)
 
-                if line == "lattice:":
-                    read_lattice = True
-                    counter = 0
-                    continue
-
-                if read_lattice and len(data) == 8:
-                    unit_cell[counter, :] = [float(data[x]) for x in range(2, 5)]
-                    counter += 1
-                    if counter == 3:
-                        counter = 0
-                        read_lattice = False
-
-                if line == "points:":
-                    read_coord = True
-                    atoms = []
-                    coords = []
-                    continue
+        for key in ("unit_cell", "supercell"):
+            if key not in cells:
+                raise IOError("Error, no %s block found in %s" % (key, yaml_filename))
 
-                if read_coord:
-                    if "symbol" in line:
-                        atoms.append(data[2])
-                    if "coordinates" in line:
-                        vector = np.array([float(data[x]) for x in range(2, 5)])
-                        coords.append(Methods.cryst_to_cart(unit_cell, vector))
-                    if "mass" in line:
-                        if not atoms[-1] in masses:
-                            masses[atoms[-1]] = float(data[1]) / MASS_RY_TO_UMA
-                    if "reduced_to" in line:
-                        if read_primitive_cell:
-                            unit_cell_itau.append(int(data[1]) - 1)
-
-                if  "supercell" in line:
-                    if read_primitive_cell:
-                        self.structure = Structure.Structure(len(atoms))
-                        self.structure.atoms = atoms
-                        self.structure.coords[:,:] = np.array(coords) * BOHR_TO_ANGSTROM
-                        self.structure.masses = masses
-                        self.structure.has_unit_cell = True
-                        self.structure.unit_cell = unit_cell.copy() * BOHR_TO_ANGSTROM
-                    read_coord = False
-                    read_lattice = False
-                    read_primitive_cell = False
-                    read_superstruct = True
-                    continue
+        supercell_matrix = np.array(supercell_matrix, dtype = np.intc)
+        if supercell_matrix.shape != (3, 3):
+            raise IOError("Error, no supercell_matrix found in %s" % yaml_filename)
+        if np.any(supercell_matrix != np.diag(np.diag(supercell_matrix))):
+            raise NotImplementedError("Error, only diagonal supercell matrices are supported.")
+        supercell = np.diag(supercell_matrix)
 
-        # Now create the superstructure
-        if read_superstruct:
-            superstruct = Structure.Structure(len(atoms))
-            superstruct.atoms = atoms
-            superstruct.coords[:,:] = np.array(coords) * BOHR_TO_ANGSTROM
-            superstruct.masses = masses
-            superstruct.unit_cell = unit_cell.copy() * BOHR_TO_ANGSTROM
-            superstruct.has_unit_cell = True
+        length_factor = _phonopy_length_to_angstrom(units.get("length"))
+        fc_factor = _phonopy_fc_to_ry_bohr2(units.get("force_constants"))
 
-        # Get the Equivalent atoms in the unit cell
+        self.structure = _structure_from_phonopy_cell(cells["unit_cell"], length_factor)
+        superstruct = _structure_from_phonopy_cell(cells["supercell"], length_factor)
         itau = superstruct.get_itau(self.structure) - 1
 
-        # Now load the Force constant matrix
+        # Load the force constant matrix and convert it to Ry/bohr^2
         if fc_filename is None:
             fc_filename = os.path.join(os.path.dirname(yaml_filename), "FORCE_CONSTANTS")
 
-        fc = np.zeros( (superstruct.N_atoms * 3, superstruct.N_atoms * 3), dtype = np.double)
-        FC_TMP = np.zeros((3,3), dtype = np.double)
+        n_rows, fc_blocks = _read_phonopy_force_constants(fc_filename)
+        s2p = [x - 1 for x in cells["supercell"]["reduced_to"]]
+        fc = _get_phonopy_fc_supercell(fc_blocks, n_rows, superstruct, s2p)
+        fc *= fc_factor
 
-        with open(fc_filename, "r") as fp:
-
-            x = 0
-            y = 0
-            counter = 0
-            FC = np.zeros((3,3), dtype = np.double)
-            for i, line in enumerate(fp.readlines()):
-                line = line.strip()
-                data = line.split()
-
-
-                if i == 0:
-                    nat_prim = int(data[0])
-                    nat_tot = int(data[1])
-                    continue
-
-                iteration = (i - 1) // 4
-                counter = (i-1) % 4
-                x = iteration // nat_tot
-                y = iteration % nat_tot
-
-                if counter > 0:
-                    for new_x in np.arange(superstruct.N_atoms)[itau == x]:
-                        fc[3 * new_x + counter -1, 3*y: 3*y + 3] = [float(fx) for fx in data]
-                        fc[3*y: 3*y + 3, 3 * new_x + counter -1] = [float(fx) for fx in data]
-                #     counter += 1
-
-                #     if counter == 3:
-                #         # Save the FC in the correct blocks
-                #         counter = 0
-                #         for ia, ib in blocks:
-                #             fc[3*ia : 3*ia + 3, 3*ib: 3*ib + 3] = FC_TMP
-                #             fc[3*ib : 3*ib + 3, 3*ia: 3*ia + 3] = FC_TMP
-
-
-                # if len(data) == 2:
-                #     #x = int(data[0]) - 1
-                #     #y = int(data[1]) - 1
-                #     #x = itau[x]
-                #     counter = 0
-
-                #     # Get the blocks
-                #     blocks = []
-                #     #print(x, y)
-                #     DR = self.structure.coords[x, :] - superstruct.coords[y,:]
-                #     for ia in range(superstruct.N_atoms):
-                #         if unit_cell_itau[itau[ia]] != x:
-                #             continue
-                #         for ib in range(superstruct.N_atoms):
-                #             if unit_cell_itau[itau[ib]] != unit_cell_itau[itau[y]]:
-                #                 continue
-
-                #             # Check if the two ia and ib are the correct block
-                #             delta_r = superstruct.coords[ia, :] - superstruct.coords[ib, :]
-                #             dist = Methods.get_closest_vector(superstruct.unit_cell, DR - delta_r)
-                #             if np.linalg.norm(dist) < __EPSILON__:
-                #                 blocks.append((ia,ib))
-
-                # elif len(data) == 3:
-                #     FC_TMP[counter, :] = [float(fx) for fx in data]
-                #     counter += 1
-
-                #     if counter == 3:
-                #         # Save the FC in the correct blocks
-                #         counter = 0
-                #         for ia, ib in blocks:
-                #             fc[3*ia : 3*ia + 3, 3*ib: 3*ib + 3] = FC_TMP
-                #             fc[3*ib : 3*ib + 3, 3*ia: 3*ia + 3] = FC_TMP
-
-
-        # Now transform back in real space
+        # Transform to the q space (with gamma as the first point)
         q_tot = symmetries.GetQGrid(self.structure.unit_cell, supercell)
+        gamma_index = np.argmin(np.sum(np.array(q_tot)**2, axis = 1))
+        q_tot[gamma_index] = q_tot[0].copy()
+        q_tot[0] = np.zeros(3, dtype = np.double)
+
         dynq = GetDynQFromFCSupercell(fc, np.array(q_tot), self.structure, superstruct, itau)
-        self.dynmats = [None] * len(q_tot)
+        self.dynmats = [dynq[iq, :, :] for iq in range(len(q_tot))]
         self.q_tot = q_tot
         self.q_stars = [q_tot]
 
-        for iq in range(len(q_tot)):
-            self.dynmats[iq] = dynq[iq, :, :]
-
         self.AdjustQStar()
 
 
@@ -5294,3 +5151,196 @@ def compute_force(indices):
         final_dyn = correct_dyn
 
     return final_dyn
+
+
+def _read_phonopy_yaml(filename):
+    """
+    Parse the subset of a phonopy.yaml file needed by load_phonopy.
+
+    Returns
+    -------
+        units : dict
+            The physical_unit block (e.g. "length", "force_constants").
+        supercell_matrix : list
+            The rows of the supercell matrix.
+        cells : dict
+            For each cell block ("primitive_cell", "unit_cell",
+            "supercell") a dict with "lattice", "symbols", "fracs",
+            "masses" and "reduced_to".
+    """
+    units = {}
+    supercell_matrix = []
+    cells = {}
+    section = None
+    subkey = None
+
+    with open(filename, "r") as fp:
+        for raw_line in fp:
+            line = raw_line.split("#")[0].rstrip()
+            stripped = line.strip()
+            if not stripped:
+                continue
+
+            data = stripped.replace(",", " ").replace("[", " ").replace("]", " ").split()
+            indented = line[0] in " \t"
+
+            if not indented:
+                if stripped.endswith(":"):
+                    section = stripped[:-1]
+                    subkey = None
+                    if section in ("primitive_cell", "unit_cell", "supercell"):
+                        cells[section] = {"lattice": [], "symbols": [], "fracs": [],
+                                          "masses": [], "reduced_to": []}
+                elif stripped.startswith("-") and section == "supercell_matrix":
+                    supercell_matrix.append([int(x) for x in data[1:4]])
+                continue
+
+            if section == "physical_unit":
+                key, _, value = stripped.partition(":")
+                units[key.strip()] = value.strip().strip('"')
+                continue
+
+            if section not in cells:
+                continue
+
+            cell = cells[section]
+            if stripped.endswith(":"):
+                subkey = stripped[:-1] if stripped[:-1] in ("lattice", "points") else None
+            elif subkey == "lattice" and stripped.startswith("-"):
+                cell["lattice"].append([float(x) for x in data[1:4]])
+            elif subkey == "points":
+                if stripped.startswith("- symbol:"):
+                    cell["symbols"].append(data[2])
+                elif stripped.startswith("coordinates:"):
+                    cell["fracs"].append([float(x) for x in data[1:4]])
+                elif stripped.startswith("mass:"):
+                    cell["masses"].append(float(data[1]))
+                elif stripped.startswith("reduced_to:"):
+                    cell["reduced_to"].append(int(data[1]))
+
+    return units, supercell_matrix, cells
+
+
+def _phonopy_length_to_angstrom(unit):
+    """
+    Conversion factor from a phonopy length unit to Angstrom.
+    The unit string comes from the physical_unit block of phonopy.yaml;
+    None falls back to the phonopy default (Angstrom).
+    """
+    factors = {"angstrom": 1.0, "au": BOHR_TO_ANGSTROM, "bohr": BOHR_TO_ANGSTROM}
+    if unit is None:
+        unit = "angstrom"
+    key = unit.strip().lower()
+    if key not in factors:
+        raise ValueError("Error, unknown phonopy length unit '%s'" % unit)
+    return factors[key]
+
+
+def _phonopy_fc_to_ry_bohr2(unit):
+    """
+    Conversion factor from a phonopy force constant unit to Ry/bohr^2.
+    The unit string comes from the physical_unit block of phonopy.yaml;
+    None falls back to the phonopy default (eV/Angstrom^2).
+    """
+    factors = {"ev/angstrom^2": BOHR_TO_ANGSTROM**2 / RY_TO_EV,
+               "ev/angstrom.au": BOHR_TO_ANGSTROM / RY_TO_EV,
+               "ev/au^2": 1.0 / RY_TO_EV,
+               "ry/au^2": 1.0,
+               "mry/au^2": 1e-3,
+               "hartree/au^2": 2.0}
+    if unit is None:
+        unit = "eV/angstrom^2"
+    key = unit.strip().lower()
+    if key not in factors:
+        raise ValueError("Error, unknown phonopy force constant unit '%s'" % unit)
+    return factors[key]
+
+
+def _structure_from_phonopy_cell(cell, length_factor):
+    """
+    Build a Structure (Angstrom) from a cell block parsed by _read_phonopy_yaml.
+    Masses are converted from AMU to Rydberg units.
+    """
+    nat = len(cell["symbols"])
+    if nat == 0 or len(cell["lattice"]) != 3:
+        raise IOError("Error, incomplete cell block in the phonopy yaml file.")
+
+    lattice = np.array(cell["lattice"], dtype = np.double) * length_factor
+    structure = Structure.Structure(nat)
+    structure.atoms = list(cell["symbols"])
+    structure.coords[:,:] = Methods.cryst_to_cart(lattice, np.array(cell["fracs"], dtype = np.double))
+    structure.masses = {atm: mass / MASS_RY_TO_UMA
+                        for atm, mass in zip(cell["symbols"], cell["masses"])}
+    structure.unit_cell = lattice
+    structure.has_unit_cell = True
+    return structure
+
+
+def _read_phonopy_force_constants(filename):
+    """
+    Read a phonopy FORCE_CONSTANTS file (full or compact format).
+
+    Returns
+    -------
+        n_rows : int
+            The number of row atoms (equal to the number of supercell
+            atoms for the full format, to the number of primitive
+            atoms for the compact one).
+        blocks : dict
+            Maps the 0-based (row_atom, column_atom) supercell indices
+            to the 3x3 force constant blocks, in the units of the file.
+    """
+    with open(filename, "r") as fp:
+        lines = [line.strip() for line in fp if line.strip()]
+
+    header = lines[0].split()
+    n_rows = int(header[0])
+    n_cols = int(header[-1])
+
+    blocks = {}
+    index = 1
+    for _ in range(n_rows * n_cols):
+        i, j = [int(x) - 1 for x in lines[index].split()]
+        blocks[(i, j)] = np.array([[float(x) for x in lines[index + k].split()]
+                                   for k in range(1, 4)], dtype = np.double)
+        index += 4
+    return n_rows, blocks
+
+
+def _get_phonopy_fc_supercell(blocks, n_rows, superstruct, s2p):
+    """
+    Expand phonopy force constant blocks to the full 3N x 3N supercell matrix.
+
+    In the compact format the rows only span the primitive atoms; the
+    missing rows are recovered through the lattice translation that
+    connects each supercell atom to its primitive representative
+    (the reduced_to field of phonopy.yaml):
+    Phi(p + T, j) = Phi(p, j - T).
+    """
+    nat_sc = superstruct.N_atoms
+    fc = np.zeros((3 * nat_sc, 3 * nat_sc), dtype = np.double)
+    fracs = Methods.covariant_coordinates(superstruct.unit_cell, superstruct.coords)
+    lookup = {_wrapped_frac_key(f): i for i, f in enumerate(fracs)}
+
+    if n_rows == nat_sc:
+        s2p = list(range(nat_sc))
+    elif len(s2p) != nat_sc:
+        raise IOError("Error, the compact FORCE_CONSTANTS format requires the reduced_to fields in phonopy.yaml.")
+
+    for i in range(nat_sc):
+        rep = s2p[i]
+        translation = fracs[i] - fracs[rep]
+        for j in range(nat_sc):
+            jp = lookup[_wrapped_frac_key(fracs[j] - translation)]
+            fc[3*i : 3*i + 3, 3*j : 3*j + 3] = blocks[(rep, jp)]
+    return fc
+
+
+def _wrapped_frac_key(frac_coords):
+    """
+    Round fractional coordinates wrapped into [0, 1) to a hashable key,
+    so translated atoms can be matched up to numerical noise.
+    """
+    wrapped = np.array(frac_coords, dtype = np.double) % 1.0
+    wrapped[wrapped > 1.0 - 1e-5] = 0.0
+    return tuple(np.round(wrapped, 4))
diff --git a/tests/TestPhonopyEigenvectors/FORCE_CONSTANTS b/tests/TestPhonopyEigenvectors/FORCE_CONSTANTS
new file mode 100644
index 00000000..0f224a5f
--- /dev/null
+++ b/tests/TestPhonopyEigenvectors/FORCE_CONSTANTS
@@ -0,0 +1,433 @@
+   2   54
+1 1
+    -2.734937511847708     0.000000000000000    -0.000000000000002
+    -0.000000000000000    -2.734937511847708    -0.000000000000002
+     0.000000000000000     0.000000000000000    -0.834439795385673
+1 2
+     0.242014973039306     0.000000000000000     0.166264987451150
+    -0.000000000000000     0.107639721043997    -0.000000000000002
+    -0.166261605472175    -0.000000000000000     0.058411070273433
+1 3
+     0.242014973039306     0.000000000000000    -0.166264987451151
+     0.000000000000000     0.107639721043997    -0.000000000000005
+     0.166261605472175    -0.000000000000000     0.058411070273437
+1 4
+     0.107639721043997    -0.000000000000000    -0.000000000000001
+    -0.000000000000000     0.242014973039306     0.166264987451155
+    -0.000000000000000    -0.166261605472175     0.058411070273432
+1 5
+    -0.071891928785864    -0.074443077411060     0.061096418304508
+    -0.074443077411060    -0.071891928785864     0.061096418304509
+    -0.061096565179435    -0.061096565179435     0.050212038239191
+1 6
+    -0.071891928785864     0.074443077411060    -0.061096418304513
+     0.074443077411060    -0.071891928785864     0.061096418304509
+     0.061096565179435    -0.061096565179435     0.050212038239191
+1 7
+     0.107639721043997    -0.000000000000000     0.000000000000001
+     0.000000000000000     0.242014973039306    -0.166264987451154
+    -0.000000000000000     0.166261605472175     0.058411070273440
+1 8
+    -0.071891928785864     0.074443077411060     0.061096418304514
+     0.074443077411060    -0.071891928785864    -0.061096418304505
+    -0.061096565179435     0.061096565179435     0.050212038239199
+1 9
+    -0.071891928785864    -0.074443077411060    -0.061096418304512
+    -0.074443077411060    -0.071891928785864    -0.061096418304510
+     0.061096565179435     0.061096565179435     0.050212038239199
+1 10
+     0.012059674667633    -0.000000000000000     0.000000000000000
+     0.000000000000000     0.012059674667633     0.000000000000003
+     0.000000000000000    -0.000000000000000     0.061747534752208
+1 11
+    -0.007672017403121    -0.000000000000000    -0.012864389634379
+    -0.000000000000000     0.009900233835320    -0.000000000000000
+    -0.005950390449198    -0.000000000000000    -0.009225405143920
+1 12
+    -0.007672017403121    -0.000000000000000     0.012864389634380
+     0.000000000000000     0.009900233835320     0.000000000000002
+     0.005950390449198    -0.000000000000000    -0.009225405143915
+1 13
+     0.009900233835320    -0.000000000000000     0.000000000000004
+     0.000000000000000    -0.007672017403121    -0.012864389634383
+    -0.000000000000000    -0.005950390449198    -0.009225405143905
+1 14
+    -0.004969269367139    -0.004971330346596    -0.007541286791378
+    -0.004971330346596    -0.004969269367139    -0.007541286791368
+    -0.004062823192360    -0.004062823192360    -0.006154986719700
+1 15
+    -0.004969269367139     0.004971330346596     0.007541286791365
+     0.004971330346596    -0.004969269367139    -0.007541286791370
+     0.004062823192360    -0.004062823192360    -0.006154986719704
+1 16
+     0.009900233835320    -0.000000000000000    -0.000000000000005
+    -0.000000000000000    -0.007672017403121     0.012864389634383
+     0.000000000000000     0.005950390449198    -0.009225405143919
+1 17
+    -0.004969269367139     0.004971330346596    -0.007541286791368
+     0.004971330346596    -0.004969269367139     0.007541286791374
+    -0.004062823192360     0.004062823192360    -0.006154986719694
+1 18
+    -0.004969269367139    -0.004971330346596     0.007541286791367
+    -0.004971330346596    -0.004969269367139     0.007541286791372
+     0.004062823192360     0.004062823192360    -0.006154986719711
+1 19
+     0.012062153655729     0.000000000000000     0.000000000000005
+    -0.000000000000000     0.012062153655729     0.000000000000006
+    -0.000000000000000    -0.000000000000000     0.061747480137757
+1 20
+    -0.007673042976415     0.000000000000000     0.005937073638392
+     0.000000000000000     0.009901780491442     0.000000000000007
+     0.012879460301284    -0.000000000000000    -0.009226082711752
+1 21
+    -0.007673042976415     0.000000000000000    -0.005937073638385
+    -0.000000000000000     0.009901780491442     0.000000000000007
+    -0.012879460301284    -0.000000000000000    -0.009226082711748
+1 22
+     0.009901780491442     0.000000000000000    -0.000000000000003
+     0.000000000000000    -0.007673042976415     0.005937073638396
+    -0.000000000000000     0.012879460301284    -0.009226082711749
+1 23
+    -0.004970047246506    -0.004972108514219     0.004063443410872
+    -0.004972108514219    -0.004970047246506     0.004063443410876
+     0.007541865973409     0.007541865973409    -0.006155442198164
+1 24
+    -0.004970047246506     0.004972108514219    -0.004063443410873
+     0.004972108514219    -0.004970047246506     0.004063443410873
+    -0.007541865973409     0.007541865973409    -0.006155442198182
+1 25
+     0.009901780491442     0.000000000000000    -0.000000000000004
+     0.000000000000000    -0.007673042976415    -0.005937073638381
+    -0.000000000000000    -0.012879460301284    -0.009226082711741
+1 26
+    -0.004970047246506     0.004972108514219     0.004063443410870
+     0.004972108514219    -0.004970047246506    -0.004063443410875
+     0.007541865973409    -0.007541865973409    -0.006155442198175
+1 27
+    -0.004970047246506    -0.004972108514219    -0.004063443410870
+    -0.004972108514219    -0.004970047246506    -0.004063443410874
+    -0.007541865973409    -0.007541865973409    -0.006155442198171
+1 28
+     0.570468108210207     0.778597679242130     0.433135383900911
+     0.778597679242130     0.570468108210207     0.433135383900913
+     0.433167720977483     0.433167720977483    -0.081266978300192
+1 29
+    -0.062793024515072     0.000000000000000     0.000000000000006
+    -0.000000000000000     0.000423135325370    -0.000172880525678
+     0.000000000000000    -0.000173064901124     0.000002558674230
+1 30
+     0.570468108210207    -0.778597679242130    -0.433135383900917
+    -0.778597679242130     0.570468108210207     0.433135383900913
+    -0.433167720977483     0.433167720977483    -0.081266978300186
+1 31
+     0.000423135325370     0.000000000000000    -0.000172880525675
+    -0.000000000000000    -0.062793024515072    -0.000000000000002
+    -0.000173064901124     0.000000000000000     0.000002558674211
+1 32
+    -0.000845159494615    -0.000000000000000     0.000000000000001
+     0.000000000000000    -0.000845159494615     0.000000000000002
+    -0.000000000000000     0.000000000000000     0.000000000000010
+1 33
+     0.000423135325370     0.000000000000000     0.000172880525674
+     0.000000000000000    -0.062793024515072    -0.000000000000005
+     0.000173064901124     0.000000000000000     0.000002558674208
+1 34
+     0.570468108210207    -0.778597679242130     0.433135383900917
+    -0.778597679242130     0.570468108210207    -0.433135383900915
+     0.433167720977483    -0.433167720977483    -0.081266978300188
+1 35
+    -0.062793024515072     0.000000000000000    -0.000000000000006
+    -0.000000000000000     0.000423135325370     0.000172880525671
+    -0.000000000000000     0.000173064901124     0.000002558674217
+1 36
+     0.570468108210207     0.778597679242130    -0.433135383900910
+     0.778597679242130     0.570468108210207    -0.433135383900914
+    -0.433167720977483    -0.433167720977483    -0.081266978300182
+1 37
+     0.000010721709158    -0.000010721709037    -0.002240772398868
+    -0.000010721709037     0.000010721709158    -0.002240772398868
+    -0.002240927803182    -0.002240927803182    -0.004415954304947
+1 38
+    -0.000021443418355     0.000000000000000    -0.000000000000004
+     0.000000000000000     0.000000000000150    -0.000010885765098
+    -0.000000000000000    -0.000010955351926    -0.000021450610848
+1 39
+     0.000010721709158     0.000010721709037     0.002240772398864
+     0.000010721709037     0.000010721709158    -0.002240772398868
+     0.002240927803182    -0.002240927803182    -0.004415954304946
+1 40
+     0.000000000000150     0.000000000000000    -0.000010885765094
+    -0.000000000000000    -0.000021443418355    -0.000000000000000
+    -0.000010955351926     0.000000000000000    -0.000021450610842
+1 41
+    -0.000000000000295     0.000000000000000     0.000000000000001
+     0.000000000000000    -0.000000000000295     0.000000000000003
+    -0.000000000000000     0.000000000000000    -0.000000000000319
+1 42
+     0.000000000000150     0.000000000000000     0.000010885765101
+     0.000000000000000    -0.000021443418355    -0.000000000000008
+     0.000010955351926     0.000000000000000    -0.000021450610838
+1 43
+     0.000010721709158     0.000010721709037    -0.002240772398861
+     0.000010721709037     0.000010721709158     0.002240772398865
+    -0.002240927803182     0.002240927803182    -0.004415954304931
+1 44
+    -0.000021443418355    -0.000000000000000    -0.000000000000003
+    -0.000000000000000     0.000000000000150     0.000010885765097
+     0.000000000000000     0.000010955351926    -0.000021450610834
+1 45
+     0.000010721709158    -0.000010721709037     0.002240772398871
+    -0.000010721709037     0.000010721709158     0.002240772398864
+     0.002240927803182     0.002240927803182    -0.004415954304929
+1 46
+     0.045496537428695     0.071036919572828    -0.121130319227708
+     0.071036919572828     0.045496537428695    -0.121130319227714
+    -0.121120558875011    -0.121120558875011     0.185557726191145
+1 47
+    -0.004178137810083    -0.000000000000000     0.000000000000000
+     0.000000000000000     0.000027981355392     0.000032875539145
+     0.000000000000000     0.000032951024819     0.000019101735599
+1 48
+     0.045496537428695    -0.071036919572828     0.121130319227703
+    -0.071036919572828     0.045496537428695    -0.121130319227716
+     0.121120558875011    -0.121120558875011     0.185557726191153
+1 49
+     0.000027981355392    -0.000000000000000     0.000032875539154
+     0.000000000000000    -0.004178137810083    -0.000000000000004
+     0.000032951024819     0.000000000000000     0.000019101735635
+1 50
+    -0.000055962710809    -0.000000000000000    -0.000000000000001
+     0.000000000000000    -0.000055962710809    -0.000000000000010
+    -0.000000000000000     0.000000000000000     0.000000000000342
+1 51
+     0.000027981355392    -0.000000000000000    -0.000032875539156
+    -0.000000000000000    -0.004178137810083     0.000000000000005
+    -0.000032951024819     0.000000000000000     0.000019101735609
+1 52
+     0.045496537428695    -0.071036919572828    -0.121130319227713
+    -0.071036919572828     0.045496537428695     0.121130319227711
+    -0.121120558875011     0.121120558875011     0.185557726191150
+1 53
+    -0.004178137810083     0.000000000000000    -0.000000000000001
+    -0.000000000000000     0.000027981355392    -0.000032875539158
+     0.000000000000000    -0.000032951024819     0.000019101735613
+1 54
+     0.045496537428695     0.071036919572828     0.121130319227712
+     0.071036919572828     0.045496537428695     0.121130319227709
+     0.121120558875011     0.121120558875011     0.185557726191167
+28 1
+     0.570468108210207     0.778597679242130     0.433135383900913
+     0.778597679242130     0.570468108210207     0.433135383900917
+     0.433167720977483     0.433167720977483    -0.081266978300186
+28 2
+     0.570468108210207    -0.778597679242130    -0.433135383900914
+    -0.778597679242130     0.570468108210207     0.433135383900910
+    -0.433167720977483     0.433167720977483    -0.081266978300182
+28 3
+    -0.062793024515072    -0.000000000000000    -0.000000000000005
+    -0.000000000000000     0.000423135325370    -0.000172880525674
+     0.000000000000000    -0.000173064901124     0.000002558674208
+28 4
+     0.570468108210207    -0.778597679242130     0.433135383900913
+    -0.778597679242130     0.570468108210207    -0.433135383900911
+     0.433167720977483    -0.433167720977483    -0.081266978300192
+28 5
+     0.570468108210207     0.778597679242130    -0.433135383900915
+     0.778597679242130     0.570468108210207    -0.433135383900917
+    -0.433167720977483    -0.433167720977483    -0.081266978300188
+28 6
+    -0.062793024515072     0.000000000000000    -0.000000000000002
+    -0.000000000000000     0.000423135325370     0.000172880525675
+     0.000000000000000     0.000173064901124     0.000002558674211
+28 7
+     0.000423135325370     0.000000000000000    -0.000172880525678
+    -0.000000000000000    -0.062793024515072    -0.000000000000006
+    -0.000173064901124    -0.000000000000000     0.000002558674230
+28 8
+     0.000423135325370     0.000000000000000     0.000172880525671
+    -0.000000000000000    -0.062793024515072     0.000000000000006
+     0.000173064901124     0.000000000000000     0.000002558674217
+28 9
+    -0.000845159494615    -0.000000000000000     0.000000000000002
+     0.000000000000000    -0.000845159494615    -0.000000000000001
+     0.000000000000000     0.000000000000000     0.000000000000010
+28 10
+     0.045496537428695     0.071036919572828    -0.121130319227716
+     0.071036919572828     0.045496537428695    -0.121130319227703
+    -0.121120558875011    -0.121120558875011     0.185557726191153
+28 11
+     0.045496537428695    -0.071036919572828     0.121130319227709
+    -0.071036919572828     0.045496537428695    -0.121130319227712
+     0.121120558875011    -0.121120558875011     0.185557726191167
+28 12
+    -0.004178137810083     0.000000000000000     0.000000000000005
+     0.000000000000000     0.000027981355392     0.000032875539156
+     0.000000000000000     0.000032951024819     0.000019101735609
+28 13
+     0.045496537428695    -0.071036919572828    -0.121130319227714
+    -0.071036919572828     0.045496537428695     0.121130319227708
+    -0.121120558875011     0.121120558875011     0.185557726191145
+28 14
+     0.045496537428695     0.071036919572828     0.121130319227711
+     0.071036919572828     0.045496537428695     0.121130319227713
+     0.121120558875011     0.121120558875011     0.185557726191150
+28 15
+    -0.004178137810083    -0.000000000000000    -0.000000000000004
+     0.000000000000000     0.000027981355392    -0.000032875539154
+     0.000000000000000    -0.000032951024819     0.000019101735635
+28 16
+     0.000027981355392    -0.000000000000000     0.000032875539145
+     0.000000000000000    -0.004178137810083    -0.000000000000000
+     0.000032951024819    -0.000000000000000     0.000019101735599
+28 17
+     0.000027981355392     0.000000000000000    -0.000032875539158
+    -0.000000000000000    -0.004178137810083     0.000000000000001
+    -0.000032951024819    -0.000000000000000     0.000019101735613
+28 18
+    -0.000055962710809    -0.000000000000000    -0.000000000000010
+     0.000000000000000    -0.000055962710809     0.000000000000001
+     0.000000000000000     0.000000000000000     0.000000000000342
+28 19
+     0.000010721709158    -0.000010721709037    -0.002240772398868
+    -0.000010721709037     0.000010721709158    -0.002240772398864
+    -0.002240927803182    -0.002240927803182    -0.004415954304946
+28 20
+     0.000010721709158     0.000010721709037     0.002240772398864
+     0.000010721709037     0.000010721709158    -0.002240772398871
+     0.002240927803182    -0.002240927803182    -0.004415954304929
+28 21
+    -0.000021443418355    -0.000000000000000    -0.000000000000008
+    -0.000000000000000     0.000000000000150    -0.000010885765101
+     0.000000000000000    -0.000010955351926    -0.000021450610838
+28 22
+     0.000010721709158     0.000010721709037    -0.002240772398868
+     0.000010721709037     0.000010721709158     0.002240772398868
+    -0.002240927803182     0.002240927803182    -0.004415954304947
+28 23
+     0.000010721709158    -0.000010721709037     0.002240772398865
+    -0.000010721709037     0.000010721709158     0.002240772398861
+     0.002240927803182     0.002240927803182    -0.004415954304931
+28 24
+    -0.000021443418355     0.000000000000000    -0.000000000000000
+    -0.000000000000000     0.000000000000150     0.000010885765094
+     0.000000000000000     0.000010955351926    -0.000021450610842
+28 25
+     0.000000000000150    -0.000000000000000    -0.000010885765098
+    -0.000000000000000    -0.000021443418355     0.000000000000004
+    -0.000010955351926     0.000000000000000    -0.000021450610848
+28 26
+     0.000000000000150     0.000000000000000     0.000010885765097
+     0.000000000000000    -0.000021443418355     0.000000000000003
+     0.000010955351926    -0.000000000000000    -0.000021450610834
+28 27
+    -0.000000000000295    -0.000000000000000     0.000000000000003
+    -0.000000000000000    -0.000000000000295    -0.000000000000001
+     0.000000000000000     0.000000000000000    -0.000000000000319
+28 28
+    -2.734937511847708     0.000000000000000    -0.000000000000002
+    -0.000000000000000    -2.734937511847708     0.000000000000002
+     0.000000000000000    -0.000000000000000    -0.834439795385673
+28 29
+     0.242014973039306    -0.000000000000000    -0.166264987451154
+     0.000000000000000     0.107639721043997    -0.000000000000001
+     0.166261605472175     0.000000000000000     0.058411070273440
+28 30
+     0.242014973039306     0.000000000000000     0.166264987451155
+     0.000000000000000     0.107639721043997     0.000000000000001
+    -0.166261605472175     0.000000000000000     0.058411070273432
+28 31
+     0.107639721043997     0.000000000000000    -0.000000000000002
+    -0.000000000000000     0.242014973039306    -0.166264987451150
+    -0.000000000000000     0.166261605472175     0.058411070273433
+28 32
+    -0.071891928785864    -0.074443077411060    -0.061096418304505
+    -0.074443077411060    -0.071891928785864    -0.061096418304514
+     0.061096565179435     0.061096565179435     0.050212038239199
+28 33
+    -0.071891928785864     0.074443077411060     0.061096418304509
+     0.074443077411060    -0.071891928785864    -0.061096418304508
+    -0.061096565179435     0.061096565179435     0.050212038239191
+28 34
+     0.107639721043997    -0.000000000000000    -0.000000000000005
+    -0.000000000000000     0.242014973039306     0.166264987451151
+    -0.000000000000000    -0.166261605472175     0.058411070273437
+28 35
+    -0.071891928785864     0.074443077411060    -0.061096418304510
+     0.074443077411060    -0.071891928785864     0.061096418304512
+     0.061096565179435    -0.061096565179435     0.050212038239199
+28 36
+    -0.071891928785864    -0.074443077411060     0.061096418304509
+    -0.074443077411060    -0.071891928785864     0.061096418304513
+    -0.061096565179435    -0.061096565179435     0.050212038239191
+28 37
+     0.012062153655729     0.000000000000000     0.000000000000006
+    -0.000000000000000     0.012062153655729    -0.000000000000005
+    -0.000000000000000     0.000000000000000     0.061747480137757
+28 38
+    -0.007673042976415    -0.000000000000000    -0.005937073638381
+    -0.000000000000000     0.009901780491442     0.000000000000004
+    -0.012879460301284     0.000000000000000    -0.009226082711741
+28 39
+    -0.007673042976415    -0.000000000000000     0.005937073638396
+    -0.000000000000000     0.009901780491442     0.000000000000003
+     0.012879460301284     0.000000000000000    -0.009226082711749
+28 40
+     0.009901780491442    -0.000000000000000     0.000000000000007
+    -0.000000000000000    -0.007673042976415    -0.005937073638392
+    -0.000000000000000    -0.012879460301284    -0.009226082711752
+28 41
+    -0.004970047246506    -0.004972108514219    -0.004063443410875
+    -0.004972108514219    -0.004970047246506    -0.004063443410870
+    -0.007541865973409    -0.007541865973409    -0.006155442198175
+28 42
+    -0.004970047246506     0.004972108514219     0.004063443410876
+     0.004972108514219    -0.004970047246506    -0.004063443410872
+     0.007541865973409    -0.007541865973409    -0.006155442198164
+28 43
+     0.009901780491442     0.000000000000000     0.000000000000007
+    -0.000000000000000    -0.007673042976415     0.005937073638385
+    -0.000000000000000     0.012879460301284    -0.009226082711748
+28 44
+    -0.004970047246506     0.004972108514219    -0.004063443410874
+     0.004972108514219    -0.004970047246506     0.004063443410870
+    -0.007541865973409     0.007541865973409    -0.006155442198171
+28 45
+    -0.004970047246506    -0.004972108514219     0.004063443410873
+    -0.004972108514219    -0.004970047246506     0.004063443410873
+     0.007541865973409     0.007541865973409    -0.006155442198182
+28 46
+     0.012059674667633    -0.000000000000000     0.000000000000003
+     0.000000000000000     0.012059674667633    -0.000000000000000
+    -0.000000000000000    -0.000000000000000     0.061747534752208
+28 47
+    -0.007672017403121     0.000000000000000     0.012864389634383
+     0.000000000000000     0.009900233835320     0.000000000000005
+     0.005950390449198    -0.000000000000000    -0.009225405143919
+28 48
+    -0.007672017403121    -0.000000000000000    -0.012864389634383
+     0.000000000000000     0.009900233835320    -0.000000000000004
+    -0.005950390449198     0.000000000000000    -0.009225405143905
+28 49
+     0.009900233835320     0.000000000000000    -0.000000000000000
+     0.000000000000000    -0.007672017403121     0.012864389634379
+    -0.000000000000000     0.005950390449198    -0.009225405143920
+28 50
+    -0.004969269367139    -0.004971330346596     0.007541286791374
+    -0.004971330346596    -0.004969269367139     0.007541286791368
+     0.004062823192360     0.004062823192360    -0.006154986719694
+28 51
+    -0.004969269367139     0.004971330346596    -0.007541286791368
+     0.004971330346596    -0.004969269367139     0.007541286791378
+    -0.004062823192360     0.004062823192360    -0.006154986719700
+28 52
+     0.009900233835320    -0.000000000000000     0.000000000000002
+     0.000000000000000    -0.007672017403121    -0.012864389634380
+    -0.000000000000000    -0.005950390449198    -0.009225405143915
+28 53
+    -0.004969269367139     0.004971330346596     0.007541286791372
+     0.004971330346596    -0.004969269367139    -0.007541286791367
+     0.004062823192360    -0.004062823192360    -0.006154986719711
+28 54
+    -0.004969269367139    -0.004971330346596    -0.007541286791370
+    -0.004971330346596    -0.004969269367139    -0.007541286791365
+    -0.004062823192360    -0.004062823192360    -0.006154986719704
\ No newline at end of file
diff --git a/tests/TestPhonopyEigenvectors/generate_reference.py b/tests/TestPhonopyEigenvectors/generate_reference.py
new file mode 100644
index 00000000..99ed70bc
--- /dev/null
+++ b/tests/TestPhonopyEigenvectors/generate_reference.py
@@ -0,0 +1,114 @@
+"""Regenerate the phonopy fixture and the reference dynamical matrices.
+
+Run once to (re)create phonopy.yaml, FORCE_CONSTANTS, reference_dynmats.npy
+and reference_qtot.npy for test_phonopy_eigenvectors. Needs phonopy and ase;
+the test itself needs neither.
+
+Two inequivalent atoms of the same species (Cu) in a tetragonal cell make a
+wrong atom assignment observable: the second atom sits off the mirror plane
+so the two sites differ. A genuine phonopy FORCE_CONSTANTS/phonopy.yaml pair
+is produced from EMT forces. An exact harmonic force law built from those
+same force constants then drives compute_phonons_finite_displacements, whose
+atom assignment is known correct; with a linear force law the finite
+displacements recover the force constants exactly.
+"""
+import os
+
+import numpy as np
+from ase import Atoms
+from ase.calculators.emt import EMT
+from ase.calculators.calculator import Calculator, all_changes
+
+import cellconstructor as CC
+import cellconstructor.Phonons
+from cellconstructor.Structure import Structure
+
+from phonopy import Phonopy
+from phonopy.structure.atoms import PhonopyAtoms
+from phonopy.file_IO import write_FORCE_CONSTANTS
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+SC = (3, 3, 3)
+EPS = 0.01
+A, C = 3.6, 4.2
+CELL = np.array([[A, 0, 0], [0, A, 0], [0, 0, C]], dtype=float)
+SCALED = np.array([[0.0, 0.0, 0.0], [0.5, 0.5, 0.35]], dtype=float)
+SYMBOLS = ["Cu", "Cu"]
+
+
+def make_phonopy_fixture():
+    unit = PhonopyAtoms(symbols=SYMBOLS, scaled_positions=SCALED, cell=CELL)
+    phonon = Phonopy(unit, supercell_matrix=np.diag(SC))
+    phonon.generate_displacements(distance=EPS)
+    forces = []
+    for sc in phonon.supercells_with_displacements:
+        atoms = Atoms(symbols=sc.symbols, scaled_positions=sc.scaled_positions,
+                      cell=sc.cell, pbc=True)
+        atoms.calc = EMT()
+        f = atoms.get_forces()
+        forces.append(f - f.mean(axis=0))
+    phonon.forces = forces
+    phonon.produce_force_constants()
+    phonon.save(os.path.join(HERE, "phonopy.yaml"),
+                settings={"force_constants": False})
+    # Write the compact format (one block row per primitive atom): it is far
+    # smaller and forces load_phonopy through the reduced_to expansion, the
+    # path where a wrong atom assignment can hide.
+    p2s_map = phonon.primitive.p2s_map
+    write_FORCE_CONSTANTS(phonon.force_constants[p2s_map],
+                          os.path.join(HERE, "FORCE_CONSTANTS"),
+                          p2s_map=p2s_map)
+    return phonon
+
+
+class HarmonicCalc(Calculator):
+    """F = -Phi u, matching each atom to a reference site by minimum image so
+    the result is independent of the supercell atom ordering."""
+    implemented_properties = ["energy", "forces"]
+
+    def __init__(self, ref_positions, ref_cell, fc_full, **kw):
+        super().__init__(**kw)
+        self.R0 = np.asarray(ref_positions, dtype=float)
+        self.H = np.asarray(ref_cell, dtype=float)
+        self.Hinv = np.linalg.inv(self.H)
+        self.PHI = np.asarray(fc_full, dtype=float)
+        self.n = len(self.R0)
+
+    def calculate(self, atoms=None, properties=("energy",),
+                  system_changes=all_changes):
+        super().calculate(atoms, properties, system_changes)
+        idx = np.empty(self.n, dtype=int)
+        U = np.zeros((self.n, 3))
+        for k, p in enumerate(atoms.get_positions()):
+            dmin = ((p[None, :] - self.R0) @ self.Hinv.T)
+            dmin = (dmin - np.round(dmin)) @ self.H
+            r = int(np.argmin(np.linalg.norm(dmin, axis=1)))
+            idx[k], U[r] = r, dmin[r]
+        forces_ref = -(self.PHI @ U.ravel()).reshape(self.n, 3)
+        forces = np.array([forces_ref[r] for r in idx])
+        self.results = {"energy": 0.5 * float(U.ravel() @ self.PHI @ U.ravel()),
+                        "forces": forces}
+
+
+def full_fc(phonon):
+    fc = phonon.force_constants
+    n = len(phonon.supercell)
+    return fc.transpose(0, 2, 1, 3).reshape(3 * n, 3 * n)
+
+
+def main():
+    phonon = make_phonopy_fixture()
+    calc = HarmonicCalc(phonon.supercell.positions, phonon.supercell.cell,
+                        full_fc(phonon))
+    struct = Structure()
+    struct.generate_from_ase_atoms(
+        Atoms(symbols=SYMBOLS, scaled_positions=SCALED, cell=CELL, pbc=True))
+    dyn = CC.Phonons.compute_phonons_finite_displacements(
+        struct, calc, epsilon=EPS, supercell=SC, use_symmetries=False)
+    np.save(os.path.join(HERE, "reference_dynmats.npy"), np.array(dyn.dynmats))
+    np.save(os.path.join(HERE, "reference_qtot.npy"), np.array(dyn.q_tot))
+    print("regenerated fixture and reference for", len(dyn.q_tot), "q points")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/tests/TestPhonopyEigenvectors/phonopy.yaml b/tests/TestPhonopyEigenvectors/phonopy.yaml
new file mode 100644
index 00000000..e9a40a60
--- /dev/null
+++ b/tests/TestPhonopyEigenvectors/phonopy.yaml
@@ -0,0 +1,398 @@
+phonopy:
+  version: "4.2.0"
+  frequency_unit_conversion_factor: 15.633302
+  symmetry_tolerance: 1.00000e-05
+
+physical_unit:
+  atomic_mass: "AMU"
+  length: "angstrom"
+  force: "eV/angstrom"
+  force_constants: "eV/angstrom^2"
+
+space_group:
+  type: "P4/nmm"
+  number: 129
+  Hall_symbol: "P 4ab 2ab -1ab"
+
+primitive_matrix:
+- [  1.000000000000000,  0.000000000000000,  0.000000000000000 ]
+- [  0.000000000000000,  1.000000000000000,  0.000000000000000 ]
+- [  0.000000000000000,  0.000000000000000,  1.000000000000000 ]
+
+supercell_matrix:
+- [   3,   0,   0 ]
+- [   0,   3,   0 ]
+- [   0,   0,   3 ]
+
+primitive_cell:
+  lattice:
+  - [     3.600000000000000,     0.000000000000000,     0.000000000000000 ] # a
+  - [     0.000000000000000,     3.600000000000000,     0.000000000000000 ] # b
+  - [     0.000000000000000,     0.000000000000000,     4.200000000000000 ] # c
+  points:
+  - symbol: Cu # 1
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 63.546000
+  - symbol: Cu # 2
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.350000000000000 ]
+    mass: 63.546000
+  reciprocal_lattice: # without 2pi
+  - [     0.277777777777778,     0.000000000000000,     0.000000000000000 ] # a*
+  - [     0.000000000000000,     0.277777777777778,     0.000000000000000 ] # b*
+  - [     0.000000000000000,     0.000000000000000,     0.238095238095238 ] # c*
+
+unit_cell:
+  lattice:
+  - [     3.600000000000000,     0.000000000000000,     0.000000000000000 ] # a
+  - [     0.000000000000000,     3.600000000000000,     0.000000000000000 ] # b
+  - [     0.000000000000000,     0.000000000000000,     4.200000000000000 ] # c
+  points:
+  - symbol: Cu # 1
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 2
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.350000000000000 ]
+    mass: 63.546000
+    reduced_to: 2
+
+supercell:
+  lattice:
+  - [    10.800000000000001,     0.000000000000000,     0.000000000000000 ] # a
+  - [     0.000000000000000,    10.800000000000001,     0.000000000000000 ] # b
+  - [     0.000000000000000,     0.000000000000000,    12.600000000000001 ] # c
+  points:
+  - symbol: Cu # 1
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.000000000000000 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 2
+    coordinates: [  0.333333333333333,  0.000000000000000,  0.000000000000000 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 3
+    coordinates: [  0.666666666666667,  0.000000000000000,  0.000000000000000 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 4
+    coordinates: [  0.000000000000000,  0.333333333333333,  0.000000000000000 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 5
+    coordinates: [  0.333333333333333,  0.333333333333333,  0.000000000000000 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 6
+    coordinates: [  0.666666666666667,  0.333333333333333,  0.000000000000000 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 7
+    coordinates: [  0.000000000000000,  0.666666666666667,  0.000000000000000 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 8
+    coordinates: [  0.333333333333333,  0.666666666666667,  0.000000000000000 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 9
+    coordinates: [  0.666666666666667,  0.666666666666667,  0.000000000000000 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 10
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.333333333333333 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 11
+    coordinates: [  0.333333333333333,  0.000000000000000,  0.333333333333333 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 12
+    coordinates: [  0.666666666666667,  0.000000000000000,  0.333333333333333 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 13
+    coordinates: [  0.000000000000000,  0.333333333333333,  0.333333333333333 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 14
+    coordinates: [  0.333333333333333,  0.333333333333333,  0.333333333333333 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 15
+    coordinates: [  0.666666666666667,  0.333333333333333,  0.333333333333333 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 16
+    coordinates: [  0.000000000000000,  0.666666666666667,  0.333333333333333 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 17
+    coordinates: [  0.333333333333333,  0.666666666666667,  0.333333333333333 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 18
+    coordinates: [  0.666666666666667,  0.666666666666667,  0.333333333333333 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 19
+    coordinates: [  0.000000000000000,  0.000000000000000,  0.666666666666667 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 20
+    coordinates: [  0.333333333333333,  0.000000000000000,  0.666666666666667 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 21
+    coordinates: [  0.666666666666667,  0.000000000000000,  0.666666666666667 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 22
+    coordinates: [  0.000000000000000,  0.333333333333333,  0.666666666666667 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 23
+    coordinates: [  0.333333333333333,  0.333333333333333,  0.666666666666667 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 24
+    coordinates: [  0.666666666666667,  0.333333333333333,  0.666666666666667 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 25
+    coordinates: [  0.000000000000000,  0.666666666666667,  0.666666666666667 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 26
+    coordinates: [  0.333333333333333,  0.666666666666667,  0.666666666666667 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 27
+    coordinates: [  0.666666666666667,  0.666666666666667,  0.666666666666667 ]
+    mass: 63.546000
+    reduced_to: 1
+  - symbol: Cu # 28
+    coordinates: [  0.166666666666667,  0.166666666666667,  0.116666666666667 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 29
+    coordinates: [  0.500000000000000,  0.166666666666667,  0.116666666666667 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 30
+    coordinates: [  0.833333333333333,  0.166666666666667,  0.116666666666667 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 31
+    coordinates: [  0.166666666666667,  0.500000000000000,  0.116666666666667 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 32
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.116666666666667 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 33
+    coordinates: [  0.833333333333333,  0.500000000000000,  0.116666666666667 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 34
+    coordinates: [  0.166666666666667,  0.833333333333333,  0.116666666666667 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 35
+    coordinates: [  0.500000000000000,  0.833333333333333,  0.116666666666667 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 36
+    coordinates: [  0.833333333333333,  0.833333333333333,  0.116666666666667 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 37
+    coordinates: [  0.166666666666667,  0.166666666666667,  0.450000000000000 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 38
+    coordinates: [  0.500000000000000,  0.166666666666667,  0.450000000000000 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 39
+    coordinates: [  0.833333333333333,  0.166666666666667,  0.450000000000000 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 40
+    coordinates: [  0.166666666666667,  0.500000000000000,  0.450000000000000 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 41
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.450000000000000 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 42
+    coordinates: [  0.833333333333333,  0.500000000000000,  0.450000000000000 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 43
+    coordinates: [  0.166666666666667,  0.833333333333333,  0.450000000000000 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 44
+    coordinates: [  0.500000000000000,  0.833333333333333,  0.450000000000000 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 45
+    coordinates: [  0.833333333333333,  0.833333333333333,  0.450000000000000 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 46
+    coordinates: [  0.166666666666667,  0.166666666666667,  0.783333333333333 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 47
+    coordinates: [  0.500000000000000,  0.166666666666667,  0.783333333333333 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 48
+    coordinates: [  0.833333333333333,  0.166666666666667,  0.783333333333333 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 49
+    coordinates: [  0.166666666666667,  0.500000000000000,  0.783333333333333 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 50
+    coordinates: [  0.500000000000000,  0.500000000000000,  0.783333333333333 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 51
+    coordinates: [  0.833333333333333,  0.500000000000000,  0.783333333333333 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 52
+    coordinates: [  0.166666666666667,  0.833333333333333,  0.783333333333333 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 53
+    coordinates: [  0.500000000000000,  0.833333333333333,  0.783333333333333 ]
+    mass: 63.546000
+    reduced_to: 28
+  - symbol: Cu # 54
+    coordinates: [  0.833333333333333,  0.833333333333333,  0.783333333333333 ]
+    mass: 63.546000
+    reduced_to: 28
+
+displacements:
+- atom:    1
+  displacement:
+    [   0.0065079137345597,  0.0000000000000000,  0.0075925660236530 ]
+  forces:
+  - [   0.0177992802437177,  0.0000000000000012,  1.2726427064583887 ]
+  - [  -0.0003074008741038,  0.0000000000000006,  1.2648896833610854 ]
+  - [  -0.0028114986418550,  0.0000000000000026,  1.2670563217165121 ]
+  - [  -0.0007097949791350,  0.0012751004749335,  1.2659760274916623 ]
+  - [   0.0009382286726692,  0.0009550135100801,  1.2656390215704321 ]
+  - [   0.0000062343569357, -0.0000226630822965,  1.2664413934936138 ]
+  - [  -0.0007097949791325, -0.0012751004749340,  1.2659760274916627 ]
+  - [   0.0009382286726681, -0.0009550135100804,  1.2656390215704321 ]
+  - [   0.0000062343569355,  0.0000226630822951,  1.2664413934936147 ]
+  - [  -0.0000786167295348,  0.0000000000000003,  1.2659503717642606 ]
+  - [   0.0000917299112541,  0.0000000000000005,  1.2665733614397625 ]
+  - [   0.0000051104251328,  0.0000000000000022,  1.2664090877894825 ]
+  - [  -0.0000649727867540,  0.0000441714334616,  1.2664922351536838 ]
+  - [   0.0000636294182696,  0.0000636443447055,  1.2665197599914098 ]
+  - [   0.0000016291074018, -0.0000016430148601,  1.2664207252082187 ]
+  - [  -0.0000649727867522, -0.0000441714334615,  1.2664922351536847 ]
+  - [   0.0000636294182681, -0.0000636443447061,  1.2665197599914115 ]
+  - [   0.0000016291074020,  0.0000016430148581,  1.2664207252082174 ]
+  - [  -0.0000773178616611,  0.0000000000000011,  1.2659589719115896 ]
+  - [  -0.0000479241320289,  0.0000000000000013,  1.2664547423205359 ]
+  - [   0.0001496504767554, -0.0000000000000002,  1.2665350087765599 ]
+  - [  -0.0000632925471690, -0.0000981264117626,  1.2664941152088776 ]
+  - [  -0.0000249554928807, -0.0000249426859263,  1.2664436025158825 ]
+  - [   0.0000884605332043, -0.0000884725430623,  1.2664955244777703 ]
+  - [  -0.0000632925471704,  0.0000981264117616,  1.2664941152088776 ]
+  - [  -0.0000249554928805,  0.0000249426859262,  1.2664436025158821 ]
+  - [   0.0000884605332039,  0.0000884725430622,  1.2664955244777703 ]
+  - [  -0.0069920286026477, -0.0083704649500550, -1.2685868389891468 ]
+  - [   0.0004104812661250,  0.0000030274050901, -1.2664218951941015 ]
+  - [  -0.0004616619314017,  0.0018096966028191, -1.2629754830975843 ]
+  - [  -0.0000018941490614,  0.0000000000000014, -1.2664225298712795 ]
+  - [   0.0000055248645796,  0.0000000000000005, -1.2664232683216274 ]
+  - [  -0.0000036451781353,  0.0000000000000001, -1.2664247740552230 ]
+  - [  -0.0069920286026484,  0.0083704649500534, -1.2685868389891473 ]
+  - [   0.0004104812661235, -0.0000030274050912, -1.2664218951941011 ]
+  - [  -0.0004616619314000, -0.0018096966028197, -1.2629754830975830 ]
+  - [   0.0000171531820055,  0.0000173067682227, -1.2663747989802683 ]
+  - [   0.0000001492271725,  0.0000000942477838, -1.2664231091317946 ]
+  - [  -0.0000173024091781,  0.0000171575410568, -1.2664042906195452 ]
+  - [   0.0000000898887326,  0.0000000000000015, -1.2664230405966956 ]
+  - [   0.0000000000000055,  0.0000000000000016, -1.2664232797799395 ]
+  - [  -0.0000000898887379,  0.0000000000000000, -1.2664231838279305 ]
+  - [   0.0000171531820059, -0.0000173067682257, -1.2663747989802709 ]
+  - [   0.0000001492271705, -0.0000000942477850, -1.2664231091317946 ]
+  - [  -0.0000173024091766, -0.0000171575410577, -1.2664042906195443 ]
+  - [   0.0006226078591432,  0.0004580246309040, -1.2670467014792355 ]
+  - [   0.0000269340202168, -0.0000001253635651, -1.2664235904510159 ]
+  - [  -0.0012029980715063,  0.0013645142440410, -1.2685957692798950 ]
+  - [  -0.0000004514157910, -0.0000000000000008, -1.2664235914684558 ]
+  - [   0.0000003608666794,  0.0000000000000002, -1.2664232811889005 ]
+  - [   0.0000000905491104,  0.0000000000000008, -1.2664231622060578 ]
+  - [   0.0006226078591442, -0.0004580246309027, -1.2670467014792350 ]
+  - [   0.0000269340202187,  0.0000001253635642, -1.2664235904510146 ]
+  - [  -0.0012029980715092, -0.0013645142440421, -1.2685957692798933 ]
+- atom:    1
+  displacement:
+    [  -0.0065079137345597,  0.0000000000000000, -0.0075925660236530 ]
+  forces:
+  - [  -0.0177981945493146,  0.0000000000000012,  1.2599716279799305 ]
+  - [   0.0003178430209524, -0.0000000000000008,  1.2679406809158587 ]
+  - [   0.0028633077311351,  0.0000000000000008,  1.2657791664896463 ]
+  - [   0.0006912250587982, -0.0012495271530835,  1.2668630659579616 ]
+  - [  -0.0009252636438326, -0.0009416841865433,  1.2671967184410908 ]
+  - [  -0.0000017392095608,  0.0000185157247861,  1.2664086494850069 ]
+  - [   0.0006912250587999,  0.0012495271530834,  1.2668630659579629 ]
+  - [  -0.0009252636438330,  0.0009416841865431,  1.2671967184410917 ]
+  - [  -0.0000017392095612, -0.0000185157247862,  1.2664086494850095 ]
+  - [   0.0000783499152729,  0.0000000000000008,  1.2668880162330685 ]
+  - [  -0.0000983853530259,  0.0000000000000008,  1.2662659517701162 ]
+  - [  -0.0000044896203055,  0.0000000000000025,  1.2664365594717895 ]
+  - [   0.0000638869487507, -0.0000462858865257,  1.2663520261567527 ]
+  - [  -0.0000627441875358, -0.0000627561936243,  1.2663281396175501 ]
+  - [  -0.0000013555922535,  0.0000013684029766,  1.2664254170099025 ]
+  - [   0.0000638869487521,  0.0000462858865257,  1.2663520261567545 ]
+  - [  -0.0000627441875372,  0.0000627561936239,  1.2663281396175519 ]
+  - [  -0.0000013555922532, -0.0000013684029786,  1.2664254170099007 ]
+  - [   0.0000796810492279,  0.0000000000000005,  1.2668966155510697 ]
+  - [   0.0000476800609686,  0.0000000000000014,  1.2663920389687640 ]
+  - [  -0.0001456957233304,  0.0000000000000004,  1.2663177535724905 ]
+  - [   0.0000655873193445,  0.0000975510024422,  1.2663538959182827 ]
+  - [   0.0000248794050766,  0.0000248654934424,  1.2664030203916603 ]
+  - [  -0.0000907529197009,  0.0000907678494641,  1.2663491641968148 ]
+  - [   0.0000655873193434, -0.0000975510024424,  1.2663538959182841 ]
+  - [   0.0000248794050767, -0.0000248654934424,  1.2664030203916594 ]
+  - [  -0.0000907529197010, -0.0000907678494639,  1.2663491641968152 ]
+  - [   0.0070107242168528,  0.0083414078477622, -1.2641832733580958 ]
+  - [  -0.0004068219072277,  0.0000003961297511, -1.2664218589169600 ]
+  - [   0.0003858141952959, -0.0017466167211086, -1.2698471483217828 ]
+  - [   0.0000009885559152,  0.0000000000000016, -1.2664247386239031 ]
+  - [  -0.0000054755855862,  0.0000000000000013, -1.2664232683216272 ]
+  - [   0.0000044870296704,  0.0000000000000006, -1.2664224824416568 ]
+  - [   0.0070107242168533, -0.0083414078477645, -1.2641832733580958 ]
+  - [  -0.0004068219072285, -0.0000003961297539, -1.2664218589169609 ]
+  - [   0.0003858141952961,  0.0017466167211067, -1.2698471483217824 ]
+  - [  -0.0000167355604412, -0.0000168620585285, -1.2664710213364452 ]
+  - [  -0.0000001298766612, -0.0000000726009184, -1.2664234353432167 ]
+  - [   0.0000168654371016, -0.0000167321818647, -1.2664421819618394 ]
+  - [  -0.0000000759794946,  0.0000000000000019, -1.2664235075332306 ]
+  - [   0.0000000000000016,  0.0000000000000017, -1.2664232797799446 ]
+  - [   0.0000000759794932,  0.0000000000000005, -1.2664233673899847 ]
+  - [  -0.0000167355604409,  0.0000168620585256, -1.2664710213364474 ]
+  - [  -0.0000001298766631,  0.0000000726009170, -1.2664234353432158 ]
+  - [   0.0000168654371034,  0.0000167321818634, -1.2664421819618386 ]
+  - [  -0.0006244162008103, -0.0004566352987515, -1.2658055942409718 ]
+  - [  -0.0000274479006615,  0.0000003755270957, -1.2664232999602945 ]
+  - [   0.0012283761516761, -0.0013993542638094, -1.2642014393689367 ]
+  - [   0.0000004126253685, -0.0000000000000000, -1.2664228739320760 ]
+  - [  -0.0000003675343092,  0.0000000000000004, -1.2664232811888954 ]
+  - [  -0.0000000450910608,  0.0000000000000014, -1.2664233004743690 ]
+  - [  -0.0006244162008105,  0.0004566352987519, -1.2658055942409723 ]
+  - [  -0.0000274479006595, -0.0000003755270963, -1.2664232999602940 ]
+  - [   0.0012283761516736,  0.0013993542638083, -1.2642014393689363 ]
diff --git a/tests/TestPhonopyEigenvectors/reference_dynmats.npy b/tests/TestPhonopyEigenvectors/reference_dynmats.npy
new file mode 100644
index 00000000..5fbe5bc7
Binary files /dev/null and b/tests/TestPhonopyEigenvectors/reference_dynmats.npy differ
diff --git a/tests/TestPhonopyEigenvectors/reference_qtot.npy b/tests/TestPhonopyEigenvectors/reference_qtot.npy
new file mode 100644
index 00000000..fbb68e76
Binary files /dev/null and b/tests/TestPhonopyEigenvectors/reference_qtot.npy differ
diff --git a/tests/TestPhonopyEigenvectors/test_phonopy_eigenvectors.py b/tests/TestPhonopyEigenvectors/test_phonopy_eigenvectors.py
new file mode 100644
index 00000000..dc3e223b
--- /dev/null
+++ b/tests/TestPhonopyEigenvectors/test_phonopy_eigenvectors.py
@@ -0,0 +1,59 @@
+import os
+
+import numpy as np
+
+import cellconstructor as CC
+import cellconstructor.Phonons
+
+# Two inequivalent atoms of the same species (Cu) in a tetragonal cell,
+# 3x3x3 supercell. The reference dynamical matrices were produced by
+# compute_phonons_finite_displacements driven by an exact harmonic force
+# law built from this fixture's own FORCE_CONSTANTS (see
+# generate_reference.py); that path has the known-correct atom assignment.
+#
+# Comparing the full matrices, not just the frequencies, is what makes this
+# test sensitive to a wrong atom assignment: swapping the two atoms leaves
+# every eigenvalue unchanged but shifts the matrix by ~1e-1, so a
+# frequency-only check cannot see it.
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+
+def _match_q(q_load, q_ref, unit_cell):
+    """Index map so q_load[map[i]] equals q_ref[i] modulo a reciprocal vector."""
+    bg = np.linalg.inv(unit_cell)
+    mapping = []
+    for q in q_ref:
+        errs = [np.linalg.norm((d := (q - q2) @ bg) - np.round(d)) for q2 in q_load]
+        j = int(np.argmin(errs))
+        assert errs[j] < 1e-6, "q point %s not found in loaded grid" % (q,)
+        mapping.append(j)
+    return mapping
+
+
+def test_load_phonopy_eigenvectors():
+    dyn = CC.Phonons.Phonons()
+    dyn.load_phonopy(os.path.join(HERE, "phonopy.yaml"))
+
+    ref = np.load(os.path.join(HERE, "reference_dynmats.npy"))
+    q_ref = np.load(os.path.join(HERE, "reference_qtot.npy"))
+    mapping = _match_q(np.array(dyn.q_tot), q_ref, dyn.structure.unit_cell)
+
+    max_diff = max(np.max(np.abs(ref[i] - dyn.dynmats[mapping[i]]))
+                   for i in range(len(q_ref)))
+    assert max_diff < 1e-3, \
+        "load_phonopy dynamical matrix deviates from the reference by %.2e" % max_diff
+
+    # The comparison must fail if the two atoms are swapped, otherwise it is
+    # blind to the atom-assignment bug it is meant to catch.
+    perm = np.arange(3 * dyn.structure.N_atoms)
+    perm[0:3], perm[3:6] = [3, 4, 5], [0, 1, 2]
+    swap_diff = max(np.max(np.abs(ref[i] - dyn.dynmats[mapping[i]][np.ix_(perm, perm)]))
+                    for i in range(len(q_ref)))
+    assert swap_diff > 1e-2, \
+        "swapping the two atoms is not detected (%.2e); the test lacks teeth" % swap_diff
+
+
+if __name__ == "__main__":
+    test_load_phonopy_eigenvectors()
+    print("ok")
diff --git a/tests/TestPhonopyInput/test_phonopy_input.py b/tests/TestPhonopyInput/test_phonopy_input.py
index f65de0d2..e4751070 100644
--- a/tests/TestPhonopyInput/test_phonopy_input.py
+++ b/tests/TestPhonopyInput/test_phonopy_input.py
@@ -1,22 +1,98 @@
-import sys, os
+import os
+
 import cellconstructor as CC
 import cellconstructor.Phonons
 
-import pytest
 import numpy as np
 
-@pytest.mark.skip(reason="Bug in read phonopy not catched by the test. Deactivated the function")
-def test_phonopy_input():
+# Frequencies (cm^-1) of the bcc Nb fixture on the commensurate q grid,
+# computed independently with phonopy 2.x from the same
+# phonopy.yaml / FORCE_CONSTANTS pair.
+REFERENCE_FREQS_CM = np.sort(np.concatenate([
+    [0.0] * 3,
+    [91.427] * 12,
+    [138.125] * 6,
+    [168.055] * 6,
+    [172.073] * 6,
+    [176.968] * 6,
+    [195.129] * 6,
+    [217.426] * 3,
+]))
+
+# Lattice parameter of the conventional cell, given as 6.19929787 au
+# in the yaml file.
+A_CONV_ANGSTROM = 6.199297870 * CC.Units.BOHR_TO_ANGSTROM
+
+
+def check_dyn(dyn):
+    assert np.allclose(dyn.structure.unit_cell, np.eye(3) * A_CONV_ANGSTROM,
+                       atol=1e-6)
+    assert np.allclose(dyn.structure.coords[1, :],
+                       np.ones(3) * A_CONV_ANGSTROM / 2, atol=1e-6)
 
-    # Go to the current directory
+    w, pols = dyn.DiagonalizeSupercell()
+    freqs_cm = np.sort(w * CC.Units.RY_TO_CM)
+    assert np.allclose(freqs_cm, REFERENCE_FREQS_CM, atol=0.01)
+
+
+def test_phonopy_input():
     total_path = os.path.dirname(os.path.abspath(__file__))
     os.chdir(total_path)
 
     dyn = CC.Phonons.Phonons()
-    dyn.load_phonopy()
+    dyn.load_phonopy("phonopy.yaml")
+
+    check_dyn(dyn)
 
     dyn.save_qe("prova")
 
 
+def test_phonopy_input_ev_angstrom(tmp_path):
+    # Convert the qe-unit fixture (au, Ry/au^2) to the phonopy default
+    # units (Angstrom, eV/Angstrom^2): the loaded dynamical matrix must
+    # be identical.
+    total_path = os.path.dirname(os.path.abspath(__file__))
+    os.chdir(total_path)
+
+    in_lattice = False
+    yaml_lines = []
+    with open("phonopy.yaml", "r") as fp:
+        for line in fp:
+            stripped = line.strip()
+            if stripped.endswith(":"):
+                in_lattice = stripped == "lattice:"
+            elif in_lattice and stripped.startswith("-"):
+                values = [float(x) for x in
+                          stripped.split("[")[1].split("]")[0].split(",")]
+                values = [v * CC.Units.BOHR_TO_ANGSTROM for v in values]
+                line = "  - [ %22.15f, %22.15f, %22.15f ]\n" % tuple(values)
+            elif stripped.startswith("length:"):
+                line = '  length: "angstrom"\n'
+            elif stripped.startswith("force_constants:"):
+                line = '  force_constants: "eV/angstrom^2"\n'
+            yaml_lines.append(line)
+
+    fc_lines = []
+    fc_factor = CC.Units.RY_TO_EV / CC.Units.BOHR_TO_ANGSTROM**2
+    with open("FORCE_CONSTANTS", "r") as fp:
+        for line in fp:
+            data = line.split()
+            if len(data) == 3:
+                values = [float(x) * fc_factor for x in data]
+                line = " %20.12f %20.12f %20.12f\n" % tuple(values)
+            fc_lines.append(line)
+
+    yaml_filename = str(tmp_path / "phonopy.yaml")
+    with open(yaml_filename, "w") as fp:
+        fp.writelines(yaml_lines)
+    with open(str(tmp_path / "FORCE_CONSTANTS"), "w") as fp:
+        fp.writelines(fc_lines)
+
+    dyn = CC.Phonons.Phonons()
+    dyn.load_phonopy(yaml_filename)
+
+    check_dyn(dyn)
+
+
 if __name__ == "__main__":
     test_phonopy_input()