diff --git a/phonopy/structure/symmetry.py b/phonopy/structure/symmetry.py index 93e069d1..53418f77 100644 --- a/phonopy/structure/symmetry.py +++ b/phonopy/structure/symmetry.py @@ -403,6 +403,7 @@ def symmetrize_borns_and_epsilon(borns, epsilon, ucell, primitive_matrix=None, + primitive=None, supercell_matrix=None, symprec=1e-5, is_symmetry=True): @@ -426,6 +427,12 @@ def symmetrize_borns_and_epsilon(borns, are returned. shape=(3, 3) dtype='double' + primitive : PhonopyAtoms + This is an alternative of giving primitive_matrix (Mp). Mp is given as + Mp = (a_u, b_u, c_u)^-1 * (a_p, b_p, c_p). + In addition, the order of atoms is alined to those of atoms in this + primitive cell for Born effective charges. No rigid rotation of + crystal structure is assumed. supercell_matrix: array_like, optional Supercell matrix. This is used to select Born effective charges in **primitive cell**. Supercell matrix is needed because primitive @@ -476,16 +483,45 @@ def symmetrize_borns_and_epsilon(borns, import warnings warnings.warn("\n".join(lines)) - if primitive_matrix is None: + if primitive_matrix is None and primitive is None: return borns_, epsilon_ else: + if primitive is not None: + pmat = np.dot(np.linalg.inv(ucell.cell.T), primitive.cell.T) + else: + pmat = primitive_matrix + scell, pcell = _get_supercell_and_primitive( ucell, - primitive_matrix=primitive_matrix, + primitive_matrix=pmat, supercell_matrix=supercell_matrix, symprec=symprec) + idx = [scell.u2u_map[i] for i in scell.s2u_map[pcell.p2s_map]] - return borns_[idx], epsilon_ + borns_in_prim = borns_[idx].copy() + + if primitive is None: + return borns_in_prim, epsilon_ + else: + idx2 = _get_mapping_between_cells(pcell, primitive) + return borns_in_prim[idx2].copy(), epsilon_ + + +def _get_mapping_between_cells(cell_from, cell_to, symprec=1e-5): + indices = [] + lattice = cell_from.cell + pos_from = cell_from.scaled_positions + for i, p_to in enumerate(cell_to.scaled_positions): + diff = pos_from - p_to + diff -= np.rint(diff) + dist = np.sqrt(np.sum(np.dot(diff, lattice) ** 2, axis=1)) + ids = np.nonzero(dist < symprec)[0] + if len(ids) == 1: + indices.append(ids[0]) + else: + msg = "Index matching didn't go well." + raise RuntimeError(msg) + return indices def _symmetrize_2nd_rank_tensor(tensor, symmetry_operations, lattice): diff --git a/test/structure/test_symmetry.py b/test/structure/test_symmetry.py index b1849e04..31cd9066 100644 --- a/test/structure/test_symmetry.py +++ b/test/structure/test_symmetry.py @@ -3,7 +3,9 @@ import unittest import numpy as np import time -from phonopy.structure.symmetry import Symmetry +from phonopy.structure.symmetry import (Symmetry, symmetrize_borns_and_epsilon, + _get_supercell_and_primitive, + _get_mapping_between_cells) from phonopy.structure.cells import get_supercell from phonopy.interface.phonopy_yaml import read_cell_yaml import os @@ -70,7 +72,86 @@ class TestSymmetry(unittest.TestCase): self.assertFalse(len_sym_nonspin == len_sym_brokenspin) +class TestSymmetrizeBornsAndEpsilon(unittest.TestCase): + def setUp(self): + self.pmat = [[0, 0.5, 0.5], [0.5, 0, 0.5], [0.5, 0.5, 0]] + self.smat = np.eye(3, dtype='int') * 2 + self.cell = read_cell_yaml(os.path.join(data_dir, "..", "NaCl.yaml")) + self.scell, self.pcell = _get_supercell_and_primitive( + self.cell, + primitive_matrix=self.pmat, + supercell_matrix=self.smat) + self.idx = [self.scell.u2u_map[i] + for i in self.scell.s2u_map[self.pcell.p2s_map]] + self.epsilon = np.array([[2.43533967, 0, 0], + [0, 2.43533967, 0], + [0, 0, 2.43533967]]) + self.borns = np.array([[[1.08875538, 0, 0], + [0, 1.08875538, 0], + [0, 0, 1.08875538]], + [[1.08875538, 0, 0], + [0, 1.08875538, 0], + [0, 0, 1.08875538]], + [[1.08875538, 0, 0], + [0, 1.08875538, 0], + [0, 0, 1.08875538]], + [[1.08875538, 0, 0], + [0, 1.08875538, 0], + [0, 0, 1.08875538]], + [[-1.08875538, 0, 0], + [0, -1.08875538, 0], + [0, 0, -1.08875538]], + [[-1.08875538, 0, 0], + [0, -1.08875538, 0], + [0, 0, -1.08875538]], + [[-1.08875538, 0, 0], + [0, -1.08875538, 0], + [0, 0, -1.08875538]], + [[-1.08875538, 0, 0], + [0, -1.08875538, 0], + [0, 0, -1.08875538]]]) + + def tearDown(self): + pass + + def test_only_symmetrize(self): + borns, epsilon = symmetrize_borns_and_epsilon( + self.borns, + self.epsilon, + self.cell) + np.testing.assert_allclose(borns, self.borns) + np.testing.assert_allclose(epsilon, self.epsilon) + # primitive_matrix=[[0, 0.5, 0.5], [0.5, 0, 0.5], [0.5, 0.5, 0]]) + + def test_with_pmat_and_smat(self): + borns, epsilon = symmetrize_borns_and_epsilon( + self.borns, + self.epsilon, + self.cell, + primitive_matrix=self.pmat, + supercell_matrix=self.smat) + np.testing.assert_allclose(borns, self.borns[self.idx]) + np.testing.assert_allclose(epsilon, self.epsilon) + + def test_with_pcell(self): + idx2 = _get_mapping_between_cells(self.pcell, self.pcell) + + np.testing.assert_array_equal( + idx2, np.arange(self.pcell.get_number_of_atoms())) + + borns, epsilon = symmetrize_borns_and_epsilon( + self.borns, + self.epsilon, + self.cell, + primitive=self.pcell) + np.testing.assert_allclose(borns, self.borns[self.idx][idx2]) + np.testing.assert_allclose(epsilon, self.epsilon) + + if __name__ == '__main__': suite = unittest.TestLoader().loadTestsFromTestCase(TestSymmetry) unittest.TextTestRunner(verbosity=2).run(suite) + suite = unittest.TestLoader().loadTestsFromTestCase( + TestSymmetrizeBornsAndEpsilon) + unittest.TextTestRunner(verbosity=2).run(suite) # unittest.main()