text: # coding: utf-8
  # Copyright (c) Max-Planck-Institut für Eisenforschung GmbH - Computational Materials Design (CM) Department
  # Distributed under the terms of "New BSD License", see the LICENSE file.
  
  from __future__ import print_function
  
  from copy import copy
  import numpy as np
  
  __author__ = "Joerg Neugebauer, Jan Janssen"
  __copyright__ = (
      "Copyright 2021, Max-Planck-Institut für Eisenforschung GmbH - "
      "Computational Materials Design (CM) Department"
  )
  __version__ = "1.0"
  __maintainer__ = "Jan Janssen"
  __email__ = "janssen@mpie.de"
  __status__ = "development"
  __date__ = "Sep 1, 2017"
  
  
  class ThermoBulk(object):
      """
      Class should provide all tools to compute bulk thermodynamic quantities. Central quantity is the Free Energy F(V,T).
      ToDo: Make it a (light weight) pyiron object (introduce a new tool rather than job object).
  
      Args:
          project:
          name:
  
      """
  
      eV_to_J_per_mol = 1.60217662e-19 * 6.022e23
      kB = 1 / 8.6173303e-5
  
      def __init__(self, project=None, name=None):
          # only for compatibility with pyiron objects
          self._project = project
          self._name = name
  
          self._volumes = None
          self._temperatures = None
          self._energies = None
          self._entropy = None
          self._pressure = None
          self._num_atoms = None
  
          self._fit_order = 3
  
      def copy(self):
          """
  
          Returns:
  
          """
          cls = self.__class__
          result = cls.__new__(cls)
          result.__init__()
          result.__dict__["_volumes"] = copy(self._volumes)
          result.__dict__["_temperatures"] = copy(self._temperatures)
          result.__dict__["_energies"] = copy(self._energies)
          result.__dict__["_fit_order"] = self._fit_order
          return result
  
      def _reset_energy(self):
          """
  
          Returns:
  
          """
          if self._volumes is not None:
              if self._temperatures is not None:
                  self._energies = np.zeros((len(self._temperatures), len(self._volumes)))
                  # self.energies = 0
  
      @property
      def num_atoms(self):
          """
  
          Returns:
  
          """
          if self._num_atoms is None:
              return 1  # normalize per cell if number of atoms unknown
          return self._num_atoms
  
      @num_atoms.setter
      def num_atoms(self, num):
          """
  
          Args:
              num:
  
          Returns:
  
          """
          self._num_atoms = num
  
      @property
      def _coeff(self):
          """
  
          Returns:
  
          """
          return np.polyfit(self._volumes, self._energies.T, deg=self._fit_order)
  
      @property
      def temperatures(self):
          """
  
          Returns:
  
          """
          return self._temperatures
  
      @property
      def _d_temp(self):
          """
  
          Returns:
  
          """
          return self.temperatures[1] - self.temperatures[0]
  
      @property
      def _d_vol(self):
          """
  
          Returns:
  
          """
          return self.volumes[1] - self.volumes[0]
  
      @temperatures.setter
      def temperatures(self, temp_lst):
          """
  
          Args:
              temp_lst:
  
          Returns:
  
          """
          if not hasattr(temp_lst, "__len__"):
              raise ValueError("Requires list as input parameter")
          len_temp = -1
          if self._temperatures is not None:
              len_temp = len(self._temperatures)
          self._temperatures = np.array(temp_lst)
          if len(temp_lst) != len_temp:
              self._reset_energy()
  
      @property
      def volumes(self):
          """
  
          Returns:
  
          """
          return self._volumes
  
      @volumes.setter
      def volumes(self, volume_lst):
          """
  
          Args:
              volume_lst:
  
          Returns:
  
          """
          if not hasattr(volume_lst, "__len__"):
              raise ValueError("Requires list as input parameter")
          len_vol = -1
          if self._volumes is not None:
              len_vol = len(self._volumes)
          self._volumes = np.array(volume_lst)
          if len(volume_lst) != len_vol:
              self._reset_energy()
  
      @property
      def entropy(self):
          """
  
          Returns:
  
          """
          if self._entropy is None:
              self._compute_thermo()
          return self._entropy
  
      @property
      def pressure(self):
          """
  
          Returns:
  
          """
          if self._pressure is None:
              self._compute_thermo()
          return self._pressure
  
      @property
      def energies(self):
          """
  
          Returns:
  
          """
          return self._energies
  
      @energies.setter
      def energies(self, erg_lst):
          """
  
          Args:
              erg_lst:
  
          Returns:
  
          """
          if np.ndim(erg_lst) == 2:
              self._energies = erg_lst
          elif np.ndim(erg_lst) == 1:
              if len(erg_lst) == len(self.volumes):
                  self._energies = np.tile(erg_lst, (len(self.temperatures), 1))
              else:
                  raise ValueError()
          else:
              self._energies = (
                  np.ones((len(self.volumes), len(self.temperatures))) * erg_lst
              )
  
      def set_temperatures(
          self, temperature_min=0, temperature_max=1500, temperature_steps=50
      ):
          """
  
          Args:
              temperature_min:
              temperature_max:
              temperature_steps:
  
          Returns:
  
          """
          self.temperatures = np.linspace(
              temperature_min, temperature_max, temperature_steps
          )
  
      def set_volumes(self, volume_min, volume_max=None, volume_steps=10):
          """
  
          Args:
              volume_min:
              volume_max:
              volume_steps:
  
          Returns:
  
          """
          if volume_max is None:
              volume_max = 1.1 * volume_min
          self.volumes = np.linspace(volume_min, volume_max, volume_steps)
  
      def meshgrid(self):
          """
  
          Returns:
  
          """
          return np.meshgrid(self.volumes, self.temperatures)
  
      def get_minimum_energy_path(self, pressure=None):
          """
  
          Args:
              pressure:
  
          Returns:
  
          """
          if pressure is not None:
              raise NotImplemented()
          v_min_lst = []
          for c in self._coeff.T:
              v_min = np.roots(np.polyder(c, 1))
              p_der2 = np.polyder(c, 2)
              p_val2 = np.polyval(p_der2, v_min)
              v_m_lst = v_min[p_val2 > 0]
              if len(v_m_lst) > 0:
                  v_min_lst.append(v_m_lst[0])
              else:
                  v_min_lst.append(np.nan)
          return np.array(v_min_lst)
  
      def get_free_energy(self, vol, pressure=None):
          """
  
          Args:
              vol:
              pressure:
  
          Returns:
  
          """
          if not pressure:
              return np.polyval(self._coeff, vol)
          else:
              raise NotImplementedError()
  
      def interpolate_volume(self, volumes, fit_order=None):
          """
  
          Args:
              volumes:
              fit_order:
  
          Returns:
  
          """
          if fit_order is not None:
              self._fit_order = fit_order
          new = self.copy()
          new.volumes = volumes
          new.energies = np.array([np.polyval(self._coeff, v) for v in volumes]).T
          return new
  
      def _compute_thermo(self):
          """
  
          Returns:
  
          """
          self._entropy, self._pressure = np.gradient(
              -self.energies, self._d_temp, self._d_vol
          )
  
      def get_free_energy_p(self):
          """
  
          Returns:
  
          """
          coeff = np.polyfit(self._volumes, self.energies.T, deg=self._fit_order)
          return np.polyval(coeff, self.get_minimum_energy_path())
  
      def get_entropy_p(self):
          """
  
          Returns:
  
          """
          s_coeff = np.polyfit(self._volumes, self.entropy.T, deg=self._fit_order)
          return np.polyval(s_coeff, self.get_minimum_energy_path())
  
      def get_entropy_v(self):
          """
  
          Returns:
  
          """
          eq_volume = self.volumes[0]
          s_coeff = np.polyfit(self.volumes, self.entropy.T, deg=self._fit_order)
          const_v = eq_volume * np.ones(len(s_coeff.T))
          return np.polyval(s_coeff, const_v)
  
      def plot_free_energy(self):
          """
  
          Returns:
  
          """
          try:
              import pylab as plt
          except ImportError:
              import matplotlib.pyplot as plt
          plt.plot(self.temperatures, self.get_free_energy_p() / self.num_atoms)
          plt.xlabel("Temperature [K]")
          plt.ylabel("Free energy [eV]")
  
      def plot_entropy(self):
          """
  
          Returns:
  
          """
          try:
              import pylab as plt
          except ImportError:
              import matplotlib.pyplot as plt
          plt.plot(
              self.temperatures,
              self.eV_to_J_per_mol / self.num_atoms * self.get_entropy_p(),
              label="S$_p$",
          )
          plt.plot(
              self.temperatures,
              self.eV_to_J_per_mol / self.num_atoms * self.get_entropy_v(),
              label="S$_V$",
          )
          plt.legend()
          plt.xlabel("Temperature [K]")
          plt.ylabel("Entropy [J K$^{-1}$ mol-atoms$^{-1}$]")
  
      def plot_heat_capacity(self, to_kB=True):
          """
  
          Args:
              to_kB:
  
          Returns:
  
          """
          try:
              import pylab as plt
          except ImportError:
              import matplotlib.pyplot as plt
          if to_kB:
              units = self.kB / self.num_atoms
              plt.ylabel("Heat capacity [kB]")
          else:
              units = self.eV_to_J_per_mol
              plt.ylabel("Heat capacity [J K$^{-1}$ mol-atoms$^{-1}$]")
          temps = self.temperatures[:-2]
          c_p = temps * np.gradient(self.get_entropy_p(), self._d_temp)[:-2]
          c_v = temps * np.gradient(self.get_entropy_v(), self._d_temp)[:-2]
          plt.plot(temps, units * c_p, label="c$_p$")
          plt.plot(temps, units * c_v, label="c$_v$")
          plt.legend(loc="lower right")
          plt.xlabel("Temperature [K]")
  
      def contour_pressure(self):
          """
  
          Returns:
  
          """
          try:
              import pylab as plt
          except ImportError:
              import matplotlib.pyplot as plt
          x, y = self.meshgrid()
          p_coeff = np.polyfit(self.volumes, self.pressure.T, deg=self._fit_order)
          p_grid = np.array([np.polyval(p_coeff, v) for v in self._volumes]).T
          plt.contourf(x, y, p_grid)
          plt.plot(self.get_minimum_energy_path(), self.temperatures)
          plt.xlabel("Volume [$\AA^3$]")
          plt.ylabel("Temperature [K]")
  
      def contour_entropy(self):
          """
  
          Returns:
  
          """
          try:
              import pylab as plt
          except ImportError:
              import matplotlib.pyplot as plt
          s_coeff = np.polyfit(self.volumes, self.entropy.T, deg=self._fit_order)
          s_grid = np.array([np.polyval(s_coeff, v) for v in self.volumes]).T
          x, y = self.meshgrid()
          plt.contourf(x, y, s_grid)
          plt.plot(self.get_minimum_energy_path(), self.temperatures)
          plt.xlabel("Volume [$\AA^3$]")
          plt.ylabel("Temperature [K]")
  
      def plot_contourf(self, ax=None, show_min_erg_path=False):
          """
  
          Args:
              ax:
              show_min_erg_path:
  
          Returns:
  
          """
          try:
              import pylab as plt
          except ImportError:
              import matplotlib.pyplot as plt
          x, y = self.meshgrid()
          if ax is None:
              fig, ax = plt.subplots(1, 1)
          ax.contourf(x, y, self.energies)
          if show_min_erg_path:
              plt.plot(self.get_minimum_energy_path(), self.temperatures, "w--")
          plt.xlabel("Volume [$\AA^3$]")
          plt.ylabel("Temperature [K]")
          return ax
  
      def plot_min_energy_path(self, *args, ax=None, **qwargs):
          """
  
          Args:
              *args:
              ax:
              **qwargs:
  
          Returns:
  
          """
          try:
              import pylab as plt
          except ImportError:
              import matplotlib.pyplot as plt
          if ax is None:
              fig, ax = plt.subplots(1, 1)
              ax.xlabel("Volume [$\AA^3$]")
              ax.ylabel("Temperature [K]")
          ax.plot(self.get_minimum_energy_path(), self.temperatures, *args, **qwargs)
          return ax
  