Coverage for /builds/debichem-team/python-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 HexagonalClosedPacked as structure 

64 elif structure == 'graphite': 

65 from ase.cluster.hexagonal import Graphite as structure 

66 else: 

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

68 raise NotImplementedError(error) 

69 

70 # Check number of surfaces 

71 nsurf = len(surfaces) 

72 if len(energies) != nsurf: 

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

74 % (nsurf,)) 

75 

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

77 energies = np.array(energies) 

78 

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

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

81 # very difficult in the general case! 

82 

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

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

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

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

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

88 latticeconstant=latticeconstant) 

89 for i, s in enumerate(surfaces): 

90 d = atoms.get_layer_distance(s) 

91 energies[i] /= d 

92 

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

94 wanted_size = size ** (1.0 / 3.0) 

95 max_e = max(energies) 

96 factor = wanted_size / max_e 

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

98 latticeconstant) 

99 if len(atoms) == 0: 

100 # Probably the cluster is very flat 

101 if debug: 

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

103 factor = 1 / energies.min() 

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

105 structure, latticeconstant) 

106 if len(atoms) == 0: 

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

108 

109 # Second guess: scale to get closer. 

110 old_factor = factor 

111 old_layers = layers 

112 old_atoms = atoms 

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

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

115 structure, latticeconstant) 

116 if len(atoms) == 0: 

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

118 atoms = old_atoms 

119 factor = old_factor 

120 layers = old_layers 

121 else: 

122 del old_atoms 

123 

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

125 below = above = None 

126 if len(atoms) <= size: 

127 below = atoms 

128 if len(atoms) >= size: 

129 above = atoms 

130 

131 # Now iterate towards the right cluster 

132 iter = 0 

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

134 if len(atoms) < size: 

135 # Find a larger cluster 

136 if debug: 

137 print('Making a larger cluster.') 

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

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

140 structure, latticeconstant) 

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

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

143 layers = new_layers 

144 else: 

145 # Find a smaller cluster 

146 if debug: 

147 print('Making a smaller cluster.') 

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

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

150 structure, latticeconstant) 

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

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

153 layers = new_layers 

154 if len(atoms) <= size: 

155 below = atoms 

156 if len(atoms) >= size: 

157 above = atoms 

158 iter += 1 

159 if iter == maxiter: 

160 raise RuntimeError('Runaway iteration.') 

161 if rounding == 'below': 

162 if debug: 

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

164 return below 

165 elif rounding == 'above': 

166 if debug: 

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

168 return above 

169 else: 

170 assert rounding == 'closest' 

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

172 atoms = above 

173 else: 

174 atoms = below 

175 if debug: 

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

177 return atoms 

178 

179 

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

181 layers1 = factor * np.array(energies) 

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

183 atoms = structure(symbol, surfaces, layers, 

184 latticeconstant=latticeconstant) 

185 return (atoms, layers)