1import numpy as np 

2 

3 

4class ClusterBase: 

5 def get_layer_distance(self, miller, layers=1, tol=1e-9, new=True): 

6 """Returns the distance between planes defined by the given miller 

7 index. 

8 """ 

9 if new: 

10 # Create lattice sample 

11 size = np.zeros(3, int) 

12 for i, m in enumerate(miller): 

13 size[i] = np.abs(m) + 2 

14 

15 m = len(self.atomic_basis) 

16 p = np.zeros((size.prod() * m, 3)) 

17 for h in range(size[0]): 

18 for k in range(size[1]): 

19 for l in range(size[2]): 

20 i = h * (size[1] * size[2]) + k * size[2] + l 

21 p[m * i:m * (i + 1)] = np.dot([h, k, l] + 

22 self.atomic_basis, 

23 self.lattice_basis) 

24 

25 # Project lattice positions on the miller direction. 

26 n = self.miller_to_direction(miller) 

27 d = np.sum(n * p, axis=1) 

28 if np.all(d < tol): 

29 # All negative incl. zero 

30 d = np.sort(np.abs(d)) 

31 reverse = True 

32 else: 

33 # Some or all positive 

34 d = np.sort(d[d > -tol]) 

35 reverse = False 

36 d = d[np.concatenate((d[1:] - d[:-1] > tol, [True]))] 

37 d = d[1:] - d[:-1] 

38 

39 # Look for a pattern in the distances between layers. A pattern is 

40 # accepted if more than 50 % of the distances obeys it. 

41 pattern = None 

42 for i in range(len(d)): 

43 for n in range(1, (len(d) - i) // 2 + 1): 

44 if np.all(np.abs(d[i:i + n] - d[i + n:i + 2 * n]) < tol): 

45 counts = 2 

46 for j in range(i + 2 * n, len(d), n): 

47 if np.all(np.abs(d[j:j + n] - d[i:i + n]) < tol): 

48 counts += 1 

49 if counts * n * 1.0 / len(d) > 0.5: 

50 pattern = d[i:i + n].copy() 

51 break 

52 if pattern is not None: 

53 break 

54 

55 if pattern is None: 

56 raise RuntimeError('Could not find layer distance for the ' + 

57 '(%i,%i,%i) surface.' % miller) 

58 if reverse: 

59 pattern = pattern[::-1] 

60 

61 if layers < 0: 

62 pattern = -1 * pattern[::-1] 

63 layers *= -1 

64 

65 map = np.arange(layers - layers % 1 + 1, dtype=int) % len(pattern) 

66 return pattern[map][:-1].sum() + layers % 1 * pattern[map][-1] 

67 

68 n = self.miller_to_direction(miller) 

69 d1 = d2 = 0.0 

70 

71 d = np.abs(np.sum(n * self.lattice_basis, axis=1)) 

72 mask = np.greater(d, 1e-10) 

73 if mask.sum() > 0: 

74 d1 = np.min(d[mask]) 

75 

76 if len(self.atomic_basis) > 1: 

77 atomic_basis = np.dot(self.atomic_basis, self.lattice_basis) 

78 d = np.sum(n * atomic_basis, axis=1) 

79 s = np.sign(d) 

80 d = np.abs(d) 

81 mask = np.greater(d, 1e-10) 

82 if mask.sum() > 0: 

83 d2 = np.min(d[mask]) 

84 s2 = s[mask][np.argmin(d[mask])] 

85 

86 if d2 > 1e-10: 

87 if s2 < 0 and d1 - d2 > 1e-10: 

88 d2 = d1 - d2 

89 elif s2 < 0 and d2 - d1 > 1e-10: 

90 d2 = 2 * d1 - d2 

91 elif s2 > 0 and d2 - d1 > 1e-10: 

92 d2 = d2 - d1 

93 

94 if np.abs(d1 - d2) < 1e-10: 

95 ld = np.array([d1]) 

96 elif np.abs(d1 - 2 * d2) < 1e-10: 

97 ld = np.array([d2]) 

98 else: 

99 assert d1 > d2, 'Something is wrong with the layer distance.' 

100 ld = np.array([d2, d1 - d2]) 

101 else: 

102 ld = np.array([d1]) 

103 

104 if len(ld) > 1: 

105 if layers < 0: 

106 ld = np.array([-ld[1], -ld[0]]) 

107 layers *= -1 

108 

109 map = np.arange(layers - (layers % 1), dtype=int) % len(ld) 

110 r = ld[map].sum() + (layers % 1) * ld[np.abs(map[-1] - 1)] 

111 else: 

112 r = ld[0] * layers 

113 

114 return r 

115 

116 def miller_to_direction(self, miller, norm=True): 

117 """Returns the direction corresponding to a given Miller index.""" 

118 d = np.dot(miller, self.resiproc_basis) 

119 if norm: 

120 d = d / np.linalg.norm(d) 

121 return d