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

""" 

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') 

assert len(cluster_groups) == len(fc_list) 

 

fc_dict = {} 

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

for cluster in clusters: 

permutation = np.argsort(cluster) 

cluster = tuple(sorted(cluster)) 

assert cluster not in fc_dict.keys(), 'Duplicate cluster' 

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() 

assert natoms == len(self.atoms), \ 

'Atoms and ForceConstants have different number of atoms' 

 

self._allowed_atoms = set(range(natoms)) 

assert len(self._allowed_atoms - set(cluster_indices)) == 0, \ 

'Not all atoms are part of a force constant' 

 

def get_fc_array(self, order): 

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

 

Parameters 

---------- 

order : int 

force constants for this order will be returned 

 

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') 

 

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: 

133 ↛ 136line 133 didn't jump to line 136, because the condition on line 133 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 

 

142 ↛ 143line 142 didn't jump to line 143, because the condition on line 142 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: 

149 ↛ 153line 149 didn't jump to line 153, because the condition on line 149 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): 

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

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

raise ValueError('Cluster order not in orders') 

160 ↛ 162line 160 didn't jump to line 162, because the condition on line 160 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. 

 

Asserts 

------- 

This will yield assert error if sum rules are not enforced. 

""" 

 

# setup orders 

191 ↛ 194line 191 didn't jump to line 194, because the condition on line 191 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) 

""" 

245 ↛ 246line 245 didn't jump to line 246, because the condition on line 245 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 

251 ↛ 256line 251 didn't jump to line 256, because the condition on line 251 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) 

""" 

270 ↛ 271line 270 didn't jump to line 271, because the condition on line 270 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('') 

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

for c in self.clusters[:3]: 

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

s.append('...') 

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 

assert fc_array.shape == (Natoms,)*order + (3,)*order 

 

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

381 ↛ 380line 381 didn't jump to line 380, because the condition on line 381 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 

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

fc2 = fc2.get_fc_array(order=2) 

 

N = fc2.shape[0] 

assert len(masses) == N, '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)