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