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

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

import itertools 

import numpy as np 

import spglib as spg 

from . import atoms as atoms_module 

from ..io.logging import logger 

 

logger = logger.getChild('relate_structures') 

 

 

def align_supercell(supercell, prim, symprec=None): 

"""Rotate and translate a supercell configuration such that it is aligned 

with the target primitive cell. 

 

Parameters 

---------- 

sc : ase.Atoms 

supercell configuration 

prim : ase.Atoms 

target primitive configuration 

symprec : float 

precision parameter forwarded to spglib 

 

Returns 

------- 

tuple(ase.Atoms, numpy.ndarray, numpy.ndarray) 

aligned supercell configuration as well as rotation matrix 

(`3x3` array) and translation vector (`3x1` array) that relate 

the input to the aligned supercell configuration. 

""" 

 

# TODO: Make sure the input is what we expect 

 

# find rotation and translation 

R, T = relate_structures(supercell, prim, symprec=symprec) 

 

# Create the aligned system 

aligned_supercell = rotate_atoms(supercell, R) 

aligned_supercell.translate(T) 

aligned_supercell.wrap() 

return aligned_supercell, R, T 

 

 

def relate_structures(reference, target, symprec=1e-5): 

"""Finds rotation and translation operations that align two structures with 

periodic boundary conditions. 

 

The rotation and translation in Cartesian coordinates will map the 

reference structure onto the target 

 

Aligning reference with target can be achieved via the transformations 

>>> R, T = relate_structures(atoms_ref, atoms_target) 

>>> atoms_ref_rotated = rotate_atoms(atoms, R) 

>>> atoms_ref_rotated.translate(T) 

>>> atoms_ref_rotated.wrap() 

>>> atoms_ref_rotated == atoms_target 

 

Parameters 

---------- 

reference : ase.Atoms 

The reference structure to be mapped 

target : ase.Atoms 

The target structure 

 

Returns 

------- 

R : numpy.ndarray 

rotation matrix in Cartesian coordinates (`3x3` array) 

T : numpy.ndarray 

translation vector in Cartesian coordinates 

""" 

 

logger.debug('Reference atoms:') 

_debug_log_atoms(reference) 

 

reference_primitive_cell = get_primitive_cell(reference, symprec=symprec) 

 

logger.debug('Reference primitive cell') 

_debug_log_atoms(reference_primitive_cell) 

 

logger.debug('Target atoms:') 

_debug_log_atoms(target) 

 

target_primitive_cell = get_primitive_cell(target, symprec=symprec) 

 

logger.debug('Target primitive cell') 

_debug_log_atoms(target_primitive_cell) 

 

logger.debug('Sane check that primitive cells can match...') 

_assert_structures_match(reference_primitive_cell, target_primitive_cell) 

 

logger.debug('Finding rotations...') 

rotations = _find_rotations(reference_primitive_cell.cell, 

target_primitive_cell.cell) 

 

logger.debug('Finding transformations...') 

96 ↛ 104line 96 didn't jump to line 104, because the loop on line 96 didn't complete for R in rotations: 

rotated_reference_primitive_cell = \ 

rotate_atoms(reference_primitive_cell, R) 

T = _find_translation(rotated_reference_primitive_cell, 

target_primitive_cell) 

if T is not None: 

break 

else: 

raise Exception(('Found no translation!\n' 

'Reference primitive cell basis:\n' 

'{}\n' 

'Target primitive cell basis:\n' 

'{}') 

.format(reference_primitive_cell.basis, 

target_primitive_cell.basis)) 

 

logger.debug(('Found rotation\n' 

'{}\n' 

'and translation\n' 

'{}') 

.format(R, T)) 

 

return R, T 

 

 

def is_rotation(R, cell_metric=None): 

"""Checks if rotation matrix is orthonormal 

 

A cell metric can be passed of the rotation matrix is in scaled coordinates 

 

Parameters 

---------- 

R : numpy.ndarray 

rotation matrix (`3x3` array) 

cell_metric : numpy.ndarray 

cell metric if the rotation is in scaled coordinates 

""" 

133 ↛ 136line 133 didn't jump to line 136, because the condition on line 133 was never false if not cell_metric: 

cell_metric = np.eye(3) 

 

V = cell_metric 

V_inv = np.linalg.inv(V) 

lhs = np.linalg.multi_dot([V_inv, R.T, V, V.T, R, V_inv.T]) 

 

return np.allclose(lhs, np.eye(3), atol=1e-4) # TODO: tol 

 

 

def _find_rotations(reference_cell_metric, target_cell_metric): 

""" Generates all proper and improper rotations aligning two cell 

metrics. """ 

 

rotations = [] 

V1 = reference_cell_metric 

for perm in itertools.permutations([0, 1, 2]): 

# Make sure the improper rotations are included 

for inv in itertools.product([1, -1], repeat=3): 

V2 = np.diag(inv) @ target_cell_metric[perm, :] 

R = np.linalg.solve(V1, V2).T 

# Make sure the rotation is orthonormal 

if is_rotation(R): 

for R_tmp in rotations: 

157 ↛ 158line 157 didn't jump to line 158, because the condition on line 157 was never true if np.allclose(R, R_tmp): # TODO: tol 

break 

else: 

rotations.append(R) 

 

assert rotations, ('Found no rotations! Reference cell metric:\n' 

'{}\n' 

'Target cell metric:\n' 

'{}').format(reference_cell_metric, target_cell_metric) 

 

logger.debug('Found {} rotations'.format(len(rotations))) 

 

return rotations 

 

 

def _assert_structures_match(ref, prim): 

""" Asserts the structures are compatible with respect to number of atoms, 

atomic numbers and volume. 

 

TODO: tol 

""" 

 

179 ↛ 180line 179 didn't jump to line 180, because the condition on line 179 was never true if len(ref) != len(prim): 

raise ValueError( 

'Number of atoms in reference primitive cell {} does not match ' 

'target primtive {}'.format(len(ref), len(prim))) 

 

184 ↛ 185line 184 didn't jump to line 185, because the condition on line 184 was never true if sorted(ref.numbers) != sorted(prim.numbers): 

raise ValueError('Atomic numbers do not match\nReference: {}\nTarget:' 

' {}\n'.format(ref.numbers, prim.numbers)) 

 

188 ↛ 189line 188 didn't jump to line 189, because the condition on line 188 was never true if not np.isclose(ref.get_volume(), prim.get_volume()): 

raise ValueError( 

'Volume for reference primitive cell {} does not match target ' 

'primitive cell {}\n'.format(ref.get_volume(), prim.get_volume())) 

 

 

def get_primitive_cell(atoms, to_primitive=True, no_idealize=True, 

symprec=1e-5): 

""" Gets primitive cell from spglib. 

 

Parameters 

---------- 

atoms : ase.Atoms 

atomic structure 

to_primitive : bool 

passed to spglib 

no_idealize : bool 

passed to spglib 

""" 

207 ↛ 208line 207 didn't jump to line 208, because the condition on line 207 was never true if not all(atoms.pbc): 

raise ValueError('atoms must have pbc.') 

spg_primitive_cell = spg.standardize_cell(atoms, to_primitive=True, 

no_idealize=True, 

symprec=symprec) 

primitive_cell = atoms_module.Atoms(cell=spg_primitive_cell[0], 

scaled_positions=spg_primitive_cell[1], 

numbers=spg_primitive_cell[2], 

pbc=True) 

return primitive_cell 

 

 

def _debug_log_atoms(atoms): 

logger.debug('cell:\n{}'.format(atoms.cell)) 

logger.debug('spos:\n{}'.format(atoms.get_scaled_positions())) 

logger.debug('pos:\n{}'.format(atoms.positions)) 

logger.debug('numbers:\n{}'.format(atoms.numbers)) 

 

 

def rotate_atoms(atoms, rotation): 

"""Rotates the cell and positions of Atoms and returns a copy 

 

Parameters 

---------- 

atoms : ase.Atoms 

atomic structure 

rotation : numpy.ndarray 

rotation matrix (`3x3` array) 

""" 

cell = np.dot(rotation, atoms.cell.T).T 

positions = np.dot(rotation, atoms.positions.T).T 

return atoms_module.Atoms(cell=cell, positions=positions, 

numbers=atoms.numbers, pbc=atoms.pbc) 

 

 

def _find_translation(reference, target): 

"""Returns the translation between two compatible atomic structures. 

 

The two structures must describe the same structure when infinitely 

repeated but differ by a translation. 

 

Parameters 

---------- 

reference : ase.Atoms 

target : ase.Atoms 

 

Returns 

------- 

numpy.ndarray or None 

translation vector or `None` if structures are incompatible 

""" 

 

atoms = atoms_module.Atoms(cell=target.cell, 

positions=reference.positions, 

numbers=reference.numbers, 

pbc=True) 

atoms.wrap() 

 

atoms_atom_0 = atoms[0] 

for atom in target: 

if atoms_atom_0.symbol != atom.symbol: 

continue 

T = atom.position - atoms_atom_0.position 

atoms_copy = atoms.copy() 

atoms_copy.positions += T 

if are_nonpaired_configurations_equal(atoms_copy, target): 

return T 

return None 

 

 

def are_nonpaired_configurations_equal(atoms1, atoms2): 

""" Checks whether two configurations are identical. To be considered 

equal the structures must have the same cell metric, elemental 

occupation, scaled positions (modulo one), and periodic boundary 

conditions. 

 

Unlike the ``__eq__`` operator of :class:`ase.Atoms` the order of the 

atoms does not matter. 

 

Parameters 

---------- 

atoms1 : ase.Atoms 

atoms2 : ase.Atoms 

 

Returns 

------- 

bool 

True if atoms are equal, False otherwise 

 

TODO: tol 

""" 

n_atoms = len(atoms1) 

299 ↛ 303line 299 didn't jump to line 303, because the condition on line 299 was never true if not (np.allclose(atoms1.cell, atoms2.cell, atol=1e-4) and 

n_atoms == len(atoms2) and 

sorted(atoms1.numbers) == sorted(atoms2.numbers) and 

all(atoms1.pbc == atoms2.pbc)): 

return False 

new_cell = (atoms1.cell + atoms2.cell) / 2 

pos = [a.position for a in atoms1] + [a.position for a in atoms2] 

num = [a.number for a in atoms1] + [a.number for a in atoms2] 

s3 = atoms_module.Atoms(cell=new_cell, positions=pos, numbers=num, 

pbc=True) 

for i in range(n_atoms): 

for j in range(n_atoms, len(s3)): 

d = s3.get_distance(i, j, mic=True) 

if abs(d) < 1e-4: # TODO: tol 

313 ↛ 314line 313 didn't jump to line 314, because the condition on line 313 was never true if s3[i].number != s3[j].number: 

return False 

break 

else: 

return False 

return True