1# Copyright (C) 2011 Atsushi Togo
2# All rights reserved.
3#
4# This file is part of phonopy.
5#
6# Redistribution and use in source and binary forms, with or without
7# modification, are permitted provided that the following conditions
8# are met:
9#
10# * Redistributions of source code must retain the above copyright
11#   notice, this list of conditions and the following disclaimer.
12#
13# * Redistributions in binary form must reproduce the above copyright
14#   notice, this list of conditions and the following disclaimer in
15#   the documentation and/or other materials provided with the
16#   distribution.
17#
18# * Neither the name of the phonopy project nor the names of its
19#   contributors may be used to endorse or promote products derived
20#   from this software without specific prior written permission.
21#
22# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
23# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
24# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
25# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
26# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
27# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
28# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
29# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
30# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
32# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
33# POSSIBILITY OF SUCH DAMAGE.
34
35import warnings
36import numpy as np
37from phonopy.units import Kb, THzToEv, EvTokJmol
38
39
40def mode_cv(temp, freqs):  # freqs (eV)
41    x = freqs / Kb / temp
42    expVal = np.exp(x)
43    return Kb * x ** 2 * expVal / (expVal - 1.0) ** 2
44
45
46def mode_F(temp, freqs):
47    return (Kb * temp * np.log(1.0 - np.exp((- freqs) / (Kb * temp)))
48            + freqs / 2)
49
50
51def mode_S(temp, freqs):
52    val = freqs / (2 * Kb * temp)
53    return (1 / (2 * temp) * freqs * np.cosh(val) / np.sinh(val)
54            - Kb * np.log(2 * np.sinh(val)))
55
56
57def mode_ZPE(temp, freqs):
58    return freqs / 2
59
60
61def mode_zero(temp, freqs):
62    return np.zeros_like(freqs)
63
64
65class ThermalPropertiesBase(object):
66    def __init__(self,
67                 mesh,
68                 is_projection=False,
69                 band_indices=None,
70                 cutoff_frequency=None,
71                 pretend_real=False):
72        self._is_projection = is_projection
73        self._band_indices = None
74
75        if cutoff_frequency is None or cutoff_frequency < 0:
76            self._cutoff_frequency = 0.0
77        else:
78            self._cutoff_frequency = cutoff_frequency
79
80        if band_indices is not None:
81            bi = np.hstack(band_indices).astype('intc')
82            self._band_indices = bi
83            self._frequencies = np.array(mesh.frequencies[:, bi],
84                                         dtype='double', order='C')
85            if mesh.eigenvectors is not None:
86                self._eigenvectors = np.array(mesh.eigenvectors[:, :, bi],
87                                              dtype='double', order='C')
88        else:
89            self._frequencies = mesh.frequencies
90            self._eigenvectors = mesh.eigenvectors
91
92        if pretend_real:
93            self._frequencies = abs(self._frequencies)
94        self._frequencies = np.array(self._frequencies,
95                                     dtype='double', order='C') * THzToEv
96        self._weights = mesh.weights
97        self._num_modes = self._frequencies.shape[1] * self._weights.sum()
98        self._num_integrated_modes = np.sum(
99            self._weights * (self._frequencies >
100                             self._cutoff_frequency).sum(axis=1))
101
102    def run_free_energy(self, t):
103        if t > 0:
104            free_energy = self._calculate_thermal_property(mode_F, t)
105        else:
106            free_energy = self._calculate_thermal_property(mode_ZPE, None)
107        return free_energy / np.sum(self._weights) * EvTokJmol
108
109    def run_heat_capacity(self, t):
110        if t > 0:
111            cv = self._calculate_thermal_property(mode_cv, t)
112        else:
113            cv = self._calculate_thermal_property(mode_zero, None)
114        return cv / np.sum(self._weights) * EvTokJmol
115
116    def run_entropy(self, t):
117        if t > 0:
118            entropy = self._calculate_thermal_property(mode_S, t)
119        else:
120            entropy = self._calculate_thermal_property(mode_zero, None)
121        return entropy / np.sum(self._weights) * EvTokJmol
122
123    def _calculate_thermal_property(self, func, t):
124        if not self._is_projection:
125            t_property = 0.0
126            for freqs, w in zip(self._frequencies, self._weights):
127                cond = freqs > self._cutoff_frequency
128                t_property += np.sum(func(t, freqs[cond])) * w
129            return t_property
130        else:
131            t_property = np.zeros(len(self._frequencies[0]), dtype='double')
132            for freqs, eigvecs2, w in zip(self._frequencies,
133                                          np.abs(self._eigenvectors) ** 2,
134                                          self._weights):
135                cond = freqs > self._cutoff_frequency
136                t_property += np.dot(eigvecs2[:, cond],
137                                     func(t, freqs[cond])) * w
138            return t_property
139
140
141class ThermalProperties(ThermalPropertiesBase):
142    def __init__(self,
143                 mesh,
144                 is_projection=False,
145                 band_indices=None,
146                 cutoff_frequency=None,
147                 pretend_real=False):
148        ThermalPropertiesBase.__init__(self,
149                                       mesh,
150                                       is_projection=is_projection,
151                                       band_indices=band_indices,
152                                       cutoff_frequency=cutoff_frequency,
153                                       pretend_real=pretend_real)
154        self._thermal_properties = None
155        self._temperatures = None
156        self._high_T_entropy = None
157        self._zero_point_energy = None
158        self._projected_thermal_properties = None
159
160        self._set_high_T_entropy_and_zero_point_energy()
161
162    @property
163    def temperatures(self):
164        return self._temperatures
165
166    def get_temperatures(self):
167        warnings.warn("ThermalProperties.get_temperatures is deprecated."
168                      "Use temperatures attribute.",
169                      DeprecationWarning)
170        return self.temperatures
171
172    @temperatures.setter
173    def temperatures(self, temperatures):
174        t_array = np.array(temperatures, dtype='double')
175        self._temperatures = np.array(
176            np.extract(np.invert(t_array < 0), t_array), dtype='double')
177
178    def set_temperatures(self, temperatures):
179        warnings.warn("ThermalProperties.set_temperatures is deprecated."
180                      "Use temperatures attribute.",
181                      DeprecationWarning)
182        self.temperatures = temperatures
183
184    @property
185    def thermal_properties(self):
186        return self._thermal_properties
187
188    def get_thermal_properties(self):
189        warnings.warn("ThermalProperties.get_thermal_properties is deprecated."
190                      "Use thermal_properties attribute.",
191                      DeprecationWarning)
192        return self.thermal_properties
193
194    @property
195    def zero_point_energy(self):
196        return self._zero_point_energy
197
198    def get_zero_point_energy(self):
199        warnings.warn("ThermalProperties.get_zero_point_energy is deprecated."
200                      "Use zero_point_energy attribute.",
201                      DeprecationWarning)
202        return self.zero_point_energy
203
204    @property
205    def high_T_entropy(self):
206        return self._high_T_entropy
207
208    def get_high_T_entropy(self):
209        warnings.warn("ThermalProperties.get_high_T_entropy is deprecated."
210                      "Use high_T_entropy attribute.",
211                      DeprecationWarning)
212        return self.high_T_entropy
213
214    @property
215    def number_of_integrated_modes(self):
216        """Number of phonon modes used for integration on sampling mesh"""
217        return self._num_integrated_modes
218
219    def get_number_of_integrated_modes(self):
220        warnings.warn("ThermalProperties.get_number_of_integrated_modes is "
221                      "deprecated. Use number_of_integrated_modes attribute.",
222                      DeprecationWarning)
223        return self.number_of_integrated_modes
224
225    @property
226    def number_of_modes(self):
227        """Number of phonon modes on sampling mesh"""
228        return self._num_modes
229
230    def get_number_of_modes(self):
231        warnings.warn("ThermalProperties.get_number_of_modes is "
232                      "deprecated. Use number_of_modes attribute.",
233                      DeprecationWarning)
234        return self.number_of_modes
235
236    def set_temperature_range(self, t_min=None, t_max=None, t_step=None):
237        if t_min is None:
238            _t_min = 10
239        elif t_min < 0:
240            _t_min = 0
241        else:
242            _t_min = t_min
243
244        if t_max is None:
245            _t_max = 1000
246        elif t_max > _t_min:
247            _t_max = t_max
248        else:
249            _t_max = _t_min
250
251        if t_step is None:
252            _t_step = 10
253        elif t_step > 0:
254            _t_step = t_step
255        else:
256            _t_step = 10
257
258        self._temperatures = np.arange(_t_min, _t_max + _t_step / 2.0, _t_step,
259                                       dtype='double')
260
261    def plot(self, plt):
262        temps, fe, entropy, cv = self._thermal_properties
263
264        plt.plot(temps, fe, 'r-')
265        plt.plot(temps, entropy, 'b-')
266        plt.plot(temps, cv, 'g-')
267        plt.legend(('Free energy [kJ/mol]', 'Entropy [J/K/mol]',
268                    r'C$_\mathrm{V}$ [J/K/mol]'),
269                   loc='best')
270        plt.grid(True)
271        plt.xlabel('Temperature [K]')
272
273    def run(self, t_step=None, t_max=None, t_min=None, lang='C'):
274        import warnings
275        if (t_step is not None or t_max is not None or t_min is not None):
276            warnings.warn("keywords for this method are depreciated. "
277                          "Use \'set_temperature_range\' or "
278                          "\'set_temperature_range\' method instead.",
279                          DeprecationWarning)
280            self.set_temperature_range(t_min=t_min, t_max=t_max, t_step=t_step)
281
282        if lang == 'C':
283            import phonopy._phonopy as phonoc
284            self._run_c_thermal_properties()
285        else:
286            self._run_py_thermal_properties()
287
288        if self._is_projection:
289            fe = []
290            entropy = []
291            cv = []
292            for t in self._temperatures:
293                fe.append(self.run_free_energy(t))
294                entropy.append(self.run_entropy(t) * 1000,)
295                cv.append(self.run_heat_capacity(t) * 1000)
296
297            self._projected_thermal_properties = [
298                self._temperatures,
299                np.array(fe, dtype='double'),
300                np.array(entropy, dtype='double'),
301                np.array(cv, dtype='double')]
302
303    def write_yaml(self, filename='thermal_properties.yaml', volume=None):
304        lines = self._get_tp_yaml_lines(volume=volume)
305        if self._is_projection:
306            lines += self._get_projected_tp_yaml_lines()
307        with open(filename, 'w') as f:
308            f.write("\n".join(lines))
309
310    def _run_c_thermal_properties(self):
311        import phonopy._phonopy as phonoc
312
313        props = np.zeros((len(self._temperatures), 3),
314                         dtype='double', order='C')
315        phonoc.thermal_properties(props,
316                                  self._temperatures,
317                                  self._frequencies,
318                                  self._weights,
319                                  self._cutoff_frequency)
320        # for f, w in zip(self._frequencies, self._weights):
321        #     phonoc.thermal_properties(
322        #         props,
323        #         self._temperatures,
324        #         np.array(f, dtype='double', order='C')[None, :],
325        #         np.array([w], dtype='intc'),
326        #         cutoff_frequency)
327
328        props /= np.sum(self._weights)
329        fe = props[:, 0] * EvTokJmol + self._zero_point_energy
330        entropy = props[:, 1] * EvTokJmol * 1000
331        cv = props[:, 2] * EvTokJmol * 1000
332        self._thermal_properties = [self._temperatures, fe, entropy, cv]
333
334    def _run_py_thermal_properties(self):
335        fe = []
336        entropy = []
337        cv = []
338        for t in self._temperatures:
339            props = self._get_py_thermal_properties(t)
340            fe.append(props[0])
341            entropy.append(props[1] * 1000)
342            cv.append(props[2] * 1000)
343        self._thermal_properties = [
344            self._temperatures,
345            np.array(fe, dtype='double'),
346            np.array(entropy, dtype='double'),
347            np.array(cv, dtype='double')]
348
349    def _get_tp_yaml_lines(self, volume=None):
350        lines = []
351        lines.append("# Thermal properties / unit cell (natom)")
352        lines.append("")
353        lines.append("unit:")
354        lines.append("  temperature:   K")
355        lines.append("  free_energy:   kJ/mol")
356        lines.append("  entropy:       J/K/mol")
357        lines.append("  heat_capacity: J/K/mol")
358        lines.append("")
359        lines.append("natom: %-5d" % (self._frequencies[0].shape[0] // 3))
360        if volume is not None:
361            lines.append("volume: %-20.10f" % volume)
362        lines.append("cutoff_frequency: %8.3f" % self._cutoff_frequency)
363        lines.append("num_modes: %d" % self._num_modes)
364        lines.append("num_integrated_modes: %d" % self._num_integrated_modes)
365        if self._band_indices is not None:
366            bi = self._band_indices + 1
367            lines.append("band_index: [ " + ("%d, " * (len(bi) - 1)) %
368                         tuple(bi[:-1]) + ("%d ]" % bi[-1]))
369        lines.append("")
370        lines.append("zero_point_energy: %15.7f" % self._zero_point_energy)
371        lines.append("high_T_entropy:    %15.7f" %
372                     (self._high_T_entropy * 1000))
373        lines.append("")
374        lines.append("thermal_properties:")
375        temperatures, fe, entropy, cv = self._thermal_properties
376        for i, t in enumerate(temperatures):
377            lines.append("- temperature:   %15.7f" % t)
378            lines.append("  free_energy:   %15.7f" % fe[i])
379            lines.append("  entropy:       %15.7f" % entropy[i])
380            # Sometimes 'nan' of C_V is returned at low temperature.
381            if np.isnan(cv[i]):
382                lines.append("  heat_capacity: %15.7f" % 0)
383            else:
384                lines.append("  heat_capacity: %15.7f" % cv[i])
385            lines.append("  energy:        %15.7f" %
386                         (fe[i] + entropy[i] * t / 1000))
387            lines.append("")
388        return lines
389
390    def _get_projected_tp_yaml_lines(self):
391        lines = []
392        lines.append("projected_thermal_properties:")
393        temperatures, fe, entropy, cv = self._projected_thermal_properties
394        for i, t in enumerate(temperatures):
395            lines.append("- temperature:   %13.7f" % t)
396            line = "  free_energy:   [ "
397            line += ", ".join(["%13.7f" % x for x in fe[i]])
398            line += " ] # %13.7f" % np.sum(fe[i])
399            lines.append(line)
400            line = "  entropy:       [ "
401            line += ", ".join(["%13.7f" % x for x in entropy[i]])
402            line += " ] # %13.7f" % np.sum(entropy[i])
403            lines.append(line)
404            # Sometimes 'nan' of C_V is returned at low temperature.
405            line = "  heat_capacity: [ "
406            sum_cv = 0.0
407            for j, cv_i in enumerate(cv[i]):
408                if np.isnan(cv_i):
409                    line += "%13.7f" % 0
410                else:
411                    sum_cv += cv_i
412                    line += "%13.7f" % cv_i
413                if j < len(cv[i]) - 1:
414                    line += ", "
415                else:
416                    line += " ]"
417            line += " # %13.7f" % sum_cv
418            lines.append(line)
419            energy = fe[i] + entropy[i] * t / 1000
420            line = "  energy:        [ "
421            line += ", ".join(["%13.7f" % x for x in energy])
422            line += " ] # %13.7f" % np.sum(energy)
423            lines.append(line)
424        return lines
425
426    def _get_py_thermal_properties(self, t):
427        return (self.run_free_energy(t),
428                self.run_entropy(t),
429                self.run_heat_capacity(t))
430
431    def _set_high_T_entropy_and_zero_point_energy(self):
432        zp_energy = 0.0
433        entropy = 0.0
434        for freqs, w in zip(self._frequencies, self._weights):
435            positive_fs = np.extract(freqs > 0.0, freqs)
436            entropy -= np.sum(np.log(positive_fs)) * w
437            zp_energy += np.sum(positive_fs) * w / 2
438        self._high_T_entropy = entropy * Kb / np.sum(self._weights) * EvTokJmol
439        self._zero_point_energy = zp_energy / np.sum(self._weights) * EvTokJmol
440