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: The chemical symbol (or atomic number) of the desired element. 

19 

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

21 list of integers. 

22 

23 energies: A list of surface energies for the surfaces. 

24 

25 size: The desired number of atoms. 

26 

27 structure: The desired crystal structure. One of the strings 

28 "fcc", "bcc", or "sc". 

29 

30 rounding (optional): Specifies what should be done if no Wulff 

31 construction corresponds to exactly the requested number of atoms. 

32 Should be a string, either "above", "below" or "closest" (the 

33 default), meaning that the nearest cluster above or below - or the 

34 closest one - is created instead. 

35 

36 latticeconstant (optional): The lattice constant. If not given, 

37 extracted from ase.data. 

38 

39 debug (optional): If non-zero, information about the iteration towards 

40 the right cluster size is printed. 

41 """ 

42 

43 if debug: 

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

45 (size, rounding)) 

46 

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

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

49 

50 # Interpret structure, if it is a string. 

51 if isinstance(structure, str): 

52 if structure == 'fcc': 

53 from ase.cluster.cubic import FaceCenteredCubic as structure 

54 elif structure == 'bcc': 

55 from ase.cluster.cubic import BodyCenteredCubic as structure 

56 elif structure == 'sc': 

57 from ase.cluster.cubic import SimpleCubic as structure 

58 elif structure == 'hcp': 

59 from ase.cluster.hexagonal import \ 

60 HexagonalClosedPacked as structure 

61 elif structure == 'graphite': 

62 from ase.cluster.hexagonal import Graphite as structure 

63 else: 

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

65 raise NotImplementedError(error) 

66 

67 # Check number of surfaces 

68 nsurf = len(surfaces) 

69 if len(energies) != nsurf: 

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

71 % (nsurf,)) 

72 

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

74 energies = np.array(energies) 

75 

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

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

78 # very difficult in the general case! 

79 

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

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

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

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

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

85 latticeconstant=latticeconstant) 

86 for i, s in enumerate(surfaces): 

87 d = atoms.get_layer_distance(s) 

88 energies[i] /= d 

89 

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

91 wanted_size = size ** (1.0 / 3.0) 

92 max_e = max(energies) 

93 factor = wanted_size / max_e 

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

95 latticeconstant) 

96 if len(atoms) == 0: 

97 # Probably the cluster is very flat 

98 if debug: 

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

100 factor = 1 / energies.min() 

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

102 structure, latticeconstant) 

103 if len(atoms) == 0: 

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

105 

106 # Second guess: scale to get closer. 

107 old_factor = factor 

108 old_layers = layers 

109 old_atoms = atoms 

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

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

112 structure, latticeconstant) 

113 if len(atoms) == 0: 

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

115 atoms = old_atoms 

116 factor = old_factor 

117 layers = old_layers 

118 else: 

119 del old_atoms 

120 

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

122 below = above = None 

123 if len(atoms) <= size: 

124 below = atoms 

125 if len(atoms) >= size: 

126 above = atoms 

127 

128 # Now iterate towards the right cluster 

129 iter = 0 

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

131 if len(atoms) < size: 

132 # Find a larger cluster 

133 if debug: 

134 print('Making a larger cluster.') 

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

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

137 structure, latticeconstant) 

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

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

140 layers = new_layers 

141 else: 

142 # Find a smaller cluster 

143 if debug: 

144 print('Making a smaller cluster.') 

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

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

147 structure, latticeconstant) 

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

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

150 layers = new_layers 

151 if len(atoms) <= size: 

152 below = atoms 

153 if len(atoms) >= size: 

154 above = atoms 

155 iter += 1 

156 if iter == maxiter: 

157 raise RuntimeError('Runaway iteration.') 

158 if rounding == 'below': 

159 if debug: 

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

161 return below 

162 elif rounding == 'above': 

163 if debug: 

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

165 return above 

166 else: 

167 assert rounding == 'closest' 

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

169 atoms = above 

170 else: 

171 atoms = below 

172 if debug: 

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

174 return atoms 

175 

176 

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

178 layers1 = factor * np.array(energies) 

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

180 atoms = structure(symbol, surfaces, layers, 

181 latticeconstant=latticeconstant) 

182 return (atoms, layers)