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

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 

57 """ 

58 

59 # Store atoms and calculator 

60 self.atoms = atoms 

61 self.calc = calc 

62 

63 # Displace all atoms in the unit cell by default 

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

65 self.name = name 

66 self.delta = delta 

67 self.center_refcell = center_refcell 

68 self.supercell = supercell 

69 

70 self.cache = MultiFileJSONCache(self.name) 

71 

72 def define_offset(self): # Reference cell offset 

73 

74 if not self.center_refcell: 

75 # Corner cell 

76 self.offset = 0 

77 else: 

78 # Center cell 

79 N_c = self.supercell 

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

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

82 N_c[2] // 2) 

83 return self.offset 

84 

85 @property # type: ignore 

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

87 def N_c(self): 

88 return self._supercell 

89 

90 @property 

91 def supercell(self): 

92 return self._supercell 

93 

94 @supercell.setter 

95 def supercell(self, supercell): 

96 assert len(supercell) == 3 

97 self._supercell = tuple(supercell) 

98 self.define_offset() 

99 self._lattice_vectors_array = self.compute_lattice_vectors() 

100 

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

102 ' instead of .lattice_vectors()') 

103 def lattice_vectors(self): 

104 return self.compute_lattice_vectors() 

105 

106 def compute_lattice_vectors(self): 

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

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

109 

110 # Lattice vectors relevative to the reference cell 

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

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

113 if self.offset == 0: 

114 R_cN += N_c // 2 

115 R_cN %= N_c 

116 R_cN -= N_c // 2 

117 return R_cN 

118 

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

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

121 

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

123 

124 def set_atoms(self, atoms): 

125 """Set the atoms to vibrate. 

126 

127 Parameters: 

128 

129 atoms: list 

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

131 

132 """ 

133 

134 assert isinstance(atoms, list) 

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

136 

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

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

139 sym_a = self.atoms.get_chemical_symbols() 

140 # List for atomic indices 

141 indices = [] 

142 for type in atoms: 

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

144 if atom == type]) 

145 else: 

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

147 indices = atoms 

148 

149 self.indices = indices 

150 

151 def _eq_disp(self): 

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

153 

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

155 from ase.vibrations.vibrations import Displacement as VDisplacement 

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

157 

158 def run(self): 

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

160 

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

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

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

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

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

166 

167 """ 

168 

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

170 # beginning with the last 

171 atoms_N = self.atoms * self.supercell 

172 

173 # Set calculator if provided 

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

175 atoms_N.calc = self.calc 

176 

177 # Do calculation on equilibrium structure 

178 eq_disp = self._eq_disp() 

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

180 if handle is not None: 

181 output = self.calculate(atoms_N, eq_disp) 

182 handle.save(output) 

183 

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

185 natoms = len(self.atoms) 

186 offset = natoms * self.offset 

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

188 

189 # Loop over all displacements 

190 for a in self.indices: 

191 for i in range(3): 

192 for sign in [-1, 1]: 

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

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

195 if handle is None: 

196 continue 

197 try: 

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

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

200 

201 result = self.calculate(atoms_N, disp) 

202 handle.save(result) 

203 finally: 

204 # Return to initial positions 

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

206 

207 def clean(self): 

208 """Delete generated files.""" 

209 if world.rank != 0: 

210 return 0 

211 

212 name = Path(self.name) 

213 

214 n = 0 

215 if name.is_dir(): 

216 for fname in name.iterdir(): 

217 fname.unlink() 

218 n += 1 

219 name.rmdir() 

220 return n 

221 

222 

223class Phonons(Displacement): 

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

225 

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

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

228 

229 2 nbj nbj 

230 nbj d E F- - F+ 

231 C = ------------ ~ ------------- , 

232 mai dR dR 2 * delta 

233 mai nbj 

234 

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

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

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

238 be symmetric in the three indices mai:: 

239 

240 nbj mai bj ai 

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

242 mai nbj ai n bj n 

243 

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

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

246 right hand-side. 

247 

248 The acoustic sum-rule:: 

249 

250 _ _ 

251 aj \ bj 

252 C (R ) = - ) C (R ) 

253 ai 0 /__ ai m 

254 (m, b) 

255 != 

256 (0, a) 

257 

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

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

260 

261 :: 

262 

263 m = 0 m = 1 m = -2 m = -1 

264 ----------------------------------------------------- 

265 | | | | | 

266 | * b | * | * | * | 

267 | | | | | 

268 | * a | * | * | * | 

269 | | | | | 

270 ----------------------------------------------------- 

271 

272 Example: 

273 

274 >>> from ase.build import bulk 

275 >>> from ase.phonons import Phonons 

276 >>> from gpaw import GPAW, FermiDirac 

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

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

279 h=0.2, 

280 occupations=FermiDirac(0.)) 

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

282 >>> ph.run() 

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

284 

285 """ 

286 

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

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

289 

290 if 'name' not in kwargs: 

291 kwargs['name'] = "phonon" 

292 

293 self.deprecate_refcell(kwargs) 

294 

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

296 

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

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

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

300 

301 # Attributes for born charges and static dielectric tensor 

302 self.Z_avv = None 

303 self.eps_vv = None 

304 

305 @staticmethod 

306 def deprecate_refcell(kwargs: dict): 

307 if 'refcell' in kwargs: 

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

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

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

311 kwargs.pop('refcell') 

312 

313 return kwargs 

314 

315 def __call__(self, atoms_N): 

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

317 return atoms_N.get_forces() 

318 

319 def calculate(self, atoms_N, disp): 

320 forces = self(atoms_N) 

321 return {'forces': forces} 

322 

323 def check_eq_forces(self): 

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

325 

326 eq_disp = self._eq_disp() 

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

328 

329 fmin = feq_av.min() 

330 fmax = feq_av.max() 

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

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

333 

334 return fmin, fmax, i_min, i_max 

335 

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

337 'likely incorrect, see ' 

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

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

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

341 

342 The charge neutrality sum-rule:: 

343 

344 _ _ 

345 \ a 

346 ) Z = 0 

347 /__ ij 

348 a 

349 

350 Parameters: 

351 

352 neutrality: bool 

353 Restore charge neutrality condition on calculated Born effective 

354 charges. 

355 name: str 

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

357 in the JSON cache. 

358 

359 """ 

360 

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

362 # unit cell 

363 Z_avv, eps_vv = self.cache[name] 

364 

365 # Neutrality sum-rule 

366 if neutrality: 

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

368 Z_avv -= Z_mean 

369 

370 self.Z_avv = Z_avv[self.indices] 

371 self.eps_vv = eps_vv 

372 

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

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

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

376 

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

378 

379 Parameters: 

380 

381 method: str 

382 Specify method for evaluating the atomic forces. 

383 symmetrize: int 

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

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

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

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

388 number of iterations that will be carried out. 

389 acoustic: bool 

390 Restore the acoustic sum rule on the force constants. 

391 cutoff: None or float 

392 Zero elements in the dynamical matrix between atoms with an 

393 interatomic distance larger than the cutoff. 

394 born: bool 

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

396 dielelctric tensor from file. 

397 

398 """ 

399 

400 method = method.lower() 

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

402 if cutoff is not None: 

403 cutoff = float(cutoff) 

404 

405 # Read Born effective charges and optical dielectric tensor 

406 if born: 

407 self.read_born_charges(**kwargs) 

408 

409 # Number of atoms 

410 natoms = len(self.indices) 

411 # Number of unit cells 

412 N = np.prod(self.supercell) 

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

414 # of eV / Ang**2 

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

416 

417 # Loop over all atomic displacements and calculate force constants 

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

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

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

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

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

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

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

425 

426 if method == 'frederiksen': 

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

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

429 

430 # Finite difference derivative 

431 C_av = fminus_av - fplus_av 

432 C_av /= 2 * self.delta 

433 

434 # Slice out included atoms 

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

436 index = 3 * i + j 

437 C_xNav[index] = C_Nav 

438 

439 # Make unitcell index the first and reshape 

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

441 

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

443 if cutoff is not None: 

444 self.apply_cutoff(C_N, cutoff) 

445 

446 # Symmetrize force constants 

447 if symmetrize: 

448 for i in range(symmetrize): 

449 # Symmetrize 

450 C_N = self.symmetrize(C_N) 

451 # Restore acoustic sum-rule 

452 if acoustic: 

453 self.acoustic(C_N) 

454 else: 

455 break 

456 

457 # Store force constants and dynamical matrix 

458 self.C_N = C_N 

459 self.D_N = C_N.copy() 

460 

461 # Add mass prefactor 

462 m_a = self.atoms.get_masses() 

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

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

465 for D in self.D_N: 

466 D *= M_inv 

467 

468 def symmetrize(self, C_N): 

469 """Symmetrize force constant matrix.""" 

470 

471 # Number of atoms 

472 natoms = len(self.indices) 

473 # Number of unit cells 

474 N = np.prod(self.supercell) 

475 

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

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

478 

479 # Shift reference cell to center index 

480 if self.offset == 0: 

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

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

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

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

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

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

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

488 if self.offset == 0: 

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

490 

491 # Change to single unit cell index shape 

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

493 

494 return C_N 

495 

496 def acoustic(self, C_N): 

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

498 

499 # Number of atoms 

500 natoms = len(self.indices) 

501 # Copy force constants 

502 C_N_temp = C_N.copy() 

503 

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

505 for C in C_N_temp: 

506 for a in range(natoms): 

507 for a_ in range(natoms): 

508 C_N[self.offset, 

509 3 * a: 3 * a + 3, 

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

511 3 * a_: 3 * a_ + 3] 

512 

513 def apply_cutoff(self, D_N, r_c): 

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

515 

516 Parameters: 

517 

518 D_N: ndarray 

519 Dynamical/force constant matrix. 

520 r_c: float 

521 Cutoff in Angstrom. 

522 

523 """ 

524 

525 # Number of atoms and primitive cells 

526 natoms = len(self.indices) 

527 N = np.prod(self.supercell) 

528 # Lattice vectors 

529 R_cN = self._lattice_vectors_array 

530 # Reshape matrix to individual atomic and cartesian dimensions 

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

532 

533 # Cell vectors 

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

535 # Atomic positions in reference cell 

536 pos_av = self.atoms.get_positions() 

537 

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

539 # larger than the cutoff 

540 for n in range(N): 

541 # Lattice vector to cell 

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

543 # Atomic positions in cell 

544 posn_av = pos_av + R_v 

545 # Loop over atoms and zero elements 

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

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

548 # Atoms where the distance is larger than the cufoff 

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

550 # Zero elements 

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

552 

553 def get_force_constant(self): 

554 """Return matrix of force constants.""" 

555 

556 assert self.C_N is not None 

557 return self.C_N 

558 

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

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

561 if modes: 

562 assert 0 

563 omega_kl, modes = omega_kl 

564 

565 from ase.spectrum.band_structure import BandStructure 

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

567 return bs 

568 

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

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

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

572 for a given momentum vector q. 

573 

574 q_scaled: q vector in scaled coordinates. 

575 

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

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

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

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

580 

581 Result: 

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

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

584 """ 

585 # Evaluate fourier sum 

586 R_cN = self._lattice_vectors_array 

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

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

589 return D_q 

590 

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

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

593 

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

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

596 eigenvalues (squared frequency), the corresponding negative frequency 

597 is returned. 

598 

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

600 respectively. 

601 

602 Parameters: 

603 

604 path_kc: ndarray 

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

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

607 dynamical matrix will be calculated. 

608 modes: bool 

609 Returns both frequencies and modes when True. 

610 born: bool 

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

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

613 contribution to the force constant accounts for the splitting 

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

615 verbose: bool 

616 Print warnings when imaginary frequncies are detected. 

617 

618 """ 

619 

620 assert self.D_N is not None 

621 if born: 

622 assert self.Z_avv is not None 

623 assert self.eps_vv is not None 

624 

625 # Dynamical matrix in real-space 

626 D_N = self.D_N 

627 

628 # Lists for frequencies and modes along path 

629 omega_kl = [] 

630 u_kl = [] 

631 

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

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

634 # Unit cell volume in Bohr^3 

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

636 

637 for q_c in path_kc: 

638 

639 # Add non-analytic part 

640 if born: 

641 # q-vector in cartesian coordinates 

642 q_v = np.dot(reci_vc, q_c) 

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

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

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

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

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

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

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

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

651 self.D_na = D_na 

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

653 

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

655 # 

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

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

658 # D_m += q_xx 

659 

660 # Evaluate fourier sum 

661 D_q = self.compute_dynamical_matrix(q_c, D_N) 

662 

663 if modes: 

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

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

666 # multiply with mass prefactor 

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

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

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

670 else: 

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

672 

673 # Sort eigenvalues in increasing order 

674 omega2_l.sort() 

675 # Use dtype=complex to handle negative eigenvalues 

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

677 

678 # Take care of imaginary frequencies 

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

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

681 

682 if verbose: 

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

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

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

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

687 

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

689 

690 omega_kl.append(omega_l.real) 

691 

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

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

694 omega_kl = s * np.asarray(omega_kl) 

695 

696 if modes: 

697 return omega_kl, np.asarray(u_kl) 

698 

699 return omega_kl 

700 

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

702 from ase.spectrum.dosdata import RawDOSData 

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

704 kpts_kc = monkhorst_pack(kpts) 

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

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

707 return dos 

708 

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

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

711 

712 Parameters: 

713 

714 qpts: tuple 

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

716 npts: int 

717 Number of energy points. 

718 delta: float 

719 Broadening of Lorentzian line-shape in eV. 

720 indices: list 

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

722 atoms will be calculated. 

723 

724 """ 

725 

726 # Monkhorst-Pack grid 

727 kpts_kc = monkhorst_pack(kpts) 

728 N = np.prod(kpts) 

729 # Get frequencies 

730 omega_kl = self.band_structure(kpts_kc) 

731 # Energy axis and dos 

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

733 dos_e = np.zeros_like(omega_e) 

734 

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

736 for omega_l in omega_kl: 

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

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

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

740 

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

742 

743 return omega_e, dos_e 

744 

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

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

747 """Write modes to trajectory file. 

748 

749 Parameters: 

750 

751 q_c: ndarray 

752 q-vector of the modes. 

753 branches: int or list 

754 Branch index of modes. 

755 kT: float 

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

757 displacements in the modes. 

758 born: bool 

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

760 repeat: tuple 

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

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

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

764 nimages: int 

765 Number of images in an oscillation. 

766 center: bool 

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

768 

769 """ 

770 

771 if isinstance(branches, int): 

772 branch_l = [branches] 

773 else: 

774 branch_l = list(branches) 

775 

776 # Calculate modes 

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

778 # Repeat atoms 

779 atoms = self.atoms * repeat 

780 # Center 

781 if center: 

782 atoms.center() 

783 

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

785 pos_Nav = atoms.get_positions() 

786 # Total number of unit cells 

787 N = np.prod(repeat) 

788 

789 # Corresponding lattice vectors R_m 

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

791 # Bloch phase 

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

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

794 

795 for lval in branch_l: 

796 

797 omega = omega_l[0, lval] 

798 u_av = u_l[0, lval] 

799 

800 # Mean displacement of a classical oscillator at temperature T 

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

802 

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

804 # Insert slice with atomic displacements for the included atoms 

805 mode_av[self.indices] = u_av 

806 # Repeat and multiply by Bloch phase factor 

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

808 

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

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

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

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

813 mode_Nav).real) 

814 traj.write(atoms)