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"""Module for calculating phonons of periodic systems.""" 

2 

3from math import pi, sqrt 

4import warnings 

5from pathlib import Path 

6 

7import numpy as np 

8import numpy.linalg as la 

9import numpy.fft as fft 

10 

11import ase 

12import ase.units as units 

13from ase.parallel import world 

14from ase.dft import monkhorst_pack 

15from ase.io.trajectory import Trajectory 

16from ase.utils.filecache import MultiFileJSONCache 

17from ase.utils import deprecated 

18 

19 

20class Displacement: 

21 """Abstract base class for phonon and el-ph supercell calculations. 

22 

23 Both phonons and the electron-phonon interaction in periodic systems can be 

24 calculated with the so-called finite-displacement method where the 

25 derivatives of the total energy and effective potential are obtained from 

26 finite-difference approximations, i.e. by displacing the atoms. This class 

27 provides the required functionality for carrying out the calculations for 

28 the different displacements in its ``run`` member function. 

29 

30 Derived classes must overwrite the ``__call__`` member function which is 

31 called for each atomic displacement. 

32 

33 """ 

34 

35 def __init__(self, atoms, calc=None, supercell=(1, 1, 1), name=None, 

36 delta=0.01, center_refcell=False, comm=None): 

37 """Init with an instance of class ``Atoms`` and a calculator. 

38 

39 Parameters: 

40 

41 atoms: Atoms object 

42 The atoms to work on. 

43 calc: Calculator 

44 Calculator for the supercell calculation. 

45 supercell: tuple 

46 Size of supercell given by the number of repetitions (l, m, n) of 

47 the small unit cell in each direction. 

48 name: str 

49 Base name to use for files. 

50 delta: float 

51 Magnitude of displacement in Ang. 

52 center_refcell: bool 

53 Reference cell in which the atoms will be displaced. If False, then 

54 corner cell in supercell is used. If True, then cell in the center 

55 of the supercell is used. 

56 comm: communicator 

57 MPI communicator for the phonon calculation. 

58 Default is to use world. 

59 """ 

60 

61 # Store atoms and calculator 

62 self.atoms = atoms 

63 self.calc = calc 

64 

65 # Displace all atoms in the unit cell by default 

66 self.indices = np.arange(len(atoms)) 

67 self.name = name 

68 self.delta = delta 

69 self.center_refcell = center_refcell 

70 self.supercell = supercell 

71 

72 if comm is None: 

73 comm = world 

74 self.comm = comm 

75 

76 self.cache = MultiFileJSONCache(self.name) 

77 

78 def define_offset(self): # Reference cell offset 

79 

80 if not self.center_refcell: 

81 # Corner cell 

82 self.offset = 0 

83 else: 

84 # Center cell 

85 N_c = self.supercell 

86 self.offset = (N_c[0] // 2 * (N_c[1] * N_c[2]) + 

87 N_c[1] // 2 * N_c[2] + 

88 N_c[2] // 2) 

89 return self.offset 

90 

91 @property # type: ignore 

92 @ase.utils.deprecated('Please use phonons.supercell instead of .N_c') 

93 def N_c(self): 

94 return self._supercell 

95 

96 @property 

97 def supercell(self): 

98 return self._supercell 

99 

100 @supercell.setter 

101 def supercell(self, supercell): 

102 assert len(supercell) == 3 

103 self._supercell = tuple(supercell) 

104 self.define_offset() 

105 self._lattice_vectors_array = self.compute_lattice_vectors() 

106 

107 @ase.utils.deprecated('Please use phonons.compute_lattice_vectors()' 

108 ' instead of .lattice_vectors()') 

109 def lattice_vectors(self): 

110 return self.compute_lattice_vectors() 

111 

112 def compute_lattice_vectors(self): 

113 """Return lattice vectors for cells in the supercell.""" 

114 # Lattice vectors -- ordered as illustrated in class docstring 

115 

116 # Lattice vectors relevative to the reference cell 

117 R_cN = np.indices(self.supercell).reshape(3, -1) 

118 N_c = np.array(self.supercell)[:, np.newaxis] 

119 if self.offset == 0: 

120 R_cN += N_c // 2 

121 R_cN %= N_c 

122 R_cN -= N_c // 2 

123 return R_cN 

124 

125 def __call__(self, *args, **kwargs): 

126 """Member function called in the ``run`` function.""" 

127 

128 raise NotImplementedError("Implement in derived classes!.") 

129 

130 def set_atoms(self, atoms): 

131 """Set the atoms to vibrate. 

132 

133 Parameters: 

134 

135 atoms: list 

136 Can be either a list of strings, ints or ... 

137 

138 """ 

139 

140 assert isinstance(atoms, list) 

141 assert len(atoms) <= len(self.atoms) 

142 

143 if isinstance(atoms[0], str): 

144 assert np.all([isinstance(atom, str) for atom in atoms]) 

145 sym_a = self.atoms.get_chemical_symbols() 

146 # List for atomic indices 

147 indices = [] 

148 for type in atoms: 

149 indices.extend([a for a, atom in enumerate(sym_a) 

150 if atom == type]) 

151 else: 

152 assert np.all([isinstance(atom, int) for atom in atoms]) 

153 indices = atoms 

154 

155 self.indices = indices 

156 

157 def _eq_disp(self): 

158 return self._disp(0, 0, 0) 

159 

160 def _disp(self, a, i, step): 

161 from ase.vibrations.vibrations import Displacement as VDisplacement 

162 return VDisplacement(a, i, np.sign(step), abs(step), self) 

163 

164 def run(self): 

165 """Run the calculations for the required displacements. 

166 

167 This will do a calculation for 6 displacements per atom, +-x, +-y, and 

168 +-z. Only those calculations that are not already done will be 

169 started. Be aware that an interrupted calculation may produce an empty 

170 file (ending with .json), which must be deleted before restarting the 

171 job. Otherwise the calculation for that displacement will not be done. 

172 

173 """ 

174 

175 # Atoms in the supercell -- repeated in the lattice vector directions 

176 # beginning with the last 

177 atoms_N = self.atoms * self.supercell 

178 

179 # Set calculator if provided 

180 assert self.calc is not None, "Provide calculator in __init__ method" 

181 atoms_N.calc = self.calc 

182 

183 # Do calculation on equilibrium structure 

184 eq_disp = self._eq_disp() 

185 with self.cache.lock(eq_disp.name) as handle: 

186 if handle is not None: 

187 output = self.calculate(atoms_N, eq_disp) 

188 handle.save(output) 

189 

190 # Positions of atoms to be displaced in the reference cell 

191 natoms = len(self.atoms) 

192 offset = natoms * self.offset 

193 pos = atoms_N.positions[offset: offset + natoms].copy() 

194 

195 # Loop over all displacements 

196 for a in self.indices: 

197 for i in range(3): 

198 for sign in [-1, 1]: 

199 disp = self._disp(a, i, sign) 

200 with self.cache.lock(disp.name) as handle: 

201 if handle is None: 

202 continue 

203 try: 

204 atoms_N.positions[offset + a, i] = \ 

205 pos[a, i] + sign * self.delta 

206 

207 result = self.calculate(atoms_N, disp) 

208 handle.save(result) 

209 finally: 

210 # Return to initial positions 

211 atoms_N.positions[offset + a, i] = pos[a, i] 

212 

213 self.comm.barrier() 

214 

215 def clean(self): 

216 """Delete generated files.""" 

217 if self.comm.rank == 0: 

218 nfiles = self._clean() 

219 else: 

220 nfiles = 0 

221 self.comm.barrier() 

222 return nfiles 

223 

224 def _clean(self): 

225 name = Path(self.name) 

226 

227 nfiles = 0 

228 if name.is_dir(): 

229 for fname in name.iterdir(): 

230 fname.unlink() 

231 nfiles += 1 

232 name.rmdir() 

233 return nfiles 

234 

235 

236class Phonons(Displacement): 

237 r"""Class for calculating phonon modes using the finite displacement method. 

238 

239 The matrix of force constants is calculated from the finite difference 

240 approximation to the first-order derivative of the atomic forces as:: 

241 

242 2 nbj nbj 

243 nbj d E F- - F+ 

244 C = ------------ ~ ------------- , 

245 mai dR dR 2 * delta 

246 mai nbj 

247 

248 where F+/F- denotes the force in direction j on atom nb when atom ma is 

249 displaced in direction +i/-i. The force constants are related by various 

250 symmetry relations. From the definition of the force constants it must 

251 be symmetric in the three indices mai:: 

252 

253 nbj mai bj ai 

254 C = C -> C (R ) = C (-R ) . 

255 mai nbj ai n bj n 

256 

257 As the force constants can only depend on the difference between the m and 

258 n indices, this symmetry is more conveniently expressed as shown on the 

259 right hand-side. 

260 

261 The acoustic sum-rule:: 

262 

263 _ _ 

264 aj \ bj 

265 C (R ) = - ) C (R ) 

266 ai 0 /__ ai m 

267 (m, b) 

268 != 

269 (0, a) 

270 

271 Ordering of the unit cells illustrated here for a 1-dimensional system (in 

272 case ``refcell=None`` in constructor!): 

273 

274 :: 

275 

276 m = 0 m = 1 m = -2 m = -1 

277 ----------------------------------------------------- 

278 | | | | | 

279 | * b | * | * | * | 

280 | | | | | 

281 | * a | * | * | * | 

282 | | | | | 

283 ----------------------------------------------------- 

284 

285 Example: 

286 

287 >>> from ase.build import bulk 

288 >>> from ase.phonons import Phonons 

289 >>> from gpaw import GPAW, FermiDirac 

290 >>> atoms = bulk('Si', 'diamond', a=5.4) 

291 >>> calc = GPAW(kpts=(5, 5, 5), 

292 h=0.2, 

293 occupations=FermiDirac(0.)) 

294 >>> ph = Phonons(atoms, calc, supercell=(5, 5, 5)) 

295 >>> ph.run() 

296 >>> ph.read(method='frederiksen', acoustic=True) 

297 

298 """ 

299 

300 def __init__(self, *args, **kwargs): 

301 """Initialize with base class args and kwargs.""" 

302 

303 if 'name' not in kwargs: 

304 kwargs['name'] = "phonon" 

305 

306 self.deprecate_refcell(kwargs) 

307 

308 Displacement.__init__(self, *args, **kwargs) 

309 

310 # Attributes for force constants and dynamical matrix in real space 

311 self.C_N = None # in units of eV / Ang**2 

312 self.D_N = None # in units of eV / Ang**2 / amu 

313 

314 # Attributes for born charges and static dielectric tensor 

315 self.Z_avv = None 

316 self.eps_vv = None 

317 

318 @staticmethod 

319 def deprecate_refcell(kwargs: dict): 

320 if 'refcell' in kwargs: 

321 warnings.warn('Keyword refcell of Phonons is deprecated.' 

322 'Please use center_refcell (bool)', FutureWarning) 

323 kwargs['center_refcell'] = bool(kwargs['refcell']) 

324 kwargs.pop('refcell') 

325 

326 return kwargs 

327 

328 def __call__(self, atoms_N): 

329 """Calculate forces on atoms in supercell.""" 

330 return atoms_N.get_forces() 

331 

332 def calculate(self, atoms_N, disp): 

333 forces = self(atoms_N) 

334 return {'forces': forces} 

335 

336 def check_eq_forces(self): 

337 """Check maximum size of forces in the equilibrium structure.""" 

338 

339 eq_disp = self._eq_disp() 

340 feq_av = self.cache[eq_disp.name]['forces'] 

341 

342 fmin = feq_av.min() 

343 fmax = feq_av.max() 

344 i_min = np.where(feq_av == fmin) 

345 i_max = np.where(feq_av == fmax) 

346 

347 return fmin, fmax, i_min, i_max 

348 

349 @deprecated('Current implementation of non-analytical correction is ' 

350 'likely incorrect, see ' 

351 'https://gitlab.com/ase/ase/-/issues/941') 

352 def read_born_charges(self, name='born', neutrality=True): 

353 r"""Read Born charges and dieletric tensor from JSON file. 

354 

355 The charge neutrality sum-rule:: 

356 

357 _ _ 

358 \ a 

359 ) Z = 0 

360 /__ ij 

361 a 

362 

363 Parameters: 

364 

365 neutrality: bool 

366 Restore charge neutrality condition on calculated Born effective 

367 charges. 

368 name: str 

369 Key used to identify the file with Born charges for the unit cell 

370 in the JSON cache. 

371 

372 """ 

373 

374 # Load file with Born charges and dielectric tensor for atoms in the 

375 # unit cell 

376 Z_avv, eps_vv = self.cache[name] 

377 

378 # Neutrality sum-rule 

379 if neutrality: 

380 Z_mean = Z_avv.sum(0) / len(Z_avv) 

381 Z_avv -= Z_mean 

382 

383 self.Z_avv = Z_avv[self.indices] 

384 self.eps_vv = eps_vv 

385 

386 def read(self, method='Frederiksen', symmetrize=3, acoustic=True, 

387 cutoff=None, born=False, **kwargs): 

388 """Read forces from json files and calculate force constants. 

389 

390 Extra keyword arguments will be passed to ``read_born_charges``. 

391 

392 Parameters: 

393 

394 method: str 

395 Specify method for evaluating the atomic forces. 

396 symmetrize: int 

397 Symmetrize force constants (see doc string at top) when 

398 ``symmetrize != 0`` (default: 3). Since restoring the acoustic sum 

399 rule breaks the symmetry, the symmetrization must be repeated a few 

400 times until the changes a insignificant. The integer gives the 

401 number of iterations that will be carried out. 

402 acoustic: bool 

403 Restore the acoustic sum rule on the force constants. 

404 cutoff: None or float 

405 Zero elements in the dynamical matrix between atoms with an 

406 interatomic distance larger than the cutoff. 

407 born: bool 

408 Read in Born effective charge tensor and high-frequency static 

409 dielelctric tensor from file. 

410 

411 """ 

412 

413 method = method.lower() 

414 assert method in ['standard', 'frederiksen'] 

415 if cutoff is not None: 

416 cutoff = float(cutoff) 

417 

418 # Read Born effective charges and optical dielectric tensor 

419 if born: 

420 self.read_born_charges(**kwargs) 

421 

422 # Number of atoms 

423 natoms = len(self.indices) 

424 # Number of unit cells 

425 N = np.prod(self.supercell) 

426 # Matrix of force constants as a function of unit cell index in units 

427 # of eV / Ang**2 

428 C_xNav = np.empty((natoms * 3, N, natoms, 3), dtype=float) 

429 

430 # Loop over all atomic displacements and calculate force constants 

431 for i, a in enumerate(self.indices): 

432 for j, v in enumerate('xyz'): 

433 # Atomic forces for a displacement of atom a in direction v 

434 # basename = '%s.%d%s' % (self.name, a, v) 

435 basename = '%d%s' % (a, v) 

436 fminus_av = self.cache[basename + '-']['forces'] 

437 fplus_av = self.cache[basename + '+']['forces'] 

438 

439 if method == 'frederiksen': 

440 fminus_av[a] -= fminus_av.sum(0) 

441 fplus_av[a] -= fplus_av.sum(0) 

442 

443 # Finite difference derivative 

444 C_av = fminus_av - fplus_av 

445 C_av /= 2 * self.delta 

446 

447 # Slice out included atoms 

448 C_Nav = C_av.reshape((N, len(self.atoms), 3))[:, self.indices] 

449 index = 3 * i + j 

450 C_xNav[index] = C_Nav 

451 

452 # Make unitcell index the first and reshape 

453 C_N = C_xNav.swapaxes(0, 1).reshape((N,) + (3 * natoms, 3 * natoms)) 

454 

455 # Cut off before symmetry and acoustic sum rule are imposed 

456 if cutoff is not None: 

457 self.apply_cutoff(C_N, cutoff) 

458 

459 # Symmetrize force constants 

460 if symmetrize: 

461 for i in range(symmetrize): 

462 # Symmetrize 

463 C_N = self.symmetrize(C_N) 

464 # Restore acoustic sum-rule 

465 if acoustic: 

466 self.acoustic(C_N) 

467 else: 

468 break 

469 

470 # Store force constants and dynamical matrix 

471 self.C_N = C_N 

472 self.D_N = C_N.copy() 

473 

474 # Add mass prefactor 

475 m_a = self.atoms.get_masses() 

476 self.m_inv_x = np.repeat(m_a[self.indices]**-0.5, 3) 

477 M_inv = np.outer(self.m_inv_x, self.m_inv_x) 

478 for D in self.D_N: 

479 D *= M_inv 

480 

481 def symmetrize(self, C_N): 

482 """Symmetrize force constant matrix.""" 

483 

484 # Number of atoms 

485 natoms = len(self.indices) 

486 # Number of unit cells 

487 N = np.prod(self.supercell) 

488 

489 # Reshape force constants to (l, m, n) cell indices 

490 C_lmn = C_N.reshape(self.supercell + (3 * natoms, 3 * natoms)) 

491 

492 # Shift reference cell to center index 

493 if self.offset == 0: 

494 C_lmn = fft.fftshift(C_lmn, axes=(0, 1, 2)).copy() 

495 # Make force constants symmetric in indices -- in case of an even 

496 # number of unit cells don't include the first cell 

497 i, j, k = 1 - np.asarray(self.supercell) % 2 

498 C_lmn[i:, j:, k:] *= 0.5 

499 C_lmn[i:, j:, k:] += \ 

500 C_lmn[i:, j:, k:][::-1, ::-1, ::-1].transpose(0, 1, 2, 4, 3).copy() 

501 if self.offset == 0: 

502 C_lmn = fft.ifftshift(C_lmn, axes=(0, 1, 2)).copy() 

503 

504 # Change to single unit cell index shape 

505 C_N = C_lmn.reshape((N, 3 * natoms, 3 * natoms)) 

506 

507 return C_N 

508 

509 def acoustic(self, C_N): 

510 """Restore acoustic sumrule on force constants.""" 

511 

512 # Number of atoms 

513 natoms = len(self.indices) 

514 # Copy force constants 

515 C_N_temp = C_N.copy() 

516 

517 # Correct atomic diagonals of R_m = (0, 0, 0) matrix 

518 for C in C_N_temp: 

519 for a in range(natoms): 

520 for a_ in range(natoms): 

521 C_N[self.offset, 

522 3 * a: 3 * a + 3, 

523 3 * a: 3 * a + 3] -= C[3 * a: 3 * a + 3, 

524 3 * a_: 3 * a_ + 3] 

525 

526 def apply_cutoff(self, D_N, r_c): 

527 """Zero elements for interatomic distances larger than the cutoff. 

528 

529 Parameters: 

530 

531 D_N: ndarray 

532 Dynamical/force constant matrix. 

533 r_c: float 

534 Cutoff in Angstrom. 

535 

536 """ 

537 

538 # Number of atoms and primitive cells 

539 natoms = len(self.indices) 

540 N = np.prod(self.supercell) 

541 # Lattice vectors 

542 R_cN = self._lattice_vectors_array 

543 # Reshape matrix to individual atomic and cartesian dimensions 

544 D_Navav = D_N.reshape((N, natoms, 3, natoms, 3)) 

545 

546 # Cell vectors 

547 cell_vc = self.atoms.cell.transpose() 

548 # Atomic positions in reference cell 

549 pos_av = self.atoms.get_positions() 

550 

551 # Zero elements with a distance to atoms in the reference cell 

552 # larger than the cutoff 

553 for n in range(N): 

554 # Lattice vector to cell 

555 R_v = np.dot(cell_vc, R_cN[:, n]) 

556 # Atomic positions in cell 

557 posn_av = pos_av + R_v 

558 # Loop over atoms and zero elements 

559 for i, a in enumerate(self.indices): 

560 dist_a = np.sqrt(np.sum((pos_av[a] - posn_av)**2, axis=-1)) 

561 # Atoms where the distance is larger than the cufoff 

562 i_a = dist_a > r_c # np.where(dist_a > r_c) 

563 # Zero elements 

564 D_Navav[n, i, :, i_a, :] = 0.0 

565 

566 def get_force_constant(self): 

567 """Return matrix of force constants.""" 

568 

569 assert self.C_N is not None 

570 return self.C_N 

571 

572 def get_band_structure(self, path, modes=False, born=False, verbose=True): 

573 omega_kl = self.band_structure(path.kpts, modes, born, verbose) 

574 if modes: 

575 assert 0 

576 omega_kl, modes = omega_kl 

577 

578 from ase.spectrum.band_structure import BandStructure 

579 bs = BandStructure(path, energies=omega_kl[None]) 

580 return bs 

581 

582 def compute_dynamical_matrix(self, q_scaled: np.ndarray, D_N: np.ndarray): 

583 """ Computation of the dynamical matrix in momentum space D_ab(q). 

584 This is a Fourier transform from real-space dynamical matrix D_N 

585 for a given momentum vector q. 

586 

587 q_scaled: q vector in scaled coordinates. 

588 

589 D_N: the dynamical matrix in real-space. It is necessary, at least 

590 currently, to provide this matrix explicitly (rather than use 

591 self.D_N) because this matrix is modified by the Born charges 

592 contributions and these modifications are momentum (q) dependent. 

593 

594 Result: 

595 D(q): two-dimensional, complex-valued array of 

596 shape=(3 * natoms, 3 * natoms). 

597 """ 

598 # Evaluate fourier sum 

599 R_cN = self._lattice_vectors_array 

600 phase_N = np.exp(-2.j * pi * np.dot(q_scaled, R_cN)) 

601 D_q = np.sum(phase_N[:, np.newaxis, np.newaxis] * D_N, axis=0) 

602 return D_q 

603 

604 def band_structure(self, path_kc, modes=False, born=False, verbose=True): 

605 """Calculate phonon dispersion along a path in the Brillouin zone. 

606 

607 The dynamical matrix at arbitrary q-vectors is obtained by Fourier 

608 transforming the real-space force constants. In case of negative 

609 eigenvalues (squared frequency), the corresponding negative frequency 

610 is returned. 

611 

612 Frequencies and modes are in units of eV and Ang/sqrt(amu), 

613 respectively. 

614 

615 Parameters: 

616 

617 path_kc: ndarray 

618 List of k-point coordinates (in units of the reciprocal lattice 

619 vectors) specifying the path in the Brillouin zone for which the 

620 dynamical matrix will be calculated. 

621 modes: bool 

622 Returns both frequencies and modes when True. 

623 born: bool 

624 Include non-analytic part given by the Born effective charges and 

625 the static part of the high-frequency dielectric tensor. This 

626 contribution to the force constant accounts for the splitting 

627 between the LO and TO branches for q -> 0. 

628 verbose: bool 

629 Print warnings when imaginary frequncies are detected. 

630 

631 """ 

632 

633 assert self.D_N is not None 

634 if born: 

635 assert self.Z_avv is not None 

636 assert self.eps_vv is not None 

637 

638 # Dynamical matrix in real-space 

639 D_N = self.D_N 

640 

641 # Lists for frequencies and modes along path 

642 omega_kl = [] 

643 u_kl = [] 

644 

645 # Reciprocal basis vectors for use in non-analytic contribution 

646 reci_vc = 2 * pi * la.inv(self.atoms.cell) 

647 # Unit cell volume in Bohr^3 

648 vol = abs(la.det(self.atoms.cell)) / units.Bohr**3 

649 

650 for q_c in path_kc: 

651 

652 # Add non-analytic part 

653 if born: 

654 # q-vector in cartesian coordinates 

655 q_v = np.dot(reci_vc, q_c) 

656 # Non-analytic contribution to force constants in atomic units 

657 qdotZ_av = np.dot(q_v, self.Z_avv).ravel() 

658 C_na = (4 * pi * np.outer(qdotZ_av, qdotZ_av) / 

659 np.dot(q_v, np.dot(self.eps_vv, q_v)) / vol) 

660 self.C_na = C_na / units.Bohr**2 * units.Hartree 

661 # Add mass prefactor and convert to eV / (Ang^2 * amu) 

662 M_inv = np.outer(self.m_inv_x, self.m_inv_x) 

663 D_na = C_na * M_inv / units.Bohr**2 * units.Hartree 

664 self.D_na = D_na 

665 D_N = self.D_N + D_na / np.prod(self.supercell) 

666 

667 # if np.prod(self.N_c) == 1: 

668 # 

669 # q_av = np.tile(q_v, len(self.indices)) 

670 # q_xx = np.vstack([q_av]*len(self.indices)*3) 

671 # D_m += q_xx 

672 

673 # Evaluate fourier sum 

674 D_q = self.compute_dynamical_matrix(q_c, D_N) 

675 

676 if modes: 

677 omega2_l, u_xl = la.eigh(D_q, UPLO='U') 

678 # Sort eigenmodes according to eigenvalues (see below) and 

679 # multiply with mass prefactor 

680 u_lx = (self.m_inv_x[:, np.newaxis] * 

681 u_xl[:, omega2_l.argsort()]).T.copy() 

682 u_kl.append(u_lx.reshape((-1, len(self.indices), 3))) 

683 else: 

684 omega2_l = la.eigvalsh(D_q, UPLO='U') 

685 

686 # Sort eigenvalues in increasing order 

687 omega2_l.sort() 

688 # Use dtype=complex to handle negative eigenvalues 

689 omega_l = np.sqrt(omega2_l.astype(complex)) 

690 

691 # Take care of imaginary frequencies 

692 if not np.all(omega2_l >= 0.): 

693 indices = np.where(omega2_l < 0)[0] 

694 

695 if verbose: 

696 print('WARNING, %i imaginary frequencies at ' 

697 'q = (% 5.2f, % 5.2f, % 5.2f) ; (omega_q =% 5.3e*i)' 

698 % (len(indices), q_c[0], q_c[1], q_c[2], 

699 omega_l[indices][0].imag)) 

700 

701 omega_l[indices] = -1 * np.sqrt(np.abs(omega2_l[indices].real)) 

702 

703 omega_kl.append(omega_l.real) 

704 

705 # Conversion factor: sqrt(eV / Ang^2 / amu) -> eV 

706 s = units._hbar * 1e10 / sqrt(units._e * units._amu) 

707 omega_kl = s * np.asarray(omega_kl) 

708 

709 if modes: 

710 return omega_kl, np.asarray(u_kl) 

711 

712 return omega_kl 

713 

714 def get_dos(self, kpts=(10, 10, 10), npts=1000, delta=1e-3, indices=None): 

715 from ase.spectrum.dosdata import RawDOSData 

716 # dos = self.dos(kpts, npts, delta, indices) 

717 kpts_kc = monkhorst_pack(kpts) 

718 omega_w = self.band_structure(kpts_kc).ravel() 

719 dos = RawDOSData(omega_w, np.ones_like(omega_w)) 

720 return dos 

721 

722 def dos(self, kpts=(10, 10, 10), npts=1000, delta=1e-3, indices=None): 

723 """Calculate phonon dos as a function of energy. 

724 

725 Parameters: 

726 

727 qpts: tuple 

728 Shape of Monkhorst-Pack grid for sampling the Brillouin zone. 

729 npts: int 

730 Number of energy points. 

731 delta: float 

732 Broadening of Lorentzian line-shape in eV. 

733 indices: list 

734 If indices is not None, the atomic-partial dos for the specified 

735 atoms will be calculated. 

736 

737 """ 

738 

739 # Monkhorst-Pack grid 

740 kpts_kc = monkhorst_pack(kpts) 

741 N = np.prod(kpts) 

742 # Get frequencies 

743 omega_kl = self.band_structure(kpts_kc) 

744 # Energy axis and dos 

745 omega_e = np.linspace(0., np.amax(omega_kl) + 5e-3, num=npts) 

746 dos_e = np.zeros_like(omega_e) 

747 

748 # Sum up contribution from all q-points and branches 

749 for omega_l in omega_kl: 

750 diff_el = (omega_e[:, np.newaxis] - omega_l[np.newaxis, :])**2 

751 dos_el = 1. / (diff_el + (0.5 * delta)**2) 

752 dos_e += dos_el.sum(axis=1) 

753 

754 dos_e *= 1. / (N * pi) * 0.5 * delta 

755 

756 return omega_e, dos_e 

757 

758 def write_modes(self, q_c, branches=0, kT=units.kB * 300, born=False, 

759 repeat=(1, 1, 1), nimages=30, center=False): 

760 """Write modes to trajectory file. 

761 

762 Parameters: 

763 

764 q_c: ndarray 

765 q-vector of the modes. 

766 branches: int or list 

767 Branch index of modes. 

768 kT: float 

769 Temperature in units of eV. Determines the amplitude of the atomic 

770 displacements in the modes. 

771 born: bool 

772 Include non-analytic contribution to the force constants at q -> 0. 

773 repeat: tuple 

774 Repeat atoms (l, m, n) times in the directions of the lattice 

775 vectors. Displacements of atoms in repeated cells carry a Bloch 

776 phase factor given by the q-vector and the cell lattice vector R_m. 

777 nimages: int 

778 Number of images in an oscillation. 

779 center: bool 

780 Center atoms in unit cell if True (default: False). 

781 

782 """ 

783 

784 if isinstance(branches, int): 

785 branch_l = [branches] 

786 else: 

787 branch_l = list(branches) 

788 

789 # Calculate modes 

790 omega_l, u_l = self.band_structure([q_c], modes=True, born=born) 

791 # Repeat atoms 

792 atoms = self.atoms * repeat 

793 # Center 

794 if center: 

795 atoms.center() 

796 

797 # Here ``Na`` refers to a composite unit cell/atom dimension 

798 pos_Nav = atoms.get_positions() 

799 # Total number of unit cells 

800 N = np.prod(repeat) 

801 

802 # Corresponding lattice vectors R_m 

803 R_cN = np.indices(repeat).reshape(3, -1) 

804 # Bloch phase 

805 phase_N = np.exp(2.j * pi * np.dot(q_c, R_cN)) 

806 phase_Na = phase_N.repeat(len(self.atoms)) 

807 

808 for lval in branch_l: 

809 

810 omega = omega_l[0, lval] 

811 u_av = u_l[0, lval] 

812 

813 # Mean displacement of a classical oscillator at temperature T 

814 u_av *= sqrt(kT) / abs(omega) 

815 

816 mode_av = np.zeros((len(self.atoms), 3), dtype=complex) 

817 # Insert slice with atomic displacements for the included atoms 

818 mode_av[self.indices] = u_av 

819 # Repeat and multiply by Bloch phase factor 

820 mode_Nav = np.vstack(N * [mode_av]) * phase_Na[:, np.newaxis] 

821 

822 with Trajectory('%s.mode.%d.traj' 

823 % (self.name, lval), 'w') as traj: 

824 for x in np.linspace(0, 2 * pi, nimages, endpoint=False): 

825 atoms.set_positions((pos_Nav + np.exp(1.j * x) * 

826 mode_Nav).real) 

827 traj.write(atoms)