Coverage for /builds/debichem-team/python-ase/ase/phonons.py : 68.25%
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
19class Displacement:
20 """Abstract base class for phonon and el-ph supercell calculations.
22 Both phonons and the electron-phonon interaction in periodic systems can be
23 calculated with the so-called finite-displacement method where the
24 derivatives of the total energy and effective potential are obtained from
25 finite-difference approximations, i.e. by displacing the atoms. This class
26 provides the required functionality for carrying out the calculations for
27 the different displacements in its ``run`` member function.
29 Derived classes must overwrite the ``__call__`` member function which is
30 called for each atomic displacement.
32 """
34 def __init__(self, atoms, calc=None, supercell=(1, 1, 1), name=None,
35 delta=0.01, center_refcell=False):
36 """Init with an instance of class ``Atoms`` and a calculator.
38 Parameters:
40 atoms: Atoms object
41 The atoms to work on.
42 calc: Calculator
43 Calculator for the supercell calculation.
44 supercell: tuple
45 Size of supercell given by the number of repetitions (l, m, n) of
46 the small unit cell in each direction.
47 name: str
48 Base name to use for files.
49 delta: float
50 Magnitude of displacement in Ang.
51 center_refcell: bool
52 Reference cell in which the atoms will be displaced. If False, then
53 corner cell in supercell is used. If True, then cell in the center
54 of the supercell is used.
56 """
58 # Store atoms and calculator
59 self.atoms = atoms
60 self.calc = calc
62 # Displace all atoms in the unit cell by default
63 self.indices = np.arange(len(atoms))
64 self.name = name
65 self.delta = delta
66 self.center_refcell = center_refcell
67 self.supercell = supercell
69 self.cache = MultiFileJSONCache(self.name)
71 def define_offset(self): # Reference cell offset
73 if not self.center_refcell:
74 # Corner cell
75 self.offset = 0
76 else:
77 # Center cell
78 N_c = self.supercell
79 self.offset = (N_c[0] // 2 * (N_c[1] * N_c[2]) +
80 N_c[1] // 2 * N_c[2] +
81 N_c[2] // 2)
82 return self.offset
84 @property # type: ignore
85 @ase.utils.deprecated('Please use phonons.supercell instead of .N_c')
86 def N_c(self):
87 return self._supercell
89 @property
90 def supercell(self):
91 return self._supercell
93 @supercell.setter
94 def supercell(self, supercell):
95 assert len(supercell) == 3
96 self._supercell = tuple(supercell)
97 self.define_offset()
98 self._lattice_vectors_array = self.compute_lattice_vectors()
100 @ase.utils.deprecated('Please use phonons.compute_lattice_vectors()'
101 ' instead of .lattice_vectors()')
102 def lattice_vectors(self):
103 return self.compute_lattice_vectors()
105 def compute_lattice_vectors(self):
106 """Return lattice vectors for cells in the supercell."""
107 # Lattice vectors -- ordered as illustrated in class docstring
109 # Lattice vectors relevative to the reference cell
110 R_cN = np.indices(self.supercell).reshape(3, -1)
111 N_c = np.array(self.supercell)[:, np.newaxis]
112 if self.offset == 0:
113 R_cN += N_c // 2
114 R_cN %= N_c
115 R_cN -= N_c // 2
116 return R_cN
118 def __call__(self, *args, **kwargs):
119 """Member function called in the ``run`` function."""
121 raise NotImplementedError("Implement in derived classes!.")
123 def set_atoms(self, atoms):
124 """Set the atoms to vibrate.
126 Parameters:
128 atoms: list
129 Can be either a list of strings, ints or ...
131 """
133 assert isinstance(atoms, list)
134 assert len(atoms) <= len(self.atoms)
136 if isinstance(atoms[0], str):
137 assert np.all([isinstance(atom, str) for atom in atoms])
138 sym_a = self.atoms.get_chemical_symbols()
139 # List for atomic indices
140 indices = []
141 for type in atoms:
142 indices.extend([a for a, atom in enumerate(sym_a)
143 if atom == type])
144 else:
145 assert np.all([isinstance(atom, int) for atom in atoms])
146 indices = atoms
148 self.indices = indices
150 def _disp(self, a, i, step):
151 from ase.vibrations.vibrations import Displacement as VDisplacement
152 return VDisplacement(a, i, np.sign(step), abs(step), self)
154 def run(self):
155 """Run the calculations for the required displacements.
157 This will do a calculation for 6 displacements per atom, +-x, +-y, and
158 +-z. Only those calculations that are not already done will be
159 started. Be aware that an interrupted calculation may produce an empty
160 file (ending with .json), which must be deleted before restarting the
161 job. Otherwise the calculation for that displacement will not be done.
163 """
165 # Atoms in the supercell -- repeated in the lattice vector directions
166 # beginning with the last
167 atoms_N = self.atoms * self.supercell
169 # Set calculator if provided
170 assert self.calc is not None, "Provide calculator in __init__ method"
171 atoms_N.calc = self.calc
173 # Do calculation on equilibrium structure
174 eq_disp = self._disp(0, 0, 0)
175 # with self.cache.lock(f'{self.name}.eq') as handle:
176 with self.cache.lock(eq_disp.name) as handle:
177 if handle is not None:
178 output = self(atoms_N)
179 # Write output to file
180 if world.rank == 0:
181 handle.save(output)
183 # Positions of atoms to be displaced in the reference cell
184 natoms = len(self.atoms)
185 offset = natoms * self.offset
186 pos = atoms_N.positions[offset: offset + natoms].copy()
188 # Loop over all displacements
189 for a in self.indices:
190 for i in range(3):
191 for sign in [-1, 1]:
192 disp = self._disp(a, i, sign)
193 # key = '%s.%d%s%s' % (self.name, a, 'xyz'[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 name = f'{self.name}.eq'
327 feq_av = self.cache[name]['forces']
329 fmin = feq_av.max()
330 fmax = feq_av.min()
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 def read_born_charges(self, name=None, neutrality=True):
337 r"""Read Born charges and dieletric tensor from JSON file.
339 The charge neutrality sum-rule::
341 _ _
342 \ a
343 ) Z = 0
344 /__ ij
345 a
347 Parameters:
349 neutrality: bool
350 Restore charge neutrality condition on calculated Born effective
351 charges.
353 """
355 # Load file with Born charges and dielectric tensor for atoms in the
356 # unit cell
357 if name is None:
358 key = '%s.born' % self.name
359 else:
360 key = name
362 Z_avv, eps_vv = self.cache[key]
364 # Neutrality sum-rule
365 if neutrality:
366 Z_mean = Z_avv.sum(0) / len(Z_avv)
367 Z_avv -= Z_mean
369 self.Z_avv = Z_avv[self.indices]
370 self.eps_vv = eps_vv
372 def read(self, method='Frederiksen', symmetrize=3, acoustic=True,
373 cutoff=None, born=False, **kwargs):
374 """Read forces from json files and calculate force constants.
376 Extra keyword arguments will be passed to ``read_born_charges``.
378 Parameters:
380 method: str
381 Specify method for evaluating the atomic forces.
382 symmetrize: int
383 Symmetrize force constants (see doc string at top) when
384 ``symmetrize != 0`` (default: 3). Since restoring the acoustic sum
385 rule breaks the symmetry, the symmetrization must be repeated a few
386 times until the changes a insignificant. The integer gives the
387 number of iterations that will be carried out.
388 acoustic: bool
389 Restore the acoustic sum rule on the force constants.
390 cutoff: None or float
391 Zero elements in the dynamical matrix between atoms with an
392 interatomic distance larger than the cutoff.
393 born: bool
394 Read in Born effective charge tensor and high-frequency static
395 dielelctric tensor from file.
397 """
399 method = method.lower()
400 assert method in ['standard', 'frederiksen']
401 if cutoff is not None:
402 cutoff = float(cutoff)
404 # Read Born effective charges and optical dielectric tensor
405 if born:
406 self.read_born_charges(**kwargs)
408 # Number of atoms
409 natoms = len(self.indices)
410 # Number of unit cells
411 N = np.prod(self.supercell)
412 # Matrix of force constants as a function of unit cell index in units
413 # of eV / Ang**2
414 C_xNav = np.empty((natoms * 3, N, natoms, 3), dtype=float)
416 # Loop over all atomic displacements and calculate force constants
417 for i, a in enumerate(self.indices):
418 for j, v in enumerate('xyz'):
419 # Atomic forces for a displacement of atom a in direction v
420 # basename = '%s.%d%s' % (self.name, a, v)
421 basename = '%d%s' % (a, v)
422 fminus_av = self.cache[basename + '-']['forces']
423 fplus_av = self.cache[basename + '+']['forces']
425 if method == 'frederiksen':
426 fminus_av[a] -= fminus_av.sum(0)
427 fplus_av[a] -= fplus_av.sum(0)
429 # Finite difference derivative
430 C_av = fminus_av - fplus_av
431 C_av /= 2 * self.delta
433 # Slice out included atoms
434 C_Nav = C_av.reshape((N, len(self.atoms), 3))[:, self.indices]
435 index = 3 * i + j
436 C_xNav[index] = C_Nav
438 # Make unitcell index the first and reshape
439 C_N = C_xNav.swapaxes(0, 1).reshape((N,) + (3 * natoms, 3 * natoms))
441 # Cut off before symmetry and acoustic sum rule are imposed
442 if cutoff is not None:
443 self.apply_cutoff(C_N, cutoff)
445 # Symmetrize force constants
446 if symmetrize:
447 for i in range(symmetrize):
448 # Symmetrize
449 C_N = self.symmetrize(C_N)
450 # Restore acoustic sum-rule
451 if acoustic:
452 self.acoustic(C_N)
453 else:
454 break
456 # Store force constants and dynamical matrix
457 self.C_N = C_N
458 self.D_N = C_N.copy()
460 # Add mass prefactor
461 m_a = self.atoms.get_masses()
462 self.m_inv_x = np.repeat(m_a[self.indices]**-0.5, 3)
463 M_inv = np.outer(self.m_inv_x, self.m_inv_x)
464 for D in self.D_N:
465 D *= M_inv
467 def symmetrize(self, C_N):
468 """Symmetrize force constant matrix."""
470 # Number of atoms
471 natoms = len(self.indices)
472 # Number of unit cells
473 N = np.prod(self.supercell)
475 # Reshape force constants to (l, m, n) cell indices
476 C_lmn = C_N.reshape(self.supercell + (3 * natoms, 3 * natoms))
478 # Shift reference cell to center index
479 if self.offset == 0:
480 C_lmn = fft.fftshift(C_lmn, axes=(0, 1, 2)).copy()
481 # Make force constants symmetric in indices -- in case of an even
482 # number of unit cells don't include the first cell
483 i, j, k = 1 - np.asarray(self.supercell) % 2
484 C_lmn[i:, j:, k:] *= 0.5
485 C_lmn[i:, j:, k:] += \
486 C_lmn[i:, j:, k:][::-1, ::-1, ::-1].transpose(0, 1, 2, 4, 3).copy()
487 if self.offset == 0:
488 C_lmn = fft.ifftshift(C_lmn, axes=(0, 1, 2)).copy()
490 # Change to single unit cell index shape
491 C_N = C_lmn.reshape((N, 3 * natoms, 3 * natoms))
493 return C_N
495 def acoustic(self, C_N):
496 """Restore acoustic sumrule on force constants."""
498 # Number of atoms
499 natoms = len(self.indices)
500 # Copy force constants
501 C_N_temp = C_N.copy()
503 # Correct atomic diagonals of R_m = (0, 0, 0) matrix
504 for C in C_N_temp:
505 for a in range(natoms):
506 for a_ in range(natoms):
507 C_N[self.offset,
508 3 * a: 3 * a + 3,
509 3 * a: 3 * a + 3] -= C[3 * a: 3 * a + 3,
510 3 * a_: 3 * a_ + 3]
512 def apply_cutoff(self, D_N, r_c):
513 """Zero elements for interatomic distances larger than the cutoff.
515 Parameters:
517 D_N: ndarray
518 Dynamical/force constant matrix.
519 r_c: float
520 Cutoff in Angstrom.
522 """
524 # Number of atoms and primitive cells
525 natoms = len(self.indices)
526 N = np.prod(self.supercell)
527 # Lattice vectors
528 R_cN = self._lattice_vectors_array
529 # Reshape matrix to individual atomic and cartesian dimensions
530 D_Navav = D_N.reshape((N, natoms, 3, natoms, 3))
532 # Cell vectors
533 cell_vc = self.atoms.cell.transpose()
534 # Atomic positions in reference cell
535 pos_av = self.atoms.get_positions()
537 # Zero elements with a distance to atoms in the reference cell
538 # larger than the cutoff
539 for n in range(N):
540 # Lattice vector to cell
541 R_v = np.dot(cell_vc, R_cN[:, n])
542 # Atomic positions in cell
543 posn_av = pos_av + R_v
544 # Loop over atoms and zero elements
545 for i, a in enumerate(self.indices):
546 dist_a = np.sqrt(np.sum((pos_av[a] - posn_av)**2, axis=-1))
547 # Atoms where the distance is larger than the cufoff
548 i_a = dist_a > r_c # np.where(dist_a > r_c)
549 # Zero elements
550 D_Navav[n, i, :, i_a, :] = 0.0
552 def get_force_constant(self):
553 """Return matrix of force constants."""
555 assert self.C_N is not None
556 return self.C_N
558 def get_band_structure(self, path, modes=False, born=False, verbose=True):
559 omega_kl = self.band_structure(path.kpts, modes, born, verbose)
560 if modes:
561 assert 0
562 omega_kl, modes = omega_kl
564 from ase.spectrum.band_structure import BandStructure
565 bs = BandStructure(path, energies=omega_kl[None])
566 return bs
568 def compute_dynamical_matrix(self, q_scaled: np.ndarray, D_N: np.ndarray):
569 """ Computation of the dynamical matrix in momentum space D_ab(q).
570 This is a Fourier transform from real-space dynamical matrix D_N
571 for a given momentum vector q.
573 q_scaled: q vector in scaled coordinates.
575 D_N: the dynamical matrix in real-space. It is necessary, at least
576 currently, to provide this matrix explicitly (rather than use
577 self.D_N) because this matrix is modified by the Born charges
578 contributions and these modifications are momentum (q) dependent.
580 Result:
581 D(q): two-dimensional, complex-valued array of
582 shape=(3 * natoms, 3 * natoms).
583 """
584 # Evaluate fourier sum
585 R_cN = self._lattice_vectors_array
586 phase_N = np.exp(-2.j * pi * np.dot(q_scaled, R_cN))
587 D_q = np.sum(phase_N[:, np.newaxis, np.newaxis] * D_N, axis=0)
588 return D_q
590 def band_structure(self, path_kc, modes=False, born=False, verbose=True):
591 """Calculate phonon dispersion along a path in the Brillouin zone.
593 The dynamical matrix at arbitrary q-vectors is obtained by Fourier
594 transforming the real-space force constants. In case of negative
595 eigenvalues (squared frequency), the corresponding negative frequency
596 is returned.
598 Frequencies and modes are in units of eV and Ang/sqrt(amu),
599 respectively.
601 Parameters:
603 path_kc: ndarray
604 List of k-point coordinates (in units of the reciprocal lattice
605 vectors) specifying the path in the Brillouin zone for which the
606 dynamical matrix will be calculated.
607 modes: bool
608 Returns both frequencies and modes when True.
609 born: bool
610 Include non-analytic part given by the Born effective charges and
611 the static part of the high-frequency dielectric tensor. This
612 contribution to the force constant accounts for the splitting
613 between the LO and TO branches for q -> 0.
614 verbose: bool
615 Print warnings when imaginary frequncies are detected.
617 """
619 assert self.D_N is not None
620 if born:
621 assert self.Z_avv is not None
622 assert self.eps_vv is not None
624 # Dynamical matrix in real-space
625 D_N = self.D_N
627 # Lists for frequencies and modes along path
628 omega_kl = []
629 u_kl = []
631 # Reciprocal basis vectors for use in non-analytic contribution
632 reci_vc = 2 * pi * la.inv(self.atoms.cell)
633 # Unit cell volume in Bohr^3
634 vol = abs(la.det(self.atoms.cell)) / units.Bohr**3
636 for q_c in path_kc:
638 # Add non-analytic part
639 if born:
640 # q-vector in cartesian coordinates
641 q_v = np.dot(reci_vc, q_c)
642 # Non-analytic contribution to force constants in atomic units
643 qdotZ_av = np.dot(q_v, self.Z_avv).ravel()
644 C_na = (4 * pi * np.outer(qdotZ_av, qdotZ_av) /
645 np.dot(q_v, np.dot(self.eps_vv, q_v)) / vol)
646 self.C_na = C_na / units.Bohr**2 * units.Hartree
647 # Add mass prefactor and convert to eV / (Ang^2 * amu)
648 M_inv = np.outer(self.m_inv_x, self.m_inv_x)
649 D_na = C_na * M_inv / units.Bohr**2 * units.Hartree
650 self.D_na = D_na
651 D_N = self.D_N + D_na / np.prod(self.supercell)
653 # if np.prod(self.N_c) == 1:
654 #
655 # q_av = np.tile(q_v, len(self.indices))
656 # q_xx = np.vstack([q_av]*len(self.indices)*3)
657 # D_m += q_xx
659 # Evaluate fourier sum
660 D_q = self.compute_dynamical_matrix(q_c, D_N)
662 if modes:
663 omega2_l, u_xl = la.eigh(D_q, UPLO='U')
664 # Sort eigenmodes according to eigenvalues (see below) and
665 # multiply with mass prefactor
666 u_lx = (self.m_inv_x[:, np.newaxis] *
667 u_xl[:, omega2_l.argsort()]).T.copy()
668 u_kl.append(u_lx.reshape((-1, len(self.indices), 3)))
669 else:
670 omega2_l = la.eigvalsh(D_q, UPLO='U')
672 # Sort eigenvalues in increasing order
673 omega2_l.sort()
674 # Use dtype=complex to handle negative eigenvalues
675 omega_l = np.sqrt(omega2_l.astype(complex))
677 # Take care of imaginary frequencies
678 if not np.all(omega2_l >= 0.):
679 indices = np.where(omega2_l < 0)[0]
681 if verbose:
682 print('WARNING, %i imaginary frequencies at '
683 'q = (% 5.2f, % 5.2f, % 5.2f) ; (omega_q =% 5.3e*i)'
684 % (len(indices), q_c[0], q_c[1], q_c[2],
685 omega_l[indices][0].imag))
687 omega_l[indices] = -1 * np.sqrt(np.abs(omega2_l[indices].real))
689 omega_kl.append(omega_l.real)
691 # Conversion factor: sqrt(eV / Ang^2 / amu) -> eV
692 s = units._hbar * 1e10 / sqrt(units._e * units._amu)
693 omega_kl = s * np.asarray(omega_kl)
695 if modes:
696 return omega_kl, np.asarray(u_kl)
698 return omega_kl
700 def get_dos(self, kpts=(10, 10, 10), npts=1000, delta=1e-3, indices=None):
701 from ase.spectrum.dosdata import RawDOSData
702 # dos = self.dos(kpts, npts, delta, indices)
703 kpts_kc = monkhorst_pack(kpts)
704 omega_w = self.band_structure(kpts_kc).ravel()
705 dos = RawDOSData(omega_w, np.ones_like(omega_w))
706 return dos
708 def dos(self, kpts=(10, 10, 10), npts=1000, delta=1e-3, indices=None):
709 """Calculate phonon dos as a function of energy.
711 Parameters:
713 qpts: tuple
714 Shape of Monkhorst-Pack grid for sampling the Brillouin zone.
715 npts: int
716 Number of energy points.
717 delta: float
718 Broadening of Lorentzian line-shape in eV.
719 indices: list
720 If indices is not None, the atomic-partial dos for the specified
721 atoms will be calculated.
723 """
725 # Monkhorst-Pack grid
726 kpts_kc = monkhorst_pack(kpts)
727 N = np.prod(kpts)
728 # Get frequencies
729 omega_kl = self.band_structure(kpts_kc)
730 # Energy axis and dos
731 omega_e = np.linspace(0., np.amax(omega_kl) + 5e-3, num=npts)
732 dos_e = np.zeros_like(omega_e)
734 # Sum up contribution from all q-points and branches
735 for omega_l in omega_kl:
736 diff_el = (omega_e[:, np.newaxis] - omega_l[np.newaxis, :])**2
737 dos_el = 1. / (diff_el + (0.5 * delta)**2)
738 dos_e += dos_el.sum(axis=1)
740 dos_e *= 1. / (N * pi) * 0.5 * delta
742 return omega_e, dos_e
744 def write_modes(self, q_c, branches=0, kT=units.kB * 300, born=False,
745 repeat=(1, 1, 1), nimages=30, center=False):
746 """Write modes to trajectory file.
748 Parameters:
750 q_c: ndarray
751 q-vector of the modes.
752 branches: int or list
753 Branch index of modes.
754 kT: float
755 Temperature in units of eV. Determines the amplitude of the atomic
756 displacements in the modes.
757 born: bool
758 Include non-analytic contribution to the force constants at q -> 0.
759 repeat: tuple
760 Repeat atoms (l, m, n) times in the directions of the lattice
761 vectors. Displacements of atoms in repeated cells carry a Bloch
762 phase factor given by the q-vector and the cell lattice vector R_m.
763 nimages: int
764 Number of images in an oscillation.
765 center: bool
766 Center atoms in unit cell if True (default: False).
768 """
770 if isinstance(branches, int):
771 branch_l = [branches]
772 else:
773 branch_l = list(branches)
775 # Calculate modes
776 omega_l, u_l = self.band_structure([q_c], modes=True, born=born)
777 # Repeat atoms
778 atoms = self.atoms * repeat
779 # Center
780 if center:
781 atoms.center()
783 # Here ``Na`` refers to a composite unit cell/atom dimension
784 pos_Nav = atoms.get_positions()
785 # Total number of unit cells
786 N = np.prod(repeat)
788 # Corresponding lattice vectors R_m
789 R_cN = np.indices(repeat).reshape(3, -1)
790 # Bloch phase
791 phase_N = np.exp(2.j * pi * np.dot(q_c, R_cN))
792 phase_Na = phase_N.repeat(len(self.atoms))
794 for l in branch_l:
796 omega = omega_l[0, l]
797 u_av = u_l[0, l]
799 # Mean displacement of a classical oscillator at temperature T
800 u_av *= sqrt(kT) / abs(omega)
802 mode_av = np.zeros((len(self.atoms), 3), dtype=complex)
803 # Insert slice with atomic displacements for the included atoms
804 mode_av[self.indices] = u_av
805 # Repeat and multiply by Bloch phase factor
806 mode_Nav = np.vstack(N * [mode_av]) * phase_Na[:, np.newaxis]
808 with Trajectory('%s.mode.%d.traj' % (self.name, l), 'w') as traj:
809 for x in np.linspace(0, 2 * pi, nimages, endpoint=False):
810 atoms.set_positions((pos_Nav + np.exp(1.j * x) *
811 mode_Nav).real)
812 traj.write(atoms)