Coverage for local_installation/dynasor/modes/tools.py: 62%

1import numbers

2import numba

3from ase.geometry import get_distances

4import numpy as np

5from fractions import Fraction

6from ase.geometry import find_mic

9# Geometry

10def get_displacements(supercell, ideal, check_mic=True):

11 """Gets displacement from a supercell that might have wrapped positions"""

12 u = supercell.positions - ideal.positions

13 return get_wrapped_displacements(u, ideal.cell, check_mic=True)

16def get_wrapped_displacements(u, cell, check_mic=True):

17 """wraps displacements using mic"""

18 if check_mic: 18 ↛ 20line 18 didn't jump to line 20 because the condition on line 18 was always true

19 u, _ = find_mic(u, cell)

20 return u

23def find_permutation(atoms, atoms_ref):

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

25 configuration onto another.

27 Parameters

28 ----------

29 atoms

30 configuration to be permuted

31 atoms_ref

32 configuration onto which to map

34 Examples

35 --------

36 After obtaining the permutation via ``p = find_permutation(atoms1, atoms2)``

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

38 to ``atoms2``.

39 """

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

41 permutation = []

42 for i in range(len(atoms_ref)):

43 dist_row = get_distances(

44 atoms.positions, atoms_ref.positions[i], cell=atoms_ref.cell, pbc=True)[1][:, 0]

45 permutation.append(np.argmin(dist_row))

47 if len(set(permutation)) != len(permutation): 47 ↛ 48line 47 didn't jump to line 48 because the condition on line 47 was never true

48 raise Exception('Duplicates in permutation')

49 for i, p in enumerate(permutation):

50 if atoms[p].symbol != atoms_ref[i].symbol: 50 ↛ 51line 50 didn't jump to line 51 because the condition on line 50 was never true

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

52 return permutation

55# Linalg

56def trace(A):

57 return A[0, 0] + A[1, 1] + A[2, 2]

60def det(A):

61 if len(A) == 2:

62 return A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0]

63 d = 0

64 for i, B in enumerate(A[0]):

65 minor = np.hstack([A[1:, :i], A[1:, i+1:]])

66 d += (-1)**i * B * det(minor)

67 return d

70def inv(A, as_fraction=True):

71 """

72 Inverts A matrix

74 Parameters

75 ----------

76 A

77 array to be inverted

78 as_fraction

79 boolean which determines if inverted matrix is returned as Fraction matrix or real matrix

80 """

82 detx2 = det(A) * 2

84 rest = (trace(A)**2 - trace(A @ A)) * np.diag([1, 1, 1]) - 2 * A * trace(A) + 2 * A @ A

86 if detx2 < 0:

87 detx2 = -detx2

88 rest = -rest

90 A_inv = np.array([Fraction(n, detx2) for n in rest.flat]).reshape(3, 3)

92 assert np.all(A @ A_inv == np.eye(len(A)))

94 # return either real matrix or as fraction matrix

95 A_inv_real = np.reshape([float(x) for x in A_inv.flat], A_inv.shape)

96 assert np.allclose(A_inv_real, np.linalg.inv(A))

98 if as_fraction: 98 ↛ 101line 98 didn't jump to line 101 because the condition on line 98 was always true

99 return A_inv

100 else:

101 return A_inv_real

102

103

104def symmetrize_eigenvectors(eigenvectors, cell=None, max_iter=1000, method='varimax'):

105 """Takes a set of vectors and tries to make them nice

106

107 Parameters

108 ----------

109 eigenvectors

110 If there are n degenerate eigenvectors and m atoms in the basis the

111 array should be (n,m,3)

112 cell

113 If default None nothing is done to the cartesian directions but a cell

114 can be provided so the directions are in scaled coordinates instead.

115 max_iter

116 The number of iterations in the symmetrization procedure

117 method

118 Can be 'varimax' or 'quartimax' or a parameter between 0: qartimax and

119 1: varimax. Depending on the choice one gets e.g. Equamax, Parsimax, etc.

120

121 """

122

123 if cell is None: 123 ↛ 127line 123 didn't jump to line 127 because the condition on line 123 was always true

124 cell = np.eye(3)

125

126 # s = band, i = basis, a = axis

127 eigenvectors = np.einsum('sia,ab->sib', eigenvectors, np.linalg.inv(cell))

128

129 components = eigenvectors.reshape(len(eigenvectors), -1).T

130

131 rotation_matrix = factor_analysis(components, iterations=max_iter, method=method)

132

133 new_eigenvectors = np.dot(components, rotation_matrix).T

134

135 new_eigenvectors = new_eigenvectors.reshape(len(new_eigenvectors), -1, 3)

136

137 new_eigenvectors = np.einsum('sia,ab->sib', new_eigenvectors, cell)

138

139 return new_eigenvectors

140

141

142def factor_analysis(L, iterations=1000, method='varimax'):

143 """Performs factor analysis on L finding rotation matrix R such that L @ R = L' is simple

144

145 In the future consider using the scikit learn methods directly but beware

146 the changes need to accomodate complex numbers

147

148 References:

149 * Sparse Modeling of Landmark and Texture Variability using the Orthomax Criterion

150 Mikkel B. Stegmann, Karl Sj¨ostrand, Rasmus Larsen

151 http://www2.imm.dtu.dk/pubdb/edoc/imm4041.pdf

152

153 * http://www.cs.ucl.ac.uk/staff/d.barber/brml,

154

155 * https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.FactorAnalysis.html # noqa

156

157 * https://stats.stackexchange.com/questions/185216/factor-rotation-methods-varimax-quartimax-oblimin-etc-what-do-the-names # noqa

158

159 """

160

161 nrow, ncol = L.shape

162 R = np.eye(ncol)

163

164 if method == 'varimax': 164 ↛ 166line 164 didn't jump to line 166 because the condition on line 164 was always true

165 gamma = 1

166 elif method == 'quartimax':

167 gamma = 0

168 else:

169 gamma = method

170

171 for _ in range(iterations):

172 LR = np.dot(L, R)

173 grad = LR * (np.abs(LR)**2 - gamma * np.mean(np.abs(LR)**2, axis=0))

174 G = L.T.conj() @ grad

175 u, s, vh = np.linalg.svd(G)

176 R = u @ vh

177

178 return R

179

180

181def group_eigvals(vals, tol=1e-6):

182 assert sorted(vals) == list(vals), vals

183

184 groups = [[0]]

185 for i in range(1, len(vals)):

186 if np.abs(vals[i] - vals[i-1]) < tol:

187 groups[-1].append(i)

188 else:

189 groups.append([i])

190 return groups

191

192

193# misc

194def as_fraction(not_fraction):

195

196 if isinstance(not_fraction, numbers.Number):

197 return Fraction(not_fraction)

198

199 if isinstance(not_fraction, np.ndarray):

200 arr = np.array([Fraction(n) for n in not_fraction.flat])

201 return arr.reshape(not_fraction.shape)

202

203 if isinstance(not_fraction, tuple):

204 arr = (Fraction(n) for n in not_fraction)

205 return arr

206

207 if isinstance(not_fraction, list):

208 arr = [Fraction(n) for n in not_fraction]

209 return arr

210

211

212@numba.njit

213def get_dynamical_matrix(fc, offsets, indices, q):

214

215 n = indices.max() + 1

216 N = len(fc)

217 D = np.zeros(shape=(n, n, 3, 3), dtype=np.complex128)

218 for ia in range(n):

219 for a in range(N):

220 if ia != indices[a]:

221 continue

222 na = offsets[a]

223 for b in range(N):

224 ib = indices[b]

225 nb = offsets[b]

226 dn = (nb - na).astype(np.float64)

227 D[ia, ib] += fc[a, b] * np.exp(2j*np.pi * np.dot(dn, q))

228 break

229 return D

230

231

232# For debug, this is a slower but perhaps more accurate variant, if the fc

233# obeys translational invariance this should give the same result

234# @numba.njit

235# def get_dynamical_matrix_full(fc, offsets, indices, q):

236#

237# n = indices.max() + 1

238# N = len(fc)

239# D = np.zeros(shape=(n, n, 3, 3), dtype=np.complex128)

240# for I in range(N):

241# i = indices[I]

242# m = offsets[I]

243# for J in range(N):

244# j = indices[J]

245# n = offsets[J]

246#

247# off = (n - m).astype(np.float64)

248#

249# phase = np.exp(1j * 2*np.pi * np.dot(off, q))

250#

251# D[i, j] += fc[I, J] * phase

252#

253# D /= (N / D.shape[0])

254#

255# return D

256

257

258def make_table(M):

259 rows = []

260 for r in M:

261 rows.append(''.join(f'{e:<20.2f}' for e in r))

262 return '\n'.join(rows)