Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

""" 

Functionality for enforcing translational sum rules 

""" 

 

import sympy 

from scipy.sparse import coo_matrix 

import numpy as np 

 

from ..input_output.logging_tools import logger 

from .utilities import SparseMatrix, BiMap 

from .orbits import get_orbits 

 

 

logger = logger.getChild('sum_rules') 

 

 

def create_constraint_map(cs): 

 

# First we need to know the global index of a eigentensor given which 

# orbit it belongs to. 

params = [] 

n = 0 

for orbit_index, orbit in enumerate(cs.orbits): 

nParams_in_orbit = len(orbit.eigentensors) 

params.append(list(range(n, n + nParams_in_orbit))) 

n += nParams_in_orbit 

 

# If we dont want the sum rules an eye matrix is created instead. 

if not cs.acoustic_sum_rules: 

row = list(range(n)) 

col = list(range(n)) 

data = [1.0]*n 

cs._cvs = coo_matrix((data, (row, col)), shape=(n, n)) 

return None 

 

M = get_translational_constraint_matrix(cs) 

 

nullspace = M.nullspace() 

logger.debug('Assembling solutions...') 

row = [] 

col = [] 

data = [] 

for c, vec in enumerate(nullspace): 

for r, _, v in vec.row_list(): 

row.append(r) 

col.append(c) 

data.append(np.float64(v)) 

# This is the final product, the constraint map (or constraint vectors) 

cs._cvs = coo_matrix((data, (row, col)), shape=(len(vec), len(nullspace))) 

 

logger.debug('{} degrees of freedom'.format(cs._cvs.shape[1])) 

logger.debug('Finished constructing constraint vectors') 

 

 

def get_translational_constraint_matrix(cs): 

logger.debug('Starting constructing constraint matrices') 

 

# First we need to know the global index of a eigentensor given which 

# orbit it belongs to. 

params = [] 

n = 0 

for orbit_index, orbit in enumerate(cs.orbits): 

nParams_in_orbit = len(orbit.eigentensors) 

params.append(list(range(n, n + nParams_in_orbit))) 

n += nParams_in_orbit 

 

# The lookup is a fast way to get the eigentensors given a cluster Also the 

# orbit index is bundled to quickly find the correct parameter indices 

lookup = {} 

for orbit_index, orbit in enumerate(cs.orbits): 

for of in orbit.orientation_families: 

for cluster_index, perm_index in zip(of.cluster_indices, 

of.permutation_indices): 

cluster = cs._cluster_list[cluster_index] 

perm = cs._permutations[perm_index] 

lookup[tuple(cluster)] = [et.transpose(perm) for et in 

of.eigentensors], orbit_index 

 

# Init just one row to be able to stack. n is the total number of 

# eigentensors 

M = SparseMatrix(1, n, 0) # This is the constraint matrix 

for order in cs.cutoffs.orders: 

logger.debug('Order {}: '.format(order)) 

 

# Find which minimal amount of prefixes that is needed 

# Essentially they are all the prototypes of order one less still 

# compatible with the cutoff. 

prefix_list = [] 

if order == 2: 

for i in np.unique(cs.wyckoff_sites): 

prefix_list.append((i,)) 

else: 

tmp_cluster_set = set() 

for cluster in cs._cluster_list: 

if len(cluster) != order: 

continue 

for i in range(order): 

c = tuple(cluster[:i]) + tuple(cluster[i+1:]) 

if c[0] < len(cs._prim): 

tmp_cluster_set.add(c) 

tmp_cluster_list = BiMap() 

for c in tmp_cluster_set: 

tmp_cluster_list.append(c) 

orbits = get_orbits(tmp_cluster_list, cs._atom_list, 

cs.rotation_matrices, cs.translation_vectors, 

cs.permutations, cs._prim, cs.symprec) 

for orb in orbits: 

prefix_list.append(tmp_cluster_list[orb.prototype_index]) 

 

m = SparseMatrix(3**order, M.shape[1], 0) # intermediate matrix 

# Now loop over the prefix and add the last atom index 

for prefix in prefix_list: 

m *= 0 

cluster = list(prefix) + [None] 

for i in range(len(cs._atom_list)): 

cluster[-1] = i 

tmp = lookup.get(tuple(sorted(cluster))) 

if tmp is None: 

continue 

ets, orbit_index = tmp 

inv_argsort = np.argsort(np.argsort(cluster)) 

ets = [et.transpose(inv_argsort) for et in ets] 

for et, col in zip(ets, params[orbit_index]): 

et_flat = et.flatten() 

for j in et_flat.nonzero()[0]: 

m[j, col] += sympy.nsimplify(et_flat[j]) 

M = SparseMatrix.vstack(M, m) 

# only compress if there are very many elements 

129 ↛ 130line 129 didn't jump to line 130, because the condition on line 129 was never true if M.nnz() > cs._config['max_number_constraint_elements']: 

M = SparseMatrix.rref(M)[0] 

 

logger.debug('Finished constructing constraint matrices') 

logger.debug('Starting constructing constraint vectors') 

return M