# Coverage for hiphive/core/eigentensors.py: 100%

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

## 67 statements

1import numpy as np

2import itertools

3import sympy

4from .tensors import rotate_tensor

5from ..input_output.logging_tools import logger

6from .utilities import SparseMatrix

9logger = logger.getChild('eigensymmetries')

12def create_eigentensors(cs):

14 for orbit in cs.orbits:

15 prototype_cluster = cs._cluster_list[orbit.prototype_index]

17 # The first round of eigentensors are created from the atom label

18 # permutation condition

19 ets = init_ets_from_label_symmetry(prototype_cluster)

21 # Loop over the symmetries and apply them iteratively to reduce the

22 # number of eigentensors

23 for rotation_index, permutation_index in orbit.eigensymmetries:

25 # Rotation amtrix in scaled coordinates

26 R = cs.rotation_matrices[rotation_index]

27 p = cs._permutations[permutation_index] # permutation vector

29 # Skip if the symmetry is the trivial identity operation

30 if np.allclose(R, np.eye(3)) and p == tuple(range(orbit.order)):

31 continue

33 p_inv = np.argsort(p)

35 # Create the intermediate constraint matrix and populate it

36 M = SparseMatrix.zeros(3**orbit.order, len(ets))

37 populate_constraint_matrix(M, ets, R, p_inv)

39 nullspace = M.nullspace()

41 new_ets = []

42 for solution in nullspace:

43 # Normalize the solution vector so it is inly integers

44 # Done by multiplying with the greatest common multiple of the

45 # denominators present

46 solution = renormalize_to_integer(solution)

47 # From the solution vector and the old eigentensors a new et is

48 # created

49 new_ets.append(assemble_new_eigentensor(ets, solution))

50 ets = new_ets

52 # If the force constant is forbidden by symmetry, break early

53 if len(ets) == 0:

54 break

56 orbit.eigentensors = ets

59def populate_constraint_matrix(M, ets, R, p_inv):

60 for et_index, et in enumerate(ets):

61 tmp = rotate_tensor(et.transpose(p_inv), R)

62 tmp = tmp - et

63 for element_index, v in enumerate(tmp.flat):

64 if v != 0:

65 M[element_index, et_index] = v

68def init_ets_from_label_symmetry(cluster):

69 """ This function creates eigentensor which fulfill the atom label

70 symmetry. The only information needed is the corresponding cluster

71 """

72 order = len(cluster)

73 assert sorted(cluster) == list(cluster)

74 multiplicities = np.unique(cluster, return_counts=True)[1]

75 ets = {}

76 # Loop over all elements, represented by multi indices

77 for multi_index in itertools.product([0, 1, 2], repeat=order):

78 # Sort the multi index based on the multiplicities

79 # e.g. cluster [1 1 2]

80 # -> [z y x] == [y z x]

81 sorted_multi_index = []

82 j = 0

83 for m in multiplicities:

84 sorted_multi_index.extend(sorted(multi_index[j:j + m]))

85 j += m

86 sorted_multi_index = tuple(sorted_multi_index)

87 if sorted_multi_index not in ets:

88 ets[sorted_multi_index] = np.zeros([3]*order, dtype=np.int64)

89 ets[sorted_multi_index][multi_index] = 1

90 return list(ets.values())

93def assemble_new_eigentensor(eigentensors, solution):

94 new_eigentensor = np.zeros(eigentensors[0].shape, dtype=np.int64)

95 for parameter, eigentensor in zip(solution, eigentensors):

96 new_eigentensor += np.int64(parameter) * eigentensor

97 return new_eigentensor

100def renormalize_to_integer(vector):

101 gcm = 1 # greates common multiple

102 for v in vector:

103 if v.is_Integer:

104 continue

105 den = v.q # Denomonator

106 gcd = sympy.gcd(gcm, den) # greates common divisor

107 tmp = den / gcd

108 gcm = gcm * tmp

109 # Now all elements should be smallest possible integers

110 vector = vector * gcm

111 return vector