Source code for pyiron.atomistics.thermodynamics.thermo_bulk

# 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 2020, Max-Planck-Institut für Eisenforschung GmbH - "
    "Computational Materials Design (CM) Department"
__version__ = "1.0"
__maintainer__ = "Jan Janssen"
__email__ = ""
__status__ = "development"
__date__ = "Sep 1, 2017"

[docs]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
[docs] 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 )
[docs] 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 )
[docs] 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)
[docs] def meshgrid(self): """ Returns: """ return np.meshgrid(self.volumes, self.temperatures)
[docs] 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)
[docs] def get_free_energy(self, vol, pressure=None): """ Args: vol: pressure: Returns: """ if not pressure: return np.polyval(self._coeff, vol) else: raise NotImplementedError()
[docs] 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 )
[docs] 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())
[docs] 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())
[docs] 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)
[docs] 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]")
[docs] 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}$]")
[docs] 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]")
[docs] 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]")
[docs] 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]")
[docs] 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
[docs] 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