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

1import numpy as np 

2import itertools 

3import sympy 

4from .tensors import rotate_tensor 

5from ..input_output.logging_tools import logger 

6from .utilities import SparseMatrix 

7 

8 

9logger = logger.getChild('eigensymmetries') 

10 

11 

12def create_eigentensors(cs): 

13 

14 for orbit in cs.orbits: 

15 prototype_cluster = cs._cluster_list[orbit.prototype_index] 

16 

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

18 # permutation condition 

19 ets = init_ets_from_label_symmetry(prototype_cluster) 

20 

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: 

24 

25 # Rotation amtrix in scaled coordinates 

26 R = cs.rotation_matrices[rotation_index] 

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

28 

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 

32 

33 p_inv = np.argsort(p) 

34 

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) 

38 

39 nullspace = M.nullspace() 

40 

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 

51 

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

53 if len(ets) == 0: 

54 break 

55 

56 orbit.eigentensors = ets 

57 

58 

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 

66 

67 

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

91 

92 

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 

98 

99 

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