Coverage for /builds/ase/ase/ase/cluster/base.py: 54.02%

1# fmt: off 

2 

3import numpy as np 

4 

5 

6class ClusterBase: 

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

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

9 index. 

10 """ 

11 if new: 

12 # Create lattice sample 

13 size = np.zeros(3, int) 

14 for i, m in enumerate(miller): 

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

16 

17 m = len(self.atomic_basis) 

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

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

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

21 for l_ in range(size[2]): 

22 i = h * (size[1] * size[2]) + k * size[2] + l_ 

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

24 self.atomic_basis, 

25 self.lattice_basis) 

26 

27 # Project lattice positions on the miller direction. 

28 n = self.miller_to_direction(miller) 

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

30 if np.all(d < tol): 

31 # All negative incl. zero 

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

33 reverse = True 

34 else: 

35 # Some or all positive 

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

37 reverse = False 

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

39 d = d[1:] - d[:-1] 

40 

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

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

43 pattern = None 

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

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

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

47 counts = 2 

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

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

50 counts += 1 

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

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

53 break 

54 if pattern is not None: 

55 break 

56 

57 if pattern is None: 

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

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

60 if reverse: 

61 pattern = pattern[::-1] 

62 

63 if layers < 0: 

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

65 layers *= -1 

66 

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

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

69 

70 n = self.miller_to_direction(miller) 

71 d1 = d2 = 0.0 

72 

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

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

75 if mask.sum() > 0: 

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

77 

78 if len(self.atomic_basis) > 1: 

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

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

81 s = np.sign(d) 

82 d = np.abs(d) 

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

84 if mask.sum() > 0: 

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

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

87 

88 if d2 > 1e-10: 

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

90 d2 = d1 - d2 

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

92 d2 = 2 * d1 - d2 

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

94 d2 = d2 - d1 

95 

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

97 ld = np.array([d1]) 

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

99 ld = np.array([d2]) 

100 else: 

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

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

103 else: 

104 ld = np.array([d1]) 

105 

106 if len(ld) > 1: 

107 if layers < 0: 

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

109 layers *= -1 

110 

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

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

113 else: 

114 r = ld[0] * layers 

115 

116 return r 

117 

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

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

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

121 if norm: 

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

123 return d