Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

import pickle 

import numpy as np 

from ase.neighborlist import NeighborList 

 

 

# TODO: Use another base class to hide internals from user 

class Cutoffs: 

""" Contains information about cutoff config, meaning which clusters 

shall be allowed (inside cutoff) and which are not. 

 

Contains conventient functionality to be used in the algorithms 

 

Note that n-body refers to number of atoms in a cluster, 

e.g. the cluster (0011) is a two-body cluster of fourth order and 

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

 

Parameters 

---------- 

cutoff_matrix 

cutoff_matrix[i][j] corresponds to the cutoff for order `i+2` and 

nbody `j+2`. Elements where j>i will be ignored. 

""" 

 

def __init__(self, cutoff_matrix): 

 

self._cutoff_matrix = np.array(cutoff_matrix) 

 

for i in range(len(self._cutoff_matrix) - 1): 

assert max(self._cutoff_matrix[i, :]) >= \ 

max(self._cutoff_matrix[i+1, :]),\ 

'Make sure cutoffs are in decreasing order' 

 

@property 

def orders(self): 

""" list : list of all orders in Cutoffs """ 

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

 

@property 

def nbodies(self): 

""" list : list of all orders in Cutoffs """ 

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

 

def get_cutoff(self, order, nbody): 

""" Returns cutoff for a given body and order. """ 

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

raise ValueError('order not in orders') 

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

raise ValueError('nbody not in nbodies') 

49 ↛ 50line 49 didn't jump to line 50, because the condition on line 49 was never true if nbody > order: 

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

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

 

@property 

def max_cutoff(self): 

""" float : maximum cutoff """ 

max_cutoff = 0 

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

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

return max_cutoff 

 

@property 

def max_nbody(self): 

""" int : maximum body """ 

return self._cutoff_matrix.shape[0] + 1 

 

@property 

def max_order(self): 

""" int : maximum order """ 

return self._cutoff_matrix.shape[1] + 1 

 

def max_nbody_cutoff(self, nbody): 

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

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

raise ValueError('nbody not in nbodies') 

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

 

def max_nbody_order(self, nbody): 

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

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

raise ValueError('nbody not in nbodies') 

arr = self._cutoff_matrix[nbody - 2, max(0, nbody - 2):] 

return np.count_nonzero(arr) + (nbody - 1) 

 

def write(self, fileobj): 

""" Write Cutoffs to file in pickle format. 

 

Parameters 

---------- 

fileobj : file-like object 

file-like object to which the cutoffs will be written to 

""" 

pickle.dump(self._cutoff_matrix, fileobj) 

 

def read(fileobj): 

""" Reads a Cutoffs instance from file. 

 

Parameters 

---------- 

fileobj : file-like object 

input file to read from 

""" 

data = pickle.load(fileobj) 

103 ↛ 106line 103 didn't jump to line 106, because the condition on line 103 was never false if type(data) is np.ndarray: 

return Cutoffs(data) 

else: 

cutoffs = data['cutoffs'] 

return CutoffMaximumBody(cutoffs, len(cutoffs) + 1) 

 

def __str__(self): 

return str(self._cutoff_matrix) 

 

 

class CutoffMaximumBody(Cutoffs): 

""" Specify cutoff-list plus maximum body 

 

Usefull when creating e.g. 6th order expansions but with only 3-body 

interactions. 

 

Parameters 

---------- 

cutoff_list : list 

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

max_nbody : int 

No clusters containing more than max_nbody atoms will be generated 

""" 

 

def __init__(self, cutoff_list, max_nbody): 

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

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

cutoff_matrix[:, order - 2] = cutoff 

super().__init__(cutoff_matrix) 

 

 

def is_cutoff_allowed(atoms, cutoff): 

""" Checks if atoms is compatible with cutoff 

 

Parameters 

---------- 

atoms : ASE Atoms object 

structure used for checking compatibility with cutoff 

cutoff : float 

cutoff to be tested 

 

Returns 

------- 

bool 

True if cutoff compatible with atoms object, else False 

""" 

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

self_interaction=False, bothways=True) 

nbrlist.update(atoms) 

 

for i in range(len(atoms)): 

neighbors, _ = nbrlist.get_neighbors(i) 

if i in neighbors: 

return False 

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

return False 

return True 

 

 

def estimate_maximum_cutoff(atoms, max_iter=11): 

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

 

Parameters 

---------- 

atoms : ASE Atoms object 

structure used for checking compatibility with cutoff 

max_iter : int 

number of iterations in binary search 

""" 

 

# First upper boundary of cutoff 

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

 

# generate all possible offsets given upper_cutoff 

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

self_interaction=False, bothways=True) 

nbrlist.update(atoms) 

all_offsets = [] 

for i in range(len(atoms)): 

_, offsets = nbrlist.get_neighbors(i) 

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

 

# find lower boundary and new upper boundary 

unique_offsets = set(all_offsets) 

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

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

for offset in unique_offsets) 

lower = upper / 2.0 

 

# run binary search between the upper and lower bounds 

for _ in range(max_iter): 

cutoff = (upper + lower) / 2 

if is_cutoff_allowed(atoms, cutoff): 

lower = cutoff 

else: 

upper = cutoff 

return lower