Coverage for /builds/debichem-team/python-ase/ase/phonons.py: 82.88%

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 

12import ase.units as 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__(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 

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 

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

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

293 ... kpts=(5, 5, 5), 

294 ... h=0.2, 

295 ... occupations=FermiDirac(0.)) 

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

297 >>> ph.run() 

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

299 

300 """ 

301 

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

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

304 

305 if 'name' not in kwargs: 

306 kwargs['name'] = "phonon" 

307 

308 self.deprecate_refcell(kwargs) 

309 

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

311 

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

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

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

315 

316 # Attributes for born charges and static dielectric tensor 

317 self.Z_avv = None 

318 self.eps_vv = None 

319 

320 @staticmethod 

321 def deprecate_refcell(kwargs: dict): 

322 if 'refcell' in kwargs: 

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

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

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

326 kwargs.pop('refcell') 

327 

328 return kwargs 

329 

330 def __call__(self, atoms_N): 

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

332 return atoms_N.get_forces() 

333 

334 def calculate(self, atoms_N, disp): 

335 forces = self(atoms_N) 

336 return {'forces': forces} 

337 

338 def check_eq_forces(self): 

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

340 

341 eq_disp = self._eq_disp() 

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

343 

344 fmin = feq_av.min() 

345 fmax = feq_av.max() 

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

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

348 

349 return fmin, fmax, i_min, i_max 

350 

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

352 'likely incorrect, see ' 

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

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

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

356 

357 The charge neutrality sum-rule:: 

358 

359 _ _ 

360 \ a 

361 ) Z = 0 

362 /__ ij 

363 a 

364 

365 Parameters: 

366 

367 neutrality: bool 

368 Restore charge neutrality condition on calculated Born effective 

369 charges. 

370 name: str 

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

372 in the JSON cache. 

373 

374 .. deprecated:: 3.22.1 

375 Current implementation of non-analytical correction is likely 

376 incorrect, see :issue:`941` 

377 """ 

378 

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

380 # unit cell 

381 Z_avv, eps_vv = self.cache[name] 

382 

383 # Neutrality sum-rule 

384 if neutrality: 

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

386 Z_avv -= Z_mean 

387 

388 self.Z_avv = Z_avv[self.indices] 

389 self.eps_vv = eps_vv 

390 

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

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

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

394 

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

396 

397 Parameters: 

398 

399 method: str 

400 Specify method for evaluating the atomic forces. 

401 symmetrize: int 

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

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

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

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

406 number of iterations that will be carried out. 

407 acoustic: bool 

408 Restore the acoustic sum rule on the force constants. 

409 cutoff: None or float 

410 Zero elements in the dynamical matrix between atoms with an 

411 interatomic distance larger than the cutoff. 

412 born: bool 

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

414 dielelctric tensor from file. 

415 

416 """ 

417 

418 method = method.lower() 

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

420 if cutoff is not None: 

421 cutoff = float(cutoff) 

422 

423 # Read Born effective charges and optical dielectric tensor 

424 if born: 

425 self.read_born_charges(**kwargs) 

426 

427 # Number of atoms 

428 natoms = len(self.indices) 

429 # Number of unit cells 

430 N = np.prod(self.supercell) 

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

432 # of eV / Ang**2 

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

434 

435 # Loop over all atomic displacements and calculate force constants 

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

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

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

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

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

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

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

443 

444 if method == 'frederiksen': 

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

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

447 

448 # Finite difference derivative 

449 C_av = fminus_av - fplus_av 

450 C_av /= 2 * self.delta 

451 

452 # Slice out included atoms 

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

454 index = 3 * i + j 

455 C_xNav[index] = C_Nav 

456 

457 # Make unitcell index the first and reshape 

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

459 

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

461 if cutoff is not None: 

462 self.apply_cutoff(C_N, cutoff) 

463 

464 # Symmetrize force constants 

465 if symmetrize: 

466 for _ in range(symmetrize): 

467 # Symmetrize 

468 C_N = self.symmetrize(C_N) 

469 # Restore acoustic sum-rule 

470 if acoustic: 

471 self.acoustic(C_N) 

472 else: 

473 break 

474 

475 # Store force constants and dynamical matrix 

476 self.C_N = C_N 

477 self.D_N = C_N.copy() 

478 

479 # Add mass prefactor 

480 m_a = self.atoms.get_masses() 

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

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

483 for D in self.D_N: 

484 D *= M_inv 

485 

486 def symmetrize(self, C_N): 

487 """Symmetrize force constant matrix.""" 

488 

489 # Number of atoms 

490 natoms = len(self.indices) 

491 # Number of unit cells 

492 N = np.prod(self.supercell) 

493 

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

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

496 

497 # Shift reference cell to center index 

498 if self.offset == 0: 

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

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

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

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

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

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

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

506 if self.offset == 0: 

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

508 

509 # Change to single unit cell index shape 

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

511 

512 return C_N 

513 

514 def acoustic(self, C_N): 

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

516 

517 # Number of atoms 

518 natoms = len(self.indices) 

519 # Copy force constants 

520 C_N_temp = C_N.copy() 

521 

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

523 for C in C_N_temp: 

524 for a in range(natoms): 

525 for a_ in range(natoms): 

526 C_N[self.offset, 

527 3 * a: 3 * a + 3, 

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

529 3 * a_: 3 * a_ + 3] 

530 

531 def apply_cutoff(self, D_N, r_c): 

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

533 

534 Parameters: 

535 

536 D_N: ndarray 

537 Dynamical/force constant matrix. 

538 r_c: float 

539 Cutoff in Angstrom. 

540 

541 """ 

542 

543 # Number of atoms and primitive cells 

544 natoms = len(self.indices) 

545 N = np.prod(self.supercell) 

546 # Lattice vectors 

547 R_cN = self._lattice_vectors_array 

548 # Reshape matrix to individual atomic and cartesian dimensions 

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

550 

551 # Cell vectors 

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

553 # Atomic positions in reference cell 

554 pos_av = self.atoms.get_positions() 

555 

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

557 # larger than the cutoff 

558 for n in range(N): 

559 # Lattice vector to cell 

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

561 # Atomic positions in cell 

562 posn_av = pos_av + R_v 

563 # Loop over atoms and zero elements 

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

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

566 # Atoms where the distance is larger than the cufoff 

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

568 # Zero elements 

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

570 

571 def get_force_constant(self): 

572 """Return matrix of force constants.""" 

573 

574 assert self.C_N is not None 

575 return self.C_N 

576 

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

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

579 

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

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

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

583 optionally calculate and return phonon modes. 

584 

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

586 respectively. 

587 

588 Parameters: 

589 

590 path : BandPath object 

591 The BandPath object defining the path in the reciprocal 

592 space over which the phonon band structure is calculated. 

593 modes : bool, optional 

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

595 Defaults to False. 

596 born : bool, optional 

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

598 the phonon calculations. 

599 Defaults to False. 

600 verbose : bool, optional 

601 If True, enables verbose output during the calculation. 

602 Defaults to True. 

603 

604 Returns: 

605 

606 BandStructure or tuple of (BandStructure, ndarray) 

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

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

609 returns a tuple, where the first element is the 

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

611 of phonon modes. 

612 

613 If modes are returned, the array is of shape 

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

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

616 mode. 

617 

618 Example: 

619 

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

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

622 >>> phonons = Phonons(...) 

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

624 """ 

625 result = self.band_structure(path.kpts, 

626 modes=modes, 

627 born=born, 

628 verbose=verbose) 

629 if modes: 

630 omega_kl, omega_modes = result 

631 else: 

632 omega_kl = result 

633 

634 from ase.spectrum.band_structure import BandStructure 

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

636 

637 # Return based on the modes flag 

638 return (bs, omega_modes) if modes else bs 

639 

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

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

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

643 for a given momentum vector q. 

644 

645 q_scaled: q vector in scaled coordinates. 

646 

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

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

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

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

651 

652 Result: 

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

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

655 """ 

656 # Evaluate fourier sum 

657 R_cN = self._lattice_vectors_array 

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

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

660 return D_q 

661 

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

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

664 

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

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

667 eigenvalues (squared frequency), the corresponding negative frequency 

668 is returned. 

669 

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

671 respectively. 

672 

673 Parameters: 

674 

675 path_kc: ndarray 

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

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

678 dynamical matrix will be calculated. 

679 modes: bool 

680 Returns both frequencies and modes when True. 

681 born: bool 

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

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

684 contribution to the force constant accounts for the splitting 

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

686 verbose: bool 

687 Print warnings when imaginary frequncies are detected. 

688 

689 Returns: 

690 

691 If modes=False: Array of energies 

692 

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

694 """ 

695 

696 assert self.D_N is not None 

697 if born: 

698 assert self.Z_avv is not None 

699 assert self.eps_vv is not None 

700 

701 # Dynamical matrix in real-space 

702 D_N = self.D_N 

703 

704 # Lists for frequencies and modes along path 

705 omega_kl = [] 

706 u_kl = [] 

707 

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

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

710 # Unit cell volume in Bohr^3 

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

712 

713 for q_c in path_kc: 

714 # Add non-analytic part 

715 if born: 

716 # q-vector in cartesian coordinates 

717 q_v = np.dot(reci_vc, q_c) 

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

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

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

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

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

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

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

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

726 self.D_na = D_na 

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

728 

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

730 # 

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

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

733 # D_m += q_xx 

734 

735 # Evaluate fourier sum 

736 D_q = self.compute_dynamical_matrix(q_c, D_N) 

737 

738 if modes: 

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

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

741 # multiply with mass prefactor. This gives the eigenmode 

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

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

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

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

746 else: 

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

748 

749 # Sort eigenvalues in increasing order 

750 omega2_l.sort() 

751 # Use dtype=complex to handle negative eigenvalues 

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

753 

754 # Take care of imaginary frequencies 

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

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

757 

758 if verbose: 

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

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

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

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

763 

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

765 

766 omega_kl.append(omega_l.real) 

767 

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

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

770 omega_kl = s * np.asarray(omega_kl) 

771 

772 if modes: 

773 return omega_kl, np.asarray(u_kl) 

774 

775 return omega_kl 

776 

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

778 """Return a phonon density of states. 

779 

780 Parameters: 

781 

782 kpts: tuple 

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

784 indices: list 

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

786 DOS for the specified atoms will be calculated. 

787 verbose: bool 

788 Print warnings when imaginary frequncies are detected. 

789 

790 Returns: 

791 A RawDOSData object containing the density of states. 

792 """ 

793 from ase.spectrum.dosdata import RawDOSData 

794 

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

796 kpts_kc = monkhorst_pack(kpts) 

797 if indices is None: 

798 # Return the total DOS 

799 omega_w = self.band_structure(kpts_kc, verbose=verbose).ravel() 

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

801 else: 

802 # Return a partial DOS 

803 omegas, amplitudes = self.band_structure(kpts_kc, 

804 modes=True, 

805 verbose=verbose) 

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

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

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

809 assert ampl_sq.ndim == 3 

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

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

812 dos = RawDOSData(omegas.ravel(), weights.ravel()) 

813 return dos 

814 

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

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

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

818 

819 Parameters: 

820 

821 kpts: tuple 

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

823 npts: int 

824 Number of energy points. 

825 delta: float 

826 Broadening of Lorentzian line-shape in eV. 

827 

828 Returns: 

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

830 

831 .. deprecated:: 3.23.1 

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

833 RawDOSData object. 

834 """ 

835 

836 # Monkhorst-Pack grid 

837 kpts_kc = monkhorst_pack(kpts) 

838 N = np.prod(kpts) 

839 # Get frequencies 

840 omega_kl = self.band_structure(kpts_kc) 

841 # Energy axis and dos 

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

843 dos_e = np.zeros_like(omega_e) 

844 

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

846 for omega_l in omega_kl: 

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

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

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

850 

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

852 

853 return omega_e, dos_e 

854 

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

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

857 """Write modes to trajectory file. 

858 

859 Parameters: 

860 

861 q_c: ndarray of shape (3,) 

862 q-vector of the modes. 

863 branches: int or list 

864 Branch index of modes. 

865 kT: float 

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

867 displacements in the modes. 

868 born: bool 

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

870 repeat: tuple 

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

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

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

874 nimages: int 

875 Number of images in an oscillation. 

876 center: bool 

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

878 

879 To exaggerate the amplitudes for better visualization, multiply 

880 kT by the square of the desired factor. 

881 """ 

882 

883 if isinstance(branches, int): 

884 branch_l = [branches] 

885 else: 

886 branch_l = list(branches) 

887 

888 # Calculate modes 

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

890 # Repeat atoms 

891 atoms = self.atoms * repeat 

892 # Center 

893 if center: 

894 atoms.center() 

895 

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

897 pos_Nav = atoms.get_positions() 

898 # Total number of unit cells 

899 N = np.prod(repeat)