1# Copyright (C) 2015 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
36from phonopy.qha import BulkModulus, QHA
37
38
39class PhonopyQHA(object):
40    """PhonopyQHA API."""
41
42    def __init__(self,
43                 volumes=None,
44                 electronic_energies=None,
45                 temperatures=None,
46                 free_energy=None,
47                 cv=None,
48                 entropy=None,
49                 eos='vinet',
50                 t_max=None,
51                 energy_plot_factor=None,
52                 verbose=False):
53        """Init method.
54
55        Notes
56        -----
57        The first two parameters have to be in this order for the backward
58        compatibility.
59
60        Parameters
61        ----------
62        volumes: array_like
63            Unit cell volumes (V) in angstrom^3.
64            dtype='double'
65            shape=(volumes,)
66        electronic_energies: array_like
67            Electronic energies (U_el) or electronic free energies (F_el) in eV.
68            It is assumed as formar if ndim==1 and latter if ndim==2.
69            dtype='double'
70            shape=(volumes,) or (temperatuers, volumes)
71        temperatures: array_like
72            Temperatures ascending order (T) in K.
73            dtype='double'
74            shape=(temperatures,)
75        free_energy: array_like
76            Phonon Helmholtz free energy (F_ph) in kJ/mol.
77            dtype='double'
78            shape=(temperatuers, volumes)
79        cv: array_like
80            Phonon heat capacity at constant volume in J/K/mol.
81            dtype='double'
82            shape=(temperatuers, volumes)
83        entropy: array_like
84            Phonon entropy at constant volume (S_ph) in J/K/mol.
85            dtype='double'
86            shape=(temperatuers, volumes)
87        eos: str
88            Equation of state used for fitting F vs V.
89            'vinet', 'murnaghan' or 'birch_murnaghan'.
90        t_max: float
91            Maximum temperature to be calculated. This has to be not
92            greater than the temperature of the third element from the
93            end of 'temperatre' elements. If max_t=None, the temperature
94            of the third element from the end is used.
95        energy_plot_factor: float
96            This value is multiplied to energy like values only in plotting.
97        verbose: boolean
98            Show log or not.
99
100        """
101        self._bulk_modulus = BulkModulus(volumes,
102                                         electronic_energies,
103                                         eos=eos)
104
105        if temperatures is not None:
106            self._qha = QHA(volumes,
107                            electronic_energies,
108                            temperatures,
109                            cv,
110                            entropy,
111                            free_energy,
112                            eos=eos,
113                            t_max=t_max,
114                            energy_plot_factor=energy_plot_factor)
115            self._qha.run(verbose=verbose)
116
117    @property
118    def bulk_modulus(self):
119        """Return bulk modulus computed without phonon contribution.
120
121        Returns
122        -------
123        float
124            Bulk modulus calculated without phonon contribution.
125
126        """
127        return self._bulk_modulus.bulk_modulus
128
129    @property
130    def thermal_expansion(self):
131        """Return thermal expansion coefficients at temperatures.
132
133        Returns
134        -------
135        list
136            Thermal expansion coefficients at temperatues.
137            shape=(temperatures, )
138
139        """
140        return self._qha.thermal_expansion
141
142    @property
143    def helmholtz_volume(self):
144        """Return free_energies at volumes.
145
146        Returns
147        -------
148        ndarray
149            Helmholtz free energies calculated at temperatures and volumes.
150            shape=(temperatures, volumes)
151
152        """
153        return self._qha.helmholtz_volume
154
155    @property
156    def volume_temperature(self):
157        """Return volumes at temperatures.
158
159        Returns
160        -------
161        ndarray
162            Equilibrium volumes at temperatures
163            shape=(temperatures,), dtype=float
164
165        """
166        return self._qha.volume_temperature
167
168    @property
169    def gibbs_temperature(self):
170        """Return Gibbs free energies at temperatures.
171
172        Returns
173        -------
174        ndarray
175            Gibbs free energies at temperatures.
176            shape=(temperatures, ), dtype=float
177
178        """
179        return self._qha.gibbs_temperature
180
181    @property
182    def bulk_modulus_temperature(self):
183        """Return bulk modulus at temperatures.
184
185        Returns
186        -------
187        ndarray
188            Bulk modulus at constant pressure and temperatures
189            shape=(temperatures, ), dtype=float
190
191        """
192        return self._qha.bulk_modulus_temperature
193
194    @property
195    def heat_capacity_P_numerical(self):
196        """Return heat capacities at constant pressure at temperatures.
197
198        These values are calculated by -T*d^2G/dT^2.
199
200        Returns
201        -------
202        list
203            Heat capacity at constant pressure and temperatures.
204            shape=(temperatures, )
205
206        """
207        return self._qha.get_heat_capacity_P_numerical()
208
209    @property
210    def heat_capacity_P_polyfit(self):
211        """Return heat capacities at constant pressure at temperatures.
212
213        Note
214        ----
215        This does not work when temperature dependent electronic_energies
216        is supplied.
217
218        Returns
219        -------
220        list
221            Heat capacities at constant pressure at temperatures, which are
222            calculated from the values obtained by polynomial fittings of
223            Cv and S.
224
225            shape=(temperatures, )
226
227        """
228        return self._qha.heat_capacity_P_polyfit
229
230    @property
231    def gruneisen_temperature(self):
232        """Return Gruneisen parameters at temperatures.
233
234        Returns
235        -------
236        list
237            Thermodynamic Gruneisen parameters at temperatures.
238            shape=(temperatures, )
239
240        """
241        return self._qha.gruneisen_temperature
242
243    def get_bulk_modulus_parameters(self):
244        """Return temperature independent bulk modulus EOS fitting parameters.
245
246        These values are those computed without phonon free energy.
247
248        (lowest energy, bulk modulus, b_prime, equilibrium volume)
249
250        """
251        return self._bulk_modulus.get_parameters()
252
253    def write_helmholtz_volume(self, filename='helmholtz-volume.dat'):
254        self._qha.write_helmholtz_volume(filename=filename)
255
256    def write_helmholtz_volume_fitted(self,
257                                      thin_number,
258                                      filename='helmholtz-volume_fitted.dat'):
259        self._qha.write_helmholtz_volume_fitted(thin_number, filename=filename)
260
261    def write_volume_temperature(self, filename='volume-temperature.dat'):
262        self._qha.write_volume_temperature(filename=filename)
263
264    def write_thermal_expansion(self,
265                                filename='thermal_expansion.dat'):
266        self._qha.write_thermal_expansion(filename=filename)
267
268    def write_gibbs_temperature(self, filename='gibbs-temperature.dat'):
269        self._qha.write_gibbs_temperature(filename=filename)
270
271    def write_bulk_modulus_temperature(self,
272                                       filename='bulk_modulus-temperature.dat'):
273        self._qha.write_bulk_modulus_temperature(filename=filename)
274
275    def plot_bulk_modulus(self):
276        """Returns matplotlib.pyplot of bulk modulus fitting curve"""
277        return self._bulk_modulus.plot()
278
279    def plot_qha(self, thin_number=10, volume_temp_exp=None):
280        """Returns matplotlib.pyplot of QHA fitting curves at temperatures"""
281        return self._qha.plot(thin_number=thin_number,
282                              volume_temp_exp=volume_temp_exp)
283
284    def plot_helmholtz_volume(self,
285                              thin_number=10,
286                              xlabel=r'Volume $(\AA^3)$',
287                              ylabel='Free energy'):
288        return self._qha.plot_helmholtz_volume(thin_number=thin_number,
289                                               xlabel=xlabel,
290                                               ylabel=ylabel)
291
292    def plot_pdf_helmholtz_volume(self,
293                                  thin_number=10,
294                                  filename='helmholtz-volume.pdf'):
295        self._qha.plot_pdf_helmholtz_volume(thin_number=thin_number,
296                                            filename=filename)
297
298    def plot_volume_temperature(self, exp_data=None):
299        return self._qha.plot_volume_temperature(exp_data=exp_data)
300
301    def plot_pdf_volume_temperature(self,
302                                    exp_data=None,
303                                    filename='volume-temperature.pdf'):
304        self._qha.plot_pdf_volume_temperature(exp_data=exp_data,
305                                              filename=filename)
306
307    def plot_thermal_expansion(self):
308        return self._qha.plot_thermal_expansion()
309
310    def plot_pdf_thermal_expansion(self,
311                                   filename='thermal_expansion.pdf'):
312        self._qha.plot_pdf_thermal_expansion(filename=filename)
313
314    def plot_gibbs_temperature(self,
315                               xlabel='Temperature (K)',
316                               ylabel='Gibbs free energy'):
317        return self._qha.plot_gibbs_temperature(xlabel=xlabel,
318                                                ylabel=ylabel)
319
320    def plot_pdf_gibbs_temperature(self, filename='gibbs-temperature.pdf'):
321        self._qha.plot_pdf_gibbs_temperature(filename=filename)
322
323    def plot_bulk_modulus_temperature(self,
324                                      xlabel='Temperature (K)',
325                                      ylabel='Bulk modulus'):
326        return self._qha.plot_bulk_modulus_temperature(xlabel=xlabel,
327                                                       ylabel=ylabel)
328
329    def plot_pdf_bulk_modulus_temperature(self,
330                                          filename='bulk_modulus-temperature.pdf'):
331        self._qha.plot_pdf_bulk_modulus_temperature(filename=filename)
332
333    def plot_heat_capacity_P_numerical(self, Z=1, exp_data=None):
334        return self._qha.plot_heat_capacity_P_numerical(Z=Z,
335                                                        exp_data=exp_data)
336
337    def plot_pdf_heat_capacity_P_numerical(self,
338                                           exp_data=None,
339                                           filename='Cp-temperature.pdf'):
340        self._qha.plot_pdf_heat_capacity_P_numerical(exp_data=exp_data,
341                                                     filename=filename)
342
343    def write_heat_capacity_P_numerical(self, filename='Cp-temperature.dat'):
344        self._qha.write_heat_capacity_P_numerical(filename=filename)
345
346    def plot_heat_capacity_P_polyfit(self, exp_data=None, Z=1):
347        return self._qha.plot_heat_capacity_P_polyfit(
348            Z=Z, exp_data=exp_data)
349
350    def plot_pdf_heat_capacity_P_polyfit(self,
351                                         exp_data=None,
352                                         filename='Cp-temperature_polyfit.pdf'):
353        self._qha.plot_pdf_heat_capacity_P_polyfit(exp_data=exp_data,
354                                                   filename=filename)
355
356    def write_heat_capacity_P_polyfit(self,
357                                      filename='Cp-temperature_polyfit.dat',
358                                      filename_ev='entropy-volume.dat',
359                                      filename_cvv='Cv-volume.dat',
360                                      filename_dsdvt='dsdv-temperature.dat'):
361        self._qha.write_heat_capacity_P_polyfit(filename=filename,
362                                                filename_ev=filename_ev,
363                                                filename_cvv=filename_cvv,
364                                                filename_dsdvt=filename_dsdvt)
365
366    def plot_gruneisen_temperature(self):
367        return self._qha.plot_gruneisen_temperature()
368
369    def plot_pdf_gruneisen_temperature(self,
370                                       filename='gruneisen-temperature.pdf'):
371        self._qha.plot_pdf_gruneisen_temperature(filename=filename)
372
373    def write_gruneisen_temperature(self, filename='gruneisen-temperature.dat'):
374        self._qha.write_gruneisen_temperature(filename=filename)
375
376    def get_bulk_modulus(self):
377        warnings.warn("PhonopyQHA.get_bulk_modulus() is deprecated."
378                      "Use bulk_modulus attribute.",
379                      DeprecationWarning)
380        return self.bulk_modulus
381
382    def get_helmholtz_volume(self):
383        warnings.warn("PhonopyQHA.get_helmholtz_volume() is deprecated."
384                      "Use helmholtz_volume attribute.",
385                      DeprecationWarning)
386        return self.helmholtz_volume
387
388    def get_volume_temperature(self):
389        warnings.warn("PhonopyQHA.get_volume_temperature() is deprecated."
390                      "Use volume_temperature attribute.",
391                      DeprecationWarning)
392        return self.volume_temperature
393
394    def get_thermal_expansion(self):
395        warnings.warn("PhonopyQHA.get_thermal_expansion() is deprecated."
396                      "Use thermal_expansion attribute.",
397                      DeprecationWarning)
398        return self.thermal_expansion
399
400    def get_gibbs_temperature(self):
401        warnings.warn("PhonopyQHA.get_gibbs_temperature() is deprecated."
402                      "Use gibbs_temperature attribute.",
403                      DeprecationWarning)
404        return self.gibbs_temperature
405
406    def get_bulk_modulus_temperature(self):
407        warnings.warn("PhonopyQHA.get_bulk_modulus_temperature() is deprecated."
408                      "Use bulk_modulus_temperature attribute.",
409                      DeprecationWarning)
410        return self.bulk_modulus_temperature
411
412    def get_heat_capacity_P_numerical(self):
413        warnings.warn("PhonopyQHA.get_heat_capacity_P_numerical() is deprecated."
414                      "Use heat_capacity_P_numerical attribute.",
415                      DeprecationWarning)
416        return self.heat_capacity_P_numerical
417
418    def get_heat_capacity_P_polyfit(self):
419        warnings.warn("PhonopyQHA.get_heat_capacity_P_polyfit() is deprecated."
420                      "Use heat_capacity_P_polyfit attribute.",
421                      DeprecationWarning)
422        return self.heat_capacity_P_polyfit
423
424    def get_gruneisen_temperature(self):
425        warnings.warn("PhonopyQHA.get_gruneisen_temperature() is deprecated."
426                      "Use gruneisen_temperature attribute.",
427                      DeprecationWarning)
428        return self.gruneisen_temperature
429