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
|
import numpy as np
delta = 1e-10
def wulff_construction(symbol, surfaces, energies, size, structure,
rounding='closest', latticeconstant=None,
debug=False, maxiter=100):
"""Create a cluster using the Wulff construction.
A cluster is created with approximately the number of atoms
specified, following the Wulff construction, i.e. minimizing the
surface energy of the cluster.
Parameters:
symbol: The chemical symbol (or atomic number) of the desired element.
surfaces: A list of surfaces. Each surface is an (h, k, l) tuple or
list of integers.
energies: A list of surface energies for the surfaces.
size: The desired number of atoms.
structure: The desired crystal structure. One of the strings
"fcc", "bcc", or "sc".
rounding (optional): Specifies what should be done if no Wulff
construction corresponds to exactly the requested number of atoms.
Should be a string, either "above", "below" or "closest" (the
default), meaning that the nearest cluster above or below - or the
closest one - is created instead.
latticeconstant (optional): The lattice constant. If not given,
extracted from ase.data.
debug (optional): If non-zero, information about the iteration towards
the right cluster size is printed.
"""
if debug:
print('Wulff: Aiming for cluster with %i atoms (%s)' %
(size, rounding))
if rounding not in ['above', 'below', 'closest']:
raise ValueError('Invalid rounding: %s' % rounding)
# Interpret structure, if it is a string.
if isinstance(structure, str):
if structure == 'fcc':
from ase.cluster.cubic import FaceCenteredCubic as structure
elif structure == 'bcc':
from ase.cluster.cubic import BodyCenteredCubic as structure
elif structure == 'sc':
from ase.cluster.cubic import SimpleCubic as structure
elif structure == 'hcp':
from ase.cluster.hexagonal import \
HexagonalClosedPacked as structure
elif structure == 'graphite':
from ase.cluster.hexagonal import Graphite as structure
else:
error = 'Crystal structure %s is not supported.' % structure
raise NotImplementedError(error)
# Check number of surfaces
nsurf = len(surfaces)
if len(energies) != nsurf:
raise ValueError('The energies array should contain %d values.'
% (nsurf,))
# Copy energies array so it is safe to modify it
energies = np.array(energies)
# We should check that for each direction, the surface energy plus
# the energy in the opposite direction is positive. But this is
# very difficult in the general case!
# Before starting, make a fake cluster just to extract the
# interlayer distances in the relevant directions, and use these
# to "renormalize" the surface energies such that they can be used
# to convert to number of layers instead of to distances.
atoms = structure(symbol, surfaces, 5 * np.ones(len(surfaces), int),
latticeconstant=latticeconstant)
for i, s in enumerate(surfaces):
d = atoms.get_layer_distance(s)
energies[i] /= d
# First guess a size that is not too large.
wanted_size = size ** (1.0 / 3.0)
max_e = max(energies)
factor = wanted_size / max_e
atoms, layers = make_atoms(symbol, surfaces, energies, factor, structure,
latticeconstant)
if len(atoms) == 0:
# Probably the cluster is very flat
if debug:
print('First try made an empty cluster, trying again.')
factor = 1 / energies.min()
atoms, layers = make_atoms(symbol, surfaces, energies, factor,
structure, latticeconstant)
if len(atoms) == 0:
raise RuntimeError('Failed to create a finite cluster.')
# Second guess: scale to get closer.
old_factor = factor
old_layers = layers
old_atoms = atoms
factor *= (size / len(atoms))**(1.0 / 3.0)
atoms, layers = make_atoms(symbol, surfaces, energies, factor,
structure, latticeconstant)
if len(atoms) == 0:
print('Second guess gave an empty cluster, discarding it.')
atoms = old_atoms
factor = old_factor
layers = old_layers
else:
del old_atoms
# Find if the cluster is too small or too large (both means perfect!)
below = above = None
if len(atoms) <= size:
below = atoms
if len(atoms) >= size:
above = atoms
# Now iterate towards the right cluster
iter = 0
while (below is None or above is None):
if len(atoms) < size:
# Find a larger cluster
if debug:
print('Making a larger cluster.')
factor = ((layers + 0.5 + delta) / energies).min()
atoms, new_layers = make_atoms(symbol, surfaces, energies, factor,
structure, latticeconstant)
assert (new_layers - layers).max() == 1
assert (new_layers - layers).min() >= 0
layers = new_layers
else:
# Find a smaller cluster
if debug:
print('Making a smaller cluster.')
factor = ((layers - 0.5 - delta) / energies).max()
atoms, new_layers = make_atoms(symbol, surfaces, energies, factor,
structure, latticeconstant)
assert (new_layers - layers).max() <= 0
assert (new_layers - layers).min() == -1
layers = new_layers
if len(atoms) <= size:
below = atoms
if len(atoms) >= size:
above = atoms
iter += 1
if iter == maxiter:
raise RuntimeError('Runaway iteration.')
if rounding == 'below':
if debug:
print('Choosing smaller cluster with %i atoms' % len(below))
return below
elif rounding == 'above':
if debug:
print('Choosing larger cluster with %i atoms' % len(above))
return above
else:
assert rounding == 'closest'
if (len(above) - size) < (size - len(below)):
atoms = above
else:
atoms = below
if debug:
print('Choosing closest cluster with %i atoms' % len(atoms))
return atoms
def make_atoms(symbol, surfaces, energies, factor, structure, latticeconstant):
layers1 = factor * np.array(energies)
layers = np.round(layers1).astype(int)
atoms = structure(symbol, surfaces, layers,
latticeconstant=latticeconstant)
return (atoms, layers)
|
