Coverage for /builds/ase/ase/ase/phonons.py: 83.04%

1# fmt: off 

2 

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

4 

5import warnings 

6from math import pi, sqrt 

7from pathlib import Path 

8 

9import numpy as np 

10import numpy.fft as fft 

11import numpy.linalg as la 

12 

13import ase 

14import ase.units as units 

15from ase.dft import monkhorst_pack 

16from ase.io.trajectory import Trajectory 

17from ase.parallel import world 

18from ase.utils import deprecated 

19from ase.utils.filecache import MultiFileJSONCache 

20 

21 

22class Displacement: 

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

24 

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

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

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

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

29 provides the required functionality for carrying out the calculations for 

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

31 

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

33 called for each atomic displacement. 

34 

35 """ 

36 

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

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

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

40 

41 Parameters: 

42 

43 atoms: Atoms object 

44 The atoms to work on. 

45 calc: Calculator 

46 Calculator for the supercell calculation. 

47 supercell: tuple 

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

49 the small unit cell in each direction. 

50 name: str 

51 Base name to use for files. 

52 delta: float 

53 Magnitude of displacement in Ang. 

54 center_refcell: bool 

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

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

57 of the supercell is used. 

58 comm: communicator 

59 MPI communicator for the phonon calculation. 

60 Default is to use world. 

61 """ 

62 

63 # Store atoms and calculator 

64 self.atoms = atoms 

65 self.calc = calc 

66 

67 # Displace all atoms in the unit cell by default 

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

69 self.name = name 

70 self.delta = delta 

71 self.center_refcell = center_refcell 

72 self.supercell = supercell 

73 

74 if comm is None: 

75 comm = world 

76 self.comm = comm 

77 

78 self.cache = MultiFileJSONCache(self.name) 

79 

80 def define_offset(self): # Reference cell offset 

81 

82 if not self.center_refcell: 

83 # Corner cell 

84 self.offset = 0 

85 else: 

86 # Center cell 

87 N_c = self.supercell 

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

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

90 N_c[2] // 2) 

91 return self.offset 

92 

93 @property 

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

95 def N_c(self): 

96 return self._supercell 

97 

98 @property 

99 def supercell(self): 

100 return self._supercell 

101 

102 @supercell.setter 

103 def supercell(self, supercell): 

104 assert len(supercell) == 3 

105 self._supercell = tuple(supercell) 

106 self.define_offset() 

107 self._lattice_vectors_array = self.compute_lattice_vectors() 

108 

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

110 ' instead of .lattice_vectors()') 

111 def lattice_vectors(self): 

112 return self.compute_lattice_vectors() 

113 

114 def compute_lattice_vectors(self): 

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

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

117 

118 # Lattice vectors relevative to the reference cell 

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

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

121 if self.offset == 0: 

122 R_cN += N_c // 2 

123 R_cN %= N_c 

124 R_cN -= N_c // 2 

125 return R_cN 

126 

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

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

129 

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

131 

132 def set_atoms(self, atoms): 

133 """Set the atoms to vibrate. 

134 

135 Parameters: 

136 

137 atoms: list 

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

139 

140 """ 

141 

142 assert isinstance(atoms, list) 

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

144 

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

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

147 sym_a = self.atoms.get_chemical_symbols() 

148 # List for atomic indices 

149 indices = [] 

150 for type in atoms: 

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

152 if atom == type]) 

153 else: 

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

155 indices = atoms 

156 

157 self.indices = indices 

158 

159 def _eq_disp(self): 

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

161 

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

163 from ase.vibrations.vibrations import Displacement as VDisplacement 

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

165 

166 def run(self): 

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

168 

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

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

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

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

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

174 

175 """ 

176 

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

178 # beginning with the last 

179 atoms_N = self.atoms * self.supercell 

180 

181 # Set calculator if provided 

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

183 atoms_N.calc = self.calc 

184 

185 # Do calculation on equilibrium structure 

186 eq_disp = self._eq_disp() 

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

188 if handle is not None: 

189 output = self.calculate(atoms_N, eq_disp) 

190 handle.save(output) 

191 

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

193 natoms = len(self.atoms) 

194 offset = natoms * self.offset 

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

196 

197 # Loop over all displacements 

198 for a in self.indices: 

199 for i in range(3): 

200 for sign in [-1, 1]: 

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

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

203 if handle is None: 

204 continue 

205 try: 

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

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

208 

209 result = self.calculate(atoms_N, disp) 

210 handle.save(result) 

211 finally: 

212 # Return to initial positions 

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

214 

215 self.comm.barrier() 

216 

217 def clean(self): 

218 """Delete generated files.""" 

219 if self.comm.rank == 0: 

220 nfiles = self._clean() 

221 else: 

222 nfiles = 0 

223 self.comm.barrier() 

224 return nfiles 

225 

226 def _clean(self): 

227 name = Path(self.name) 

228 

229 nfiles = 0 

230 if name.is_dir(): 

231 for fname in name.iterdir(): 

232 fname.unlink() 

233 nfiles += 1 

234 name.rmdir() 

235 return nfiles 

236 

237 

238class Phonons(Displacement): 

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

240 

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

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

243 

244 2 nbj nbj 

245 nbj d E F- - F+ 

246 C = ------------ ~ ------------- , 

247 mai dR dR 2 * delta 

248 mai nbj 

249 

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

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

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

253 be symmetric in the three indices mai:: 

254 

255 nbj mai bj ai 

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

257 mai nbj ai n bj n 

258 

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

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

261 right hand-side. 

262 

263 The acoustic sum-rule:: 

264 

265 _ _ 

266 aj \ bj 

267 C (R ) = - ) C (R ) 

268 ai 0 /__ ai m 

269 (m, b) 

270 != 

271 (0, a) 

272 

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

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

275 

276 :: 

277 

278 m = 0 m = 1 m = -2 m = -1 

279 ----------------------------------------------------- 

280 | | | | | 

281 | * b | * | * | * | 

282 | | | | | 

283 | * a | * | * | * | 

284 | | | | | 

285 ----------------------------------------------------- 

286 

287 Example: 

288 

289 >>> from ase.build import bulk 

290 >>> from ase.phonons import Phonons 

291 >>> from gpaw import GPAW, FermiDirac 

292 

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

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

295 ... kpts=(5, 5, 5), 

296 ... h=0.2, 

297 ... occupations=FermiDirac(0.)) 

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

299 >>> ph.run() 

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

301 

302 """ 

303 

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

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

306 

307 if 'name' not in kwargs: 

308 kwargs['name'] = "phonon" 

309 

310 self.deprecate_refcell(kwargs) 

311 

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

313 

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

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

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

317 

318 # Attributes for born charges and static dielectric tensor 

319 self.Z_avv = None 

320 self.eps_vv = None 

321 

322 @staticmethod 

323 def deprecate_refcell(kwargs: dict): 

324 if 'refcell' in kwargs: 

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

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

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

328 kwargs.pop('refcell') 

329 

330 return kwargs 

331 

332 def __call__(self, atoms_N): 

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

334 return atoms_N.get_forces() 

335 

336 def calculate(self, atoms_N, disp): 

337 forces = self(atoms_N) 

338 return {'forces': forces} 

339 

340 def check_eq_forces(self): 

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

342 

343 eq_disp = self._eq_disp() 

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

345 

346 fmin = feq_av.min() 

347 fmax = feq_av.max() 

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

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

350 

351 return fmin, fmax, i_min, i_max 

352 

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

354 'likely incorrect, see ' 

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

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

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

358 

359 The charge neutrality sum-rule:: 

360 

361 _ _ 

362 \ a 

363 ) Z = 0 

364 /__ ij 

365 a 

366 

367 Parameters: 

368 

369 neutrality: bool 

370 Restore charge neutrality condition on calculated Born effective 

371 charges. 

372 name: str 

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

374 in the JSON cache. 

375 

376 .. deprecated:: 3.22.1 

377 Current implementation of non-analytical correction is likely 

378 incorrect, see :issue:`941` 

379 """ 

380 

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

382 # unit cell 

383 Z_avv, eps_vv = self.cache[name] 

384 

385 # Neutrality sum-rule 

386 if neutrality: 

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

388 Z_avv -= Z_mean 

389 

390 self.Z_avv = Z_avv[self.indices] 

391 self.eps_vv = eps_vv 

392 

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

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

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

396 

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

398 

399 Parameters: 

400 

401 method: str 

402 Specify method for evaluating the atomic forces. 

403 symmetrize: int 

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

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

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

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

408 number of iterations that will be carried out. 

409 acoustic: bool 

410 Restore the acoustic sum rule on the force constants. 

411 cutoff: None or float 

412 Zero elements in the dynamical matrix between atoms with an 

413 interatomic distance larger than the cutoff. 

414 born: bool 

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

416 dielelctric tensor from file. 

417 

418 """ 

419 

420 method = method.lower() 

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

422 if cutoff is not None: 

423 cutoff = float(cutoff) 

424 

425 # Read Born effective charges and optical dielectric tensor 

426 if born: 

427 self.read_born_charges(**kwargs) 

428 

429 # Number of atoms 

430 natoms = len(self.indices) 

431 # Number of unit cells 

432 N = np.prod(self.supercell) 

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

434 # of eV / Ang**2 

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

436 

437 # Loop over all atomic displacements and calculate force constants 

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

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

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

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

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

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

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

445 

446 if method == 'frederiksen': 

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

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

449 

450 # Finite difference derivative 

451 C_av = fminus_av - fplus_av 

452 C_av /= 2 * self.delta 

453 

454 # Slice out included atoms 

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

456 index = 3 * i + j 

457 C_xNav[index] = C_Nav 

458 

459 # Make unitcell index the first and reshape 

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

461 

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

463 if cutoff is not None: 

464 self.apply_cutoff(C_N, cutoff) 

465 

466 # Symmetrize force constants 

467 if symmetrize: 

468 for _ in range(symmetrize): 

469 # Symmetrize 

470 C_N = self.symmetrize(C_N) 

471 # Restore acoustic sum-rule 

472 if acoustic: 

473 self.acoustic(C_N) 

474 else: 

475 break 

476 

477 # Store force constants and dynamical matrix 

478 self.C_N = C_N 

479 self.D_N = C_N.copy() 

480 

481 # Add mass prefactor 

482 m_a = self.atoms.get_masses() 

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

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

485 for D in self.D_N: 

486 D *= M_inv 

487 

488 def symmetrize(self, C_N): 

489 """Symmetrize force constant matrix.""" 

490 

491 # Number of atoms 

492 natoms = len(self.indices) 

493 # Number of unit cells 

494 N = np.prod(self.supercell) 

495 

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

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

498 

499 # Shift reference cell to center index 

500 if self.offset == 0: 

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

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

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

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

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

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

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

508 if self.offset == 0: 

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

510 

511 # Change to single unit cell index shape 

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

513 

514 return C_N 

515 

516 def acoustic(self, C_N): 

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

518 

519 # Number of atoms 

520 natoms = len(self.indices) 

521 # Copy force constants 

522 C_N_temp = C_N.copy() 

523 

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

525 for C in C_N_temp: 

526 for a in range(natoms): 

527 for a_ in range(natoms): 

528 C_N[self.offset, 

529 3 * a: 3 * a + 3, 

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

531 3 * a_: 3 * a_ + 3] 

532 

533 def apply_cutoff(self, D_N, r_c): 

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

535 

536 Parameters: 

537 

538 D_N: ndarray 

539 Dynamical/force constant matrix. 

540 r_c: float 

541 Cutoff in Angstrom. 

542 

543 """ 

544 

545 # Number of atoms and primitive cells 

546 natoms = len(self.indices) 

547 N = np.prod(self.supercell) 

548 # Lattice vectors 

549 R_cN = self._lattice_vectors_array 

550 # Reshape matrix to individual atomic and cartesian dimensions 

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

552 

553 # Cell vectors 

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

555 # Atomic positions in reference cell 

556 pos_av = self.atoms.get_positions() 

557 

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

559 # larger than the cutoff 

560 for n in range(N): 

561 # Lattice vector to cell 

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

563 # Atomic positions in cell 

564 posn_av = pos_av + R_v 

565 # Loop over atoms and zero elements 

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

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

568 # Atoms where the distance is larger than the cufoff 

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

570 # Zero elements 

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

572 

573 def get_force_constant(self): 

574 """Return matrix of force constants.""" 

575 

576 assert self.C_N is not None 

577 return self.C_N 

578 

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

580 """Calculate and return the phonon band structure. 

581 

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

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

584 `band_structure` method of the `Phonons` class. The method can 

585 optionally calculate and return phonon modes. 

586 

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

588 respectively. 

589 

590 Parameters: 

591 

592 path : BandPath object 

593 The BandPath object defining the path in the reciprocal 

594 space over which the phonon band structure is calculated. 

595 modes : bool, optional 

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

597 Defaults to False. 

598 born : bool, optional 

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

600 the phonon calculations. 

601 Defaults to False. 

602 verbose : bool, optional 

603 If True, enables verbose output during the calculation. 

604 Defaults to True. 

605 

606 Returns: 

607 

608 BandStructure or tuple of (BandStructure, ndarray) 

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

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

611 returns a tuple, where the first element is the 

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

613 of phonon modes. 

614 

615 If modes are returned, the array is of shape 

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

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

618 mode. 

619 

620 Example: 

621 

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

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

624 >>> phonons = Phonons(...) 

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

626 """ 

627 result = self.band_structure(path.kpts, 

628 modes=modes, 

629 born=born, 

630 verbose=verbose) 

631 if modes: 

632 omega_kl, omega_modes = result 

633 else: 

634 omega_kl = result 

635 

636 from ase.spectrum.band_structure import BandStructure 

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

638 

639 # Return based on the modes flag 

640 return (bs, omega_modes) if modes else bs 

641 

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

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

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

645 for a given momentum vector q. 

646 

647 q_scaled: q vector in scaled coordinates. 

648 

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

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

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

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

653 

654 Result: 

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

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

657 """ 

658 # Evaluate fourier sum 

659 R_cN = self._lattice_vectors_array 

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

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

662 return D_q 

663 

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

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

666 

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

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

669 eigenvalues (squared frequency), the corresponding negative frequency 

670 is returned. 

671 

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

673 respectively. 

674 

675 Parameters: 

676 

677 path_kc: ndarray 

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

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

680 dynamical matrix will be calculated. 

681 modes: bool 

682 Returns both frequencies and modes when True. 

683 born: bool 

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

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

686 contribution to the force constant accounts for the splitting 

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

688 verbose: bool 

689 Print warnings when imaginary frequncies are detected. 

690 

691 Returns: 

692 

693 If modes=False: Array of energies 

694 

695 If modes=True: Tuple of two arrays with energies and modes. 

696 """ 

697 

698 assert self.D_N is not None 

699 if born: 

700 assert self.Z_avv is not None 

701 assert self.eps_vv is not None 

702 

703 # Dynamical matrix in real-space 

704 D_N = self.D_N 

705 

706 # Lists for frequencies and modes along path 

707 omega_kl = [] 

708 u_kl = [] 

709 

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

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

712 # Unit cell volume in Bohr^3 

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

714 

715 for q_c in path_kc: 

716 # Add non-analytic part 

717 if born: 

718 # q-vector in cartesian coordinates 

719 q_v = np.dot(reci_vc, q_c) 

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

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

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

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

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

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

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

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

728 self.D_na = D_na 

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

730 

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

732 # 

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

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

735 # D_m += q_xx 

736 

737 # Evaluate fourier sum 

738 D_q = self.compute_dynamical_matrix(q_c, D_N) 

739 

740 if modes: 

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

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

743 # multiply with mass prefactor. This gives the eigenmode 

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

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

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

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

748 else: 

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

750 

751 # Sort eigenvalues in increasing order 

752 omega2_l.sort() 

753 # Use dtype=complex to handle negative eigenvalues 

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

755 

756 # Take care of imaginary frequencies 

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

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

759 

760 if verbose: 

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

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

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

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

765 

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

767 

768 omega_kl.append(omega_l.real) 

769 

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

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

772 omega_kl = s * np.asarray(omega_kl) 

773 

774 if modes: 

775 return omega_kl, np.asarray(u_kl) 

776 

777 return omega_kl 

778 

779 def get_dos(self, kpts=(10, 10, 10), indices=None, verbose=True): 

780 """Return a phonon density of states. 

781 

782 Parameters: 

783 

784 kpts: tuple 

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

786 indices: list 

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

788 DOS for the specified atoms will be calculated. 

789 verbose: bool 

790 Print warnings when imaginary frequncies are detected. 

791 

792 Returns: 

793 A RawDOSData object containing the density of states. 

794 """ 

795 from ase.spectrum.dosdata import RawDOSData 

796 

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

798 kpts_kc = monkhorst_pack(kpts) 

799 if indices is None: 

800 # Return the total DOS 

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

802 assert omega_w.ndim == 2 

803 n_kpt = omega_w.shape[0] 

804 omega_w = omega_w.ravel() 

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

806 else: 

807 # Return a partial DOS 

808 omegas, amplitudes = self.band_structure(kpts_kc, 

809 modes=True, 

810 verbose=verbose) 

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

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

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

814 assert ampl_sq.ndim == 3 

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

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

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

818 return dos 

819 

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

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

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

823 

824 Parameters: 

825 

826 kpts: tuple 

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

828 npts: int 

829 Number of energy points. 

830 delta: float 

831 Broadening of Lorentzian line-shape in eV. 

832 

833 Returns: 

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

835 

836 .. deprecated:: 3.23.1 

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

838 RawDOSData object. 

839 """ 

840 

841 # Monkhorst-Pack grid 

842 kpts_kc = monkhorst_pack(kpts) 

843 N = np.prod(kpts) 

844 # Get frequencies 

845 omega_kl = self.band_structure(kpts_kc) 

846 # Energy axis and dos 

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

848 dos_e = np.zeros_like(omega_e) 

849 

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

851 for omega_l in omega_kl: 

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

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

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

855 

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

857 

858 return omega_e, dos_e 

859 

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

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

862 """Write modes to trajectory file. 

863 

864 Parameters: 

865 

866 q_c: ndarray of shape (3,) 

867 q-vector of the modes. 

868 branches: int or list 

869 Branch index of modes. 

870 kT: float 

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

872 displacements in the modes. 

873 born: bool 

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

875 repeat: tuple 

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

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

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

879 nimages: int 

880 Number of images in an oscillation. 

881 center: bool 

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

883 

884 To exaggerate the amplitudes for better visualization, multiply 

885 kT by the square of the desired factor. 

886 """ 

887 

888 if isinstance(branches, int): 

889 branch_l = [branches] 

890 else: 

891 branch_l = list(branches) 

892 

893 # Calculate modes 

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

895 # Repeat atoms 

896 atoms = self.atoms * repeat 

897 # Center 

898 if center: