Coverage for /builds/kinetik161/ase/ase/phonons.py: 75.78%
322 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +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 i 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 omega_kl = self.band_structure(path.kpts, modes, born, verbose)
579 if modes:
580 assert 0
581 omega_kl, modes = omega_kl
583 from ase.spectrum.band_structure import BandStructure
584 bs = BandStructure(path, energies=omega_kl[None])
585 return bs
587 def compute_dynamical_matrix(self, q_scaled: np.ndarray, D_N: np.ndarray):
588 """ Computation of the dynamical matrix in momentum space D_ab(q).
589 This is a Fourier transform from real-space dynamical matrix D_N
590 for a given momentum vector q.
592 q_scaled: q vector in scaled coordinates.
594 D_N: the dynamical matrix in real-space. It is necessary, at least
595 currently, to provide this matrix explicitly (rather than use
596 self.D_N) because this matrix is modified by the Born charges
597 contributions and these modifications are momentum (q) dependent.
599 Result:
600 D(q): two-dimensional, complex-valued array of
601 shape=(3 * natoms, 3 * natoms).
602 """
603 # Evaluate fourier sum
604 R_cN = self._lattice_vectors_array
605 phase_N = np.exp(-2.j * pi * np.dot(q_scaled, R_cN))
606 D_q = np.sum(phase_N[:, np.newaxis, np.newaxis] * D_N, axis=0)
607 return D_q
609 def band_structure(self, path_kc, modes=False, born=False, verbose=True):
610 """Calculate phonon dispersion along a path in the Brillouin zone.
612 The dynamical matrix at arbitrary q-vectors is obtained by Fourier
613 transforming the real-space force constants. In case of negative
614 eigenvalues (squared frequency), the corresponding negative frequency
615 is returned.
617 Frequencies and modes are in units of eV and Ang/sqrt(amu),
618 respectively.
620 Parameters:
622 path_kc: ndarray
623 List of k-point coordinates (in units of the reciprocal lattice
624 vectors) specifying the path in the Brillouin zone for which the
625 dynamical matrix will be calculated.
626 modes: bool
627 Returns both frequencies and modes when True.
628 born: bool
629 Include non-analytic part given by the Born effective charges and
630 the static part of the high-frequency dielectric tensor. This
631 contribution to the force constant accounts for the splitting
632 between the LO and TO branches for q -> 0.
633 verbose: bool
634 Print warnings when imaginary frequncies are detected.
636 """
638 assert self.D_N is not None
639 if born:
640 assert self.Z_avv is not None
641 assert self.eps_vv is not None
643 # Dynamical matrix in real-space
644 D_N = self.D_N
646 # Lists for frequencies and modes along path
647 omega_kl = []
648 u_kl = []
650 # Reciprocal basis vectors for use in non-analytic contribution
651 reci_vc = 2 * pi * la.inv(self.atoms.cell)
652 # Unit cell volume in Bohr^3
653 vol = abs(la.det(self.atoms.cell)) / units.Bohr**3
655 for q_c in path_kc:
657 # Add non-analytic part
658 if born:
659 # q-vector in cartesian coordinates
660 q_v = np.dot(reci_vc, q_c)
661 # Non-analytic contribution to force constants in atomic units
662 qdotZ_av = np.dot(q_v, self.Z_avv).ravel()
663 C_na = (4 * pi * np.outer(qdotZ_av, qdotZ_av) /
664 np.dot(q_v, np.dot(self.eps_vv, q_v)) / vol)
665 self.C_na = C_na / units.Bohr**2 * units.Hartree
666 # Add mass prefactor and convert to eV / (Ang^2 * amu)
667 M_inv = np.outer(self.m_inv_x, self.m_inv_x)
668 D_na = C_na * M_inv / units.Bohr**2 * units.Hartree
669 self.D_na = D_na
670 D_N = self.D_N + D_na / np.prod(self.supercell)
672 # if np.prod(self.N_c) == 1:
673 #
674 # q_av = np.tile(q_v, len(self.indices))
675 # q_xx = np.vstack([q_av]*len(self.indices)*3)
676 # D_m += q_xx
678 # Evaluate fourier sum
679 D_q = self.compute_dynamical_matrix(q_c, D_N)
681 if modes:
682 omega2_l, u_xl = la.eigh(D_q, UPLO='U')
683 # Sort eigenmodes according to eigenvalues (see below) and
684 # multiply with mass prefactor
685 u_lx = (self.m_inv_x[:, np.newaxis] *
686 u_xl[:, omega2_l.argsort()]).T.copy()
687 u_kl.append(u_lx.reshape((-1, len(self.indices), 3)))
688 else:
689 omega2_l = la.eigvalsh(D_q, UPLO='U')
691 # Sort eigenvalues in increasing order
692 omega2_l.sort()
693 # Use dtype=complex to handle negative eigenvalues
694 omega_l = np.sqrt(omega2_l.astype(complex))
696 # Take care of imaginary frequencies
697 if not np.all(omega2_l >= 0.):
698 indices = np.where(omega2_l < 0)[0]
700 if verbose:
701 print('WARNING, %i imaginary frequencies at '
702 'q = (% 5.2f, % 5.2f, % 5.2f) ; (omega_q =% 5.3e*i)'
703 % (len(indices), q_c[0], q_c[1], q_c[2],
704 omega_l[indices][0].imag))
706 omega_l[indices] = -1 * np.sqrt(np.abs(omega2_l[indices].real))
708 omega_kl.append(omega_l.real)
710 # Conversion factor: sqrt(eV / Ang^2 / amu) -> eV
711 s = units._hbar * 1e10 / sqrt(units._e * units._amu)
712 omega_kl = s * np.asarray(omega_kl)
714 if modes:
715 return omega_kl, np.asarray(u_kl)
717 return omega_kl
719 def get_dos(self, kpts=(10, 10, 10), npts=1000, delta=1e-3, indices=None):
720 from ase.spectrum.dosdata import RawDOSData
722 # dos = self.dos(kpts, npts, delta, indices)
723 kpts_kc = monkhorst_pack(kpts)
724 omega_w = self.band_structure(kpts_kc).ravel()
725 dos = RawDOSData(omega_w, np.ones_like(omega_w))
726 return dos
728 def dos(self, kpts=(10, 10, 10), npts=1000, delta=1e-3, indices=None):
729 """Calculate phonon dos as a function of energy.
731 Parameters:
733 qpts: tuple
734 Shape of Monkhorst-Pack grid for sampling the Brillouin zone.
735 npts: int
736 Number of energy points.
737 delta: float
738 Broadening of Lorentzian line-shape in eV.
739 indices: list
740 If indices is not None, the atomic-partial dos for the specified
741 atoms will be calculated.
743 """
745 # Monkhorst-Pack grid
746 kpts_kc = monkhorst_pack(kpts)
747 N = np.prod(kpts)
748 # Get frequencies
749 omega_kl = self.band_structure(kpts_kc)
750 # Energy axis and dos
751 omega_e = np.linspace(0., np.amax(omega_kl) + 5e-3, num=npts)
752 dos_e = np.zeros_like(omega_e)
754 # Sum up contribution from all q-points and branches
755 for omega_l in omega_kl:
756 diff_el = (omega_e[:, np.newaxis] - omega_l[np.newaxis, :])**2
757 dos_el = 1. / (diff_el + (0.5 * delta)**2)
758 dos_e += dos_el.sum(axis=1)
760 dos_e *= 1. / (N * pi) * 0.5 * delta
762 return omega_e, dos_e
764 def write_modes(self, q_c, branches=0, kT=units.kB * 300, born=False,
765 repeat=(1, 1, 1), nimages=30, center=False):
766 """Write modes to trajectory file.
768 Parameters:
770 q_c: ndarray
771 q-vector of the modes.
772 branches: int or list
773 Branch index of modes.
774 kT: float
775 Temperature in units of eV. Determines the amplitude of the atomic
776 displacements in the modes.
777 born: bool
778 Include non-analytic contribution to the force constants at q -> 0.
779 repeat: tuple
780 Repeat atoms (l, m, n) times in the directions of the lattice
781 vectors. Displacements of atoms in repeated cells carry a Bloch
782 phase factor given by the q-vector and the cell lattice vector R_m.
783 nimages: int
784 Number of images in an oscillation.
785 center: bool
786 Center atoms in unit cell if True (default: False).
788 """
790 if isinstance(branches, int):
791 branch_l = [branches]
792 else:
793 branch_l = list(branches)
795 # Calculate modes
796 omega_l, u_l = self.band_structure([q_c], modes=True, born=born)
797 # Repeat atoms
798 atoms = self.atoms * repeat
799 # Center
800 if center:
801 atoms.center()
803 # Here ``Na`` refers to a composite unit cell/atom dimension
804 pos_Nav = atoms.get_positions()
805 # Total number of unit cells
806 N = np.prod(repeat)
808 # Corresponding lattice vectors R_m
809 R_cN = np.indices(repeat).reshape(3, -1)
810 # Bloch phase
811 phase_N = np.exp(2.j * pi * np.dot(q_c, R_cN))
812 phase_Na = phase_N.repeat(len(self.atoms))
814 for lval in branch_l:
816 omega = omega_l[0, lval]
817 u_av = u_l[0, lval]
819 # Mean displacement of a classical oscillator at temperature T
820 u_av *= sqrt(kT) / abs(omega)
822 mode_av = np.zeros((len(self.atoms), 3), dtype=complex)
823 # Insert slice with atomic displacements for the included atoms
824 mode_av[self.indices] = u_av
825 # Repeat and multiply by Bloch phase factor
826 mode_Nav = np.vstack(N * [mode_av]) * phase_Na[:, np.newaxis]
828 with Trajectory('%s.mode.%d.traj'
829 % (self.name, lval), 'w') as traj:
830 for x in np.linspace(0, 2 * pi, nimages, endpoint=False):
831 atoms.set_positions((pos_Nav + np.exp(1.j * x) *
832 mode_Nav).real)
833 traj.write(atoms)