 r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1"""

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

3component in its own unit cell.

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

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

24def orthogonal_basis(X, Y=None):

26 is_1d = Y is None

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

28 b = X

29 if not is_1d:

30 b = Y

31 b = complete_cell(b)

33 Q = np.linalg.qr(b.T).T

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

35 Q = -Q

37 if np.dot(b, Q) < 0:

38 Q = -Q

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

41 Q = -Q

43 if is_1d:

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

46 return Q

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

59def traverse_graph(atoms, kcutoff):

60 if kcutoff is None:

61 kcutoff = select_cutoff(atoms)

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)

69 rda.check()

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

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

74 # build supercomponent by mapping components into visited cells

75 positions = []

76 numbers = []

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

79 indices = np.where(components == c)

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)

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

93def select_chain_rotation(scaled):

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

96 for s in scaled:

97 vhat = np.array([s, s, 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)

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)

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

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]

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

118 index = np.argsort(sizes)

119 basis = basis_points[index]

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

121 norm = np.linalg.norm(vector)

122 vhat = vector / norm

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)

128 # move atoms into new cell

129 scaled[:, 2] %= 1.0

131 # subtract barycentre in x and y directions

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

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

135 scaled = select_chain_rotation(scaled)

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

144 # construct a new atoms object containing the isolated chain

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

148def construct_inplane_basis(atoms, k, v):

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]

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

155 indices = np.argsort(sizes)

156 basis_points = basis_points[indices]

158 # identify primitive basis

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

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

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

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

164 continue

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)

170 _, a, b = best

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

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

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

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)

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

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

186 # move atoms into new cell

187 scaled[:, :2] %= 1.0

189 # subtract barycentre in z-direction

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

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 *= 4 * zmax

198 # construct a new atoms object containing the isolated chain

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

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

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)

218def isolate_components(atoms, kcutoff=None):

219 """Isolates components by dimensionality type.

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:

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

227 The cell has no physical meaning.

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.

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

233 The z cell direction has no physical meaning.

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

237 Parameters:

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.

245 Returns:

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)

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

255 v = sorted(list(v))

257 # identify the components which constitute the component

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

259 dim = ranks[k]

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)

270 result = collections.defaultdict(list)

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

272 result[dim].append(atoms)

273 return result