Coverage for ase / phonons.py: 83.04%

1"""Module for calculating phonons of periodic systems.""" 

2 

3import warnings 

4from math import pi, sqrt 

5from pathlib import Path 

6 

7import numpy as np 

8import numpy.fft as fft 

9import numpy.linalg as la 

10 

11import ase 

12from ase import Atoms, units 

13from ase.dft import monkhorst_pack 

14from ase.io.trajectory import Trajectory 

15from ase.parallel import world 

16from ase.utils import deprecated 

17from ase.utils.filecache import MultiFileJSONCache 

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

36 self, 

37 atoms: Atoms, 

38 calc=None, 

39 supercell: tuple[int, int, int] = (1, 1, 1), 

40 name: str | None = None, 

41 delta: float = 0.01, 

42 center_refcell: bool = False, 

43 comm=None, 

44 ): 

45 """Init with an instance of :class:`~ase.Atoms` and a calculator. 

46 

47 Parameters 

48 ---------- 

49 atoms: :class:`~ase.Atoms` 

50 The atoms to work on. 

51 calc: Calculator 

52 Calculator for the supercell calculation. 

53 supercell: tuple[int, int, int] 

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

55 the small unit cell in each direction. 

56 name: str 

57 Base name to use for files. 

58 delta: float 

59 Magnitude of displacement in Ang. 

60 center_refcell: bool 

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

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

63 of the supercell is used. 

64 comm: communicator 

65 MPI communicator for the phonon calculation. 

66 Default is to use world. 

67 """ 

68 

69 # Store atoms and calculator 

70 self.atoms = atoms 

71 self.calc = calc 

72 

73 # Displace all atoms in the unit cell by default 

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

75 self.name = name 

76 self.delta = delta 

77 self.center_refcell = center_refcell 

78 self.supercell = supercell 

79 

80 if comm is None: 

81 comm = world 

82 self.comm = comm 

83 

84 self.cache = MultiFileJSONCache(self.name) 

85 

86 def define_offset(self): # Reference cell offset 

87 if not self.center_refcell: 

88 # Corner cell 

89 self.offset = 0 

90 else: 

91 # Center cell 

92 N_c = self.supercell 

93 self.offset = ( 

94 N_c[0] // 2 * (N_c[1] * N_c[2]) 

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

96 + N_c[2] // 2 

97 ) 

98 return self.offset 

99 

100 @property 

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

102 def N_c(self): 

103 return self._supercell 

104 

105 @property 

106 def supercell(self): 

107 return self._supercell 

108 

109 @supercell.setter 

110 def supercell(self, supercell): 

111 assert len(supercell) == 3 

112 self._supercell = tuple(supercell) 

113 self.define_offset() 

114 self._lattice_vectors_array = self.compute_lattice_vectors() 

115 

116 @ase.utils.deprecated( 

117 'Please use phonons.compute_lattice_vectors()' 

118 ' instead of .lattice_vectors()' 

119 ) 

120 def lattice_vectors(self): 

121 return self.compute_lattice_vectors() 

122 

123 def compute_lattice_vectors(self): 

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

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

126 

127 # Lattice vectors relevative to the reference cell 

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

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

130 if self.offset == 0: 

131 R_cN += N_c // 2 

132 R_cN %= N_c 

133 R_cN -= N_c // 2 

134 return R_cN 

135 

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

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

138 

139 raise NotImplementedError('Implement in derived classes!.') 

140 

141 def set_atoms(self, atoms): 

142 """Set the atoms to vibrate. 

143 

144 Parameters 

145 ---------- 

146 atoms: list 

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

148 

149 """ 

150 

151 assert isinstance(atoms, list) 

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

153 

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

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

156 sym_a = self.atoms.get_chemical_symbols() 

157 # List for atomic indices 

158 indices = [] 

159 for type in atoms: 

160 indices.extend( 

161 [a for a, atom in enumerate(sym_a) if atom == type] 

162 ) 

163 else: 

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

165 indices = atoms 

166 

167 self.indices = indices 

168 

169 def _eq_disp(self): 

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

171 

172 def _disp(self, a, i: int, step): 

173 from ase.vibrations.vibrations import Displacement as VDisplacement 

174 

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

176 

177 def run(self): 

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

179 

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

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

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

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

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

185 

186 """ 

187 

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

189 # beginning with the last 

190 atoms_N = self.atoms * self.supercell 

191 

192 # Set calculator if provided 

193 assert self.calc is not None, 'Provide calculator in __init__ method' 

194 atoms_N.calc = self.calc 

195 

196 # Do calculation on equilibrium structure 

197 eq_disp = self._eq_disp() 

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

199 if handle is not None: 

200 output = self.calculate(atoms_N, eq_disp) 

201 handle.save(output) 

202 

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

204 natoms = len(self.atoms) 

205 offset = natoms * self.offset 

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

207 

208 # Loop over all displacements 

209 for a in self.indices: 

210 for i in range(3): 

211 for sign in [-1, 1]: 

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

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

214 if handle is None: 

215 continue 

216 try: 

217 atoms_N.positions[offset + a, i] = ( 

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

219 ) 

220 

221 result = self.calculate(atoms_N, disp) 

222 handle.save(result) 

223 finally: 

224 # Return to initial positions 

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

226 

227 self.comm.barrier() 

228 

229 def clean(self): 

230 """Delete generated files.""" 

231 if self.comm.rank == 0: 

232 nfiles = self._clean() 

233 else: 

234 nfiles = 0 

235 self.comm.barrier() 

236 return nfiles 

237 

238 def _clean(self): 

239 name = Path(self.name) 

240 

241 nfiles = 0 

242 if name.is_dir(): 

243 for fname in name.iterdir(): 

244 fname.unlink() 

245 nfiles += 1 

246 name.rmdir() 

247 return nfiles 

248 

249 

250class Phonons(Displacement): 

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

252 

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

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

255 

256 2 nbj nbj 

257 nbj d E F- - F+ 

258 C = ------------ ~ ------------- , 

259 mai dR dR 2 * delta 

260 mai nbj 

261 

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

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

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

265 be symmetric in the three indices mai:: 

266 

267 nbj mai bj ai 

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

269 mai nbj ai n bj n 

270 

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

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

273 right hand-side. 

274 

275 The acoustic sum-rule:: 

276 

277 _ _ 

278 aj \ bj 

279 C (R ) = - ) C (R ) 

280 ai 0 /__ ai m 

281 (m, b) 

282 != 

283 (0, a) 

284 

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

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

287 

288 :: 

289 

290 m = 0 m = 1 m = -2 m = -1 

291 ----------------------------------------------------- 

292 | | | | | 

293 | * b | * | * | * | 

294 | | | | | 

295 | * a | * | * | * | 

296 | | | | | 

297 ----------------------------------------------------- 

298 

299 Examples 

300 -------- 

301 >>> from ase.build import bulk 

302 >>> from ase.phonons import Phonons 

303 >>> from gpaw import GPAW, FermiDirac 

304 

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

306 >>> calc = GPAW(mode='fd', 

307 ... kpts=(5, 5, 5), 

308 ... h=0.2, 

309 ... occupations=FermiDirac(0.)) 

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

311 >>> ph.run() 

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

313 

314 """ 

315 

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

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

318 

319 if 'name' not in kwargs: 

320 kwargs['name'] = 'phonon' 

321 

322 self.deprecate_refcell(kwargs) 

323 

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

325 

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

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

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

329 

330 # Attributes for born charges and static dielectric tensor 

331 self.Z_avv = None 

332 self.eps_vv = None 

333 

334 @staticmethod 

335 def deprecate_refcell(kwargs: dict): 

336 if 'refcell' in kwargs: 

337 warnings.warn( 

338 'Keyword refcell of Phonons is deprecated.' 

339 'Please use center_refcell (bool)', 

340 FutureWarning, 

341 ) 

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

343 kwargs.pop('refcell') 

344 

345 return kwargs 

346 

347 def __call__(self, atoms_N: Atoms): 

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

349 return atoms_N.get_forces() 

350 

351 def calculate(self, atoms_N: Atoms, disp) -> dict[str, np.ndarray]: 

352 forces = self(atoms_N) 

353 return {'forces': forces} 

354 

355 def check_eq_forces(self): 

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

357 

358 eq_disp = self._eq_disp() 

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

360 

361 fmin = feq_av.min() 

362 fmax = feq_av.max() 

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

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

365 

366 return fmin, fmax, i_min, i_max 

367 

368 @deprecated( 

369 'Current implementation of non-analytical correction is ' 

370 'likely incorrect, see ' 

371 'https://gitlab.com/ase/ase/-/issues/941' 

372 ) 

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

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

375 

376 The charge neutrality sum-rule:: 

377 

378 _ _ 

379 \ a 

380 ) Z = 0 

381 /__ ij 

382 a 

383 

384 Parameters 

385 ---------- 

386 

387 neutrality: bool 

388 Restore charge neutrality condition on calculated Born effective 

389 charges. 

390 name: str 

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

392 in the JSON cache. 

393 

394 .. deprecated:: 3.22.1 

395 Current implementation of non-analytical correction is likely 

396 incorrect, see :issue:`941` 

397 """ 

398 

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

400 # unit cell 

401 Z_avv, eps_vv = self.cache[name] 

402 

403 # Neutrality sum-rule 

404 if neutrality: 

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

406 Z_avv -= Z_mean 

407 

408 self.Z_avv = Z_avv[self.indices] 

409 self.eps_vv = eps_vv 

410 

411 def read( 

412 self, 

413 method: str = 'Frederiksen', 

414 symmetrize: int = 3, 

415 acoustic: bool = True, 

416 cutoff: float | None = None, 

417 born: bool = False, 

418 **kwargs, 

419 ): 

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

421 

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

423 

424 Parameters 

425 ---------- 

426 method: str 

427 Specify method for evaluating the atomic forces. 

428 symmetrize: int 

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

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

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

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

433 number of iterations that will be carried out. 

434 acoustic: bool 

435 Restore the acoustic sum rule on the force constants. 

436 cutoff: None or float 

437 Zero elements in the dynamical matrix between atoms with an 

438 interatomic distance larger than the cutoff. 

439 born: bool 

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

441 dielelctric tensor from file. 

442 

443 """ 

444 

445 method = method.lower() 

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

447 if cutoff is not None: 

448 cutoff = float(cutoff) 

449 

450 # Read Born effective charges and optical dielectric tensor 

451 if born: 

452 self.read_born_charges(**kwargs) 

453 

454 # Number of atoms 

455 natoms = len(self.indices) 

456 # Number of unit cells 

457 N = np.prod(self.supercell) 

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

459 # of eV / Ang**2 

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

461 

462 # Loop over all atomic displacements and calculate force constants 

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

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

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

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

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

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

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

470 

471 if method == 'frederiksen': 

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

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

474 

475 # Finite difference derivative 

476 C_av = fminus_av - fplus_av 

477 C_av /= 2 * self.delta 

478 

479 # Slice out included atoms 

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

481 index = 3 * i + j 

482 C_xNav[index] = C_Nav 

483 

484 # Make unitcell index the first and reshape 

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

486 

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

488 if cutoff is not None: 

489 self.apply_cutoff(C_N, cutoff) 

490 

491 # Symmetrize force constants 

492 if symmetrize: 

493 for _ in range(symmetrize): 

494 # Symmetrize 

495 C_N = self.symmetrize(C_N) 

496 # Restore acoustic sum-rule 

497 if acoustic: 

498 self.acoustic(C_N) 

499 else: 

500 break 

501 

502 # Store force constants and dynamical matrix 

503 self.C_N = C_N 

504 self.D_N = C_N.copy() 

505 

506 # Add mass prefactor 

507 m_a = self.atoms.get_masses() 

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

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

510 for D in self.D_N: 

511 D *= M_inv 

512 

513 def symmetrize(self, C_N: np.ndarray) -> np.ndarray: 

514 """Symmetrize force constant matrix.""" 

515 

516 # Number of atoms 

517 natoms = len(self.indices) 

518 # Number of unit cells 

519 N = np.prod(self.supercell) 

520 

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

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

523 

524 # Shift reference cell to center index 

525 if self.offset == 0: 

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

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

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

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

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

531 C_lmn[i:, j:, k:] += ( 

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

533 ) 

534 if self.offset == 0: 

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

536 

537 # Change to single unit cell index shape 

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

539 

540 return C_N 

541 

542 def acoustic(self, C_N: np.ndarray) -> None: 

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

544 

545 # Number of atoms 

546 natoms = len(self.indices) 

547 # Copy force constants 

548 C_N_temp = C_N.copy() 

549 

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

551 for C in C_N_temp: 

552 for a in range(natoms): 

553 for a_ in range(natoms): 

554 C_N[self.offset, 3 * a : 3 * a + 3, 3 * a : 3 * a + 3] -= C[ 

555 3 * a : 3 * a + 3, 3 * a_ : 3 * a_ + 3 

556 ] 

557 

558 def apply_cutoff(self, D_N: np.ndarray, r_c: float) -> None: 

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

560 

561 Parameters 

562 ---------- 

563 D_N: ndarray 

564 Dynamical/force constant matrix. 

565 r_c: float 

566 Cutoff in Angstrom. 

567 

568 """ 

569 

570 # Number of atoms and primitive cells 

571 natoms = len(self.indices) 

572 N = np.prod(self.supercell) 

573 # Lattice vectors 

574 R_cN = self._lattice_vectors_array 

575 # Reshape matrix to individual atomic and cartesian dimensions 

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

577 

578 # Cell vectors 

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

580 # Atomic positions in reference cell 

581 pos_av = self.atoms.get_positions() 

582 

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

584 # larger than the cutoff 

585 for n in range(N): 

586 # Lattice vector to cell 

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

588 # Atomic positions in cell 

589 posn_av = pos_av + R_v 

590 # Loop over atoms and zero elements 

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

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

593 # Atoms where the distance is larger than the cufoff 

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

595 # Zero elements 

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

597 

598 def get_force_constant(self) -> np.ndarray: 

599 """Return matrix of force constants.""" 

600 

601 assert self.C_N is not None 

602 return self.C_N 

603 

604 def get_band_structure( 

605 self, 

606 path, 

607 modes: bool = False, 

608 born: bool = False, 

609 verbose: bool = True, 

610 ): 

611 r"""Calculate and return the phonon band structure. 

612 

613 This method computes the phonon band structure for a given path 

614 in reciprocal space. It is a wrapper around the internal 

615 :meth:`~Phonons.band_structure` method of the :class:`Phonons` class. 

616 The method can optionally calculate and return phonon modes. 

617 

618 Frequencies and modes are in units of eV and 

619 :math:`1/\sqrt{\mathrm{amu}}`, respectively. 

620 

621 Parameters 

622 ---------- 

623 path : BandPath object 

624 The BandPath object defining the path in the reciprocal 

625 space over which the phonon band structure is calculated. 

626 modes : bool, optional 

627 If True, phonon modes will also be calculated and returned. 

628 Defaults to False. 

629 born : bool, optional 

630 If True, includes the effect of Born effective charges in 

631 the phonon calculations. 

632 Defaults to False. 

633 verbose : bool, optional 

634 If True, enables verbose output during the calculation. 

635 Defaults to True. 

636 

637 Returns 

638 ------- 

639 BandStructure or tuple of (BandStructure, ndarray) 

640 If ``modes`` is False, returns a ``BandStructure`` object 

641 containing the phonon band structure. If ``modes`` is True, 

642 returns a tuple, where the first element is the 

643 ``BandStructure`` object and the second element is an ndarray 

644 of phonon modes. 

645 

646 If modes are returned, the array is of shape 

647 (k-point, bands, atoms, 3) and the norm-squared of the mode 

648 is `1 / m_{eff}`, where `m_{eff}` is the effective mass of the 

649 mode. 

650 

651 Examples 

652 -------- 

653 >>> from ase.dft.kpoints import BandPath 

654 >>> path = BandPath(...) # Define the band path 

655 >>> phonons = Phonons(...) 

656 >>> bs, modes = phonons.get_band_structure(path, modes=True) 

657 """ 

658 result = self.band_structure( 

659 path.kpts, modes=modes, born=born, verbose=verbose 

660 ) 

661 if modes: 

662 omega_kl, omega_modes = result 

663 else: 

664 omega_kl = result 

665 

666 from ase.spectrum.band_structure import BandStructure 

667 

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

669 

670 # Return based on the modes flag 

671 return (bs, omega_modes) if modes else bs 

672 

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

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

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

676 for a given momentum vector q. 

677 

678 Parameters 

679 ---------- 

680 q_scaled : np.ndarray 

681 q vector in scaled coordinates. 

682 D_N : np.ndarray 

683 Dynamical matrix in real-space. 

684 It is necessary, at least currently, to provide this matrix 

685 explicitly (rather than use self.D_N) because this matrix is 

686 modified by the Born charges contributions and these modifications 

687 are momentum (q) dependent. 

688 

689 Returns 

690 ------- 

691 D_q : np.ndarray 

692 2D complex-valued array D(q) with shape=(3 * natoms, 3 * natoms). 

693 

694 """ 

695 # Evaluate fourier sum 

696 R_cN = self._lattice_vectors_array 

697 phase_N = np.exp(-2.0j * pi * np.dot(q_scaled, R_cN)) 

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

699 return D_q 

700 

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

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

703 

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

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

706 eigenvalues (squared frequency), the corresponding negative frequency 

707 is returned. 

708 

709 Frequencies and modes are in units of eV and 1/sqrt(amu), 

710 respectively. 

711 

712 Parameters 

713 ---------- 

714 path_kc: ndarray 

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

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

717 dynamical matrix will be calculated. 

718 modes: bool 

719 Returns both frequencies and modes when True. 

720 born: bool 

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

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

723 contribution to the force constant accounts for the splitting 

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

725 verbose: bool 

726 Print warnings when imaginary frequncies are detected. 

727 

728 Returns 

729 ------- 

730 omega_kl : np.ndarray 

731 Energies. 

732 eigenvectors : np.ndarray 

733 Eigenvectors (only when ``modes`` is ``True``). 

734 

735 """ 

736 

737 assert self.D_N is not None 

738 if born: 

739 assert self.Z_avv is not None 

740 assert self.eps_vv is not None 

741 

742 # Dynamical matrix in real-space 

743 D_N = self.D_N 

744 

745 # Lists for frequencies and modes along path 

746 omega_kl = [] 

747 u_kl = [] 

748 

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

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

751 # Unit cell volume in Bohr^3 

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

753 

754 for q_c in path_kc: 

755 # Add non-analytic part 

756 if born: 

757 # q-vector in cartesian coordinates 

758 q_v = np.dot(reci_vc, q_c) 

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

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

761 C_na = ( 

762 4 

763 * pi 

764 * np.outer(qdotZ_av, qdotZ_av) 

765 / np.dot(q_v, np.dot(self.eps_vv, q_v)) 

766 / vol 

767 ) 

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

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

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

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

772 self.D_na = D_na 

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

774 

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

776 # 

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

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

779 # D_m += q_xx 

780 

781 # Evaluate fourier sum 

782 D_q = self.compute_dynamical_matrix(q_c, D_N) 

783 

784 if modes: 

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

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

787 # multiply with mass prefactor. This gives the eigenmode 

788 # (which is now not normalized!) with units 1/sqrt(amu). 

789 u_lx = ( 

790 self.m_inv_x[:, np.newaxis] * u_xl[:, omega2_l.argsort()] 

791 ).T.copy() 

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

793 else: 

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

795 

796 # Sort eigenvalues in increasing order 

797 omega2_l.sort() 

798 # Use dtype=complex to handle negative eigenvalues 

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

800 

801 # Take care of imaginary frequencies 

802 if not np.all(omega2_l >= 0.0): 

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

804 

805 if verbose: 

806 print( 

807 'WARNING, %i imaginary frequencies at ' 

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

809 % ( 

810 len(indices), 

811 q_c[0], 

812 q_c[1], 

813 q_c[2], 

814 omega_l[indices][0].imag, 

815 ) 

816 ) 

817 

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

819 

820 omega_kl.append(omega_l.real) 

821 

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

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

824 omega_kl = s * np.asarray(omega_kl) 

825 

826 if modes: 

827 return omega_kl, np.asarray(u_kl) 

828 

829 return omega_kl 

830 

831 def get_dos( 

832 self, 

833 kpts: tuple[int, int, int] = (10, 10, 10), 

834 indices: list | None = None, 

835 verbose: bool = True, 

836 ): 

837 """Return a phonon density of states. 

838 

839 Parameters 

840 ---------- 

841 kpts: tuple 

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

843 indices: list 

844 If indices is not None, the amplitude-weighted atomic-partial 

845 DOS for the specified atoms will be calculated. 

846 verbose: bool 

847 Print warnings when imaginary frequncies are detected. 

848 

849 Returns 

850 ------- 

851 RawDOSData 

852 Density of states. 

853 

854 """ 

855 from ase.spectrum.dosdata import RawDOSData 

856 

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

858 kpts_kc = monkhorst_pack(kpts) 

859 if indices is None: 

860 # Return the total DOS 

861 omega_w = self.band_structure(kpts_kc, verbose=verbose) 

862 assert omega_w.ndim == 2 

863 n_kpt = omega_w.shape[0] 

864 omega_w = omega_w.ravel() 

865 dos = RawDOSData(omega_w, np.ones_like(omega_w) / n_kpt) 

866 else: 

867 # Return a partial DOS 

868 omegas, amplitudes = self.band_structure( 

869 kpts_kc, modes=True, verbose=verbose 

870 ) 

871 # omegas.shape = (k-points, bands) 

872 # amplitudes.shape = (k-points, bands, atoms, 3) 

873 ampl_sq = (np.abs(amplitudes) ** 2).sum(axis=3) 

874 assert ampl_sq.ndim == 3 

875 assert ampl_sq.shape == omegas.shape + (len(self.indices),) 

876 weights = ampl_sq[:, :, indices].sum(axis=2) / ampl_sq.sum(axis=2) 

877 dos = RawDOSData(omegas.ravel(), weights.ravel() / omegas.shape[0]) 

878 return dos 

879 

880 @deprecated('Please use Phonons.get_dos() instead of Phonons.dos().') 

881 def dos( 

882 self, 

883 kpts: tuple[int, int, int] = (10, 10, 10), 

884 npts: int = 1000, 

885 delta: float = 1e-3, 

886 ): 

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

888 

889 Parameters 

890 ---------- 

891 kpts: tuple 

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

893 npts: int 

894 Number of energy points. 

895 delta: float 

896 Broadening of Lorentzian line-shape in eV. 

897 

898 Returns 

899 ------- 

900 Tuple of (frequencies, dos). The frequencies are in units of eV. 

901 

902 .. deprecated:: 3.23.1 

903 Please use the ``.get_dos()`` method instead, it returns a proper 

904 RawDOSData object. 

905 """ 

906 

907 # Monkhorst-Pack grid 

908 kpts_kc = monkhorst_pack(kpts) 

909 N = np.prod(kpts) 

910 # Get frequencies 

911 omega_kl = self.band_structure(kpts_kc) 

912 # Energy axis and dos 

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