Coverage for /builds/ase/ase/ase/phonons.py : 75.54%
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, 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 # type: ignore
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
290 >>> atoms = bulk('Si', 'diamond', a=5.4)
291 >>> calc = GPAW(kpts=(5, 5, 5),
292 h=0.2,
293 occupations=FermiDirac(0.))
294 >>> ph = Phonons(atoms, calc, supercell=(5, 5, 5))
295 >>> ph.run()
296 >>> ph.read(method='frederiksen', acoustic=True)
298 """
300 def __init__(self, *args, **kwargs):
301 """Initialize with base class args and kwargs."""
303 if 'name' not in kwargs:
304 kwargs['name'] = "phonon"
306 self.deprecate_refcell(kwargs)
308 Displacement.__init__(self, *args, **kwargs)
310 # Attributes for force constants and dynamical matrix in real space
311 self.C_N = None # in units of eV / Ang**2
312 self.D_N = None # in units of eV / Ang**2 / amu
314 # Attributes for born charges and static dielectric tensor
315 self.Z_avv = None
316 self.eps_vv = None
318 @staticmethod
319 def deprecate_refcell(kwargs: dict):
320 if 'refcell' in kwargs:
321 warnings.warn('Keyword refcell of Phonons is deprecated.'
322 'Please use center_refcell (bool)', FutureWarning)
323 kwargs['center_refcell'] = bool(kwargs['refcell'])
324 kwargs.pop('refcell')
326 return kwargs
328 def __call__(self, atoms_N):
329 """Calculate forces on atoms in supercell."""
330 return atoms_N.get_forces()
332 def calculate(self, atoms_N, disp):
333 forces = self(atoms_N)
334 return {'forces': forces}
336 def check_eq_forces(self):
337 """Check maximum size of forces in the equilibrium structure."""
339 eq_disp = self._eq_disp()
340 feq_av = self.cache[eq_disp.name]['forces']
342 fmin = feq_av.min()
343 fmax = feq_av.max()
344 i_min = np.where(feq_av == fmin)
345 i_max = np.where(feq_av == fmax)
347 return fmin, fmax, i_min, i_max
349 @deprecated('Current implementation of non-analytical correction is '
350 'likely incorrect, see '
351 'https://gitlab.com/ase/ase/-/issues/941')
352 def read_born_charges(self, name='born', neutrality=True):
353 r"""Read Born charges and dieletric tensor from JSON file.
355 The charge neutrality sum-rule::
357 _ _
358 \ a
359 ) Z = 0
360 /__ ij
361 a
363 Parameters:
365 neutrality: bool
366 Restore charge neutrality condition on calculated Born effective
367 charges.
368 name: str
369 Key used to identify the file with Born charges for the unit cell
370 in the JSON cache.
372 """
374 # Load file with Born charges and dielectric tensor for atoms in the
375 # unit cell
376 Z_avv, eps_vv = self.cache[name]
378 # Neutrality sum-rule
379 if neutrality:
380 Z_mean = Z_avv.sum(0) / len(Z_avv)
381 Z_avv -= Z_mean
383 self.Z_avv = Z_avv[self.indices]
384 self.eps_vv = eps_vv
386 def read(self, method='Frederiksen', symmetrize=3, acoustic=True,
387 cutoff=None, born=False, **kwargs):
388 """Read forces from json files and calculate force constants.
390 Extra keyword arguments will be passed to ``read_born_charges``.
392 Parameters:
394 method: str
395 Specify method for evaluating the atomic forces.
396 symmetrize: int
397 Symmetrize force constants (see doc string at top) when
398 ``symmetrize != 0`` (default: 3). Since restoring the acoustic sum
399 rule breaks the symmetry, the symmetrization must be repeated a few
400 times until the changes a insignificant. The integer gives the
401 number of iterations that will be carried out.
402 acoustic: bool
403 Restore the acoustic sum rule on the force constants.
404 cutoff: None or float
405 Zero elements in the dynamical matrix between atoms with an
406 interatomic distance larger than the cutoff.
407 born: bool
408 Read in Born effective charge tensor and high-frequency static
409 dielelctric tensor from file.
411 """
413 method = method.lower()
414 assert method in ['standard', 'frederiksen']
415 if cutoff is not None:
416 cutoff = float(cutoff)
418 # Read Born effective charges and optical dielectric tensor
419 if born:
420 self.read_born_charges(**kwargs)
422 # Number of atoms
423 natoms = len(self.indices)
424 # Number of unit cells
425 N = np.prod(self.supercell)
426 # Matrix of force constants as a function of unit cell index in units
427 # of eV / Ang**2
428 C_xNav = np.empty((natoms * 3, N, natoms, 3), dtype=float)
430 # Loop over all atomic displacements and calculate force constants
431 for i, a in enumerate(self.indices):
432 for j, v in enumerate('xyz'):
433 # Atomic forces for a displacement of atom a in direction v
434 # basename = '%s.%d%s' % (self.name, a, v)
435 basename = '%d%s' % (a, v)
436 fminus_av = self.cache[basename + '-']['forces']
437 fplus_av = self.cache[basename + '+']['forces']
439 if method == 'frederiksen':
440 fminus_av[a] -= fminus_av.sum(0)
441 fplus_av[a] -= fplus_av.sum(0)
443 # Finite difference derivative
444 C_av = fminus_av - fplus_av
445 C_av /= 2 * self.delta
447 # Slice out included atoms
448 C_Nav = C_av.reshape((N, len(self.atoms), 3))[:, self.indices]
449 index = 3 * i + j
450 C_xNav[index] = C_Nav
452 # Make unitcell index the first and reshape
453 C_N = C_xNav.swapaxes(0, 1).reshape((N,) + (3 * natoms, 3 * natoms))
455 # Cut off before symmetry and acoustic sum rule are imposed
456 if cutoff is not None:
457 self.apply_cutoff(C_N, cutoff)
459 # Symmetrize force constants
460 if symmetrize:
461 for i in range(symmetrize):
462 # Symmetrize
463 C_N = self.symmetrize(C_N)
464 # Restore acoustic sum-rule
465 if acoustic:
466 self.acoustic(C_N)
467 else:
468 break
470 # Store force constants and dynamical matrix
471 self.C_N = C_N
472 self.D_N = C_N.copy()
474 # Add mass prefactor
475 m_a = self.atoms.get_masses()
476 self.m_inv_x = np.repeat(m_a[self.indices]**-0.5, 3)
477 M_inv = np.outer(self.m_inv_x, self.m_inv_x)
478 for D in self.D_N:
479 D *= M_inv
481 def symmetrize(self, C_N):
482 """Symmetrize force constant matrix."""
484 # Number of atoms
485 natoms = len(self.indices)
486 # Number of unit cells
487 N = np.prod(self.supercell)
489 # Reshape force constants to (l, m, n) cell indices
490 C_lmn = C_N.reshape(self.supercell + (3 * natoms, 3 * natoms))
492 # Shift reference cell to center index
493 if self.offset == 0:
494 C_lmn = fft.fftshift(C_lmn, axes=(0, 1, 2)).copy()
495 # Make force constants symmetric in indices -- in case of an even
496 # number of unit cells don't include the first cell
497 i, j, k = 1 - np.asarray(self.supercell) % 2
498 C_lmn[i:, j:, k:] *= 0.5
499 C_lmn[i:, j:, k:] += \
500 C_lmn[i:, j:, k:][::-1, ::-1, ::-1].transpose(0, 1, 2, 4, 3).copy()
501 if self.offset == 0:
502 C_lmn = fft.ifftshift(C_lmn, axes=(0, 1, 2)).copy()
504 # Change to single unit cell index shape
505 C_N = C_lmn.reshape((N, 3 * natoms, 3 * natoms))
507 return C_N
509 def acoustic(self, C_N):
510 """Restore acoustic sumrule on force constants."""
512 # Number of atoms
513 natoms = len(self.indices)
514 # Copy force constants
515 C_N_temp = C_N.copy()
517 # Correct atomic diagonals of R_m = (0, 0, 0) matrix
518 for C in C_N_temp:
519 for a in range(natoms):
520 for a_ in range(natoms):
521 C_N[self.offset,
522 3 * a: 3 * a + 3,
523 3 * a: 3 * a + 3] -= C[3 * a: 3 * a + 3,
524 3 * a_: 3 * a_ + 3]
526 def apply_cutoff(self, D_N, r_c):
527 """Zero elements for interatomic distances larger than the cutoff.
529 Parameters:
531 D_N: ndarray
532 Dynamical/force constant matrix.
533 r_c: float
534 Cutoff in Angstrom.
536 """
538 # Number of atoms and primitive cells
539 natoms = len(self.indices)
540 N = np.prod(self.supercell)
541 # Lattice vectors
542 R_cN = self._lattice_vectors_array
543 # Reshape matrix to individual atomic and cartesian dimensions
544 D_Navav = D_N.reshape((N, natoms, 3, natoms, 3))
546 # Cell vectors
547 cell_vc = self.atoms.cell.transpose()
548 # Atomic positions in reference cell
549 pos_av = self.atoms.get_positions()
551 # Zero elements with a distance to atoms in the reference cell
552 # larger than the cutoff
553 for n in range(N):
554 # Lattice vector to cell
555 R_v = np.dot(cell_vc, R_cN[:, n])
556 # Atomic positions in cell
557 posn_av = pos_av + R_v
558 # Loop over atoms and zero elements
559 for i, a in enumerate(self.indices):
560 dist_a = np.sqrt(np.sum((pos_av[a] - posn_av)**2, axis=-1))
561 # Atoms where the distance is larger than the cufoff
562 i_a = dist_a > r_c # np.where(dist_a > r_c)
563 # Zero elements
564 D_Navav[n, i, :, i_a, :] = 0.0
566 def get_force_constant(self):
567 """Return matrix of force constants."""
569 assert self.C_N is not None
570 return self.C_N
572 def get_band_structure(self, path, modes=False, born=False, verbose=True):
573 omega_kl = self.band_structure(path.kpts, modes, born, verbose)
574 if modes:
575 assert 0
576 omega_kl, modes = omega_kl
578 from ase.spectrum.band_structure import BandStructure
579 bs = BandStructure(path, energies=omega_kl[None])
580 return bs
582 def compute_dynamical_matrix(self, q_scaled: np.ndarray, D_N: np.ndarray):
583 """ Computation of the dynamical matrix in momentum space D_ab(q).
584 This is a Fourier transform from real-space dynamical matrix D_N
585 for a given momentum vector q.
587 q_scaled: q vector in scaled coordinates.
589 D_N: the dynamical matrix in real-space. It is necessary, at least
590 currently, to provide this matrix explicitly (rather than use
591 self.D_N) because this matrix is modified by the Born charges
592 contributions and these modifications are momentum (q) dependent.
594 Result:
595 D(q): two-dimensional, complex-valued array of
596 shape=(3 * natoms, 3 * natoms).
597 """
598 # Evaluate fourier sum
599 R_cN = self._lattice_vectors_array
600 phase_N = np.exp(-2.j * pi * np.dot(q_scaled, R_cN))
601 D_q = np.sum(phase_N[:, np.newaxis, np.newaxis] * D_N, axis=0)
602 return D_q
604 def band_structure(self, path_kc, modes=False, born=False, verbose=True):
605 """Calculate phonon dispersion along a path in the Brillouin zone.
607 The dynamical matrix at arbitrary q-vectors is obtained by Fourier
608 transforming the real-space force constants. In case of negative
609 eigenvalues (squared frequency), the corresponding negative frequency
610 is returned.
612 Frequencies and modes are in units of eV and Ang/sqrt(amu),
613 respectively.
615 Parameters:
617 path_kc: ndarray
618 List of k-point coordinates (in units of the reciprocal lattice
619 vectors) specifying the path in the Brillouin zone for which the
620 dynamical matrix will be calculated.
621 modes: bool
622 Returns both frequencies and modes when True.
623 born: bool
624 Include non-analytic part given by the Born effective charges and
625 the static part of the high-frequency dielectric tensor. This
626 contribution to the force constant accounts for the splitting
627 between the LO and TO branches for q -> 0.
628 verbose: bool
629 Print warnings when imaginary frequncies are detected.
631 """
633 assert self.D_N is not None
634 if born:
635 assert self.Z_avv is not None
636 assert self.eps_vv is not None
638 # Dynamical matrix in real-space
639 D_N = self.D_N
641 # Lists for frequencies and modes along path
642 omega_kl = []
643 u_kl = []
645 # Reciprocal basis vectors for use in non-analytic contribution
646 reci_vc = 2 * pi * la.inv(self.atoms.cell)
647 # Unit cell volume in Bohr^3
648 vol = abs(la.det(self.atoms.cell)) / units.Bohr**3
650 for q_c in path_kc:
652 # Add non-analytic part
653 if born:
654 # q-vector in cartesian coordinates
655 q_v = np.dot(reci_vc, q_c)
656 # Non-analytic contribution to force constants in atomic units
657 qdotZ_av = np.dot(q_v, self.Z_avv).ravel()
658 C_na = (4 * pi * np.outer(qdotZ_av, qdotZ_av) /
659 np.dot(q_v, np.dot(self.eps_vv, q_v)) / vol)
660 self.C_na = C_na / units.Bohr**2 * units.Hartree
661 # Add mass prefactor and convert to eV / (Ang^2 * amu)
662 M_inv = np.outer(self.m_inv_x, self.m_inv_x)
663 D_na = C_na * M_inv / units.Bohr**2 * units.Hartree
664 self.D_na = D_na
665 D_N = self.D_N + D_na / np.prod(self.supercell)
667 # if np.prod(self.N_c) == 1:
668 #
669 # q_av = np.tile(q_v, len(self.indices))
670 # q_xx = np.vstack([q_av]*len(self.indices)*3)
671 # D_m += q_xx
673 # Evaluate fourier sum
674 D_q = self.compute_dynamical_matrix(q_c, D_N)
676 if modes:
677 omega2_l, u_xl = la.eigh(D_q, UPLO='U')
678 # Sort eigenmodes according to eigenvalues (see below) and
679 # multiply with mass prefactor
680 u_lx = (self.m_inv_x[:, np.newaxis] *
681 u_xl[:, omega2_l.argsort()]).T.copy()
682 u_kl.append(u_lx.reshape((-1, len(self.indices), 3)))
683 else:
684 omega2_l = la.eigvalsh(D_q, UPLO='U')
686 # Sort eigenvalues in increasing order
687 omega2_l.sort()
688 # Use dtype=complex to handle negative eigenvalues
689 omega_l = np.sqrt(omega2_l.astype(complex))
691 # Take care of imaginary frequencies
692 if not np.all(omega2_l >= 0.):
693 indices = np.where(omega2_l < 0)[0]
695 if verbose:
696 print('WARNING, %i imaginary frequencies at '
697 'q = (% 5.2f, % 5.2f, % 5.2f) ; (omega_q =% 5.3e*i)'
698 % (len(indices), q_c[0], q_c[1], q_c[2],
699 omega_l[indices][0].imag))
701 omega_l[indices] = -1 * np.sqrt(np.abs(omega2_l[indices].real))
703 omega_kl.append(omega_l.real)
705 # Conversion factor: sqrt(eV / Ang^2 / amu) -> eV
706 s = units._hbar * 1e10 / sqrt(units._e * units._amu)
707 omega_kl = s * np.asarray(omega_kl)
709 if modes:
710 return omega_kl, np.asarray(u_kl)
712 return omega_kl
714 def get_dos(self, kpts=(10, 10, 10), npts=1000, delta=1e-3, indices=None):
715 from ase.spectrum.dosdata import RawDOSData
716 # dos = self.dos(kpts, npts, delta, indices)
717 kpts_kc = monkhorst_pack(kpts)
718 omega_w = self.band_structure(kpts_kc).ravel()
719 dos = RawDOSData(omega_w, np.ones_like(omega_w))
720 return dos
722 def dos(self, kpts=(10, 10, 10), npts=1000, delta=1e-3, indices=None):
723 """Calculate phonon dos as a function of energy.
725 Parameters:
727 qpts: tuple
728 Shape of Monkhorst-Pack grid for sampling the Brillouin zone.
729 npts: int
730 Number of energy points.
731 delta: float
732 Broadening of Lorentzian line-shape in eV.
733 indices: list
734 If indices is not None, the atomic-partial dos for the specified
735 atoms will be calculated.
737 """
739 # Monkhorst-Pack grid
740 kpts_kc = monkhorst_pack(kpts)
741 N = np.prod(kpts)
742 # Get frequencies
743 omega_kl = self.band_structure(kpts_kc)
744 # Energy axis and dos
745 omega_e = np.linspace(0., np.amax(omega_kl) + 5e-3, num=npts)
746 dos_e = np.zeros_like(omega_e)
748 # Sum up contribution from all q-points and branches
749 for omega_l in omega_kl:
750 diff_el = (omega_e[:, np.newaxis] - omega_l[np.newaxis, :])**2
751 dos_el = 1. / (diff_el + (0.5 * delta)**2)
752 dos_e += dos_el.sum(axis=1)
754 dos_e *= 1. / (N * pi) * 0.5 * delta
756 return omega_e, dos_e
758 def write_modes(self, q_c, branches=0, kT=units.kB * 300, born=False,
759 repeat=(1, 1, 1), nimages=30, center=False):
760 """Write modes to trajectory file.
762 Parameters:
764 q_c: ndarray
765 q-vector of the modes.
766 branches: int or list
767 Branch index of modes.
768 kT: float
769 Temperature in units of eV. Determines the amplitude of the atomic
770 displacements in the modes.
771 born: bool
772 Include non-analytic contribution to the force constants at q -> 0.
773 repeat: tuple
774 Repeat atoms (l, m, n) times in the directions of the lattice
775 vectors. Displacements of atoms in repeated cells carry a Bloch
776 phase factor given by the q-vector and the cell lattice vector R_m.
777 nimages: int
778 Number of images in an oscillation.
779 center: bool
780 Center atoms in unit cell if True (default: False).
782 """
784 if isinstance(branches, int):
785 branch_l = [branches]
786 else:
787 branch_l = list(branches)
789 # Calculate modes
790 omega_l, u_l = self.band_structure([q_c], modes=True, born=born)
791 # Repeat atoms
792 atoms = self.atoms * repeat
793 # Center
794 if center:
795 atoms.center()
797 # Here ``Na`` refers to a composite unit cell/atom dimension
798 pos_Nav = atoms.get_positions()
799 # Total number of unit cells
800 N = np.prod(repeat)
802 # Corresponding lattice vectors R_m
803 R_cN = np.indices(repeat).reshape(3, -1)
804 # Bloch phase
805 phase_N = np.exp(2.j * pi * np.dot(q_c, R_cN))
806 phase_Na = phase_N.repeat(len(self.atoms))
808 for lval in branch_l:
810 omega = omega_l[0, lval]
811 u_av = u_l[0, lval]
813 # Mean displacement of a classical oscillator at temperature T
814 u_av *= sqrt(kT) / abs(omega)
816 mode_av = np.zeros((len(self.atoms), 3), dtype=complex)
817 # Insert slice with atomic displacements for the included atoms
818 mode_av[self.indices] = u_av
819 # Repeat and multiply by Bloch phase factor
820 mode_Nav = np.vstack(N * [mode_av]) * phase_Na[:, np.newaxis]
822 with Trajectory('%s.mode.%d.traj'
823 % (self.name, lval), 'w') as traj:
824 for x in np.linspace(0, 2 * pi, nimages, endpoint=False):
825 atoms.set_positions((pos_Nav + np.exp(1.j * x) *
826 mode_Nav).real)
827 traj.write(atoms)