Coverage for /builds/kinetik161/ase/ase/cluster/wulff.py: 51.52%

1import numpy as np 

2 

3delta = 1e-10 

4 

5 

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

7 rounding='closest', latticeconstant=None, 

8 debug=False, maxiter=100): 

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

10 

11 A cluster is created with approximately the number of atoms 

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

13 surface energy of the cluster. 

14 

15 Parameters 

16 ---------- 

17 symbol : str or int 

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

19 

20 surfaces : list 

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

22 integers. 

23 

24 energies : list 

25 A list of surface energies for the surfaces. 

26 

27 size : int 

28 The desired number of atoms. 

29 

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

31 The desired crystal structure. 

32 

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

34 Specifies what should be done if no Wulff construction corresponds 

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

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

37 closest one - is created instead. 

38 

39 latticeconstant : float (optional) 

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

41 

42 debug : bool, default False 

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

44 size is printed. 

45 """ 

46 

47 if debug: 

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

49 (size, rounding)) 

50 

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

52 raise ValueError(f'Invalid rounding: {rounding}') 

53 

54 # Interpret structure, if it is a string. 

55 if isinstance(structure, str): 

56 if structure == 'fcc': 

57 from ase.cluster.cubic import FaceCenteredCubic as structure 

58 elif structure == 'bcc': 

59 from ase.cluster.cubic import BodyCenteredCubic as structure 

60 elif structure == 'sc': 

61 from ase.cluster.cubic import SimpleCubic as structure 

62 elif structure == 'hcp': 

63 from ase.cluster.hexagonal import \ 

64 HexagonalClosedPacked as structure 

65 elif structure == 'graphite': 

66 from ase.cluster.hexagonal import Graphite as structure 

67 else: 

68 error = f'Crystal structure {structure} is not supported.' 

69 raise NotImplementedError(error) 

70 

71 # Check number of surfaces 

72 nsurf = len(surfaces) 

73 if len(energies) != nsurf: 

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

75 % (nsurf,)) 

76 

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

78 energies = np.array(energies) 

79 

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

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

82 # very difficult in the general case! 

83 

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

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

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

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

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

89 latticeconstant=latticeconstant) 

90 for i, s in enumerate(surfaces): 

91 d = atoms.get_layer_distance(s) 

92 energies[i] /= d 

93 

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

95 wanted_size = size ** (1.0 / 3.0) 

96 max_e = max(energies) 

97 factor = wanted_size / max_e 

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

99 latticeconstant) 

100 if len(atoms) == 0: 

101 # Probably the cluster is very flat 

102 if debug: 

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

104 factor = 1 / energies.min() 

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

106 structure, latticeconstant) 

107 if len(atoms) == 0: 

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

109 

110 # Second guess: scale to get closer. 

111 old_factor = factor 

112 old_layers = layers 

113 old_atoms = atoms 

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

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

116 structure, latticeconstant) 

117 if len(atoms) == 0: 

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

119 atoms = old_atoms 

120 factor = old_factor 

121 layers = old_layers 

122 else: 

123 del old_atoms 

124 

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

126 below = above = None 

127 if len(atoms) <= size: 

128 below = atoms 

129 if len(atoms) >= size: 

130 above = atoms 

131 

132 # Now iterate towards the right cluster 

133 iter = 0 

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

135 if len(atoms) < size: 

136 # Find a larger cluster 

137 if debug: 

138 print('Making a larger cluster.') 

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

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

141 structure, latticeconstant) 

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

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

144 layers = new_layers 

145 else: 

146 # Find a smaller cluster 

147 if debug: 

148 print('Making a smaller cluster.') 

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

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

151 structure, latticeconstant) 

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

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

154 layers = new_layers 

155 if len(atoms) <= size: 

156 below = atoms 

157 if len(atoms) >= size: 

158 above = atoms 

159 iter += 1 

160 if iter == maxiter: 

161 raise RuntimeError('Runaway iteration.') 

162 if rounding == 'below': 

163 if debug: 

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

165 return below 

166 elif rounding == 'above': 

167 if debug: 

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

169 return above 

170 else: 

171 assert rounding == 'closest' 

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

173 atoms = above 

174 else: 

175 atoms = below 

176 if debug: 

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

178 return atoms 

179 

180 

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

182 layers1 = factor * np.array(energies) 

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

184 atoms = structure(symbol, surfaces, layers, 

185 latticeconstant=latticeconstant) 

186 return (atoms, layers)