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

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

""" 

This module contains various support/utility functions. 

""" 

 

from typing import List 

import numpy as np 

 

from ase import Atoms 

from ase.calculators.calculator import Calculator 

from ase.geometry import find_mic 

from ase.geometry import get_distances 

from hiphive import ClusterSpace, ForceConstants 

from hiphive.force_constant_model import ForceConstantModel 

from .io.logging import logger 

 

 

logger = logger.getChild('utilities') 

 

 

def get_displacements(atoms: Atoms, atoms_ideal: Atoms) -> np.ndarray: 

"""Returns the the smallest possible displacements between a 

displaced configuration relative to an ideal (reference) 

configuration. 

 

Notes 

----- 

* uses :func:`ase.geometry.find_mic` 

* assumes periodic boundary conditions in all directions 

 

Parameters 

---------- 

atoms 

configuration with displaced atoms 

atoms_ideal 

ideal configuration relative to which displacements are computed 

""" 

37 ↛ 38line 37 didn't jump to line 38, because the condition on line 37 was never true if not np.array_equal(atoms.numbers, atoms_ideal.numbers): 

raise ValueError('Atomic numbers do not match.') 

 

displacements = [] 

for pos, ideal_pos in zip(atoms.positions, atoms_ideal.positions): 

v_ij = np.array([pos - ideal_pos]) 

displacements.append(find_mic(v_ij, atoms_ideal.cell, pbc=True)[0][0]) 

return np.array(displacements) 

 

 

def prepare_structures(structures: List[Atoms], atoms_ideal: Atoms, 

calc: Calculator) -> None: 

"""Prepares a set of structures in the format suitable for a 

:class:`StructureContainer <hiphive.StructureContainer>`. This 

includes retrieving the forces using the provided calculator, 

resetting the positions to the ideal configuration, and adding 

arrays with forces and displacements. 

 

Notes 

----- 

Changes the structures in place. 

 

Parameters 

---------- 

structures 

list of input configurations 

atoms_ideal 

reference configuration relative to which displacements 

are computed 

calc 

ASE calculator used for computing forces 

""" 

for atoms in structures: 

atoms.set_calculator(calc) 

forces = atoms.get_forces() 

displacements = get_displacements(atoms, atoms_ideal) 

atoms.positions = atoms_ideal.get_positions() 

atoms.new_array('forces', forces) 

atoms.new_array('displacements', displacements) 

 

 

def find_permutation(atoms: Atoms, atoms_ref: Atoms) -> List[int]: 

""" Returns the best permutation of atoms for mapping one 

configuration onto another. 

 

Parameters 

---------- 

atoms 

configuration to be permuted 

atoms_ref 

configuration onto which to map 

 

Examples 

-------- 

After obtaining the permutation via 

``` 

p = find_permutation(atoms1, atoms2) 

``` 

the reordered structure ``atoms1[p]`` will give the closest match 

to ``atoms2``. 

""" 

assert np.linalg.norm(atoms.cell - atoms_ref.cell) < 1e-6 

dist_matrix = get_distances( 

atoms.positions, atoms_ref.positions, cell=atoms_ref.cell, pbc=True)[1] 

permutation = [] 

for i in range(len(atoms_ref)): 

dist_row = dist_matrix[:, i] 

permutation.append(np.argmin(dist_row)) 

 

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

raise Exception('Duplicates in permutation') 

for i, p in enumerate(permutation): 

if atoms[p].symbol != atoms_ref[i].symbol: 

raise Exception('Matching lattice sites have different occupation') 

return permutation 

 

 

class Shell: 

""" 

Neighbor Shell class 

 

Parameters 

---------- 

types : list or tuple 

atomic types for neighbor shell 

distance : float 

interatomic distance for neighbor shell 

count : int 

number of pairs in the neighbor shell 

""" 

 

def __init__(self, types, distance, count=0): 

self.types = types 

self.distance = distance 

self.count = count 

 

def __str__(self): 

s = '{}-{} distance: {:10.6f} count: {}'.format(*self.types, self.distance, self.count) 

return s 

 

__repr__ = __str__ 

 

 

def get_neighbor_shells(atoms: Atoms, cutoff: float, dist_tol: float = 1e-5) -> List[Shell]: 

""" Returns a list of neighbor shells. 

 

Distances are grouped into shells via the following algorithm: 

 

1. Find smallest atomic distance `d_min` 

 

2. Find all pair distances in the range `d_min + 1 * dist_tol` 

 

3. Construct a shell from these and pop them from distance list 

 

4. Go to 1. 

 

Parameters 

---------- 

atoms 

configuration used for finding shells 

cutoff 

exclude neighbor shells which have a distance larger than this value 

dist_tol 

distance tolerance 

""" 

 

# TODO: Remove this once this feature have been in ASE long enough 

try: 

from ase.neighborlist import neighbor_list 

except ImportError: 

raise ImportError('get_neighbor_shells uses new functionality from ASE' 

', please update ASE in order to use this function.') 

 

# get distances 

ijd = neighbor_list('ijd', atoms, cutoff) 

ijd = list(zip(*ijd)) 

ijd.sort(key=lambda x: x[2]) 

 

# sort into shells 

symbols = atoms.get_chemical_symbols() 

shells = [] 

for i, j, d in ijd: 

types = tuple(sorted([symbols[i], symbols[j]])) 

for shell in shells: 

if abs(d - shell.distance) < dist_tol and types == shell.types: 

shell.count += 1 

break 

else: 

shell = Shell(types, d, 1) 

shells.append(shell) 

shells.sort(key=lambda x: (x.distance, x.types, x.count)) 

 

# warning if two shells are within 2 * tol 

for i, s1 in enumerate(shells): 

for j, s2 in enumerate(shells[i+1:]): 

if s1.types != s2.types: 

continue 

if not s1.distance < s2.distance - 2 * dist_tol: 

logger.warning('Found two shells within 2 * dist_tol') 

 

return shells 

 

 

def extract_parameters(fcs: ForceConstants, cs: ClusterSpace): 

""" Extracts parameters from force constants. 

 

This function can be used to extract parameters to create a 

ForceConstantPotential from a known set of force constants. 

The return values come from NumPy's `lstsq function 

<https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html>`_. 

 

Parameters 

---------- 

fcs 

force constants 

cs 

cluster space 

""" 

fcm = ForceConstantModel(fcs.supercell, cs) 

return np.linalg.lstsq(*fcm.get_fcs_sensing(fcs), rcond=None)[0]