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

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

""" 

This module provides functionality for storing and handling of force constants 

""" 

import tarfile 

import numpy as np 

 

from ase import units 

from itertools import permutations, product 

from scipy.linalg import eig 

from string import ascii_lowercase, ascii_uppercase 

 

from .io.read_write_files import add_items_to_tarfile_pickle, read_items_pickle 

 

 

class ForceConstants: 

"""Container class for force constants. 

 

Either specify `fc_dict` or both `cluster_groups` and `fc_list`. 

 

Parameters 

---------- 

fc_dict : dict 

dict which holds all force constants with clusters as keys and the 

respective force constant as value 

cluster_groups : list 

list of groups of clusters, clusters in the same group should have 

identical force constants. 

fc_list : list 

list of force constants, one force constant for each cluster group 

atoms: ASE Atoms object 

supercell corresponding to the fcs 

 

Attributes 

---------- 

_fc_dict : dict 

dictionary that holds all force constants with clusters as keys and the 

respective force constant as value 

""" 

 

def __init__(self, fc_dict=None, cluster_groups=None, fc_list=None, 

atoms=None): 

 

# setup fc_dict 

if fc_dict is None: 

self._setup_fc_dict(cluster_groups, fc_list) 

self._sparse = True 

self._cluster_groups = cluster_groups 

self._fc_list = fc_list 

else: 

self._fc_dict = fc_dict 

self._sparse = False 

 

# find which orders are in the fc_dict 

orders = np.unique([len(c) for c in self._fc_dict.keys()]) 

self._orders = orders.tolist() 

 

# setup atoms and allowed_atoms 

self._setup_atoms(atoms) 

 

def _setup_fc_dict(self, cluster_groups, fc_list): 

""" Setup fc_dict from cluster_groups and fc_list """ 

62 ↛ 63line 62 didn't jump to line 63, because the condition on line 62 was never true if cluster_groups is None or fc_list is None: 

raise ValueError('Invalid input') 

64 ↛ 65line 64 didn't jump to line 65, because the condition on line 64 was never true if len(cluster_groups) != len(fc_list): 

raise ValueError('cluster_groups and fc_list not of same length') 

 

fc_dict = {} 

for clusters, fc in zip(cluster_groups, fc_list): 

for cluster in clusters: 

permutation = np.argsort(cluster) 

cluster = tuple(sorted(cluster)) 

72 ↛ 73line 72 didn't jump to line 73, because the condition on line 72 was never true if cluster in fc_dict.keys(): 

raise ValueError('cluster_groups contains duplicates') 

fc_dict[cluster] = fc.transpose(permutation) 

 

self._fc_dict = fc_dict 

 

def _setup_atoms(self, atoms): 

""" Setup atoms and sanity check its properties with fc_dict """ 

cluster_indices = [ind for c in self._fc_dict.keys() for ind in c] 

natoms = np.max(cluster_indices) + 1 

if atoms is None: 

self._atoms = None 

else: 

self._atoms = atoms.copy() 

86 ↛ 87line 86 didn't jump to line 87, because the condition on line 86 was never true if natoms != len(self.atoms): 

raise ValueError('Atoms and ForceConstants have different ' 

'number of atoms') 

 

self._allowed_atoms = set(range(natoms)) 

91 ↛ 92line 91 didn't jump to line 92, because the condition on line 91 was never true if len(self._allowed_atoms - set(cluster_indices)) != 0: 

raise ValueError('Not all atoms are part of a force constant') 

 

def get_fc_array(self, order, format='phonopy'): 

""" Return force constants in array format for specified order. 

 

Parameters 

---------- 

order : int 

force constants for this order will be returned 

format : string 

specify which format (shape) the numpy array should have, 

available are phonopy and ase 

Returns 

------- 

NumPy array 

NumPy array with shape `(N,)*order + (3,)*order` where `N` is 

the number of atoms 

""" 

if order not in self.orders: 

raise ValueError('Order not in orders') 

112 ↛ 113line 112 didn't jump to line 113, because the condition on line 112 was never true if format not in ['ase', 'phonopy']: 

raise ValueError('Format must be either ase or phonopy') 

 

if format == 'ase': 

116 ↛ 117line 116 didn't jump to line 117, because the condition on line 116 was never true if order != 2: 

raise ValueError('ASE format works only for order 2') 

fc_array = self.get_fc_array(order, format='phonopy') 

return fc_array.transpose([0, 2, 1, 3]).reshape( 

self.natoms * 3, self.natoms * 3) 

 

fc_array = np.zeros((self.natoms,)*order + (3,)*order) 

for c, fc in self.get_fc_dict(order=order, permutations=True).items(): 

fc_array[c] = fc 

return fc_array 

 

def get_fc_dict(self, order=None, permutations=False): 

""" Return force constant dictionary for one specific order. 

 

Parameters 

---------- 

order : int 

fcs returned for this order 

permutations : bool 

if True returns all permutations of cluster, else only force 

constants for sorted cluster 

 

Returns 

------- 

dict 

dictionary with keys corresponding to clusters and values to the 

respective force constant 

""" 

 

# Return all fcs 

if order is None: 

147 ↛ 150line 147 didn't jump to line 150, because the condition on line 147 was never false if not permutations: 

return self._fc_dict 

else: 

fc_dict = {} 

for c, fc in self._fc_dict.items(): 

for c_perm, fc_perm in get_permuted_fcs(c, fc).items(): 

fc_dict[c_perm] = fc_perm 

return fc_dict 

 

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

raise ValueError('Order not in orders') 

 

# Return fcs for specific order 

fc_order = {} 

for c, fc in self._fc_dict.items(): 

if len(c) == order: 

163 ↛ 167line 163 didn't jump to line 167, because the condition on line 163 was never false if permutations: 

for c_perm, fc_perm in get_permuted_fcs(c, fc).items(): 

fc_order[c_perm] = fc_perm 

else: 

fc_order[c] = fc 

return fc_order 

 

def __getitem__(self, cluster): 

171 ↛ 187line 171 didn't jump to line 187, because the condition on line 171 was never false if isinstance(cluster, (tuple, list)): 

172 ↛ 173line 172 didn't jump to line 173, because the condition on line 172 was never true if len(cluster) not in self.orders: 

raise ValueError('Cluster order not in orders') 

174 ↛ 176line 174 didn't jump to line 176, because the condition on line 174 was never true if len(self._allowed_atoms.intersection(set(cluster))) != \ 

len(set(cluster)): 

raise ValueError('Cluster outside range of atomic indices') 

 

try: 

sorted_cluster = tuple(sorted(cluster)) 

if cluster == sorted_cluster: 

return self._fc_dict[cluster] 

inv_perm = np.argsort(np.argsort(cluster)) 

return self._fc_dict[sorted_cluster].transpose(inv_perm) 

except KeyError: 

return np.zeros((3,)*len(cluster)) 

else: 

raise ValueError('Cluster must be tuple or list') 

 

def assert_acoustic_sum_rules(self, order=None, tol=1e-6): 

""" Asserts that acoustic sum rules are enforced for fcs 

 

Parameters 

---------- 

order : int, Optional 

specifies which order to check, if None all are checked 

tol : float 

numeric tolerance for checking sum rules 

 

Raises 

------ 

This method raise an assertion error if sum rules are not enforced. 

""" 

 

# setup orders 

205 ↛ 208line 205 didn't jump to line 208, because the condition on line 205 was never false if order is None: 

orders = self.orders 

else: 

if order not in self.orders: 

raise ValueError('Order not available in FCS') 

orders = [order] 

 

atomic_indices = range(self.natoms) 

for order in orders: 

assert_msg = 'Acoustic sum rule order {} not enforced for atom' + \ 

' {}'*(order-1) + ' x' 

for ijk in product(atomic_indices, repeat=order-1): 

fc_sum_ijk = np.zeros((3, )*order) 

for l in atomic_indices: 

cluster = ijk + (l, ) 

fc_sum_ijk += self[cluster] 

assert np.linalg.norm(fc_sum_ijk) < tol, \ 

assert_msg.format(order, *ijk) 

 

def __len__(self): 

return len(self._fc_dict) 

 

@property 

def atoms(self): 

""" ASE Atoms object: supercell corresponding to force constants """ 

return self._atoms 

 

@property 

def natoms(self): 

""" int : Number of atoms (maximum index in a cluster +1) """ 

return len(self._allowed_atoms) 

 

@property 

def orders(self): 

""" list : List of the orders for which force constants exists """ 

return self._orders 

 

@property 

def sparse(self): 

""" bool : If True the object was initialized with sparse data """ 

return self._sparse 

 

@property 

def clusters(self): 

return sorted(self._fc_dict.keys(), key=lambda c: (len(c), c)) 

 

def write(self, f): 

"""Write a ForceConstants instance to a file. 

 

Parameters 

---------- 

f : str or file object 

name of input file (str) or stream to write to (file object) 

""" 

259 ↛ 260line 259 didn't jump to line 260, because the condition on line 259 was never true if isinstance(f, str): 

tar_file = tarfile.open(mode='w', name=f) 

else: 

tar_file = tarfile.open(mode='w', fileobj=f) 

 

# add the needed info to recreate the fcs object 

265 ↛ 270line 265 didn't jump to line 270, because the condition on line 265 was never false if self.sparse: 

items = dict(cluster_groups=self._cluster_groups, 

fc_list=self._fc_list, atoms=self.atoms) 

add_items_to_tarfile_pickle(tar_file, items, 'repr') 

else: 

items = dict(fc_dict=self._fc_dict, atoms=self.atoms) 

add_items_to_tarfile_pickle(tar_file, items, 'repr') 

 

tar_file.close() 

 

@staticmethod 

def read(f): 

"""Load a ForceConstants instance from file 

 

Parameters 

---------- 

f : string or file object 

name of input file (string) or stream to load from (file object) 

""" 

284 ↛ 285line 284 didn't jump to line 285, because the condition on line 284 was never true if isinstance(f, str): 

tar_file = tarfile.open(mode='r', name=f) 

else: 

tar_file = tarfile.open(mode='r', fileobj=f) 

 

items = read_items_pickle(tar_file, 'repr') 

return ForceConstants(**items) 

 

def print_cluster(self, cluster: tuple): 

""" 

Prints ForceConstants[cluster] in a nice format. 

 

Parameters 

---------- 

cluster : tuple 

tuple of ints indicating the sites belonging to the cluster 

""" 

print(self._repr_fc(cluster)) 

 

def __str__(self): 

s = [] 

s.append(' ForceConstants '.center(54, '=')) 

s.append('Orders: {}'.format(self.orders)) 

s.append('Atoms: {}'.format(self.atoms)) 

s.append('') 

309 ↛ 316line 309 didn't jump to line 316, because the condition on line 309 was never false if len(self) > 10: 

for c in self.clusters[:3]: 

s.append(self._repr_fc(c)+'\n') 

s.append('...\n') 

for c in self.clusters[-3:]: 

s.append(self._repr_fc(c)+'\n') 

else: 

for c in self.clusters: 

s.append(self._repr_fc(c)+'\n') 

return '\n'.join(s) 

 

def __repr__(self): 

fc_dict_str = '{{{}: {}, ..., {}: {}}}'.format( 

self.clusters[0], self[self.clusters[0]].round(5), 

self.clusters[-1], self[self.clusters[-1]].round(5)) 

return 'ForceConstants(fc_dict={}, atoms={!r})'.format( 

fc_dict_str, self.atoms) 

 

def _repr_fc(self, cluster: tuple, precision=5, suppress=True) -> str: 

""" 

Representation for single cluster and its force constant. 

 

Parameters 

---------- 

cluster : tuple 

tuple of ints indicating the sites belonging to the cluster 

""" 

s = [] 

s.append('Cluster: {}'.format(cluster)) 

for atom_index in cluster: 

s.append(str(self.atoms[atom_index])) 

s.append('force constant:') 

s.append(np.array_str(self[cluster], precision=precision, 

suppress_small=suppress)) 

return '\n'.join(s) 

 

 

def get_permuted_fcs(cluster, fc): 

"""Generate all permutations of cluster and their force constant given a 

single cluster and its force constant. 

 

Parameters 

---------- 

cluster : tuple 

tuple of ints indicating the sites belonging to the cluster 

fc : NumPy array 

force constant for cluster 

 

Returns 

------- 

dict 

dict with keys corresponding to clusters and values to their force 

constant 

""" 

fc_dict = {} 

fc = fc.transpose(np.argsort(cluster)) 

for perm_cluster in set(permutations(cluster)): 

inv_perm = np.argsort(np.argsort(perm_cluster)) 

fc_perm = fc.transpose(inv_perm) 

fc_dict[perm_cluster] = fc_perm 

return fc_dict 

 

 

def fc_array_to_dict(fc_array, cut_tol=1e-8): 

"""Convert a force constant from array to dictionary format. 

 

Parameters 

---------- 

fc_array : array 

force constant array for one order 

cut_tol : float 

tolerance for which force constants to keep in dict 

 

Returns 

------- 

dict 

sparse dictionary with keys corresponding to clusters and values to the 

respective force constant 

""" 

 

fc_dict = {} 

Natoms = fc_array.shape[0] 

order = len(fc_array.shape) // 2 

392 ↛ 393line 392 didn't jump to line 393, because the condition on line 392 was never true if fc_array.shape != (Natoms,)*order + (3,)*order: 

raise ValueError('fc_array shape does not have expected shape') 

 

for cluster in product(*[range(Natoms)]*order): 

396 ↛ 395line 396 didn't jump to line 395, because the condition on line 396 was never false if np.linalg.norm(fc_array[cluster]) > cut_tol: 

if cluster == tuple(sorted(cluster)): 

fc_dict[cluster] = fc_array[cluster] 

return fc_dict 

 

 

def compute_gamma_frequencies(fc2, masses): 

""" Computes gamma frequencies using second order force constants 

 

Assumes fc2 is in eV/A2 

 

Parameters 

---------- 

fc2 : ForceConstants object, NumPy array 

Second order force constants 

masses : list 

list of masses for each atom 

 

Returns 

------- 

gamma_frequencies : NumPy array 

array of gamma frequencies in THz 

""" 

 

# get fc2_array 

421 ↛ 424line 421 didn't jump to line 424, because the condition on line 421 was never false if isinstance(fc2, ForceConstants): 

fc2 = fc2.get_fc_array(order=2) 

 

N = fc2.shape[0] 

425 ↛ 426line 425 didn't jump to line 426, because the condition on line 425 was never true if len(masses) != N: 

raise ValueError('Length of masses not compatible with fc2') 

mass_matrix = np.sqrt(np.outer(masses, masses)) 

 

# divide with mass matrix 

fc2_tmp = np.zeros((N, N, 3, 3)) 

for pair in product(range(N), repeat=2): 

fc2_tmp[pair] = fc2[pair] / mass_matrix[pair] 

 

# reshape into matrix and solve eigenvalues 

fc2_tmp = fc2_tmp.transpose([0, 2, 1, 3]).reshape(N * 3, N * 3) 

eigen_vals, _ = eig(fc2_tmp) 

eigen_vals *= 1e20 / units.J * units.kg # [1/s**2] 

eigen_vals.sort() 

 

# if negative eigenval, set frequency to negative 

gamma_frequencies = [] 

for val in eigen_vals.real: 

if val >= 0: 

gamma_frequencies.append(np.sqrt(val)) 

else: 

gamma_frequencies.append(-np.sqrt(np.abs(val))) 

 

# Sort and convert to THz 

gamma_frequencies = np.array(gamma_frequencies) / 1e12 / (2 * np.pi) 

gamma_frequencies.sort() 

return gamma_frequencies 

 

 

def symbolize_force_constant(fc, tol=1e-10): 

"""Carry out a symbolic symmetrization of a force constant tensor. 

 

Parameters 

---------- 

fc : NumPy array 

force constant tensor 

tol : float 

tolerance used to decide whether two elements are identical 

 

Returns 

------- 

NumPy array 

symbolic representation of force constant matrix 

""" 

 

fc_int = np.round(fc / tol).astype(int) 

 

fc_chars = np.empty(fc_int.shape, dtype=object).flatten() 

all_chars = ascii_lowercase + ascii_uppercase 

lookup_chars = {} 

for i, val in enumerate(fc_int.flatten()): 

if val == 0: 

fc_chars[i] = 0 

elif val in lookup_chars.keys(): 

fc_chars[i] = lookup_chars[val] 

elif -val in lookup_chars.keys(): 

fc_chars[i] = '-{:}'.format(lookup_chars[-val]) 

else: 

lookup_chars[val] = all_chars[len(lookup_chars.keys())] 

fc_chars[i] = lookup_chars[val] 

return fc_chars.reshape(fc_int.shape)