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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
|
from __future__ import print_function
import numpy as np
from ase.data import atomic_numbers as ref_atomic_numbers
from ase.spacegroup import Spacegroup
from ase.cluster.base import ClusterBase
from ase.cluster.cluster import Cluster
class ClusterFactory(ClusterBase):
directions = [[1, 0, 0],
[0, 1, 0],
[0, 0, 1]]
atomic_basis = np.array([[0., 0., 0.]])
element_basis = None
Cluster = Cluster # Make it possible to change the class of the object returned.
def __call__(self, symbols, surfaces, layers, latticeconstant=None,
center=None, vacuum=0.0, debug=0):
self.debug = debug
# Interpret symbol
self.set_atomic_numbers(symbols)
# Interpret lattice constant
if latticeconstant is None:
if self.element_basis is None:
self.lattice_constant = self.get_lattice_constant()
else:
raise ValueError("A lattice constant must be specified for a compound")
else:
self.lattice_constant = latticeconstant
self.set_basis()
if self.debug:
print("Lattice constant(s):", self.lattice_constant)
print("Lattice basis:\n", self.lattice_basis)
print("Resiprocal basis:\n", self.resiproc_basis)
print("Atomic basis:\n", self.atomic_basis)
self.set_surfaces_layers(surfaces, layers)
self.set_lattice_size(center)
if self.debug:
print("Center position:", self.center.round(2))
print("Base lattice size:", self.size)
cluster = self.make_cluster(vacuum)
cluster.symmetry = self.xtal_name
cluster.surfaces = self.surfaces.copy()
cluster.lattice_basis = self.lattice_basis.copy()
cluster.atomic_basis = self.atomic_basis.copy()
cluster.resiproc_basis = self.resiproc_basis.copy()
return cluster
def make_cluster(self, vacuum):
# Make the base crystal by repeating the unit cell
size = np.array(self.size)
translations = np.zeros((size.prod(), 3))
for h in range(size[0]):
for k in range(size[1]):
for l in range(size[2]):
i = h * (size[1] * size[2]) + k * size[2] + l
translations[i] = np.dot([h, k, l], self.lattice_basis)
atomic_basis = np.dot(self.atomic_basis, self.lattice_basis)
positions = np.zeros((len(translations) * len(atomic_basis), 3))
numbers = np.zeros(len(positions))
n = len(atomic_basis)
for i, trans in enumerate(translations):
positions[n*i:n*(i+1)] = atomic_basis + trans
numbers[n*i:n*(i+1)] = self.atomic_numbers
# Remove all atoms that is outside the defined surfaces
for s, l in zip(self.surfaces, self.layers):
n = self.miller_to_direction(s)
rmax = self.get_layer_distance(s, l + 0.1)
r = np.dot(positions - self.center, n)
mask = np.less(r, rmax)
if self.debug > 1:
print("Cutting %s at %i layers ~ %.3f A" % (s, l, rmax))
positions = positions[mask]
numbers = numbers[mask]
# Fit the cell, so it only just consist the atoms
min = np.zeros(3)
max = np.zeros(3)
for i in range(3):
v = self.directions[i]
r = np.dot(positions, v)
min[i] = r.min()
max[i] = r.max()
cell = max - min + vacuum
positions = positions - min + vacuum / 2.0
self.center = self.center - min + vacuum / 2.0
return self.Cluster(symbols=numbers, positions=positions, cell=cell)
def set_atomic_numbers(self, symbols):
"Extract atomic number from element"
# The types that can be elements: integers and strings
atomic_numbers = []
if self.element_basis is None:
if isinstance(symbols, str):
atomic_numbers.append(ref_atomic_numbers[symbols])
elif isinstance(symbols, int):
atomic_numbers.append(symbols)
else:
raise TypeError("The symbol argument must be a " +
"string or an atomic number.")
element_basis = [0] * len(self.atomic_basis)
else:
if isinstance(symbols, (list, tuple)):
nsymbols = len(symbols)
else:
nsymbols = 0
nelement_basis = max(self.element_basis) + 1
if nsymbols != nelement_basis:
raise TypeError("The symbol argument must be a sequence " +
"of length %d" % (nelement_basis,) +
" (one for each kind of lattice position")
for s in symbols:
if isinstance(s, str):
atomic_numbers.append(ref_atomic_numbers[s])
elif isinstance(s, int):
atomic_numbers.append(s)
else:
raise TypeError("The symbol argument must be a " +
"string or an atomic number.")
element_basis = self.element_basis
self.atomic_numbers = [atomic_numbers[n] for n in element_basis]
assert len(self.atomic_numbers) == len(self.atomic_basis)
def set_lattice_size(self, center):
if center is None:
offset = np.zeros(3)
else:
offset = np.array(center)
if (offset > 1.0).any() or (offset < 0.0).any():
raise ValueError("Center offset must lie within the lattice unit \
cell.")
max = np.ones(3)
min = -np.ones(3)
v = np.linalg.inv(self.lattice_basis.T)
for s, l in zip(self.surfaces, self.layers):
n = self.miller_to_direction(s) * self.get_layer_distance(s, l)
k = np.round(np.dot(v, n), 2)
for i in range(3):
if k[i] > 0.0:
k[i] = np.ceil(k[i])
elif k[i] < 0.0:
k[i] = np.floor(k[i])
if self.debug > 1:
print("Spaning %i layers in %s in lattice basis ~ %s" % (l, s, k))
max[k > max] = k[k > max]
min[k < min] = k[k < min]
self.center = np.dot(offset - min, self.lattice_basis)
self.size = (max - min + np.ones(3)).astype(int)
def set_surfaces_layers(self, surfaces, layers):
if len(surfaces) != len(layers):
raise ValueError("Improper size of surface and layer arrays: %i != %i"
% (len(surfaces), len(layers)))
sg = Spacegroup(self.spacegroup)
surfaces = np.array(surfaces)
layers = np.array(layers)
for i, s in enumerate(surfaces):
s = reduce_miller(s)
surfaces[i] = s
surfaces_full = surfaces.copy()
layers_full = layers.copy()
for s, l in zip(surfaces, layers):
equivalent_surfaces = sg.equivalent_reflections(s.reshape(-1, 3))
for es in equivalent_surfaces:
# If the equivalent surface (es) is not in the surface list,
# then append it.
if not np.equal(es, surfaces_full).all(axis=1).any():
surfaces_full = np.append(surfaces_full, es.reshape(1, 3), axis=0)
layers_full = np.append(layers_full, l)
self.surfaces = surfaces_full.copy()
self.layers = layers_full.copy()
def get_resiproc_basis(self, basis):
"""Returns the resiprocal basis to a given lattice (crystal) basis"""
k = 1 / np.dot(basis[0], cross(basis[1], basis[2]))
# The same as the inversed basis matrix transposed
return k * np.array([cross(basis[1], basis[2]),
cross(basis[2], basis[0]),
cross(basis[0], basis[1])])
# Helping functions
def cross(a, b):
"""The cross product of two vectors."""
return np.array([a[1]*b[2] - b[1]*a[2],
a[2]*b[0] - b[2]*a[0],
a[0]*b[1] - b[0]*a[1]])
def GCD(a,b):
"""Greatest Common Divisor of a and b."""
#print "--"
while a != 0:
#print a,b,">",
a,b = b%a,a
#print a,b
return b
def reduce_miller(hkl):
"""Reduce Miller index to the lowest equivalent integers."""
hkl = np.array(hkl)
old = hkl.copy()
d = GCD(GCD(hkl[0], hkl[1]), hkl[2])
while d != 1:
hkl = hkl // d
d = GCD(GCD(hkl[0], hkl[1]), hkl[2])
if np.dot(old, hkl) > 0:
return hkl
else:
return -hkl
|
