From ca8479c981dca2d4befa144b258f1204eeb72355 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 9 Jun 2026 22:12:15 +0200 Subject: [PATCH 1/2] Fix units and force constant expansion in Phonons.load_phonopy load_phonopy assumed Bohr coordinates and Ry/au^2 force constants regardless of the calculator that produced the phonopy files, and expanded the compact FORCE_CONSTANTS format with a wrong atom mapping, so the spectrum was wrong away from Gamma. The function was disabled with a NotImplementedError and its test was skipped. Read the physical_unit block of phonopy.yaml and convert lengths and force constants to the internal Angstrom and Ry/bohr^2 conventions, covering the unit sets phonopy writes for its calculators. Expand compact force constants through the lattice translation that connects each supercell atom to its primitive representative (the reduced_to field). Re-enable the TestPhonopyInput test and check structure and frequencies against values computed independently with phonopy, for both the au/Ry fixture and a converted Angstrom/eV fixture. Addresses SSCHAcode/python-sscha#167 --- UserGuide/gettingstarted.rst | 2 +- cellconstructor/Phonons.py | 408 +++++++++++-------- tests/TestPhonopyInput/test_phonopy_input.py | 88 +++- 3 files changed, 312 insertions(+), 186 deletions(-) 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/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() From 162d429559620a276e3a6a890e43a210ff95034d Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 11 Jun 2026 23:54:45 +0200 Subject: [PATCH 2/2] Add eigenvector-sensitive test for load_phonopy atom assignment The existing TestPhonopyInput fixture has one atom per primitive cell and checks only the supercell frequencies, so it cannot detect a wrong atom assignment: permuting atoms leaves the eigenvalues unchanged and corrupts only the eigenvectors. Add TestPhonopyEigenvectors: two inequivalent same-species (Cu) atoms in a tetragonal cell, 3x3x3 supercell, compact FORCE_CONSTANTS so the reduced_to expansion is exercised. The test compares the full dynamical matrices against compute_phonons_finite_displacements, whose atom assignment is correct by construction. The reference is driven by an exact harmonic force law built from the same force constants and is committed precomputed, so the test needs only numpy; generate_reference.py reproduces it with phonopy/ase. The current loader matches the reference to 1.5e-6 Ry/bohr^2, while swapping the two atoms shifts the matrix by 9e-2 with identical frequencies. --- tests/TestPhonopyEigenvectors/FORCE_CONSTANTS | 433 ++++++++++++++++++ .../generate_reference.py | 114 +++++ tests/TestPhonopyEigenvectors/phonopy.yaml | 398 ++++++++++++++++ .../reference_dynmats.npy | Bin 0 -> 15680 bytes .../reference_qtot.npy | Bin 0 -> 776 bytes .../test_phonopy_eigenvectors.py | 59 +++ 6 files changed, 1004 insertions(+) create mode 100644 tests/TestPhonopyEigenvectors/FORCE_CONSTANTS create mode 100644 tests/TestPhonopyEigenvectors/generate_reference.py create mode 100644 tests/TestPhonopyEigenvectors/phonopy.yaml create mode 100644 tests/TestPhonopyEigenvectors/reference_dynmats.npy create mode 100644 tests/TestPhonopyEigenvectors/reference_qtot.npy create mode 100644 tests/TestPhonopyEigenvectors/test_phonopy_eigenvectors.py 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 0000000000000000000000000000000000000000..5fbe5bc7554a5e8b1897189dbdea4bdcc5fe6b15 GIT binary patch literal 15680 zcmbVT30RG3*glbc%QA*+C9*G(#P{-L8D=cm8EaV@nG7bPLM0PLLzYOx2t$bogYv!l zj=d~N3k@M55p61!?eBc={hsMN|6JF9UBAmU=lVwP{XF;c+|T>m=Nxv}PZ>LYx{+q3 zW_70pPRr-Jb{c5biOjdM>D1k<(@#rXSGdky;;__pfs^hJM$BEj+^O>aFL$21%&GFZ zbD#d*&1@=vcj<1n*6jcDtM1_o_2#dDTk}juc7B}p+ zHRhlyp1*Ccaf|aT<8~9XR=ZPO>CG61GlsC&u$YTgIA{zQJ-=$J8WQqm6Bdm+%kd0ksI^|L*sinrJH#gy+d zjjrq=N%y{XpCU)H{CGcz5$jj&9gJ^ST1;vcz6C!n8u0q`x7p1iEAz9STX^h+k8P$7 z@6sM(cWBvYv3Ycim8zasWA^-H8QLWx8ug8rp4to@@qF8_V*i1Ij5p9MN&B|hyH`DW z$Q{FbtQ*mPsTjD!>&Xqbm9qK05yQ_l-YCC%e)`MQc00wr``$#gns8XWaA~MnLBM9Q zS%--&ngsjH951jW;f__q3SqY@UY<0&h~fJJKi2@y5AiGzuMqee;Gaf*-Pn0~I`DqO z^R?`wDPID2f?tmG66eb;Sbhn3ZUVdpD@rm;#pPnHUpvNEbDJu9BYp=KUwp03PtL{p zPYh4|!+_rt)boI&ReZhAS%L%d+l>5p{r0i?$?Mg8#WZz(;$z^M8%yb%<`oJds(9ss zgnC#4KOW~BGk)SnJm(|dM&MhU@s;NwzisTicoulQ@f@qFXOg7WPl>?!RvfxxcWTFD zp$2OGn7;SZH<;HDd13z;1 zMw?Y4$FrHdU`SN8B57AuzCP9I%`!zl;>rj`pEAdb1zs1#8;*E~fj(IB$3dL0jq}IxJRHw~;A_wLiWK>+KtD=>7lQuRQPnd^2vzm78UbIqi8LyBQ{NKd zZ9cmn)+3+vsO%Pny930hzdJpyS2RZKXOliMb#1uVa8AdPDeabu=O%lnmo!}>w(sj+ zUeS0vpBLG(=lTg<27mfrd^2A`Kgh2y@JxW0ig;Ow*MEwFFZX2lawhT%!1)s3c_81% zs(H~$w(Pvv1^KOF`O3Dy`z9mn-OTdEa{JtG=L#2G#CXI{XYpmu*No*ScLiP#oLA|` zry9ZiEBEl;Tdib#oKP)?)zcpN_5@#Z;H4m5CgSn@`XgVC*ZF+2p+g)~4ep;pp?@)} z=OggHiGK7(egf)e2E2~o+lld&3t0W+2{_*ycu$b8x}L=b`l&;qZ+WN5B5$`H#S+2! ze{enx{0h-eYw&M}es4!Vcs(=NeJeLVarSUaiwXn2 zT{nEO$D%!%Qa<#}=aJ7#Tj1SBytmL#GWy>Y`EtB)%Zb|1Th9#k#|J6TV4l{XAA^t| z_>0|v*B*S^K|k5ZZzbyY4tV~kZ>m8()%z<;f8R3y^G4wQhwsx(*v}o|2j0LB{DJ)) zi~ZUh`52{J>hZ|B-Xl{zrU)^B>`-RQl02;rF?Ca@?y{p6iPA{)gxH z(g1G*;&p=GO98$K{2AY0z9w(yb=#DtzR!J(=|$9|6Zl^>IA4Tk_#c`3pGNR|lhyT; zWwrm2)&3{lK)=y?|NQCRPhsf^S<)x$KlmT9yW0QAs|@@A-`DUzdjI@$M7-~e8ZQj| zbF_3(y>In?>H_mWVn@WYgC97A{rw&GH{Z{g5BWUsJg~nTsrI#22*-I5=ed8Wi|5U1 z|D%r&|08n0!u?qo@J6%q2L8E3KaUM}_5bsb(L?D1^7{q(tyR}g=Kg0G{2tG@2JpFl zhR1h{KJ!Dl!F?Vry;t{3=6peSw@`YH`_c{h z@%x3WK&jsEp|2ZS@$*&Xc>HW_*S4_HD!yTG>>_AU2QYwqSej;d3vzHck{*9YY52)-PzJMg%FeujA6 z&@cEOxhwoxD$Z8`Z#MX5sq~Yi_dnd<-$H%if8;i(XLI-g?nf{$BG0b_^fL_p#}a-* zt)Dk~{S01Ev(Av)IYJ@y&G$8*7w{Fi-$Q-Hd&rObeeQpDEE})2sTnV+-@ivoC736^ zAAQk}!N{*S>IeTLa{ZV?KRL+H9reps`yW~Df3(8jMc-Pj*%>Xz@2-2X|HsR8gJZydU(uR z-ao;aco;MM?!fB}ykm&xgm}fk&jx79=&W4e&#EI-*rhj+=DtYmod(5J~?YhK(c zJA38Z4cl>E%ES4WIDZ1qr||5xc!7WB{qZtAWVJe?>ufm!@ts+Gu_yAgX8FmVaXyEg z7e(;f*6SP3rKRNp=YL9@aOmiPVkzm6U2N0s&SW+6>k7QSzzavbd5Cuu_K& zPBwHK{a)Wc&-GJ#{WeA`wJ(ZnU*OY;xS}2xaQ-yTm*e>ro;|_$D&s49y!m^lcV$0T zp&#Q5baNrq+)fjd%-D&kEq9|6-}4pUz|n^sVT} z!&-@6Q?J>xO#9p`J$(e-ZS{`wjhwT+h$hd70~{_1~W63w)bVK95f8 zc8r`f?46LR)6eC9=Y`|(`HNK7S3Ze+V}Q@=*qM`hABPlvf#TIq{^=h#653wrcci#|%9wD+gTvV8;vv4|uEKYN z7{~M@53@gYcfgQry?)v|bS-}JHdh#W#eT$s>hAO!=EXvVH;NuWynyn?r#>|5#qgJT z6Ay+jPwY7?$i0I;KY5f6FLxEa13funUaaL+VFwa!Ny_&RX5M5q@>8CpeaRi*?QzlF zKaS$wbu`td-A5r5@oO~83URd=On82Ot%!9PEsY?``5-TjC-Z)?dR9o#k31~0(}5-r z|CGKx;b2i}>P~NhZzu4*VYD%4lr)6wMZ9{QuUfry=*jSFuO|t}7xOJI_WtIQoSCG* z&)NQHNS=&O`&_WmiC#j#l>GdCi5C2nXA^JY3A|n4D<1ket?A-w^xCGw{p{Mg(?rzo zqYm$+6TJ*QxuBk*PlcqU7ikjL4~t))@ceqpdnecIQXo0&^fSPT#^L-y@LN?g&%Sm+ zd&2e8a{9FubuH_%{s#GyB;?EWNJi)1CA(hg`SSZ#(a-x;^g8nGfO=jynmBLmN^|lj z;xWAr(Cfz?{3q!2L-hHHihiE1()VYc!F{XLufU6VRjp@)UO!?8>bcOYSHH{~O~@aJ zU-Cdm3_aeSa{YK)o&f+#+WdLssiZ`yrA|L>-RUjxwFFa}R;V5XlOUlNaed4D|p_Gn{m9_r_j>9t5w^fS_lYSAyIZ(ni=`6rAl?Tx+o#=S)~% z?YfiKZS7uaHvAzz*tM*D+qMw#w})=ek2n=bF2lrur%XOdmSMBEd>k8&bI*i6-94@D?4`fG+8e?W>_{1kJ`b`xd*d9kMk$jke|r%lWl;v zu@Wyl!H+(ywwVUs$2pA=$wYn+k)Kk}Wvduo{!Y5@$ncfV?bgfp8(i_5;QEue0(_e@ zzM}H;n%>IKhu9Ay$!flGZ&f|BrMa_){&u6|2O+8`Ze+LFHEHCJteyjrFT>L)cuyPm zp)TM%m+@us99e!5?{<=XcP4&rDV7OWk6ajI;$TcVSFNW{fndHj$~@wWJKg?6s()5t zFnNUho-w|xehXOr{MY;I`QBoDl{$g335s5}EwjlohGJv7(3&&3Ii=ZgukbK9hOu9ruFZv)1c)n8++~JN<8gkg|C0#O|%$%?=Ze1L4I@4Z$&@%y=kk3 zm3Yl38Pu~L?NT+q=8aUn?2v7~&1<1+fq8WH%9qm5yS^3nSozSdr$1E$rc5Gtz^_!* z&pxdFKBh`tHD>6?H=*>CaO6gts5e5>`b{6anf_YpubR&c!A@OI^rt(+bL>mTsrx4; zs{19M-Picu0vjJ{K6XK`lIxQw$D4uvM62~9DD|_O&G1b(5iRi6Ku&#yQ-cULhhG*yNE5%rCyHSk*^xpq?s=ZTRdCsLr!*|eA4=) zKg&bYS$Q3r(V2_}-+It%Eb_xVimP)QG@SZtmK1-%ddTc!)yROV_4Luu!KmlBs`~t% z^&{OZlw1d%EA;dNd~ZTOQd&?K({(=N+VPpDo+tgOpvIFAtLphhQtKyAX!3qW;iZx9 zgu06t*=&52B~3NUZvEQDgHB&Ca?F;`!F0)g?GuePtApb}_8R6+l|IGpq-|$N&7U~E z6XKUPF}Xe{Pv}>buclmzz31(GPR?S|$3iJTzts&eLrJ=u+@) zgnfAh`6Xlj>@J!*WxZpTlrr*yHh5;aq|%S3P`ZtM{yq3BeF$UwFFj!=xvAzW?*-qx z*k1v8R*$zBd6Mkww)2~%_>o({!~PM|RQ1$K>izXmFuAmMj*Ht1Y38Hoc;hu%DIMpB zA|Ds#^D;;R?Ip^xo4E6rVzXaxgiiH_3*KVD9;fuiMMNrKPtKZYcHiY{<_#c*EH2hCtr9V^r zPdWTgb%Xs?CMo_WoB5x4rNQ(G_-0`~l>E8>aq%aALO+R^m*bef;pn%rzplgolm-*E z|H)MQXCEWt_Hf!5zp;5zzr@GglJYgAJMOQs;5%OJf3R=hSF{$y!2gKM{}f2`Lz;Zk zVTuto@IR8`e|~`fImqwx4$#j)dgZ_NhusvW*EyT?`wIR?R{YO*@IU+f=~0z_N(JtJ z)>icc(b%_@{=&fjh-1OmMCJFGe|A1)O$_{x+CNu*f6Dw%i7;s5h_6Pu?UlM2#n(ROu&^z2EkM|EUO|2L4Cn{zuFG0Q&i{YCkoMXHr&I#@7!1=ScuD z*k8hH#s9qH{)ar%`Jat?|ARScq2QP==|k*V0M{8D(O!%L)ourN8e{$h}-Wd4jEa52bOZEF%_$a#`#8vHoL<9dL z7wi1bWafW5R0Pvwz_Y>psQphS^FKMvuV@ne=^cZ8%l%I>_7{IYJ-_Pb2O9Q%DP#-$ zb2`JT^gl%HpQ8o)rd#XI`YQwerw+Bz`JWH)KmYL0D>^Ini?irOxId_meR&P}rC|Sr zDE=o;dSfv=^JY{HwSU%>3-Q?JKY+i|A0POicRMNfKWV`G1N#enMV0?~M2R|;W z{ZBPg;Q02)xXqFB(_MQ{CVOQ{U8mjZ``#@}igHU$DX#b^cnU7{x0cKyT{LIgU#s_< zc&*oo=C$tyijQuMbvRudg8N59jg&;~Na*J^W}4s~;`yliQS0|^5XYNeHwrY%l{^M5 z9C*j6I=yve&rCb#aQXLJg@;B}#7YO}cCt8KeFM!Xm@P~*3t~ECqOAOHclywny?^g@ zyq5i#VZn$OLcV;h{I%7U^xR+H+JBt`F=P43hjnE@JO)m~J7{#@+bXx)jfzsT|6o5%QKURuz@Q%|il zs`UBpz|%53xd`<%TQq-sk9GNincKL9HwM2DRyGJd-e9sZJ&~Sxx%O}y+6nn`s%OUS z7+AOR^R5i=?ZNo6c>X?AM>inVoAmbEIIjFqhEQT^J-t;#sc;+hbRFO|`3EsqQog=u z6hyW}?b*%60dd$(e>+`yDJR>-x9`Gq~A65ygslWwM4&_zdzqNh^&~n zXj$LNdU~OrQfR|+^SfDsp?-EAd%354;0I|&ZbbO&lIN15pTiEMGx$oV=R;LJ&CvIO zbPw|V8~Jkm6nZAE2)$J#aQ$$;ihe%&Df)TPHHh6`ct4`i^dae}XYZ=@EZ$A3=?W3% zN)`S5>OjoFSJ6*B=%-5%;re-VJ*xqo`%v1f)(_L8JL}Kaliq~uCm;2^0={ol^U{Jw zDf;O{J%D!=dNkC}whQ5B=C>;l-YWer6PSLi^!-u#sp&wNeg={!YW=YMD)po5_XtxT zy65th@J|!6xPBzApOMeCQ)W7T5fuGY_FJc)+pUA>R^-R$S&3J+o9Og-Cd8ElKtBbT z7p9-yv@80p(vN8n@rHh475!uj*P&nbebjRy)!4t-W#fIpazy*T?E)W46AoLSPcjN3 z?m@@)`*#Q;S>V^bYCkJKzpywKVt;9PKi2;wA1a*hUwf?JO`dKljeA@0LYSD`t1Q_w zL3oS#d=i#c&%Egi$=dqug5fLsQuLqItMc=JoFfA_m;6(M@czAkUOm52shSQeWUjL+_hV&qC-Y4ScOs z_n!%I>0!Dqxk*RHcSQ&ZbG#qvHp!JFpq~8unK^$?>)6_>L~7gP>_gk?MzsEelZ{tb z&dZsUx7*iu{rA4zTD{@D;FVV7GW41Sy>38$r=ZURFB?sPbC4P6#m1TDfq6Ydw zzs<;!$TL6gpFND-L-9WJnyx>^UUb=as~uLl7YNZkl4|wvyDq!|U;h2OqMuS5+8*}@ zn@4NX=k%u6&Bhc<2Kte?e)#@zfqtGoc3HByV-QjFb54K%%3MFU?ff0)*j#1$SrbAA z==8nWm8>)953lD(%wsb2(?>NggjltzlQ5|A`zj^hO8q!i_NT#QR~rAH`gt$4QS@U> z6#abgqjNgvO%3Vt?|pl5(ebFEGTtC!MKf~utF>V{?!2Jg6 zPo@w3t5V;$^1bwYp&xlG=A{htGEb+UcYgFjsMmrOmO-Sh4o~Yqux~}KpWfBYy*&O2 z6CB1*TGiW?|NFrh>CD)=-dgFMaI)o@wl2YEr9-RAj0C^!)NE63PIB0G8r*t~`OU<` z;(%`(e3>+35%e>Z&7n^dR!=iG_FA!pSe3*!-H`iAa5A28{ri~Hf{QAiR(SPjX@C2h zzO;Gu*qK}JE@$U7(4nRt4Mn^ms(6E0Jnb_2z-ah~1fT8H-*4HhVf(j}Y#m;rFKt+f zH!fv4?RIgOMaL)n-w$e;el)FF9OGXZj@>HuoV|re(?0jeEDMyIW*pygBW|<2AMuB> z_-vjvLzSP`uBgY(X__gob|03Jc`AMy&3j|I<(A_!Y03$wd&7-WrE}|I=3I^RrC*PK z>|cCmIc4+0{@%ixkxdNwvUtcZI(ZAVjUm#yF)IdL?fqrCmUaueItJe?xKlCFQ=!eDE==EdBm&N;6Kg2*k gl7W6C1O2f4{?!lS>GeYm^dlSShZ^XI80bg%A4uPAZU6uP literal 0 HcmV?d00001 diff --git a/tests/TestPhonopyEigenvectors/reference_qtot.npy b/tests/TestPhonopyEigenvectors/reference_qtot.npy new file mode 100644 index 0000000000000000000000000000000000000000..fbb68e76e6f9b1c618864618e4dc446acb49a607 GIT binary patch literal 776 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I$7<~j<-nmP)#3giMV1{}b-n)9vdmhJX1TEsOFh_>v9^5Ha^I86P1G;x@E zdmQ30^)UT#8mb<~M>hv14x@3J1Mz0V9(4DkiNno76NlT2CT 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")