Coverage for /builds/hweiske/ase/ase/phonons.py: 77.40%
323 statements
« prev ^ index » next coverage.py v7.2.7, created at 2024-04-22 11:22 +0000
1"""Module for calculating phonons of periodic systems."""
3import warnings
4from math import pi, sqrt
5from pathlib import Path
7import numpy as np
8import numpy.fft as fft
9import numpy.linalg as la
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
20class Displacement:
21 """Abstract base class for phonon and el-ph supercell calculations.
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.
30 Derived classes must overwrite the ``__call__`` member function which is
31 called for each atomic displacement.
33 """
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.
39 Parameters:
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 """
61 # Store atoms and calculator
62 self.atoms = atoms
63 self.calc = calc
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
72 if comm is None:
73 comm = world
74 self.comm = comm
76 self.cache = MultiFileJSONCache(self.name)
78 def define_offset(self): # Reference cell offset
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
91 @property
92 @ase.utils.deprecated('Please use phonons.supercell instead of .N_c')
93 def N_c(self):
94 return self._supercell
96 @property
97 def supercell(self):
98 return self._supercell
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()
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()
112 def compute_lattice_vectors(self):
113 """Return lattice vectors for cells in the supercell."""
114 # Lattice vectors -- ordered as illustrated in class docstring
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
125 def __call__(self, *args, **kwargs):
126 """Member function called in the ``run`` function."""
128 raise NotImplementedError("Implement in derived classes!.")
130 def set_atoms(self, atoms):
131 """Set the atoms to vibrate.
133 Parameters:
135 atoms: list
136 Can be either a list of strings, ints or ...
138 """
140 assert isinstance(atoms, list)
141 assert len(atoms) <= len(self.atoms)
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
155 self.indices = indices
157 def _eq_disp(self):
158 return self._disp(0, 0, 0)
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)
164 def run(self):
165 """Run the calculations for the required displacements.
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.
173 """
175 # Atoms in the supercell -- repeated in the lattice vector directions
176 # beginning with the last
177 atoms_N = self.atoms * self.supercell
179 # Set calculator if provided
180 assert self.calc is not None, "Provide calculator in __init__ method"
181 atoms_N.calc = self.calc
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)
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()
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
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]
213 self.comm.barrier()
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
224 def _clean(self):
225 name = Path(self.name)
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
236class Phonons(Displacement):
237 r"""Class for calculating phonon modes using the finite displacement method.
239 The matrix of force constants is calculated from the finite difference
240 approximation to the first-order derivative of the atomic forces as::
242 2 nbj nbj
243 nbj d E F- - F+
244 C = ------------ ~ ------------- ,
245 mai dR dR 2 * delta
246 mai nbj
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::
253 nbj mai bj ai
254 C = C -> C (R ) = C (-R ) .
255 mai nbj ai n bj n
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.
261 The acoustic sum-rule::
263 _ _
264 aj \ bj
265 C (R ) = - ) C (R )
266 ai 0 /__ ai m
267 (m, b)
268 !=
269 (0, a)
271 Ordering of the unit cells illustrated here for a 1-dimensional system (in
272 case ``refcell=None`` in constructor!):
274 ::
276 m = 0 m = 1 m = -2 m = -1
277 -----------------------------------------------------
278 | | | | |
279 | * b | * | * | * |
280 | | | | |
281 | * a | * | * | * |
282 | | | | |
283 -----------------------------------------------------
285 Example:
287 >>> from ase.build import bulk
288 >>> from ase.phonons import Phonons
289 >>> from gpaw import GPAW, FermiDirac
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)
300 """
302 def __init__(self, *args, **kwargs):
303 """Initialize with base class args and kwargs."""
305 if 'name' not in kwargs:
306 kwargs['name'] = "phonon"
308 self.deprecate_refcell(kwargs)
310 Displacement.__init__(self, *args, **kwargs)
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
316 # Attributes for born charges and static dielectric tensor
317 self.Z_avv = None
318 self.eps_vv = None
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')
328 return kwargs
330 def __call__(self, atoms_N):
331 """Calculate forces on atoms in supercell."""
332 return atoms_N.get_forces()
334 def calculate(self, atoms_N, disp):
335 forces = self(atoms_N)
336 return {'forces': forces}
338 def check_eq_forces(self):
339 """Check maximum size of forces in the equilibrium structure."""
341 eq_disp = self._eq_disp()
342 feq_av = self.cache[eq_disp.name]['forces']
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)
349 return fmin, fmax, i_min, i_max
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.
357 The charge neutrality sum-rule::
359 _ _
360 \ a
361 ) Z = 0
362 /__ ij
363 a
365 Parameters:
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.
374 .. deprecated:: 3.22.1
375 Current implementation of non-analytical correction is likely
376 incorrect, see :issue:`941`
377 """
379 # Load file with Born charges and dielectric tensor for atoms in the
380 # unit cell
381 Z_avv, eps_vv = self.cache[name]
383 # Neutrality sum-rule
384 if neutrality:
385 Z_mean = Z_avv.sum(0) / len(Z_avv)
386 Z_avv -= Z_mean
388 self.Z_avv = Z_avv[self.indices]
389 self.eps_vv = eps_vv
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.
395 Extra keyword arguments will be passed to ``read_born_charges``.
397 Parameters:
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.
416 """
418 method = method.lower()
419 assert method in ['standard', 'frederiksen']
420 if cutoff is not None:
421 cutoff = float(cutoff)
423 # Read Born effective charges and optical dielectric tensor
424 if born:
425 self.read_born_charges(**kwargs)
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)
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']
444 if method == 'frederiksen':
445 fminus_av[a] -= fminus_av.sum(0)
446 fplus_av[a] -= fplus_av.sum(0)
448 # Finite difference derivative
449 C_av = fminus_av - fplus_av
450 C_av /= 2 * self.delta
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
457 # Make unitcell index the first and reshape
458 C_N = C_xNav.swapaxes(0, 1).reshape((N,) + (3 * natoms, 3 * natoms))
460 # Cut off before symmetry and acoustic sum rule are imposed
461 if cutoff is not None:
462 self.apply_cutoff(C_N, cutoff)
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
475 # Store force constants and dynamical matrix
476 self.C_N = C_N
477 self.D_N = C_N.copy()
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
486 def symmetrize(self, C_N):
487 """Symmetrize force constant matrix."""
489 # Number of atoms
490 natoms = len(self.indices)
491 # Number of unit cells
492 N = np.prod(self.supercell)
494 # Reshape force constants to (l, m, n) cell indices
495 C_lmn = C_N.reshape(self.supercell + (3 * natoms, 3 * natoms))
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()
509 # Change to single unit cell index shape
510 C_N = C_lmn.reshape((N, 3 * natoms, 3 * natoms))
512 return C_N
514 def acoustic(self, C_N):
515 """Restore acoustic sumrule on force constants."""
517 # Number of atoms
518 natoms = len(self.indices)
519 # Copy force constants
520 C_N_temp = C_N.copy()
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]
531 def apply_cutoff(self, D_N, r_c):
532 """Zero elements for interatomic distances larger than the cutoff.
534 Parameters:
536 D_N: ndarray
537 Dynamical/force constant matrix.
538 r_c: float
539 Cutoff in Angstrom.
541 """
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))
551 # Cell vectors
552 cell_vc = self.atoms.cell.transpose()
553 # Atomic positions in reference cell
554 pos_av = self.atoms.get_positions()
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
571 def get_force_constant(self):
572 """Return matrix of force constants."""
574 assert self.C_N is not None
575 return self.C_N
577 def get_band_structure(self, path, modes=False, born=False, verbose=True):
578 """Calculate and return the phonon band structure.
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.
585 Parameters:
587 path : BandPath object
588 The BandPath object defining the path in the reciprocal
589 space over which the phonon band structure is calculated.
590 modes : bool, optional
591 If True, phonon modes will also be calculated and returned.
592 Defaults to False.
593 born : bool, optional
594 If True, includes the effect of Born effective charges in
595 the phonon calculations.
596 Defaults to False.
597 verbose : bool, optional
598 If True, enables verbose output during the calculation.
599 Defaults to True.
601 Returns:
603 BandStructure or tuple of (BandStructure, ndarray)
604 If `modes` is False, returns a `BandStructure` object
605 containing the phonon band structure. If `modes` is True,
606 returns a tuple, where the first element is the
607 `BandStructure` object and the second element is an ndarray
608 of phonon modes.
610 Example:
612 >>> from ase.dft.kpoints import BandPath
613 >>> path = BandPath(...) # Define the band path
614 >>> phonons = Phonons(...)
615 >>> bs, modes = phonons.get_band_structure(path, modes=True)
616 """
617 result = self.band_structure(path.kpts,
618 modes=modes,
619 born=born,
620 verbose=verbose)
621 if modes:
622 omega_kl, omega_modes = result
623 else:
624 omega_kl = result
626 from ase.spectrum.band_structure import BandStructure
627 bs = BandStructure(path, energies=omega_kl[None])
629 # Return based on the modes flag
630 return (bs, omega_modes) if modes else bs
632 def compute_dynamical_matrix(self, q_scaled: np.ndarray, D_N: np.ndarray):
633 """ Computation of the dynamical matrix in momentum space D_ab(q).
634 This is a Fourier transform from real-space dynamical matrix D_N
635 for a given momentum vector q.
637 q_scaled: q vector in scaled coordinates.
639 D_N: the dynamical matrix in real-space. It is necessary, at least
640 currently, to provide this matrix explicitly (rather than use
641 self.D_N) because this matrix is modified by the Born charges
642 contributions and these modifications are momentum (q) dependent.
644 Result:
645 D(q): two-dimensional, complex-valued array of
646 shape=(3 * natoms, 3 * natoms).
647 """
648 # Evaluate fourier sum
649 R_cN = self._lattice_vectors_array
650 phase_N = np.exp(-2.j * pi * np.dot(q_scaled, R_cN))
651 D_q = np.sum(phase_N[:, np.newaxis, np.newaxis] * D_N, axis=0)
652 return D_q
654 def band_structure(self, path_kc, modes=False, born=False, verbose=True):
655 """Calculate phonon dispersion along a path in the Brillouin zone.
657 The dynamical matrix at arbitrary q-vectors is obtained by Fourier
658 transforming the real-space force constants. In case of negative
659 eigenvalues (squared frequency), the corresponding negative frequency
660 is returned.
662 Frequencies and modes are in units of eV and Ang/sqrt(amu),
663 respectively.
665 Parameters:
667 path_kc: ndarray
668 List of k-point coordinates (in units of the reciprocal lattice
669 vectors) specifying the path in the Brillouin zone for which the
670 dynamical matrix will be calculated.
671 modes: bool
672 Returns both frequencies and modes when True.
673 born: bool
674 Include non-analytic part given by the Born effective charges and
675 the static part of the high-frequency dielectric tensor. This
676 contribution to the force constant accounts for the splitting
677 between the LO and TO branches for q -> 0.
678 verbose: bool
679 Print warnings when imaginary frequncies are detected.
681 """
683 assert self.D_N is not None
684 if born:
685 assert self.Z_avv is not None
686 assert self.eps_vv is not None
688 # Dynamical matrix in real-space
689 D_N = self.D_N
691 # Lists for frequencies and modes along path
692 omega_kl = []
693 u_kl = []
695 # Reciprocal basis vectors for use in non-analytic contribution
696 reci_vc = 2 * pi * la.inv(self.atoms.cell)
697 # Unit cell volume in Bohr^3
698 vol = abs(la.det(self.atoms.cell)) / units.Bohr**3
700 for q_c in path_kc:
702 # Add non-analytic part
703 if born:
704 # q-vector in cartesian coordinates
705 q_v = np.dot(reci_vc, q_c)
706 # Non-analytic contribution to force constants in atomic units
707 qdotZ_av = np.dot(q_v, self.Z_avv).ravel()
708 C_na = (4 * pi * np.outer(qdotZ_av, qdotZ_av) /
709 np.dot(q_v, np.dot(self.eps_vv, q_v)) / vol)
710 self.C_na = C_na / units.Bohr**2 * units.Hartree
711 # Add mass prefactor and convert to eV / (Ang^2 * amu)
712 M_inv = np.outer(self.m_inv_x, self.m_inv_x)
713 D_na = C_na * M_inv / units.Bohr**2 * units.Hartree
714 self.D_na = D_na
715 D_N = self.D_N + D_na / np.prod(self.supercell)
717 # if np.prod(self.N_c) == 1:
718 #
719 # q_av = np.tile(q_v, len(self.indices))
720 # q_xx = np.vstack([q_av]*len(self.indices)*3)
721 # D_m += q_xx
723 # Evaluate fourier sum
724 D_q = self.compute_dynamical_matrix(q_c, D_N)
726 if modes:
727 omega2_l, u_xl = la.eigh(D_q, UPLO='U')
728 # Sort eigenmodes according to eigenvalues (see below) and
729 # multiply with mass prefactor
730 u_lx = (self.m_inv_x[:, np.newaxis] *
731 u_xl[:, omega2_l.argsort()]).T.copy()
732 u_kl.append(u_lx.reshape((-1, len(self.indices), 3)))
733 else:
734 omega2_l = la.eigvalsh(D_q, UPLO='U')
736 # Sort eigenvalues in increasing order
737 omega2_l.sort()
738 # Use dtype=complex to handle negative eigenvalues
739 omega_l = np.sqrt(omega2_l.astype(complex))
741 # Take care of imaginary frequencies
742 if not np.all(omega2_l >= 0.):
743 indices = np.where(omega2_l < 0)[0]
745 if verbose:
746 print('WARNING, %i imaginary frequencies at '
747 'q = (% 5.2f, % 5.2f, % 5.2f) ; (omega_q =% 5.3e*i)'
748 % (len(indices), q_c[0], q_c[1], q_c[2],
749 omega_l[indices][0].imag))
751 omega_l[indices] = -1 * np.sqrt(np.abs(omega2_l[indices].real))
753 omega_kl.append(omega_l.real)
755 # Conversion factor: sqrt(eV / Ang^2 / amu) -> eV
756 s = units._hbar * 1e10 / sqrt(units._e * units._amu)
757 omega_kl = s * np.asarray(omega_kl)
759 if modes:
760 return omega_kl, np.asarray(u_kl)
762 return omega_kl
764 def get_dos(self, kpts=(10, 10, 10), npts=1000, delta=1e-3, indices=None):
765 from ase.spectrum.dosdata import RawDOSData
767 # dos = self.dos(kpts, npts, delta, indices)
768 kpts_kc = monkhorst_pack(kpts)
769 omega_w = self.band_structure(kpts_kc).ravel()
770 dos = RawDOSData(omega_w, np.ones_like(omega_w))
771 return dos
773 def dos(self, kpts=(10, 10, 10), npts=1000, delta=1e-3, indices=None):
774 """Calculate phonon dos as a function of energy.
776 Parameters:
778 qpts: tuple
779 Shape of Monkhorst-Pack grid for sampling the Brillouin zone.
780 npts: int
781 Number of energy points.
782 delta: float
783 Broadening of Lorentzian line-shape in eV.
784 indices: list
785 If indices is not None, the atomic-partial dos for the specified
786 atoms will be calculated.
788 """
790 # Monkhorst-Pack grid
791 kpts_kc = monkhorst_pack(kpts)
792 N = np.prod(kpts)
793 # Get frequencies
794 omega_kl = self.band_structure(kpts_kc)
795 # Energy axis and dos
796 omega_e = np.linspace(0., np.amax(omega_kl) + 5e-3, num=npts)
797 dos_e = np.zeros_like(omega_e)
799 # Sum up contribution from all q-points and branches
800 for omega_l in omega_kl:
801 diff_el = (omega_e[:, np.newaxis] - omega_l[np.newaxis, :])**2
802 dos_el = 1. / (diff_el + (0.5 * delta)**2)
803 dos_e += dos_el.sum(axis=1)
805 dos_e *= 1. / (N * pi) * 0.5 * delta
807 return omega_e, dos_e
809 def write_modes(self, q_c, branches=0, kT=units.kB * 300, born=False,
810 repeat=(1, 1, 1), nimages=30, center=False):
811 """Write modes to trajectory file.
813 Parameters:
815 q_c: ndarray
816 q-vector of the modes.
817 branches: int or list
818 Branch index of modes.
819 kT: float
820 Temperature in units of eV. Determines the amplitude of the atomic
821 displacements in the modes.
822 born: bool
823 Include non-analytic contribution to the force constants at q -> 0.
824 repeat: tuple
825 Repeat atoms (l, m, n) times in the directions of the lattice
826 vectors. Displacements of atoms in repeated cells carry a Bloch
827 phase factor given by the q-vector and the cell lattice vector R_m.
828 nimages: int
829 Number of images in an oscillation.
830 center: bool
831 Center atoms in unit cell if True (default: False).
833 """
835 if isinstance(branches, int):
836 branch_l = [branches]
837 else:
838 branch_l = list(branches)
840 # Calculate modes
841 omega_l, u_l = self.band_structure([q_c], modes=True, born=born)
842 # Repeat atoms
843 atoms = self.atoms * repeat
844 # Center
845 if center:
846 atoms.center()
848 # Here ``Na`` refers to a composite unit cell/atom dimension
849 pos_Nav = atoms.get_positions()
850 # Total number of unit cells
851 N = np.prod(repeat)
853 # Corresponding lattice vectors R_m
854 R_cN = np.indices(repeat).reshape(3, -1)
855 # Bloch phase
856 phase_N = np.exp(2.j * pi * np.dot(q_c, R_cN))
857 phase_Na = phase_N.repeat(len(self.atoms))
859 for lval in branch_l:
861 omega = omega_l[0, lval]
862 u_av = u_l[0, lval]
864 # Mean displacement of a classical oscillator at temperature T
865 u_av *= sqrt(kT) / abs(omega)
867 mode_av = np.zeros((len(self.atoms), 3), dtype=complex)
868 # Insert slice with atomic displacements for the included atoms
869 mode_av[self.indices] = u_av
870 # Repeat and multiply by Bloch phase factor
871 mode_Nav = np.vstack(N * [mode_av]) * phase_Na[:, np.newaxis]
873 with Trajectory('%s.mode.%d.traj'
874 % (self.name, lval), 'w') as traj:
875 for x in np.linspace(0, 2 * pi, nimages, endpoint=False):
876 atoms.set_positions((pos_Nav + np.exp(1.j * x) *
877 mode_Nav).real)
878 traj.write(atoms)