1import numpy as np 

2 

3 

4delta = 1e-10 

5 

6 

7def wulff_construction(symbol, surfaces, energies, size, structure, 

8 rounding='closest', latticeconstant=None, 

9 debug=False, maxiter=100): 

10 """Create a cluster using the Wulff construction. 

11 

12 A cluster is created with approximately the number of atoms 

13 specified, following the Wulff construction, i.e. minimizing the 

14 surface energy of the cluster. 

15 

16 Parameters 

17 ---------- 

18 symbol : str or int 

19 The chemical symbol (or atomic number) of the desired element. 

20 

21 surfaces : list 

22 A list of surfaces. Each surface is an (h, k, l) tuple or list of 

23 integers. 

24 

25 energies : list 

26 A list of surface energies for the surfaces. 

27 

28 size : int 

29 The desired number of atoms. 

30 

31 structure : {'fcc', bcc', 'sc'} 

32 The desired crystal structure. 

33 

34 rounding : {'closest', 'above', 'below'} 

35 Specifies what should be done if no Wulff construction corresponds 

36 to exactly the requested number of atoms. 'above', 'below', and 

37 'closest' mean that the nearest cluster above or below - or the 

38 closest one - is created instead. 

39 

40 latticeconstant : float (optional) 

41 The lattice constant. If not given, extracted from `ase.data`. 

42 

43 debug : bool, default False 

44 If True, information about the iteration towards the right cluster 

45 size is printed. 

46 """ 

47 

48 if debug: 

49 print('Wulff: Aiming for cluster with %i atoms (%s)' % 

50 (size, rounding)) 

51 

52 if rounding not in ['above', 'below', 'closest']: 

53 raise ValueError('Invalid rounding: %s' % rounding) 

54 

55 # Interpret structure, if it is a string. 

56 if isinstance(structure, str): 

57 if structure == 'fcc': 

58 from ase.cluster.cubic import FaceCenteredCubic as structure 

59 elif structure == 'bcc': 

60 from ase.cluster.cubic import BodyCenteredCubic as structure 

61 elif structure == 'sc': 

62 from ase.cluster.cubic import SimpleCubic as structure 

63 elif structure == 'hcp': 

64 from ase.cluster.hexagonal import \ 

65 HexagonalClosedPacked as structure 

66 elif structure == 'graphite': 

67 from ase.cluster.hexagonal import Graphite as structure 

68 else: 

69 error = 'Crystal structure %s is not supported.' % structure 

70 raise NotImplementedError(error) 

71 

72 # Check number of surfaces 

73 nsurf = len(surfaces) 

74 if len(energies) != nsurf: 

75 raise ValueError('The energies array should contain %d values.' 

76 % (nsurf,)) 

77 

78 # Copy energies array so it is safe to modify it 

79 energies = np.array(energies) 

80 

81 # We should check that for each direction, the surface energy plus 

82 # the energy in the opposite direction is positive. But this is 

83 # very difficult in the general case! 

84 

85 # Before starting, make a fake cluster just to extract the 

86 # interlayer distances in the relevant directions, and use these 

87 # to "renormalize" the surface energies such that they can be used 

88 # to convert to number of layers instead of to distances. 

89 atoms = structure(symbol, surfaces, 5 * np.ones(len(surfaces), int), 

90 latticeconstant=latticeconstant) 

91 for i, s in enumerate(surfaces): 

92 d = atoms.get_layer_distance(s) 

93 energies[i] /= d 

94 

95 # First guess a size that is not too large. 

96 wanted_size = size ** (1.0 / 3.0) 

97 max_e = max(energies) 

98 factor = wanted_size / max_e 

99 atoms, layers = make_atoms(symbol, surfaces, energies, factor, structure, 

100 latticeconstant) 

101 if len(atoms) == 0: 

102 # Probably the cluster is very flat 

103 if debug: 

104 print('First try made an empty cluster, trying again.') 

105 factor = 1 / energies.min() 

106 atoms, layers = make_atoms(symbol, surfaces, energies, factor, 

107 structure, latticeconstant) 

108 if len(atoms) == 0: 

109 raise RuntimeError('Failed to create a finite cluster.') 

110 

111 # Second guess: scale to get closer. 

112 old_factor = factor 

113 old_layers = layers 

114 old_atoms = atoms 

115 factor *= (size / len(atoms))**(1.0 / 3.0) 

116 atoms, layers = make_atoms(symbol, surfaces, energies, factor, 

117 structure, latticeconstant) 

118 if len(atoms) == 0: 

119 print('Second guess gave an empty cluster, discarding it.') 

120 atoms = old_atoms 

121 factor = old_factor 

122 layers = old_layers 

123 else: 

124 del old_atoms 

125 

126 # Find if the cluster is too small or too large (both means perfect!) 

127 below = above = None 

128 if len(atoms) <= size: 

129 below = atoms 

130 if len(atoms) >= size: 

131 above = atoms 

132 

133 # Now iterate towards the right cluster 

134 iter = 0 

135 while (below is None or above is None): 

136 if len(atoms) < size: 

137 # Find a larger cluster 

138 if debug: 

139 print('Making a larger cluster.') 

140 factor = ((layers + 0.5 + delta) / energies).min() 

141 atoms, new_layers = make_atoms(symbol, surfaces, energies, factor, 

142 structure, latticeconstant) 

143 assert (new_layers - layers).max() == 1 

144 assert (new_layers - layers).min() >= 0 

145 layers = new_layers 

146 else: 

147 # Find a smaller cluster 

148 if debug: 

149 print('Making a smaller cluster.') 

150 factor = ((layers - 0.5 - delta) / energies).max() 

151 atoms, new_layers = make_atoms(symbol, surfaces, energies, factor, 

152 structure, latticeconstant) 

153 assert (new_layers - layers).max() <= 0 

154 assert (new_layers - layers).min() == -1 

155 layers = new_layers 

156 if len(atoms) <= size: 

157 below = atoms 

158 if len(atoms) >= size: 

159 above = atoms 

160 iter += 1 

161 if iter == maxiter: 

162 raise RuntimeError('Runaway iteration.') 

163 if rounding == 'below': 

164 if debug: 

165 print('Choosing smaller cluster with %i atoms' % len(below)) 

166 return below 

167 elif rounding == 'above': 

168 if debug: 

169 print('Choosing larger cluster with %i atoms' % len(above)) 

170 return above 

171 else: 

172 assert rounding == 'closest' 

173 if (len(above) - size) < (size - len(below)): 

174 atoms = above 

175 else: 

176 atoms = below 

177 if debug: 

178 print('Choosing closest cluster with %i atoms' % len(atoms)) 

179 return atoms 

180 

181 

182def make_atoms(symbol, surfaces, energies, factor, structure, latticeconstant): 

183 layers1 = factor * np.array(energies) 

184 layers = np.round(layers1).astype(int) 

185 atoms = structure(symbol, surfaces, layers, 

186 latticeconstant=latticeconstant) 

187 return (atoms, layers)