Coverage for hiphive/cutoffs.py: 90%

1import pickle 

2from typing import Union, BinaryIO, TextIO 

3import numpy as np 

4from ase import Atoms 

5from ase.neighborlist import NeighborList 

6from hiphive.input_output.pretty_table_prints import table_array_to_string 

7 

8 

9class Cutoffs: 

10 """ This class maintains information about the cutoff configuration, 

11 i.e. which clusters will be included ("inside cutoff"). It also 

12 encapsulates functionality that is used e.g., during cluster space 

13 construction. 

14 

15 Here, `n-body` refers to number of atoms in a cluster. For example 

16 the cluster (0011) is a two-body cluster of fourth order and the 

17 cluster (123) is a three-body cluster of third order. 

18 

19 Parameters 

20 ---------- 

21 cutoff_matrix : numpy.ndarray 

22 The matrix element `ij` provides to the cutoff for order `j+2` 

23 and nbody `i+2`; elements with `i>j` will be ignored. 

24 """ 

25 

26 def __init__(self, cutoff_matrix): 

27 

28 self._cutoff_matrix = np.array(cutoff_matrix, dtype=float) 

29 

30 if len(self._cutoff_matrix.shape) != 2: 30 ↛ 31line 30 didn't jump to line 31 because the condition on line 30 was never true

31 raise ValueError('Please specify cutoff matrix as a 2D array') 

32 for i, row in enumerate(self._cutoff_matrix): 

33 if np.any(row[i:] < 0): 33 ↛ 34line 33 didn't jump to line 34 because the condition on line 33 was never true

34 raise ValueError('Negative number as cutoff') 

35 row[:i] = np.nan 

36 self._cutoff_matrix = self._cutoff_matrix[:(self.max_nbody-1), :(self.max_order-1)] 

37 

38 @property 

39 def cutoff_matrix(self) -> np.ndarray: 

40 """ Copy of cutoff matrix. """ 

41 return self._cutoff_matrix.copy() 

42 

43 @property 

44 def orders(self) -> list[int]: 

45 """ Allowed orders. """ 

46 return list(range(2, self.max_order + 1)) 

47 

48 @property 

49 def nbodies(self) -> list[int]: 

50 """ Allowed bodies. """ 

51 return list(range(2, self.max_nbody + 1)) 

52 

53 def get_cutoff(self, order: int, nbody: int) -> float: 

54 """ 

55 Returns the cutoff for a given body and order. 

56 

57 Parameters 

58 ---------- 

59 order 

60 nbody 

61 

62 Raises 

63 ------ 

64 ValueError 

65 If `order` is not in orders. 

66 ValueError 

67 If `nbody` is not in nbodies. 

68 ValueError 

69 If `nbody` is larger than order. 

70 """ 

71 if order not in self.orders: 71 ↛ 72line 71 didn't jump to line 72 because the condition on line 71 was never true

72 raise ValueError('order not in orders') 

73 if nbody not in self.nbodies: 73 ↛ 74line 73 didn't jump to line 74 because the condition on line 73 was never true

74 raise ValueError('nbody not in nbodies') 

75 if nbody > order: 75 ↛ 76line 75 didn't jump to line 76 because the condition on line 75 was never true

76 raise ValueError('nbody can not be larger than order') 

77 return self._cutoff_matrix[nbody - 2, order - 2] 

78 

79 @property 

80 def max_cutoff(self) -> float: 

81 """ Maximum cutoff. """ 

82 max_cutoff = 0 

83 for i, row in enumerate(self._cutoff_matrix): 

84 max_cutoff = max(max_cutoff, np.max(row[i:])) 

85 return max_cutoff 

86 

87 @property 

88 def max_nbody(self) -> int: 

89 """ Maximum body. """ 

90 nbody = 1 

91 for i, row in enumerate(self._cutoff_matrix): 

92 if np.any(row[i:]): 92 ↛ 91line 92 didn't jump to line 91 because the condition on line 92 was always true

93 nbody = i + 2 

94 return nbody 

95 

96 @property 

97 def max_order(self) -> int: 

98 """ Maximum order. """ 

99 order = None 

100 for col in range(self._cutoff_matrix.shape[1]): 

101 if np.any(self._cutoff_matrix[:col + 1, col]): 101 ↛ 100line 101 didn't jump to line 100 because the condition on line 101 was always true

102 order = col + 2 

103 return order 

104 

105 def max_nbody_cutoff(self, nbody: int): 

106 """ Return maximum cutoff for a given body. 

107 

108 Parameters 

109 ---------- 

110 nbody 

111 """ 

112 if nbody not in self.nbodies: 112 ↛ 113line 112 didn't jump to line 113 because the condition on line 112 was never true

113 raise ValueError('nbody not in nbodies.') 

114 return np.max(self._cutoff_matrix[nbody - 2, max(0, nbody - 2):]) 

115 

116 def max_nbody_order(self, nbody: int): 

117 """ Returns maximum order for a given body. 

118 

119 Parameters 

120 ---------- 

121 nbody 

122 """ 

123 if nbody not in self.nbodies: 123 ↛ 124line 123 didn't jump to line 124 because the condition on line 123 was never true

124 raise ValueError('nbody not in nbodies.') 

125 row = self._cutoff_matrix[nbody - 2] 

126 max_order = None 

127 for order, cutoff in enumerate(row[nbody-2:], start=nbody): 

128 if cutoff: 128 ↛ 127line 128 didn't jump to line 127 because the condition on line 128 was always true

129 max_order = order 

130 return max_order 

131 

132 def write(self, fileobj: Union[BinaryIO, TextIO]): 

133 """ Writes instance to file. 

134 

135 Parameters 

136 ---------- 

137 fileobj 

138 File-like object to which the cutoffs will be written. 

139 """ 

140 pickle.dump(self._cutoff_matrix, fileobj) 

141 

142 def read(fileobj: Union[BinaryIO, TextIO]): 

143 """ Reads a :class:`Cutoffs` instance from file. 

144 

145 Parameters 

146 ---------- 

147 fileobj 

148 Input file to read from. 

149 """ 

150 data = pickle.load(fileobj) 

151 return Cutoffs(data) 

152 

153 def to_filename_tag(self) -> str: 

154 """ Simple function turning cutoffs into a string to be used in, 

155 e.g., file names. """ 

156 s = [] 

157 for i, c in enumerate(self._cutoff_matrix.tolist(), start=2): 

158 s.append('{}body-{}'.format(i, '_'.join(map(str, c)))) 

159 return '_'.join(s) 

160 

161 def __str__(self): 

162 cutoff_matrix = self._cutoff_matrix.copy() 

163 cutoff_matrix = np.vstack(([[None] * len(self.orders)], cutoff_matrix)) 

164 s = table_array_to_string(cutoff_matrix) 

165 

166 width = max(len(c) for c in s.split('\n')) 

167 header = ' Cutoffs '.center(width, '=') + '\n' 

168 bottom = '\n' + ''.center(width, '=') 

169 s = header + s + bottom 

170 return s 

171 

172 def __repr__(self): 

173 return 'Cutoffs({!r})'.format(self._cutoff_matrix) 

174 

175 

176class CutoffMaximumBody(Cutoffs): 

177 """ Class for specifying cutoff-list plus maximum body. 

178 

179 Usefull when creating, e.g., sixth order expansions but with 

180 only 3-body interactions. 

181 

182 Parameters 

183 ---------- 

184 cutoff_list : list[float] 

185 List of cutoffs for order 2, 3, etc. Must be in decresing order. 

186 max_nbody : int 

187 No clusters containing more than `max_nbody` atoms will be generated. 

188 """ 

189 

190 def __init__(self, cutoff_list, max_nbody): 

191 cutoff_matrix = np.zeros((max_nbody - 1, len(cutoff_list))) 

192 for order, cutoff in enumerate(cutoff_list, start=2): 

193 cutoff_matrix[:, order - 2] = cutoff 

194 super().__init__(cutoff_matrix) 

195 

196 

197def is_cutoff_allowed(atoms: Atoms, cutoff: float) -> bool: 

198 """ Checks if atoms is compatible with cutoff. 

199 

200 Parameters 

201 ---------- 

202 atoms 

203 Structure used for checking compatibility with cutoff. 

204 cutoff 

205 Cutoff to be tested. 

206 

207 Returns 

208 ------- 

209 True if `cutoff` is compatible with `atoms`, else False. 

210 """ 

211 nbrlist = NeighborList(cutoffs=[cutoff / 2] * len(atoms), skin=0, 

212 self_interaction=False, bothways=True) 

213 nbrlist.update(atoms) 

214 

215 for i in range(len(atoms)): 

216 neighbors, _ = nbrlist.get_neighbors(i) 

217 if i in neighbors: 

218 return False 

219 if len(neighbors) != len(set(neighbors)): 

220 return False 

221 return True 

222 

223 

224def estimate_maximum_cutoff(atoms: Atoms, max_iter: int = 11) -> float: 

225 """ Estimates the maximum possible cutoff given the atoms object. 

226 

227 Parameters 

228 ---------- 

229 atoms 

230 Structure used for checking compatibility with cutoff. 

231 max_iter 

232 Number of iterations in binary search. 

233 """ 

234 

235 # First upper boundary of cutoff 

236 upper_cutoff = min(np.linalg.norm(atoms.cell, axis=1)) 

237 

238 # generate all possible offsets given upper_cutoff 

239 nbrlist = NeighborList(cutoffs=[upper_cutoff / 2] * len(atoms), skin=0, 

240 self_interaction=False, bothways=True) 

241 nbrlist.update(atoms) 

242 all_offsets = [] 

243 for i in range(len(atoms)): 

244 _, offsets = nbrlist.get_neighbors(i) 

245 all_offsets.extend([tuple(offset) for offset in offsets]) 

246 

247 # find lower boundary and new upper boundary 

248 unique_offsets = set(all_offsets) 

249 unique_offsets.discard((0, 0, 0)) 

250 upper = min(np.linalg.norm(np.dot(offset, atoms.cell)) 

251 for offset in unique_offsets) 

252 lower = upper / 2.0 

253 

254 # run binary search between the upper and lower bounds 

255 for _ in range(max_iter): 

256 cutoff = (upper + lower) / 2 

257 if is_cutoff_allowed(atoms, cutoff): 

258 lower = cutoff 

259 else: 

260 upper = cutoff 

261 return lower 

262 

263 

264class BaseClusterFilter: 

265 """Base cluster filter class. 

266 

267 This filter simply accepts all proposed clusters. A proper 

268 subclass must implement the same methods. 

269 """ 

270 def setup(self, atoms: Atoms): 

271 """ The filter is passed the environment of the primitive cell. 

272 

273 Parameters 

274 ---------- 

275 atoms 

276 Non-pbc primitive cell plus neighboring atoms. 

277 """ 

278 self._atoms = atoms 

279 

280 def __call__(self, cluster: tuple[int]): 

281 """ Returns True or False when a cluster is proposed. 

282 

283 Parameters 

284 ---------- 

285 cluster 

286 Indices of proposed cluster referenced to the internal 

287 :class:`Atoms` object. 

288 """ 

289 return True