1""" 

2Implements functions for extracting ('isolating') a low-dimensional material 

3component in its own unit cell. 

4 

5This uses the rank-determination method described in: 

6Definition of a scoring parameter to identify low-dimensional materials 

7components 

8P.M. Larsen, M. Pandey, M. Strange, and K. W. Jacobsen 

9Phys. Rev. Materials 3 034003, 2019 

10https://doi.org/10.1103/PhysRevMaterials.3.034003 

11""" 

12import itertools 

13import collections 

14import numpy as np 

15 

16from ase import Atoms 

17from ase.geometry.cell import complete_cell 

18from ase.geometry.dimensionality import analyze_dimensionality 

19from ase.geometry.dimensionality import rank_determination 

20from ase.geometry.dimensionality.bond_generator import next_bond 

21from ase.geometry.dimensionality.interval_analysis import merge_intervals 

22 

23 

24def orthogonal_basis(X, Y=None): 

25 

26 is_1d = Y is None 

27 b = np.zeros((3, 3)) 

28 b[0] = X 

29 if not is_1d: 

30 b[1] = Y 

31 b = complete_cell(b) 

32 

33 Q = np.linalg.qr(b.T)[0].T 

34 if np.dot(b[0], Q[0]) < 0: 

35 Q[0] = -Q[0] 

36 

37 if np.dot(b[2], Q[1]) < 0: 

38 Q[1] = -Q[1] 

39 

40 if np.linalg.det(Q) < 0: 

41 Q[2] = -Q[2] 

42 

43 if is_1d: 

44 Q = Q[[1, 2, 0]] 

45 

46 return Q 

47 

48 

49def select_cutoff(atoms): 

50 intervals = analyze_dimensionality(atoms, method='RDA', merge=False) 

51 dimtype = max(merge_intervals(intervals), key=lambda x: x.score).dimtype 

52 m = next(e for e in intervals if e.dimtype == dimtype) 

53 if m.b == float("inf"): 

54 return m.a + 0.1 

55 else: 

56 return (m.a + m.b) / 2 

57 

58 

59def traverse_graph(atoms, kcutoff): 

60 if kcutoff is None: 

61 kcutoff = select_cutoff(atoms) 

62 

63 rda = rank_determination.RDA(len(atoms)) 

64 for (k, i, j, offset) in next_bond(atoms): 

65 if k > kcutoff: 

66 break 

67 rda.insert_bond(i, j, offset) 

68 

69 rda.check() 

70 return rda.graph.find_all(), rda.all_visited, rda.ranks 

71 

72 

73def build_supercomponent(atoms, components, k, v, anchor=True): 

74 # build supercomponent by mapping components into visited cells 

75 positions = [] 

76 numbers = [] 

77 

78 for c, offset in dict(v[::-1]).items(): 

79 indices = np.where(components == c)[0] 

80 ps = atoms.positions + np.dot(offset, atoms.get_cell()) 

81 positions.extend(ps[indices]) 

82 numbers.extend(atoms.numbers[indices]) 

83 positions = np.array(positions) 

84 numbers = np.array(numbers) 

85 

86 # select an 'anchor' atom, which will lie at the origin 

87 anchor_index = next((i for i in range(len(atoms)) if components[i] == k)) 

88 if anchor: 

89 positions -= atoms.positions[anchor_index] 

90 return positions, numbers 

91 

92 

93def select_chain_rotation(scaled): 

94 

95 best = (-1, [1, 0, 0]) 

96 for s in scaled: 

97 vhat = np.array([s[0], s[1], 0]) 

98 norm = np.linalg.norm(vhat) 

99 if norm < 1E-6: 

100 continue 

101 vhat /= norm 

102 obj = np.sum(np.dot(scaled, vhat)**2) 

103 best = max(best, (obj, vhat), key=lambda x: x[0]) 

104 _, vhat = best 

105 cost, sint, _ = vhat 

106 rot = np.array([[cost, -sint, 0], [sint, cost, 0], [0, 0, 1]]) 

107 return np.dot(scaled, rot) 

108 

109 

110def isolate_chain(atoms, components, k, v): 

111 

112 # identify the vector along the chain; this is the new cell vector 

113 basis_points = np.array([offset for c, offset in v if c == k]) 

114 assert len(basis_points) >= 2 

115 assert (0, 0, 0) in [tuple(e) for e in basis_points] 

116 

117 sizes = np.linalg.norm(basis_points, axis=1) 

118 index = np.argsort(sizes)[1] 

119 basis = basis_points[index] 

120 vector = np.dot(basis, atoms.get_cell()) 

121 norm = np.linalg.norm(vector) 

122 vhat = vector / norm 

123 

124 # project atoms into new basis 

125 positions, numbers = build_supercomponent(atoms, components, k, v) 

126 scaled = np.dot(positions, orthogonal_basis(vhat).T / norm) 

127 

128 # move atoms into new cell 

129 scaled[:, 2] %= 1.0 

130 

131 # subtract barycentre in x and y directions 

132 scaled[:, :2] -= np.mean(scaled, axis=0)[:2] 

133 

134 # pick a good chain rotation (i.e. non-random) 

135 scaled = select_chain_rotation(scaled) 

136 

137 # make cell large enough in x and y directions 

138 init_cell = norm * np.eye(3) 

139 pos = np.dot(scaled, init_cell) 

140 rmax = np.max(np.linalg.norm(pos[:, :2], axis=1)) 

141 rmax = max(1, rmax) 

142 cell = np.diag([4 * rmax, 4 * rmax, norm]) 

143 

144 # construct a new atoms object containing the isolated chain 

145 return Atoms(numbers=numbers, positions=pos, cell=cell, pbc=[0, 0, 1]) 

146 

147 

148def construct_inplane_basis(atoms, k, v): 

149 

150 basis_points = np.array([offset for c, offset in v if c == k]) 

151 assert len(basis_points) >= 3 

152 assert (0, 0, 0) in [tuple(e) for e in basis_points] 

153 

154 sizes = np.linalg.norm(basis_points, axis=1) 

155 indices = np.argsort(sizes) 

156 basis_points = basis_points[indices] 

157 

158 # identify primitive basis 

159 best = (float("inf"), None) 

160 for u, v in itertools.combinations(basis_points, 2): 

161 

162 basis = np.array([[0, 0, 0], u, v]) 

163 if np.linalg.matrix_rank(basis) < 2: 

164 continue 

165 

166 a = np.dot(u, atoms.get_cell()) 

167 b = np.dot(v, atoms.get_cell()) 

168 norm = np.linalg.norm(np.cross(a, b)) 

169 best = min(best, (norm, a, b), key=lambda x: x[0]) 

170 _, a, b = best 

171 return a, b, orthogonal_basis(a, b) 

172 

173 

174def isolate_monolayer(atoms, components, k, v): 

175 

176 a, b, basis = construct_inplane_basis(atoms, k, v) 

177 

178 # project atoms into new basis 

179 c = np.cross(a, b) 

180 c /= np.linalg.norm(c) 

181 init_cell = np.dot(np.array([a, b, c]), basis.T) 

182 

183 positions, numbers = build_supercomponent(atoms, components, k, v) 

184 scaled = np.linalg.solve(init_cell.T, np.dot(positions, basis.T).T).T 

185 

186 # move atoms into new cell 

187 scaled[:, :2] %= 1.0 

188 

189 # subtract barycentre in z-direction 

190 scaled[:, 2] -= np.mean(scaled, axis=0)[2] 

191 

192 # make cell large enough in z-direction 

193 pos = np.dot(scaled, init_cell) 

194 zmax = np.max(np.abs(pos[:, 2])) 

195 cell = np.copy(init_cell) 

196 cell[2] *= 4 * zmax 

197 

198 # construct a new atoms object containing the isolated chain 

199 return Atoms(numbers=numbers, positions=pos, cell=cell, pbc=[1, 1, 0]) 

200 

201 

202def isolate_bulk(atoms, components, k, v): 

203 positions, numbers = build_supercomponent(atoms, components, k, v, 

204 anchor=False) 

205 atoms = Atoms(numbers=numbers, positions=positions, cell=atoms.cell, 

206 pbc=[1, 1, 1]) 

207 atoms.wrap(eps=0) 

208 return atoms 

209 

210 

211def isolate_cluster(atoms, components, k, v): 

212 positions, numbers = build_supercomponent(atoms, components, k, v) 

213 positions -= np.min(positions, axis=0) 

214 cell = np.diag(np.max(positions, axis=0)) 

215 return Atoms(numbers=numbers, positions=positions, cell=cell, pbc=False) 

216 

217 

218def isolate_components(atoms, kcutoff=None): 

219 """Isolates components by dimensionality type. 

220 

221 Given a k-value cutoff the components (connected clusters) are 

222 identified. For each component an Atoms object is created, which contains 

223 that component only. The geometry of the resulting Atoms object depends 

224 on the component dimensionality type: 

225 

226 0D: The cell is a tight box around the atoms. pbc=[0, 0, 0]. 

227 The cell has no physical meaning. 

228 

229 1D: The chain is aligned along the z-axis. pbc=[0, 0, 1]. 

230 The x and y cell directions have no physical meaning. 

231 

232 2D: The layer is aligned in the x-y plane. pbc=[1, 1, 0]. 

233 The z cell direction has no physical meaning. 

234 

235 3D: The original cell is used. pbc=[1, 1, 1]. 

236 

237 Parameters: 

238 

239 atoms: ASE atoms object 

240 The system to analyze. 

241 kcutoff: float 

242 The k-value cutoff to use. Default=None, in which case the 

243 dimensionality scoring parameter is used to select the cutoff. 

244 

245 Returns: 

246 

247 components: dict 

248 key: the component dimenionalities. 

249 values: a list of Atoms objects for each dimensionality type. 

250 """ 

251 data = {} 

252 components, all_visited, ranks = traverse_graph(atoms, kcutoff) 

253 

254 for k, v in all_visited.items(): 

255 v = sorted(list(v)) 

256 

257 # identify the components which constitute the component 

258 key = tuple(np.unique([c for c, offset in v])) 

259 dim = ranks[k] 

260 

261 if dim == 0: 

262 data[('0D', key)] = isolate_cluster(atoms, components, k, v) 

263 elif dim == 1: 

264 data[('1D', key)] = isolate_chain(atoms, components, k, v) 

265 elif dim == 2: 

266 data[('2D', key)] = isolate_monolayer(atoms, components, k, v) 

267 elif dim == 3: 

268 data[('3D', key)] = isolate_bulk(atoms, components, k, v) 

269 

270 result = collections.defaultdict(list) 

271 for (dim, _), atoms in data.items(): 

272 result[dim].append(atoms) 

273 return result