Coverage for /builds/ase/ase/ase/phonons.py : 74.92%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""Module for calculating phonons of periodic systems."""
3from math import pi, sqrt
4import warnings
5from pathlib import Path
7import numpy as np
8import numpy.linalg as la
9import numpy.fft as fft
11import ase
12import ase.units as units
13from ase.parallel import world
14from ase.dft import monkhorst_pack
15from ase.io.trajectory import Trajectory
16from ase.utils.filecache import MultiFileJSONCache
17from ase.utils import deprecated
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):
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.
57 """
59 # Store atoms and calculator
60 self.atoms = atoms
61 self.calc = calc
63 # Displace all atoms in the unit cell by default
64 self.indices = np.arange(len(atoms))
65 self.name = name
66 self.delta = delta
67 self.center_refcell = center_refcell
68 self.supercell = supercell
70 self.cache = MultiFileJSONCache(self.name)
72 def define_offset(self): # Reference cell offset
74 if not self.center_refcell:
75 # Corner cell
76 self.offset = 0
77 else:
78 # Center cell
79 N_c = self.supercell
80 self.offset = (N_c[0] // 2 * (N_c[1] * N_c[2]) +
81 N_c[1] // 2 * N_c[2] +
82 N_c[2] // 2)
83 return self.offset
85 @property # type: ignore
86 @ase.utils.deprecated('Please use phonons.supercell instead of .N_c')
87 def N_c(self):
88 return self._supercell
90 @property
91 def supercell(self):
92 return self._supercell
94 @supercell.setter
95 def supercell(self, supercell):
96 assert len(supercell) == 3
97 self._supercell = tuple(supercell)
98 self.define_offset()
99 self._lattice_vectors_array = self.compute_lattice_vectors()
101 @ase.utils.deprecated('Please use phonons.compute_lattice_vectors()'
102 ' instead of .lattice_vectors()')
103 def lattice_vectors(self):
104 return self.compute_lattice_vectors()
106 def compute_lattice_vectors(self):
107 """Return lattice vectors for cells in the supercell."""
108 # Lattice vectors -- ordered as illustrated in class docstring
110 # Lattice vectors relevative to the reference cell
111 R_cN = np.indices(self.supercell).reshape(3, -1)
112 N_c = np.array(self.supercell)[:, np.newaxis]
113 if self.offset == 0:
114 R_cN += N_c // 2
115 R_cN %= N_c
116 R_cN -= N_c // 2
117 return R_cN
119 def __call__(self, *args, **kwargs):
120 """Member function called in the ``run`` function."""
122 raise NotImplementedError("Implement in derived classes!.")
124 def set_atoms(self, atoms):
125 """Set the atoms to vibrate.
127 Parameters:
129 atoms: list
130 Can be either a list of strings, ints or ...
132 """
134 assert isinstance(atoms, list)
135 assert len(atoms) <= len(self.atoms)
137 if isinstance(atoms[0], str):
138 assert np.all([isinstance(atom, str) for atom in atoms])
139 sym_a = self.atoms.get_chemical_symbols()
140 # List for atomic indices
141 indices = []
142 for type in atoms:
143 indices.extend([a for a, atom in enumerate(sym_a)
144 if atom == type])
145 else:
146 assert np.all([isinstance(atom, int) for atom in atoms])
147 indices = atoms
149 self.indices = indices
151 def _eq_disp(self):
152 return self._disp(0, 0, 0)
154 def _disp(self, a, i, step):
155 from ase.vibrations.vibrations import Displacement as VDisplacement
156 return VDisplacement(a, i, np.sign(step), abs(step), self)
158 def run(self):
159 """Run the calculations for the required displacements.
161 This will do a calculation for 6 displacements per atom, +-x, +-y, and
162 +-z. Only those calculations that are not already done will be
163 started. Be aware that an interrupted calculation may produce an empty
164 file (ending with .json), which must be deleted before restarting the
165 job. Otherwise the calculation for that displacement will not be done.
167 """
169 # Atoms in the supercell -- repeated in the lattice vector directions
170 # beginning with the last
171 atoms_N = self.atoms * self.supercell
173 # Set calculator if provided
174 assert self.calc is not None, "Provide calculator in __init__ method"
175 atoms_N.calc = self.calc
177 # Do calculation on equilibrium structure
178 eq_disp = self._eq_disp()
179 with self.cache.lock(eq_disp.name) as handle:
180 if handle is not None:
181 output = self.calculate(atoms_N, eq_disp)
182 handle.save(output)
184 # Positions of atoms to be displaced in the reference cell
185 natoms = len(self.atoms)
186 offset = natoms * self.offset
187 pos = atoms_N.positions[offset: offset + natoms].copy()
189 # Loop over all displacements
190 for a in self.indices:
191 for i in range(3):
192 for sign in [-1, 1]:
193 disp = self._disp(a, i, sign)
194 with self.cache.lock(disp.name) as handle:
195 if handle is None:
196 continue
197 try:
198 atoms_N.positions[offset + a, i] = \
199 pos[a, i] + sign * self.delta
201 result = self.calculate(atoms_N, disp)
202 handle.save(result)
203 finally:
204 # Return to initial positions
205 atoms_N.positions[offset + a, i] = pos[a, i]
207 def clean(self):
208 """Delete generated files."""
209 if world.rank != 0:
210 return 0
212 name = Path(self.name)
214 n = 0
215 if name.is_dir():
216 for fname in name.iterdir():
217 fname.unlink()
218 n += 1
219 name.rmdir()
220 return n
223class Phonons(Displacement):
224 r"""Class for calculating phonon modes using the finite displacement method.
226 The matrix of force constants is calculated from the finite difference
227 approximation to the first-order derivative of the atomic forces as::
229 2 nbj nbj
230 nbj d E F- - F+
231 C = ------------ ~ ------------- ,
232 mai dR dR 2 * delta
233 mai nbj
235 where F+/F- denotes the force in direction j on atom nb when atom ma is
236 displaced in direction +i/-i. The force constants are related by various
237 symmetry relations. From the definition of the force constants it must
238 be symmetric in the three indices mai::
240 nbj mai bj ai
241 C = C -> C (R ) = C (-R ) .
242 mai nbj ai n bj n
244 As the force constants can only depend on the difference between the m and
245 n indices, this symmetry is more conveniently expressed as shown on the
246 right hand-side.
248 The acoustic sum-rule::
250 _ _
251 aj \ bj
252 C (R ) = - ) C (R )
253 ai 0 /__ ai m
254 (m, b)
255 !=
256 (0, a)
258 Ordering of the unit cells illustrated here for a 1-dimensional system (in
259 case ``refcell=None`` in constructor!):
261 ::
263 m = 0 m = 1 m = -2 m = -1
264 -----------------------------------------------------
265 | | | | |
266 | * b | * | * | * |
267 | | | | |
268 | * a | * | * | * |
269 | | | | |
270 -----------------------------------------------------
272 Example:
274 >>> from ase.build import bulk
275 >>> from ase.phonons import Phonons
276 >>> from gpaw import GPAW, FermiDirac
277 >>> atoms = bulk('Si', 'diamond', a=5.4)
278 >>> calc = GPAW(kpts=(5, 5, 5),
279 h=0.2,
280 occupations=FermiDirac(0.))
281 >>> ph = Phonons(atoms, calc, supercell=(5, 5, 5))
282 >>> ph.run()
283 >>> ph.read(method='frederiksen', acoustic=True)
285 """
287 def __init__(self, *args, **kwargs):
288 """Initialize with base class args and kwargs."""
290 if 'name' not in kwargs:
291 kwargs['name'] = "phonon"
293 self.deprecate_refcell(kwargs)
295 Displacement.__init__(self, *args, **kwargs)
297 # Attributes for force constants and dynamical matrix in real space
298 self.C_N = None # in units of eV / Ang**2
299 self.D_N = None # in units of eV / Ang**2 / amu
301 # Attributes for born charges and static dielectric tensor
302 self.Z_avv = None
303 self.eps_vv = None
305 @staticmethod
306 def deprecate_refcell(kwargs: dict):
307 if 'refcell' in kwargs:
308 warnings.warn('Keyword refcell of Phonons is deprecated.'
309 'Please use center_refcell (bool)', FutureWarning)
310 kwargs['center_refcell'] = bool(kwargs['refcell'])
311 kwargs.pop('refcell')
313 return kwargs
315 def __call__(self, atoms_N):
316 """Calculate forces on atoms in supercell."""
317 return atoms_N.get_forces()
319 def calculate(self, atoms_N, disp):
320 forces = self(atoms_N)
321 return {'forces': forces}
323 def check_eq_forces(self):
324 """Check maximum size of forces in the equilibrium structure."""
326 eq_disp = self._eq_disp()
327 feq_av = self.cache[eq_disp.name]['forces']
329 fmin = feq_av.min()
330 fmax = feq_av.max()
331 i_min = np.where(feq_av == fmin)
332 i_max = np.where(feq_av == fmax)
334 return fmin, fmax, i_min, i_max
336 @deprecated('Current implementation of non-analytical correction is '
337 'likely incorrect, see '
338 'https://gitlab.com/ase/ase/-/issues/941')
339 def read_born_charges(self, name='born', neutrality=True):
340 r"""Read Born charges and dieletric tensor from JSON file.
342 The charge neutrality sum-rule::
344 _ _
345 \ a
346 ) Z = 0
347 /__ ij
348 a
350 Parameters:
352 neutrality: bool
353 Restore charge neutrality condition on calculated Born effective
354 charges.
355 name: str
356 Key used to identify the file with Born charges for the unit cell
357 in the JSON cache.
359 """
361 # Load file with Born charges and dielectric tensor for atoms in the
362 # unit cell
363 Z_avv, eps_vv = self.cache[name]
365 # Neutrality sum-rule
366 if neutrality:
367 Z_mean = Z_avv.sum(0) / len(Z_avv)
368 Z_avv -= Z_mean
370 self.Z_avv = Z_avv[self.indices]
371 self.eps_vv = eps_vv
373 def read(self, method='Frederiksen', symmetrize=3, acoustic=True,
374 cutoff=None, born=False, **kwargs):
375 """Read forces from json files and calculate force constants.
377 Extra keyword arguments will be passed to ``read_born_charges``.
379 Parameters:
381 method: str
382 Specify method for evaluating the atomic forces.
383 symmetrize: int
384 Symmetrize force constants (see doc string at top) when
385 ``symmetrize != 0`` (default: 3). Since restoring the acoustic sum
386 rule breaks the symmetry, the symmetrization must be repeated a few
387 times until the changes a insignificant. The integer gives the
388 number of iterations that will be carried out.
389 acoustic: bool
390 Restore the acoustic sum rule on the force constants.
391 cutoff: None or float
392 Zero elements in the dynamical matrix between atoms with an
393 interatomic distance larger than the cutoff.
394 born: bool
395 Read in Born effective charge tensor and high-frequency static
396 dielelctric tensor from file.
398 """
400 method = method.lower()
401 assert method in ['standard', 'frederiksen']
402 if cutoff is not None:
403 cutoff = float(cutoff)
405 # Read Born effective charges and optical dielectric tensor
406 if born:
407 self.read_born_charges(**kwargs)
409 # Number of atoms
410 natoms = len(self.indices)
411 # Number of unit cells
412 N = np.prod(self.supercell)
413 # Matrix of force constants as a function of unit cell index in units
414 # of eV / Ang**2
415 C_xNav = np.empty((natoms * 3, N, natoms, 3), dtype=float)
417 # Loop over all atomic displacements and calculate force constants
418 for i, a in enumerate(self.indices):
419 for j, v in enumerate('xyz'):
420 # Atomic forces for a displacement of atom a in direction v
421 # basename = '%s.%d%s' % (self.name, a, v)
422 basename = '%d%s' % (a, v)
423 fminus_av = self.cache[basename + '-']['forces']
424 fplus_av = self.cache[basename + '+']['forces']
426 if method == 'frederiksen':
427 fminus_av[a] -= fminus_av.sum(0)
428 fplus_av[a] -= fplus_av.sum(0)
430 # Finite difference derivative
431 C_av = fminus_av - fplus_av
432 C_av /= 2 * self.delta
434 # Slice out included atoms
435 C_Nav = C_av.reshape((N, len(self.atoms), 3))[:, self.indices]
436 index = 3 * i + j
437 C_xNav[index] = C_Nav
439 # Make unitcell index the first and reshape
440 C_N = C_xNav.swapaxes(0, 1).reshape((N,) + (3 * natoms, 3 * natoms))
442 # Cut off before symmetry and acoustic sum rule are imposed
443 if cutoff is not None:
444 self.apply_cutoff(C_N, cutoff)
446 # Symmetrize force constants
447 if symmetrize:
448 for i in range(symmetrize):
449 # Symmetrize
450 C_N = self.symmetrize(C_N)
451 # Restore acoustic sum-rule
452 if acoustic:
453 self.acoustic(C_N)
454 else:
455 break
457 # Store force constants and dynamical matrix
458 self.C_N = C_N
459 self.D_N = C_N.copy()
461 # Add mass prefactor
462 m_a = self.atoms.get_masses()
463 self.m_inv_x = np.repeat(m_a[self.indices]**-0.5, 3)
464 M_inv = np.outer(self.m_inv_x, self.m_inv_x)
465 for D in self.D_N:
466 D *= M_inv
468 def symmetrize(self, C_N):
469 """Symmetrize force constant matrix."""
471 # Number of atoms
472 natoms = len(self.indices)
473 # Number of unit cells
474 N = np.prod(self.supercell)
476 # Reshape force constants to (l, m, n) cell indices
477 C_lmn = C_N.reshape(self.supercell + (3 * natoms, 3 * natoms))
479 # Shift reference cell to center index
480 if self.offset == 0:
481 C_lmn = fft.fftshift(C_lmn, axes=(0, 1, 2)).copy()
482 # Make force constants symmetric in indices -- in case of an even
483 # number of unit cells don't include the first cell
484 i, j, k = 1 - np.asarray(self.supercell) % 2
485 C_lmn[i:, j:, k:] *= 0.5
486 C_lmn[i:, j:, k:] += \
487 C_lmn[i:, j:, k:][::-1, ::-1, ::-1].transpose(0, 1, 2, 4, 3).copy()
488 if self.offset == 0:
489 C_lmn = fft.ifftshift(C_lmn, axes=(0, 1, 2)).copy()
491 # Change to single unit cell index shape
492 C_N = C_lmn.reshape((N, 3 * natoms, 3 * natoms))
494 return C_N
496 def acoustic(self, C_N):
497 """Restore acoustic sumrule on force constants."""
499 # Number of atoms
500 natoms = len(self.indices)
501 # Copy force constants
502 C_N_temp = C_N.copy()
504 # Correct atomic diagonals of R_m = (0, 0, 0) matrix
505 for C in C_N_temp:
506 for a in range(natoms):
507 for a_ in range(natoms):
508 C_N[self.offset,
509 3 * a: 3 * a + 3,
510 3 * a: 3 * a + 3] -= C[3 * a: 3 * a + 3,
511 3 * a_: 3 * a_ + 3]
513 def apply_cutoff(self, D_N, r_c):
514 """Zero elements for interatomic distances larger than the cutoff.
516 Parameters:
518 D_N: ndarray
519 Dynamical/force constant matrix.
520 r_c: float
521 Cutoff in Angstrom.
523 """
525 # Number of atoms and primitive cells
526 natoms = len(self.indices)
527 N = np.prod(self.supercell)
528 # Lattice vectors
529 R_cN = self._lattice_vectors_array
530 # Reshape matrix to individual atomic and cartesian dimensions
531 D_Navav = D_N.reshape((N, natoms, 3, natoms, 3))
533 # Cell vectors
534 cell_vc = self.atoms.cell.transpose()
535 # Atomic positions in reference cell
536 pos_av = self.atoms.get_positions()
538 # Zero elements with a distance to atoms in the reference cell
539 # larger than the cutoff
540 for n in range(N):
541 # Lattice vector to cell
542 R_v = np.dot(cell_vc, R_cN[:, n])
543 # Atomic positions in cell
544 posn_av = pos_av + R_v
545 # Loop over atoms and zero elements
546 for i, a in enumerate(self.indices):
547 dist_a = np.sqrt(np.sum((pos_av[a] - posn_av)**2, axis=-1))
548 # Atoms where the distance is larger than the cufoff
549 i_a = dist_a > r_c # np.where(dist_a > r_c)
550 # Zero elements
551 D_Navav[n, i, :, i_a, :] = 0.0
553 def get_force_constant(self):
554 """Return matrix of force constants."""
556 assert self.C_N is not None
557 return self.C_N
559 def get_band_structure(self, path, modes=False, born=False, verbose=True):
560 omega_kl = self.band_structure(path.kpts, modes, born, verbose)
561 if modes:
562 assert 0
563 omega_kl, modes = omega_kl
565 from ase.spectrum.band_structure import BandStructure
566 bs = BandStructure(path, energies=omega_kl[None])
567 return bs
569 def compute_dynamical_matrix(self, q_scaled: np.ndarray, D_N: np.ndarray):
570 """ Computation of the dynamical matrix in momentum space D_ab(q).
571 This is a Fourier transform from real-space dynamical matrix D_N
572 for a given momentum vector q.
574 q_scaled: q vector in scaled coordinates.
576 D_N: the dynamical matrix in real-space. It is necessary, at least
577 currently, to provide this matrix explicitly (rather than use
578 self.D_N) because this matrix is modified by the Born charges
579 contributions and these modifications are momentum (q) dependent.
581 Result:
582 D(q): two-dimensional, complex-valued array of
583 shape=(3 * natoms, 3 * natoms).
584 """
585 # Evaluate fourier sum
586 R_cN = self._lattice_vectors_array
587 phase_N = np.exp(-2.j * pi * np.dot(q_scaled, R_cN))
588 D_q = np.sum(phase_N[:, np.newaxis, np.newaxis] * D_N, axis=0)
589 return D_q
591 def band_structure(self, path_kc, modes=False, born=False, verbose=True):
592 """Calculate phonon dispersion along a path in the Brillouin zone.
594 The dynamical matrix at arbitrary q-vectors is obtained by Fourier
595 transforming the real-space force constants. In case of negative
596 eigenvalues (squared frequency), the corresponding negative frequency
597 is returned.
599 Frequencies and modes are in units of eV and Ang/sqrt(amu),
600 respectively.
602 Parameters:
604 path_kc: ndarray
605 List of k-point coordinates (in units of the reciprocal lattice
606 vectors) specifying the path in the Brillouin zone for which the
607 dynamical matrix will be calculated.
608 modes: bool
609 Returns both frequencies and modes when True.
610 born: bool
611 Include non-analytic part given by the Born effective charges and
612 the static part of the high-frequency dielectric tensor. This
613 contribution to the force constant accounts for the splitting
614 between the LO and TO branches for q -> 0.
615 verbose: bool
616 Print warnings when imaginary frequncies are detected.
618 """
620 assert self.D_N is not None
621 if born:
622 assert self.Z_avv is not None
623 assert self.eps_vv is not None
625 # Dynamical matrix in real-space
626 D_N = self.D_N
628 # Lists for frequencies and modes along path
629 omega_kl = []
630 u_kl = []
632 # Reciprocal basis vectors for use in non-analytic contribution
633 reci_vc = 2 * pi * la.inv(self.atoms.cell)
634 # Unit cell volume in Bohr^3
635 vol = abs(la.det(self.atoms.cell)) / units.Bohr**3
637 for q_c in path_kc:
639 # Add non-analytic part
640 if born:
641 # q-vector in cartesian coordinates
642 q_v = np.dot(reci_vc, q_c)
643 # Non-analytic contribution to force constants in atomic units
644 qdotZ_av = np.dot(q_v, self.Z_avv).ravel()
645 C_na = (4 * pi * np.outer(qdotZ_av, qdotZ_av) /
646 np.dot(q_v, np.dot(self.eps_vv, q_v)) / vol)
647 self.C_na = C_na / units.Bohr**2 * units.Hartree
648 # Add mass prefactor and convert to eV / (Ang^2 * amu)
649 M_inv = np.outer(self.m_inv_x, self.m_inv_x)
650 D_na = C_na * M_inv / units.Bohr**2 * units.Hartree
651 self.D_na = D_na
652 D_N = self.D_N + D_na / np.prod(self.supercell)
654 # if np.prod(self.N_c) == 1:
655 #
656 # q_av = np.tile(q_v, len(self.indices))
657 # q_xx = np.vstack([q_av]*len(self.indices)*3)
658 # D_m += q_xx
660 # Evaluate fourier sum
661 D_q = self.compute_dynamical_matrix(q_c, D_N)
663 if modes:
664 omega2_l, u_xl = la.eigh(D_q, UPLO='U')
665 # Sort eigenmodes according to eigenvalues (see below) and
666 # multiply with mass prefactor
667 u_lx = (self.m_inv_x[:, np.newaxis] *
668 u_xl[:, omega2_l.argsort()]).T.copy()
669 u_kl.append(u_lx.reshape((-1, len(self.indices), 3)))
670 else:
671 omega2_l = la.eigvalsh(D_q, UPLO='U')
673 # Sort eigenvalues in increasing order
674 omega2_l.sort()
675 # Use dtype=complex to handle negative eigenvalues
676 omega_l = np.sqrt(omega2_l.astype(complex))
678 # Take care of imaginary frequencies
679 if not np.all(omega2_l >= 0.):
680 indices = np.where(omega2_l < 0)[0]
682 if verbose:
683 print('WARNING, %i imaginary frequencies at '
684 'q = (% 5.2f, % 5.2f, % 5.2f) ; (omega_q =% 5.3e*i)'
685 % (len(indices), q_c[0], q_c[1], q_c[2],
686 omega_l[indices][0].imag))
688 omega_l[indices] = -1 * np.sqrt(np.abs(omega2_l[indices].real))
690 omega_kl.append(omega_l.real)
692 # Conversion factor: sqrt(eV / Ang^2 / amu) -> eV
693 s = units._hbar * 1e10 / sqrt(units._e * units._amu)
694 omega_kl = s * np.asarray(omega_kl)
696 if modes:
697 return omega_kl, np.asarray(u_kl)
699 return omega_kl
701 def get_dos(self, kpts=(10, 10, 10), npts=1000, delta=1e-3, indices=None):
702 from ase.spectrum.dosdata import RawDOSData
703 # dos = self.dos(kpts, npts, delta, indices)
704 kpts_kc = monkhorst_pack(kpts)
705 omega_w = self.band_structure(kpts_kc).ravel()
706 dos = RawDOSData(omega_w, np.ones_like(omega_w))
707 return dos
709 def dos(self, kpts=(10, 10, 10), npts=1000, delta=1e-3, indices=None):
710 """Calculate phonon dos as a function of energy.
712 Parameters:
714 qpts: tuple
715 Shape of Monkhorst-Pack grid for sampling the Brillouin zone.
716 npts: int
717 Number of energy points.
718 delta: float
719 Broadening of Lorentzian line-shape in eV.
720 indices: list
721 If indices is not None, the atomic-partial dos for the specified
722 atoms will be calculated.
724 """
726 # Monkhorst-Pack grid
727 kpts_kc = monkhorst_pack(kpts)
728 N = np.prod(kpts)
729 # Get frequencies
730 omega_kl = self.band_structure(kpts_kc)
731 # Energy axis and dos
732 omega_e = np.linspace(0., np.amax(omega_kl) + 5e-3, num=npts)
733 dos_e = np.zeros_like(omega_e)
735 # Sum up contribution from all q-points and branches
736 for omega_l in omega_kl:
737 diff_el = (omega_e[:, np.newaxis] - omega_l[np.newaxis, :])**2
738 dos_el = 1. / (diff_el + (0.5 * delta)**2)
739 dos_e += dos_el.sum(axis=1)
741 dos_e *= 1. / (N * pi) * 0.5 * delta
743 return omega_e, dos_e
745 def write_modes(self, q_c, branches=0, kT=units.kB * 300, born=False,
746 repeat=(1, 1, 1), nimages=30, center=False):
747 """Write modes to trajectory file.
749 Parameters:
751 q_c: ndarray
752 q-vector of the modes.
753 branches: int or list
754 Branch index of modes.
755 kT: float
756 Temperature in units of eV. Determines the amplitude of the atomic
757 displacements in the modes.
758 born: bool
759 Include non-analytic contribution to the force constants at q -> 0.
760 repeat: tuple
761 Repeat atoms (l, m, n) times in the directions of the lattice
762 vectors. Displacements of atoms in repeated cells carry a Bloch
763 phase factor given by the q-vector and the cell lattice vector R_m.
764 nimages: int
765 Number of images in an oscillation.
766 center: bool
767 Center atoms in unit cell if True (default: False).
769 """
771 if isinstance(branches, int):
772 branch_l = [branches]
773 else:
774 branch_l = list(branches)
776 # Calculate modes
777 omega_l, u_l = self.band_structure([q_c], modes=True, born=born)
778 # Repeat atoms
779 atoms = self.atoms * repeat
780 # Center
781 if center:
782 atoms.center()
784 # Here ``Na`` refers to a composite unit cell/atom dimension
785 pos_Nav = atoms.get_positions()
786 # Total number of unit cells
787 N = np.prod(repeat)
789 # Corresponding lattice vectors R_m
790 R_cN = np.indices(repeat).reshape(3, -1)
791 # Bloch phase
792 phase_N = np.exp(2.j * pi * np.dot(q_c, R_cN))
793 phase_Na = phase_N.repeat(len(self.atoms))
795 for lval in branch_l:
797 omega = omega_l[0, lval]
798 u_av = u_l[0, lval]
800 # Mean displacement of a classical oscillator at temperature T
801 u_av *= sqrt(kT) / abs(omega)
803 mode_av = np.zeros((len(self.atoms), 3), dtype=complex)
804 # Insert slice with atomic displacements for the included atoms
805 mode_av[self.indices] = u_av
806 # Repeat and multiply by Bloch phase factor
807 mode_Nav = np.vstack(N * [mode_av]) * phase_Na[:, np.newaxis]
809 with Trajectory('%s.mode.%d.traj'
810 % (self.name, lval), 'w') as traj:
811 for x in np.linspace(0, 2 * pi, nimages, endpoint=False):
812 atoms.set_positions((pos_Nav + np.exp(1.j * x) *
813 mode_Nav).real)
814 traj.write(atoms)