Source code for directupsampling.thermo

import numpy as np
import pandas as pd

from directupsampling.data import GPa_to_eV_per_Angstrom3, kB
from directupsampling.eos import check_pressure_unit
from directupsampling.fesurface import FreeEnergySurface




class Thermo:

    def __init__(
        self,
        fes: FreeEnergySurface,
        tmax=1000,
        pressure=0.0001,
        temperatures=None,
        pressure_unit: str = "GPa",
    ):
        """Assumes pressure to be in GPa."""
        unit_conversion = check_pressure_unit(pressure_unit)

        self.fes = fes
        self.pressure = pressure * unit_conversion
        if temperatures is None:
            self.temperatures = np.linspace(1, tmax, tmax, dtype=int)
        else:
            self.temperatures = temperatures

        self._properties = None

    @property
    def properties(self):
        if self._properties is None:
            self.calculate_properties()

        return self._properties



    def calculate_properties(self, veps=1e-4):
        """At the moment without defects"""
        props = pd.DataFrame(index=self.temperatures)
        props.index.name = "temperature"
        prev_volume = prev_temp = prev_gibbs = prev_entropy = None
        for i, temp in enumerate(self.temperatures):
            volume = self.fes.get_volume(self.pressure, temp)
            # , ini_guess=prev_volume)
            gibbs = self.fes.get_gibbs_energy(self.pressure, temp)
            if i == 0:
                first_volume = volume
                vcoef = 0
                entropy = 0
            else:
                vcoef = (volume - prev_volume) / (temp - prev_temp) / volume
                entropy = -(gibbs - prev_gibbs) / (temp - prev_temp)

            vrel = (volume - first_volume) / first_volume
            vm = volume - veps
            vp = volume + veps
            pm = self.fes.get_pressure(vm, temp)
            pp = self.fes.get_pressure(vp, temp)
            # differs from old fortran code: pp = -(Fp2-Fp1)/Veps

            enthalpy = gibbs + temp * entropy
            compress = -2 / (vp + vm) * (vp - vm) / ((pp - pm))  # /2)
            # no /2 due to pp diff above
            bm_isot = 1 / compress

            if i > 2:
                cp = temp * (entropy - prev_entropy) / (temp - prev_temp)
            else:
                cp = 0

            tvbm_term = temp * volume * bm_isot * vcoef**2
            cv = cp - tvbm_term

            if cv > 1e-5:
                grueneisen = volume * vcoef * bm_isot / cv
                bm_adia = bm_isot * cp / cv
            else:
                grueneisen = -1
                bm_adia = bm_isot

            props.loc[temp, "volume_expansion"] = volume
            props.loc[temp, "volume_expansion_relative"] = vrel * 100
            props.loc[temp, "volume_expansion_coefficient"] = vcoef * 1e5
            props.loc[temp, "Gibbs_energy"] = gibbs
            props.loc[temp, "enthalpy"] = enthalpy
            props.loc[temp, "PV_term"] = self.pressure * volume
            props.loc[temp, "entropy"] = entropy / kB
            props.loc[temp, "heat_capacity_isobaric"] = cp / kB
            props.loc[temp, "heat_capacity_isochoric"] = cv / kB
            props.loc[temp, "TVbetaSqrBT"] = tvbm_term / kB
            props.loc[temp, "compressibility_isothermal"] = (
                compress * GPa_to_eV_per_Angstrom3
            )
            props.loc[temp, "bulk_modulus_isothermal"] = (
                bm_isot / GPa_to_eV_per_Angstrom3
            )
            props.loc[temp, "bulk_modulus_adiabatic"] = (
                bm_adia / GPa_to_eV_per_Angstrom3
            )
            props.loc[temp, "Grueneisen_parameter"] = grueneisen

            prev_volume, prev_temp = volume, temp
            prev_gibbs = gibbs
            prev_entropy = entropy

        self._properties = props



    # if defect:
    #     thermo['Gibbs_energy'] +=
    #
    # thermo['free_energy_atVeqT0K']
    #
    # thermo['free_energy_vsV_atTmax_noET0K']
    # thermo['free_energy_vsV_atTmax']

    # if some_structure:
    #     thermo['aLat_expansion']
    #     thermo['aLat_expansion_relative']
    #     thermo['aLat_expansion_coefficient']
    # if hcp:
    #     thermo['cLat_expansion']
    #     thermo['cLat_expansion_relative']
    #     thermo['cLat_expansion_coefficient']
    #     thermo['cBya']
    #
    #
    #
    #   if (Gs%nG>0) then
    #     ! Gibbs energy of defects
    #     call defectGibbsEnergy(Gs,Pm1,i,Gdm1)
    #     call defectGibbsEnergy(Gs,P  ,i,Gd0 )
    #     call defectGibbsEnergy(Gs,Pp1,i,Gdp1)
    #   else
    #     Gdm1=0; Gd0=0; Gdp1=0
    #   endif

    #   ! total Gibbs energies
    #   Gm1 = Gpm1 + Gdm1;    G = Gp0 + Gd0;    Gp1 = Gpp1 + Gdp1

    #   ! equilibrium volumes
    #   Vm1 = (G-Gm1)/(P-Pm1)
    #   Vp1 = (Gp1-G)/(Pp1-P)

    #   ! main quantities at single T
    #   VP(i) = (Vm1+Vp1)/2

    #   if (strType.NE."none") then
    #     ! lattice constant expansion
    #     if (strType.EQ."fcc") aP(i)=(4.*VP(i))**(1./3.);
    #     if (strType.EQ."bcc") aP(i)=(2.*VP(i))**(1./3.);
    #     ! for hcp also c and cBya expansion
    #     if (strType.EQ."hcp") then
    #       call cByaRatio(Fs,VP(i),cBya(i))
    #       ! hcp conversion from volume and cBya to a and c
    #       nAtoms=2
    #       aP(i) = ((nAtoms*VP(i)/sin(pi/3.))/cBya(i))**(1./3.)
    #       c(i) = cBya(i)*aP(i)
    #     endif

    #     ! relative lattice constant expansion
    #     raP(i)=(aP(i)-aP(1))/aP(1)
    #     if (strType.EQ."hcp") rc(i)=(c(i)-c(1))/c(1)
    #     ! lattice constant coefficient
    #     if (i>1) then
    #       aCoef(i) = (aP(i)-aP(i-1))/(T(i)-T(i-1))/aP(i)
    #       if (strType.EQ."hcp") cCoef(i) = (c(i)-c(i-1))/(T(i)-T(i-1))/c(i)
    #     else
    #       aCoef(i) = 0.
    #       if (strType.EQ."hcp") cCoef(i) = 0.
    #     endif
    #   endif