Update symmetrize_borns_and_epsilon to accept primitive cell to get the same order of born charges as those atoms in primitive cell

This commit is contained in:
Atsushi Togo 2019-06-17 22:39:40 +09:00
parent 90b5b68a1f
commit db63b4c9f3
2 changed files with 121 additions and 4 deletions

View File

@ -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):

View File

@ -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()