1# coding: utf-8 2# Copyright (c) Pymatgen Development Team. 3# Distributed under the terms of the MIT License. 4 5""" 6This module defines classes to represent the phonon density of states, etc. 7""" 8 9import numpy as np 10import scipy.constants as const 11from monty.functools import lazy_property 12from monty.json import MSONable 13 14from pymatgen.core.structure import Structure 15from pymatgen.util.coord import get_linear_interpolated_value 16 17BOLTZ_THZ_PER_K = const.value("Boltzmann constant in Hz/K") / const.tera # Boltzmann constant in THz/K 18THZ_TO_J = const.value("hertz-joule relationship") * const.tera 19 20 21def coth(x): 22 """ 23 Coth function. 24 25 Args: 26 x (): value 27 28 Returns: 29 coth(x) 30 """ 31 return 1.0 / np.tanh(x) 32 33 34class PhononDos(MSONable): 35 """ 36 Basic DOS object. All other DOS objects are extended versions of this 37 object. 38 """ 39 40 def __init__(self, frequencies, densities): 41 """ 42 Args: 43 frequencies: A sequences of frequencies in THz 44 densities: A list representing the density of states. 45 """ 46 self.frequencies = np.array(frequencies) 47 self.densities = np.array(densities) 48 49 def get_smeared_densities(self, sigma): 50 """ 51 Returns the densities, but with a Gaussian smearing of 52 std dev sigma applied. 53 54 Args: 55 sigma: Std dev of Gaussian smearing function. 56 57 Returns: 58 Gaussian-smeared densities. 59 """ 60 61 from scipy.ndimage.filters import gaussian_filter1d 62 63 diff = [self.frequencies[i + 1] - self.frequencies[i] for i in range(len(self.frequencies) - 1)] 64 avgdiff = sum(diff) / len(diff) 65 66 smeared_dens = gaussian_filter1d(self.densities, sigma / avgdiff) 67 return smeared_dens 68 69 def __add__(self, other): 70 """ 71 Adds two DOS together. Checks that frequency scales are the same. 72 Otherwise, a ValueError is thrown. 73 74 Args: 75 other: Another DOS object. 76 77 Returns: 78 Sum of the two DOSs. 79 """ 80 if not all(np.equal(self.frequencies, other.frequencies)): 81 raise ValueError("Frequencies of both DOS are not compatible!") 82 densities = self.densities + other.densities 83 return PhononDos(self.frequencies, densities) 84 85 def __radd__(self, other): 86 """ 87 Reflected addition of two DOS objects 88 89 Args: 90 other: Another DOS object. 91 92 Returns: 93 Sum of the two DOSs. 94 """ 95 96 return self.__add__(other) 97 98 def get_interpolated_value(self, frequency): 99 """ 100 Returns interpolated density for a particular frequency. 101 102 Args: 103 frequency: frequency to return the density for. 104 """ 105 return get_linear_interpolated_value(self.frequencies, self.densities, frequency) 106 107 def __str__(self): 108 """ 109 Returns a string which can be easily plotted (using gnuplot). 110 """ 111 stringarray = ["#{:30s} {:30s}".format("Frequency", "Density")] 112 for i, frequency in enumerate(self.frequencies): 113 stringarray.append("{:.5f} {:.5f}".format(frequency, self.densities[i])) 114 return "\n".join(stringarray) 115 116 @classmethod 117 def from_dict(cls, d): 118 """ 119 Returns PhononDos object from dict representation of PhononDos. 120 """ 121 return cls(d["frequencies"], d["densities"]) 122 123 def as_dict(self): 124 """ 125 Json-serializable dict representation of PhononDos. 126 """ 127 return { 128 "@module": self.__class__.__module__, 129 "@class": self.__class__.__name__, 130 "frequencies": list(self.frequencies), 131 "densities": list(self.densities), 132 } 133 134 @lazy_property 135 def ind_zero_freq(self): 136 """ 137 Index of the first point for which the freqencies are equal or greater than zero. 138 """ 139 ind = np.searchsorted(self.frequencies, 0) 140 if ind >= len(self.frequencies): 141 raise ValueError("No positive frequencies found") 142 return ind 143 144 @lazy_property 145 def _positive_frequencies(self): 146 """ 147 Numpy array containing the list of positive frequencies 148 """ 149 return self.frequencies[self.ind_zero_freq :] 150 151 @lazy_property 152 def _positive_densities(self): 153 """ 154 Numpy array containing the list of densities corresponding to positive frequencies 155 """ 156 return self.densities[self.ind_zero_freq :] 157 158 def cv(self, t, structure=None): 159 """ 160 Constant volume specific heat C_v at temperature T obtained from the integration of the DOS. 161 Only positive frequencies will be used. 162 Result in J/(K*mol-c). A mol-c is the abbreviation of a mole-cell, that is, the number 163 of Avogadro times the atoms in a unit cell. To compare with experimental data the result 164 should be divided by the number of unit formulas in the cell. If the structure is provided 165 the division is performed internally and the result is in J/(K*mol) 166 167 Args: 168 t: a temperature in K 169 structure: the structure of the system. If not None it will be used to determine the numer of 170 formula units 171 Returns: 172 Constant volume specific heat C_v 173 """ 174 175 if t == 0: 176 return 0 177 178 freqs = self._positive_frequencies 179 dens = self._positive_densities 180 181 def csch2(x): 182 return 1.0 / (np.sinh(x) ** 2) 183 184 wd2kt = freqs / (2 * BOLTZ_THZ_PER_K * t) 185 cv = np.trapz(wd2kt ** 2 * csch2(wd2kt) * dens, x=freqs) 186 cv *= const.Boltzmann * const.Avogadro 187 188 if structure: 189 formula_units = structure.composition.num_atoms / structure.composition.reduced_composition.num_atoms 190 cv /= formula_units 191 192 return cv 193 194 def entropy(self, t, structure=None): 195 """ 196 Vibrational entropy at temperature T obtained from the integration of the DOS. 197 Only positive frequencies will be used. 198 Result in J/(K*mol-c). A mol-c is the abbreviation of a mole-cell, that is, the number 199 of Avogadro times the atoms in a unit cell. To compare with experimental data the result 200 should be divided by the number of unit formulas in the cell. If the structure is provided 201 the division is performed internally and the result is in J/(K*mol) 202 203 Args: 204 t: a temperature in K 205 structure: the structure of the system. If not None it will be used to determine the numer of 206 formula units 207 Returns: 208 Vibrational entropy 209 """ 210 211 if t == 0: 212 return 0 213 214 freqs = self._positive_frequencies 215 dens = self._positive_densities 216 217 wd2kt = freqs / (2 * BOLTZ_THZ_PER_K * t) 218 s = np.trapz((wd2kt * coth(wd2kt) - np.log(2 * np.sinh(wd2kt))) * dens, x=freqs) 219 220 s *= const.Boltzmann * const.Avogadro 221 222 if structure: 223 formula_units = structure.composition.num_atoms / structure.composition.reduced_composition.num_atoms 224 s /= formula_units 225 226 return s 227 228 def internal_energy(self, t, structure=None): 229 """ 230 Phonon contribution to the internal energy at temperature T obtained from the integration of the DOS. 231 Only positive frequencies will be used. 232 Result in J/mol-c. A mol-c is the abbreviation of a mole-cell, that is, the number 233 of Avogadro times the atoms in a unit cell. To compare with experimental data the result 234 should be divided by the number of unit formulas in the cell. If the structure is provided 235 the division is performed internally and the result is in J/mol 236 237 Args: 238 t: a temperature in K 239 structure: the structure of the system. If not None it will be used to determine the numer of 240 formula units 241 Returns: 242 Phonon contribution to the internal energy 243 """ 244 245 if t == 0: 246 return self.zero_point_energy(structure=structure) 247 248 freqs = self._positive_frequencies 249 dens = self._positive_densities 250 251 wd2kt = freqs / (2 * BOLTZ_THZ_PER_K * t) 252 e = np.trapz(freqs * coth(wd2kt) * dens, x=freqs) / 2 253 254 e *= THZ_TO_J * const.Avogadro 255 256 if structure: 257 formula_units = structure.composition.num_atoms / structure.composition.reduced_composition.num_atoms 258 e /= formula_units 259 260 return e 261 262 def helmholtz_free_energy(self, t, structure=None): 263 """ 264 Phonon contribution to the Helmholtz free energy at temperature T obtained from the integration of the DOS. 265 Only positive frequencies will be used. 266 Result in J/mol-c. A mol-c is the abbreviation of a mole-cell, that is, the number 267 of Avogadro times the atoms in a unit cell. To compare with experimental data the result 268 should be divided by the number of unit formulas in the cell. If the structure is provided 269 the division is performed internally and the result is in J/mol 270 271 Args: 272 t: a temperature in K 273 structure: the structure of the system. If not None it will be used to determine the numer of 274 formula units 275 Returns: 276 Phonon contribution to the Helmholtz free energy 277 """ 278 279 if t == 0: 280 return self.zero_point_energy(structure=structure) 281 282 freqs = self._positive_frequencies 283 dens = self._positive_densities 284 285 wd2kt = freqs / (2 * BOLTZ_THZ_PER_K * t) 286 f = np.trapz(np.log(2 * np.sinh(wd2kt)) * dens, x=freqs) 287 288 f *= const.Boltzmann * const.Avogadro * t 289 290 if structure: 291 formula_units = structure.composition.num_atoms / structure.composition.reduced_composition.num_atoms 292 f /= formula_units 293 294 return f 295 296 def zero_point_energy(self, structure=None): 297 """ 298 Zero point energy energy of the system. Only positive frequencies will be used. 299 Result in J/mol-c. A mol-c is the abbreviation of a mole-cell, that is, the number 300 of Avogadro times the atoms in a unit cell. To compare with experimental data the result 301 should be divided by the number of unit formulas in the cell. If the structure is provided 302 the division is performed internally and the result is in J/mol 303 304 Args: 305 structure: the structure of the system. If not None it will be used to determine the numer of 306 formula units 307 Returns: 308 Phonon contribution to the internal energy 309 """ 310 311 freqs = self._positive_frequencies 312 dens = self._positive_densities 313 314 zpe = 0.5 * np.trapz(freqs * dens, x=freqs) 315 zpe *= THZ_TO_J * const.Avogadro 316 317 if structure: 318 formula_units = structure.composition.num_atoms / structure.composition.reduced_composition.num_atoms 319 zpe /= formula_units 320 321 return zpe 322 323 324class CompletePhononDos(PhononDos): 325 """ 326 This wrapper class defines a total dos, and also provides a list of PDos. 327 328 .. attribute:: pdos 329 330 Dict of partial densities of the form {Site:Densities} 331 """ 332 333 def __init__(self, structure, total_dos, pdoss): 334 """ 335 Args: 336 structure: Structure associated with this particular DOS. 337 total_dos: total Dos for structure 338 pdoss: The pdoss are supplied as an {Site: Densities} 339 """ 340 super().__init__(frequencies=total_dos.frequencies, densities=total_dos.densities) 341 self.pdos = {s: np.array(d) for s, d in pdoss.items()} 342 self.structure = structure 343 344 def get_site_dos(self, site): 345 """ 346 Get the Dos for a site. 347 348 Args: 349 site: Site in Structure associated with CompletePhononDos. 350 351 Returns: 352 PhononDos containing summed orbital densities for site. 353 """ 354 return PhononDos(self.frequencies, self.pdos[site]) 355 356 def get_element_dos(self): 357 """ 358 Get element projected Dos. 359 360 Returns: 361 dict of {Element: Dos} 362 """ 363 364 el_dos = {} 365 for site, atom_dos in self.pdos.items(): 366 el = site.specie 367 if el not in el_dos: 368 el_dos[el] = np.array(atom_dos) 369 else: 370 el_dos[el] += np.array(atom_dos) 371 return {el: PhononDos(self.frequencies, densities) for el, densities in el_dos.items()} 372 373 @classmethod 374 def from_dict(cls, d): 375 """ 376 Returns CompleteDos object from dict representation. 377 """ 378 tdos = PhononDos.from_dict(d) 379 struct = Structure.from_dict(d["structure"]) 380 pdoss = {} 381 for at, pdos in zip(struct, d["pdos"]): 382 pdoss[at] = pdos 383 384 return cls(struct, tdos, pdoss) 385 386 def as_dict(self): 387 """ 388 Json-serializable dict representation of CompletePhononDos. 389 """ 390 d = { 391 "@module": self.__class__.__module__, 392 "@class": self.__class__.__name__, 393 "structure": self.structure.as_dict(), 394 "frequencies": list(self.frequencies), 395 "densities": list(self.densities), 396 "pdos": [], 397 } 398 if len(self.pdos) > 0: 399 for at in self.structure: 400 d["pdos"].append(list(self.pdos[at])) 401 return d 402 403 def __str__(self): 404 return "Complete phonon DOS for " + str(self.structure) 405