wmaee.scopes.debye

Collection of routines for thermodynamics with the Debye model.

David Holec david.holec@unileoben.ac.at

   1"""
   2Collection of routines for thermodynamics with the Debye model.
   3
   4David Holec
   5david.holec@unileoben.ac.at
   6"""
   7
   8from typing import Optional, Union, Any
   9from copy import deepcopy
  10from math import exp
  11import numpy as np
  12import scipy.integrate as integrate
  13from scipy.optimize import curve_fit
  14from ase import Atoms
  15
  16
  17def get_rho(V: float, material: Union[str, Atoms]) -> float:
  18    """
  19    Compute mass density rho [g/cm^3] from per-atom volume V [Å^3/atom]
  20    and either an ASE Atoms object or a chemical formula.
  21
  22    Parameters
  23    ----------
  24    V : float
  25        Per-atom volume in Å^3/atom.
  26    material : str or ase.Atoms
  27        Chemical formula string (e.g., 'Fe2O3', 'Si', 'AlN') or an ASE Atoms object.
  28
  29    Returns
  30    -------
  31    float
  32        Mass density in g/cm^3.
  33
  34    Notes
  35    -----
  36    - Conversions: 1 amu = 1.66054e-24 g and 1 Å^3 = 1e-24 cm^3,
  37      so rho = 1.66054 * (average_mass_in_amu) / V.
  38    """
  39    if isinstance(material, str):
  40        atoms = Atoms(material)
  41    elif isinstance(material, Atoms):
  42        atoms = material
  43    else:
  44        raise TypeError("material must be a chemical formula string or an ase.Atoms object")
  45
  46    M_amu = atoms.get_masses().mean()  # average mass per atom (amu)
  47    return 1.66054 * M_amu / V
  48
  49
  50    
  51    
  52class Debye:
  53    """
  54    Implementation of the harmonic Debye approximation.
  55
  56    This class stores and provides accessors for thermodynamic and elastic
  57    descriptors required to evaluate properties in the harmonic Debye model,
  58    such as the Helmholtz free energy, heat capacity, and entropy.
  59
  60    Parameters
  61    ----------
  62    at_per_fu : int, optional
  63        Number of atoms in the primitive cell, $N_\\mathrm{at}$, by default 1.
  64
  65    Attributes
  66    ----------
  67    data : dict
  68        Internal storage for material fields:
  69        - 'V': float or None
  70            Volume per atom, $V$ [Å^3/atom].
  71        - 'rho': float or None
  72            Mass density, $\\rho$ [g/cm^3].
  73        - 'Nat': int
  74            Number of atoms in the primitive cell, $N_\\mathrm{at}$ [1].
  75        - 'B': float or None
  76            Bulk modulus, $B$ [GPa].
  77        - 'nu': float or None
  78            Poisson's ratio, $\\nu$ [1].
  79        - 'v': float or None
  80            Mean (Debye) sound velocity, $v$ [m/s].
  81        - 'ThetaD': float or None
  82            Debye temperature, $\\Theta_D$ [K].
  83        - 'E0': float or None
  84            Equilibrium total energy per atom, $E_0$ [eV/atom].
  85
  86    Notes
  87    -----
  88    In the harmonic Debye approximation, vibrational contributions depend
  89    primarily on the Debye temperature $\\Theta_D$ (or equivalently the
  90    Debye frequency), which can be related to the elastic constants (here
  91    represented by $B$ and $\\nu$) and the density $\\rho$ through the mean
  92    sound velocity $v$. The model is commonly used to approximate the
  93    vibrational free energy
  94    $$
  95    F_\\mathrm{vib}(T) = k_B T
  96    \\left[
  97      3 N_\\mathrm{at} \\ln\\!\\left(1 - e^{-\\Theta_D/T}\\right)
  98      - 9 N_\\mathrm{at} \\frac{T}{\\Theta_D} \\int_0^{\\Theta_D/T} \\frac{x^3}{e^x - 1}\\,dx
  99    \\right],
 100    $$
 101    and the total Helmholtz free energy can be expressed as
 102    $F(T) = E_0 + F_\\mathrm{vib}(T)$ at fixed volume in the harmonic limit.
 103
 104    Unit conventions are per-atom (where applicable) to avoid ambiguity when
 105    comparing different structures or formula units.
 106    """
 107
 108    def __init__(self, at_per_fu: int = 1):
 109        self.data = {
 110            'V': None,        # float: Å^3/atom
 111            'rho': None,      # float: g/cm^3
 112            'Nat': at_per_fu, # int: atoms in primitive cell
 113            'B': None,        # float: GPa
 114            'nu': None,       # float: dimensionless
 115            'v': None,        # float: m/s
 116            'ThetaD': None,   # float: K
 117            'E0': None        # float: eV/atom
 118        }
 119
 120    # -------------------------
 121    # Volume V [Å^3/atom]
 122    # -------------------------
 123    @property
 124    def V(self) -> Optional[float]:
 125        """
 126        Volume per atom, $V$.
 127
 128        Returns
 129        -------
 130        float or None
 131            Volume per atom $V$ [Å^3/atom], or ``None`` if unset.
 132
 133        Notes
 134        -----
 135        - Use per-atom volume to avoid dependence on the choice of the
 136          crystallographic cell or formula unit size.
 137        """
 138        return self.data['V']
 139
 140    @V.setter
 141    def V(self, V: float) -> None:
 142        """
 143        Set the volume per atom.
 144
 145        Parameters
 146        ----------
 147        V : float
 148            Volume per atom $V$ [Å^3/atom].
 149
 150        Notes
 151        -----
 152        - No validation is performed here; consider enforcing $V > 0$ upstream.
 153        """
 154        self.data['V'] = V
 155
 156    # -------------------------
 157    # Density rho [g/cm^3]
 158    # -------------------------
 159    @property
 160    def rho(self) -> Optional[float]:
 161        """
 162        Mass density, $\\rho$.
 163
 164        Returns
 165        -------
 166        float or None
 167            Mass density $\\rho$ [g/cm^3], or ``None`` if unset.
 168        """
 169        return self.data['rho']
 170
 171    @rho.setter
 172    def rho(self, rho: float) -> None:
 173        """
 174        Set the mass density.
 175
 176        Parameters
 177        ----------
 178        rho : float
 179            Mass density $\\rho$ [g/cm^3].
 180
 181        Notes
 182        -----
 183        - In Debye-model workflows, $\\rho$ and elastic properties determine
 184          the mean sound velocity $v$ (and thus $\\Theta_D$).
 185        """
 186        self.data['rho'] = rho
 187
 188    # -------------------------
 189    # Number of atoms Nat [1]
 190    # -------------------------
 191    @property
 192    def Nat(self) -> int:
 193        """
 194        Number of atoms in the primitive cell, $N_\\mathrm{at}$.
 195
 196        Returns
 197        -------
 198        int
 199            $N_\\mathrm{at}$ [1].
 200        """
 201        return self.data['Nat']
 202
 203    @Nat.setter
 204    def Nat(self, Nat: int) -> None:
 205        """
 206        Set the number of atoms in the primitive cell.
 207
 208        Parameters
 209        ----------
 210        Nat : int
 211            $N_\\mathrm{at}$ [1].
 212
 213        Notes
 214        -----
 215        - Should be a positive integer.
 216        """
 217        self.data['Nat'] = Nat
 218
 219    # -------------------------
 220    # Bulk modulus B [GPa]
 221    # -------------------------
 222    @property
 223    def B(self) -> Optional[float]:
 224        """
 225        Bulk modulus, $B$.
 226
 227        Returns
 228        -------
 229        float or None
 230            Bulk modulus $B$ [GPa], or ``None`` if unset.
 231        """
 232        return self.data['B']
 233
 234    @B.setter
 235    def B(self, B: float) -> None:
 236        """
 237        Set the bulk modulus.
 238
 239        Parameters
 240        ----------
 241        B : float
 242            Bulk modulus $B$ [GPa].
 243        """
 244        self.data['B'] = B
 245
 246    # -------------------------
 247    # Poisson's ratio nu [1]
 248    # -------------------------
 249    @property
 250    def nu(self) -> Optional[float]:
 251        """
 252        Poisson's ratio, $\\nu$.
 253
 254        Returns
 255        -------
 256        float or None
 257            Poisson's ratio $\\nu$ [dimensionless], or ``None`` if unset.
 258
 259        Notes
 260        -----
 261        - For mechanically stable, isotropic media, $-1 < \\nu < 0.5$.
 262        """
 263        return self.data['nu']
 264
 265    @nu.setter
 266    def nu(self, nu: float) -> None:
 267        """
 268        Set Poisson's ratio.
 269
 270        Parameters
 271        ----------
 272        nu : float
 273            Poisson's ratio $\\nu$ [dimensionless].
 274        """
 275        self.data['nu'] = nu
 276
 277    # -------------------------
 278    # Mean sound velocity v [m/s]
 279    # -------------------------
 280    @property
 281    def v(self) -> Optional[float]:
 282        """
 283        Mean (Debye) sound velocity, $v$.
 284
 285        Returns
 286        -------
 287        float or None
 288            Mean sound velocity $v$ [m/s], or ``None`` if unset.
 289
 290        Notes
 291        -----
 292        - Often derived from elastic constants and density
 293          (e.g., via averages of longitudinal and transverse speeds).
 294        """
 295        return self.data['v']
 296
 297    @v.setter
 298    def v(self, v: float) -> None:
 299        """
 300        Set the mean (Debye) sound velocity.
 301
 302        Parameters
 303        ----------
 304        v : float
 305            Mean sound velocity $v$ [m/s].
 306        """
 307        self.data['v'] = v
 308
 309    # -------------------------
 310    # Debye temperature ThetaD [K]
 311    # -------------------------
 312    @property
 313    def ThetaD(self) -> Optional[float]:
 314        """
 315        Debye temperature, $\\Theta_D$.
 316
 317        Returns
 318        -------
 319        float or None
 320            Debye temperature $\\Theta_D$ [K], or ``None`` if unset.
 321
 322        Notes
 323        -----
 324        - $\\Theta_D$ can be computed from $v$ and number density.
 325        - In the harmonic Debye model, $\\Theta_D$ sets the vibrational
 326          energy scale that governs $F_\\mathrm{vib}(T)$ and $C_V(T)$.
 327        """
 328        return self.data['ThetaD']
 329
 330    @ThetaD.setter
 331    def ThetaD(self, ThetaD: float) -> None:
 332        """
 333        Set the Debye temperature.
 334
 335        Parameters
 336        ----------
 337        ThetaD : float
 338            Debye temperature $\\Theta_D$ [K].
 339        """
 340        self.data['ThetaD'] = ThetaD
 341
 342    # -------------------------
 343    # Equilibrium energy E0 [eV/atom]
 344    # -------------------------
 345    @property
 346    def E0(self) -> Optional[float]:
 347        """
 348        Equilibrium total energy per atom, $E_0$.
 349
 350        Returns
 351        -------
 352        float or None
 353            $E_0$ [eV/atom], or ``None`` if unset.
 354
 355        Notes
 356        -----
 357        - In the harmonic limit at fixed volume, the Helmholtz free energy is
 358          $F(T) = E_0 + F_\\mathrm{vib}(T)$.
 359        """
 360        return self.data['E0']
 361
 362    @E0.setter
 363    def E0(self, E0: float) -> None:
 364        """
 365        Set the equilibrium total energy per atom.
 366
 367        Parameters
 368        ----------
 369        E0 : float
 370            $E_0$ [eV/atom].
 371        """
 372        self.data['E0'] = E0
 373
 374    
 375    def calculate_mean_sound_velocity(self) -> None:
 376        """
 377        Compute and store the mean (Debye) sound velocity, v.
 378    
 379        The mean sound velocity is estimated from Poisson's ratio and the
 380        bulk modulus and density via:
 381        $$
 382        v = f(\nu)\,\sqrt{B/\rho},
 383        $$
 384        where
 385        $$
 386        f(\nu) =
 387        \left[
 388            \frac{3}{
 389                2\left(\tfrac{2}{3}\tfrac{1+\nu}{1-2\nu}\right)^{3/2}
 390              + \left(\tfrac{1}{3}\tfrac{1+\nu}{1-\nu}\right)^{3/2}
 391            }
 392        \right]^{1/3}.
 393        $$
 394    
 395        Parameters
 396        ----------
 397        None
 398    
 399        Returns
 400        -------
 401        None
 402            The computed mean sound velocity is stored in ``self.data['v']`` [m/s].
 403    
 404        Notes
 405        -----
 406        - Units:
 407          - Input: ``B`` in [GPa], ``rho`` in [g/cm^3], ``nu`` dimensionless.
 408          - Output: ``v`` in [m/s].
 409          - The factor ``1e3`` accounts for the conversion
 410            $\\sqrt{\\text{GPa}/(\\text{g cm}^{-3})} \\to \\text{m s}^{-1}$,
 411            since $1\\,\\text{GPa}=10^9\\,\\text{Pa}$ and
 412            $1\\,\\text{g cm}^{-3}=10^3\\,\\text{kg m}^{-3}$.
 413        - This method updates ``self.data['v']`` in place.
 414        - Requires that ``nu``, ``B``, and ``rho`` are set beforehand.
 415        """
 416        v = self.data['nu']
 417        fnu = (3/(2*(2/3*(1+v)/(1-2*v))**(3/2) + (1/3*(1+v)/(1-v))**(3/2)))**(1/3)
 418        self.data['v'] = 1e3 * fnu * (self.data['B']/self.data['rho'])**0.5
 419    
 420    
 421    def calculate_Debye_temperature(self) -> None:
 422        """
 423        Compute and store the Debye temperature, Θ_D.
 424    
 425        The Debye temperature is evaluated as
 426        $$
 427        \\Theta_D \\,=\\, \\frac{h}{k_B}\\, v \\,\\left(\\frac{3}{4\\pi\\,\\Omega}\\right)^{1/3},
 428        $$
 429        where $\\Omega = N_\\mathrm{at}\\,V$ is the primitive-cell volume and
 430        $v$ is the mean sound velocity.
 431    
 432        Parameters
 433        ----------
 434        None
 435    
 436        Returns
 437        -------
 438        None
 439            The computed Debye temperature is stored in ``self.data['ThetaD']`` [K].
 440    
 441        Notes
 442        -----
 443        - Units:
 444          - Input: ``V`` in [Å^3/atom], ``Nat`` dimensionless (atoms/cell),
 445            so $\\Omega$ is in [Å^3/cell]; ``v`` in [m/s].
 446          - Constants: ``h`` in [J·s], ``kB`` in [J·K^-1].
 447          - Output: ``ThetaD`` in [K].
 448          - The factor ``1e10`` converts [Å^-1] to [m^-1].
 449        - This method calls ``calculate_mean_sound_velocity()`` to ensure ``v`` is up to date.
 450        - Requires that ``V``, ``Nat``, and the inputs needed for ``v`` are set.
 451        """
 452        kB = 1.38064852e-23  # J K^-1  (m2 kg s-2 K-1)
 453        h = 6.62607015e-34   # J s     (m2 kg s^-1)
 454        from math import pi
 455    
 456        self.calculate_mean_sound_velocity()
 457        Omega = self.data['Nat'] * self.data['V']  # Å^3 per cell
 458        self.data['ThetaD'] = 1e10 * h/kB * (3/(4*pi*Omega))**(1/3) * self.data['v']
 459    
 460    
 461    def DebyeIntegral(self, TD: float | list[float], order: int = 3) -> float | list[float]:
 462        """
 463        Evaluate the Debye integral of order n.
 464    
 465        This computes the dimensionless Debye function
 466        $$
 467        D_n(x) \\,=\\, \\frac{n}{x^n} \\int_0^x \\frac{t^n}{e^t - 1}\\,dt,
 468        $$
 469        where typically $n=3$ and $x=\\Theta_D/T$.
 470    
 471        Parameters
 472        ----------
 473        TD : float or array-like
 474            Upper bound $x = \\Theta_D/T$ (dimensionless).
 475        order : int, optional
 476            Debye integral order $n$ (default is 3).
 477    
 478        Returns
 479        -------
 480        float or list of float
 481            Value(s) of $D_n(x)$ at the provided bound(s). Returns a float if
 482            ``TD`` is scalar, otherwise a list of floats.
 483    
 484        Notes
 485        -----
 486        - Numerical evaluation is performed via ``scipy.integrate.quad``.
 487        - This helper is used in the vibrational free-energy expression of
 488          the harmonic Debye model.
 489        """
 490    
 491        def __f(x):
 492            return x**order/(exp(x)-1)
 493    
 494        if hasattr(TD, '__len__'):
 495            return [order/(TDD**order) * integrate.quad(__f, 0, TDD)[0] for TDD in TD]
 496        else:
 497            return order/(TD**order) * integrate.quad(__f, 0, TD)[0]
 498    
 499    
 500    def Fvib_harm(self, T: float) -> float:
 501        """
 502        Vibrational Helmholtz free energy in the harmonic Debye model, per atom.
 503    
 504        The vibrational contribution is computed as
 505        $$
 506        F_\\mathrm{vib}(T)
 507          \\,=\\, \\frac{9}{8} k_B \\Theta_D
 508          \\, - \\, k_B T
 509          \\left[
 510            D_3\\!\\left(\\frac{\\Theta_D}{T}\\right)
 511            \\, -\\, 3\\ln\\!\\left(1 - e^{-\\Theta_D/T}\\right)
 512          \\right],
 513        $$
 514        where $D_3(x)$ is the Debye function of order 3.
 515    
 516        Parameters
 517        ----------
 518        T : float
 519            Temperature $T$ [K].
 520    
 521        Returns
 522        -------
 523        float
 524            Vibrational free energy per atom, $F_\\mathrm{vib}(T)$ [eV/atom].
 525    
 526        Notes
 527        -----
 528        - Uses ``kB = 8.617333262145e-5`` eV/K.
 529        - Calls ``calculate_Debye_temperature()`` internally; ensure the material
 530          parameters needed for $\\Theta_D$ are set.
 531        - The result is per atom, consistent with the class unit conventions.
 532        """
 533        kB = 8.617333262145e-5  # eV K^-1
 534        from numpy import log, exp
 535    
 536        self.calculate_Debye_temperature()
 537        ThetaD = self.data['ThetaD']
 538        debInt = self.DebyeIntegral(ThetaD / T)
 539    
 540        return 9/8 * kB * ThetaD - kB * T * (debInt - 3 * log(1 - exp(-ThetaD / T)))
 541    
 542    
 543    def F_harm(self, T: float) -> float:
 544        """
 545        Total Helmholtz free energy in the harmonic Debye model, per atom.
 546    
 547        The total free energy at fixed volume in the harmonic limit is
 548        $$
 549        F(T) \\,=\\, E_0 \\, + \\, F_\\mathrm{vib}(T).
 550        $$
 551    
 552        Parameters
 553        ----------
 554        T : float
 555            Temperature $T$ [K].
 556    
 557        Returns
 558        -------
 559        float
 560            Total Helmholtz free energy per atom, $F(T)$ [eV/atom].
 561    
 562        Notes
 563        -----
 564        - Requires that ``E0`` is set (equilibrium energy per atom), and that the
 565          parameters needed to compute $\\Theta_D$ are available.
 566        - Internally calls ``Fvib_harm(T)``.
 567        """
 568        return self.data['E0'] + self.Fvib_harm(T)
 569
 570
 571
 572
 573class DebyeQHA:
 574    """
 575    Quasi-harmonic Debye model container for equation-of-state parameters.
 576
 577    This class stores per-atom reference quantities and material information
 578    needed for quasi-harmonic (volume-dependent) thermodynamic calculations,
 579    typically combining the Debye model for vibrational contributions with an
 580    equation of state (e.g., Birch–Murnaghan) for the static lattice energy.
 581
 582    Parameters
 583    ----------
 584    E0 : float
 585        Reference (equilibrium) energy per atom, $E_0$ [eV/atom].
 586    V0 : float
 587        Reference (equilibrium) volume per atom, $V_0$ [Å^3/atom].
 588    B0 : float
 589        Bulk modulus at $V_0$, $B_0$ [GPa].
 590    Bp : float
 591        Pressure derivative of the bulk modulus at $V_0$, $B'_0$ [dimensionless].
 592    material : str or ase.Atoms
 593        Chemical formula string (e.g., 'Al', 'Fe2O3') or an ASE ``Atoms`` object
 594        describing the material for which thermodynamic properties are computed.
 595    at_per_fu : int, optional
 596        Number of atoms in the primitive cell, $N_\\mathrm{at}$ [1], by default 1.
 597
 598    Attributes
 599    ----------
 600    data : dict
 601        Internal storage for model parameters:
 602        - 'E0': float
 603            $E_0$ [eV/atom].
 604        - 'V0': float
 605            $V_0$ [Å^3/atom].
 606        - 'B0': float
 607            $B_0$ [GPa].
 608        - 'Bp': float
 609            $B'_0$ [1].
 610        - 'material': str or ase.Atoms
 611            Material descriptor (formula or ``Atoms``).
 612        - 'Nat': int
 613            $N_\\mathrm{at}$ [1].
 614        - 'volumes': sequence of float or None
 615            Mesh of volumes per atom, $\\{V\\}$ [Å^3/atom], for QHA scans.
 616
 617    Notes
 618    -----
 619    - Units are per atom throughout to keep a consistent convention.
 620    - The pair $(B_0, B'_0)$ is commonly used with the third-order Birch–Murnaghan
 621      equation of state to describe $E(V)$ around $V_0$.
 622    """
 623
 624    def __init__(
 625        self,
 626        E0: float,
 627        V0: float,
 628        B0: float,
 629        Bp: float,
 630        material: Union[str, Atoms],
 631        at_per_fu: int = 1
 632    ) -> None:
 633        """
 634        Initialize the DebyeQHA container with equation-of-state parameters.
 635
 636        See class docstring for a detailed description of parameters and units.
 637        """
 638        self.data = {
 639            'E0': E0,           # float: eV/atom
 640            'V0': V0,           # float: Å^3/atom
 641            'B0': B0,           # float: GPa
 642            'Bp': Bp,           # float: 1 (dimensionless)
 643            'nu': None,         # float: 1 (dimensionless)
 644            'material': material,
 645            'Nat': at_per_fu,   # int: atoms in primitive cell
 646            'volumes': None     # Sequence[float]: Å^3/atom
 647        }
 648        self._deb = []          # array of Debye (HA) objects
 649
 650    # Nat [atoms in primitive cell]
 651    @property
 652    def Nat(self) -> int:
 653        """
 654        Number of atoms in the primitive cell, $N_\\mathrm{at}$.
 655    
 656        Returns
 657        -------
 658        int
 659            $N_\\mathrm{at}$ [1].
 660        """
 661        return self.data['Nat']
 662    
 663    @Nat.setter
 664    def Nat(self, Nat: int) -> None:
 665        """
 666        Set the number of atoms in the primitive cell.
 667    
 668        Parameters
 669        ----------
 670        Nat : int
 671            $N_\\mathrm{at}$ [1].
 672    
 673        Returns
 674        -------
 675        None
 676    
 677        Raises
 678        ------
 679        ValueError
 680            If ``Nat`` is not a positive integer.
 681        """
 682        if not isinstance(Nat, int) or Nat <= 0:
 683            raise ValueError("Nat must be a positive integer")
 684        self.data['Nat'] = Nat
 685
 686    # E0 [eV/atom]
 687    @property
 688    def E0(self) -> float:
 689        """
 690        Equilibrium energy per atom, $E_0$.
 691    
 692        Returns
 693        -------
 694        float
 695            $E_0$ [eV/atom].
 696        """
 697        return self.data['E0']
 698    
 699    @E0.setter
 700    def E0(self, E0: float) -> None:
 701        """
 702        Set the equilibrium energy per atom.
 703    
 704        Parameters
 705        ----------
 706        E0 : float
 707            $E_0$ [eV/atom].
 708    
 709        Returns
 710        -------
 711        None
 712        """
 713        self.data['E0'] = float(E0)
 714
 715
 716    # V0 [Å^3/atom]
 717    @property
 718    def V0(self) -> float:
 719        """
 720        Equilibrium volume per atom, $V_0$.
 721    
 722        Returns
 723        -------
 724        float
 725            $V_0$ [Å^3/atom].
 726        """
 727        return self.data['V0']
 728    
 729    @V0.setter
 730    def V0(self, V0: float) -> None:
 731        """
 732        Set the equilibrium volume per atom.
 733    
 734        Parameters
 735        ----------
 736        V0 : float
 737            $V_0$ [Å^3/atom].
 738    
 739        Returns
 740        -------
 741        None
 742        """
 743        self.data['V0'] = float(V0)
 744
 745
 746    # B0 [GPa]
 747    @property
 748    def B0(self) -> float:
 749        """
 750        Bulk modulus at $V_0$, $B_0$.
 751    
 752        Returns
 753        -------
 754        float
 755            $B_0$ [GPa].
 756        """
 757        return self.data['B0']
 758    
 759    @B0.setter
 760    def B0(self, B0: float) -> None:
 761        """
 762        Set the bulk modulus at $V_0$.
 763    
 764        Parameters
 765        ----------
 766        B0 : float
 767            $B_0$ [GPa].
 768    
 769        Returns
 770        -------
 771        None
 772        """
 773        self.data['B0'] = float(B0)
 774
 775
 776    # Bp [dimensionless]
 777    @property
 778    def Bp(self) -> float:
 779        """
 780        Pressure derivative of the bulk modulus at $V_0$, $B'_0$.
 781    
 782        Returns
 783        -------
 784        float
 785            $B'_0$ [dimensionless].
 786        """
 787        return self.data['Bp']
 788    
 789    @Bp.setter
 790    def Bp(self, Bp: float) -> None:
 791        """
 792        Set the pressure derivative of the bulk modulus at $V_0$.
 793    
 794        Parameters
 795        ----------
 796        Bp : float
 797            $B'_0$ [dimensionless].
 798    
 799        Returns
 800        -------
 801        None
 802        """
 803        self.data['Bp'] = float(Bp)
 804
 805
 806    # material [formula or ASE Atoms]
 807    @property
 808    def material(self) -> Union[str, Atoms]:
 809        """
 810        Material descriptor.
 811    
 812        Returns
 813        -------
 814        str or ase.Atoms
 815            Chemical formula string or an ASE ``Atoms`` object.
 816        """
 817        return self.data['material']
 818    
 819    @material.setter
 820    def material(self, material: Union[str, Atoms]) -> None:
 821        """
 822        Set the material descriptor.
 823    
 824        Parameters
 825        ----------
 826        material : str or ase.Atoms
 827            Chemical formula string (e.g., 'Al', 'Fe2O3') or an ASE ``Atoms`` object.
 828    
 829        Returns
 830        -------
 831        None
 832    
 833        Raises
 834        ------
 835        TypeError
 836            If ``material`` is not a string or ``ase.Atoms`` instance.
 837        """
 838        if not isinstance(material, (str, Atoms)):
 839            raise TypeError("material must be a chemical formula string or an ase.Atoms object")
 840        self.data['material'] = material
 841
 842    # Poisson's ratio
 843    @property
 844    def nu(self) -> Optional[float]:
 845        """
 846        Poisson's ratio, $\\nu$.
 847
 848        Returns
 849        -------
 850        float or None
 851            Poisson's ratio $\\nu$ [dimensionless], or ``None`` if unset.
 852
 853        Notes
 854        -----
 855        - For mechanically stable, isotropic media, $-1 < \\nu < 0.5$.
 856        """
 857        return self.data['nu']
 858
 859    @nu.setter
 860    def nu(self, nu: float) -> None:
 861        """
 862        Set Poisson's ratio.
 863
 864        Parameters
 865        ----------
 866        nu : float
 867            Poisson's ratio $\\nu$ [dimensionless].
 868        """
 869        self.data['nu'] = nu
 870
 871    # volumes [Å^3/atom] (array-like)
 872    @property
 873    def volumes(self) -> np.ndarray:
 874        """
 875        Volume mesh per atom, {V}.
 876    
 877        Returns
 878        -------
 879        numpy.ndarray or None
 880            1D array of volumes per atom [Å^3/atom], or ``None`` if unset.
 881    
 882        Notes
 883        -----
 884        - Used to sample the quasi-harmonic free energy across different volumes.
 885        - Stored internally in ``self.data['volumes']``.
 886        """
 887        return self.data['volumes']
 888    
 889    
 890    @volumes.setter
 891    def volumes(self, volumes: Any) -> None:
 892        """
 893        Set the volume mesh per atom.
 894    
 895        Parameters
 896        ----------
 897        volumes : array-like
 898            Any array-like of volumes per atom [Å^3/atom]. Examples include
 899            list, tuple, numpy.ndarray, pandas Series.
 900    
 901        Returns
 902        -------
 903        None
 904    
 905        Raises
 906        ------
 907        ValueError
 908            If the provided mesh cannot be converted to a 1D float array or is empty.
 909    
 910        Notes
 911        -----
 912        - The input is converted to a 1D numpy array of dtype float and stored in
 913          ``self.data['volumes']``.
 914        - Ensure values are per atom and, typically, positive.
 915        """
 916        arr = np.asarray(volumes, dtype=float)
 917        if arr.ndim != 1 or arr.size == 0:
 918            raise ValueError("volumes must be a non-empty 1D array-like of floats.")
 919        self.data['volumes'] = arr
 920        self._deb = []
 921        for V in self.data['volumes']:
 922            d = Debye()
 923            d.nu = self.nu
 924            d.V = V
 925            d.E0 = self.get_E0(V)
 926            d.B = self.get_bulk_modulus(V)
 927            d.rho = get_rho(V, self.material)
 928            self._deb.append(deepcopy(d))
 929
 930    
 931    def get_E0(self, V: float) -> float:
 932        """
 933        Energy at volume from stored Birch–Murnaghan parameters.
 934    
 935        Evaluates the energy per atom at a given per-atom volume using the
 936        Birch–Murnaghan equation of state with the parameters stored in the instance.
 937    
 938        Parameters
 939        ----------
 940        V : float
 941            Volume per atom, $V$ [Å^3/atom].
 942    
 943        Returns
 944        -------
 945        float
 946            Energy per atom, $E(V)$ [eV/atom].
 947    
 948        Notes
 949        -----
 950        - Uses the instance's ``E0``, ``V0``, ``B0`` (in GPa), and ``Bp`` values.
 951        - Delegates to ``BM_EoS(V, E0, V0, B0, Bp)``.
 952        """
 953        return self.BM_EoS(V, self.E0, self.V0, self.B0, self.Bp)
 954    
 955    
 956    def BM_EoS(self, V: float, E0: float, V0: float, B0: float, B0p: float) -> float:
 957        """
 958        Third-order Birch–Murnaghan equation of state, energy vs. volume per atom.
 959    
 960        Computes the per-atom energy as a function of per-atom volume using:
 961        $$
 962        E(V) \;=\; E_0 \;+\; \frac{9\,V_0\,B_0}{16}
 963        \\left\\{
 964          \\left[\\left(\\frac{V_0}{V}\\right)^{2/3}-1\\right]^3 B'_0
 965          +
 966          \\left[\\left(\\frac{V_0}{V}\\right)^{2/3}-1\\right]^2
 967          \\left[6 - 4 \\left(\\frac{V_0}{V}\\right)^{2/3}\\right]
 968        \\right\\}.
 969        $$
 970    
 971        Parameters
 972        ----------
 973        V : float
 974            Volume per atom, $V$ [Å^3/atom].
 975        E0 : float
 976            Reference energy per atom, $E_0$ [eV/atom].
 977        V0 : float
 978            Reference volume per atom, $V_0$ [Å^3/atom].
 979        B0 : float
 980            Bulk modulus at $V_0$, $B_0$ [GPa].
 981        B0p : float
 982            Pressure derivative of the bulk modulus at $V_0$, $B'_0$ [dimensionless].
 983    
 984        Returns
 985        -------
 986        float
 987            Energy per atom, $E(V)$ [eV/atom].
 988    
 989        Notes
 990        -----
 991        - If a consistent energy unit is required, convert $B_0$ to [eV/Å^3] via
 992          $1\\,\\text{eV}/\\text{Å}^3 = 160.21766208\\,\\text{GPa}$ before use.
 993        """
 994        return E0 + 9*V0*B0/16 * (((V0/V)**(2/3)-1)**3 * B0p + ((V0/V)**(2/3)-1)**2* (6-4*(V0/V)**(2/3)))
 995    
 996    
 997    def fit_BM_EoS(self, vols, enes):
 998        """
 999        Fit Birch–Murnaghan (third-order) EoS parameters to energy–volume data.
1000    
1001        Performs a nonlinear least-squares fit of the BM equation of state to the
1002        provided per-atom volume and energy data, returning fitted parameters.
1003    
1004        Parameters
1005        ----------
1006        vols : array-like
1007            Volumes per atom, $\\{V\\}$ [Å^3/atom].
1008        enes : array-like
1009            Energies per atom, $\\{E\\}$ [eV/atom].
1010    
1011        Returns
1012        -------
1013        numpy.ndarray
1014            Array of fitted parameters with shape (4,), ordered as:
1015            [E0, V0, B0, Bp], where
1016            - E0: float, E0 [eV/atom]
1017            - V0: float, V0 [Å^3/atom]
1018            - B0: float, B0 [GPa]
1019            - Bp: float, B0' [dimensionless]
1020    
1021        Notes
1022        -----
1023        - Initial guesses used for the fit are:
1024          $E_0 = \\min(E)$,
1025          $V_0 = \\tfrac{1}{2}(\\min(V)+\\max(V))$,
1026          $B_0 \\approx 1.5\\,\\text{eV}/\\text{Å}^3$,
1027          $B'_0 = 4$.
1028        - Converts the fitted bulk modulus from [eV/Å^3] to [GPa] using
1029          $1\\,\\text{eV}/\\text{Å}^3 \\approx 160.2\\,\\text{GPa}$.
1030        - Assumes a callable ``EoS(V, E0, V0, B, Bp)`` compatible with ``curve_fit``.
1031        """
1032        eosfit, _ = curve_fit(self.BM_EoS, vols, enes, [min(enes), 0.5*(min(vols)+max(vols)), 1.5, 4])
1033        eosfit[2] *= 160.2
1034        return eosfit
1035
1036        
1037    import numpy as np
1038    from numpy.typing import ArrayLike, NDArray
1039    
1040    def calc_equilibrium(self, T: ArrayLike) -> NDArray[np.floating]:
1041        """
1042        Fit BM3 EoS at each temperature and return the fitted parameters.
1043    
1044        For each temperature in T, computes the harmonic free energies across the
1045        stored volume mesh (via self._deb[i].F_harm(T_i)), fits the third-order
1046        Birch–Murnaghan EoS to E(V), and returns the fitted parameters.
1047    
1048        Parameters
1049        ----------
1050        T : array-like
1051            Temperatures, {T} [K]. Can be a scalar or array-like.
1052    
1053        Returns
1054        -------
1055        numpy.ndarray
1056            - If T is scalar: shape (4,) ordered as [E0, V0, B0, Bp]
1057            - If T is array-like: shape (N, 4) where N = len(T)
1058              ordered as [E0, V0, B0, Bp] for each temperature.
1059    
1060        Notes
1061        -----
1062        - Assumes self.volumes is the per-atom volume mesh [Å^3/atom].
1063        - Assumes self._deb is an iterable of Debye objects aligned with self.volumes,
1064          each providing F_harm(T) -> energy [eV/atom] at temperature T.
1065        """
1066        is_scalar = np.isscalar(T) or (np.ndim(T) == 0)
1067        T_arr = np.atleast_1d(T).astype(float)
1068    
1069        fits = []
1070        for t in T_arr:
1071            enes = [d.F_harm(t) for d in self._deb]
1072            eosfit = self.fit_BM_EoS(self.volumes, enes)  # returns [E0, V0, B0, Bp]
1073            fits.append(eosfit)
1074    
1075        fits = np.asarray(fits, dtype=float)
1076        return fits[0] if is_scalar else fits
1077    
1078    
1079    def calc_F(self, T: ArrayLike) -> NDArray[np.floating] | float:
1080        """
1081        Equilibrium energy E0(T) from BM3 fit at each temperature.
1082    
1083        Parameters
1084        ----------
1085        T : array-like
1086            Temperatures, {T} [K]. Can be a scalar or array-like.
1087    
1088        Returns
1089        -------
1090        float or numpy.ndarray
1091            - If T is scalar: scalar float E0(T) [eV/atom]
1092            - If T is array-like: array of E0(T) [eV/atom] with shape (N,)
1093        """
1094        eosfits = self.calc_equilibrium(T)
1095        if eosfits.ndim == 1:
1096            return float(eosfits[0])
1097        return eosfits[:, 0]
1098    
1099    
1100    def calc_V(self, T: ArrayLike) -> NDArray[np.floating] | float:
1101        """
1102        Equilibrium volume V0(T) from BM3 fit at each temperature.
1103    
1104        Parameters
1105        ----------
1106        T : array-like
1107            Temperatures, {T} [K]. Can be a scalar or array-like.
1108    
1109        Returns
1110        -------
1111        float or numpy.ndarray
1112            - If T is scalar: scalar float V0(T) [Å^3/atom]
1113            - If T is array-like: array of V0(T) [Å^3/atom] with shape (N,)
1114        """
1115        eosfits = self.calc_equilibrium(T)
1116        if eosfits.ndim == 1:
1117            return float(eosfits[1])
1118        return eosfits[:, 1]
1119    
1120    
1121    def calc_B(self, T: ArrayLike) -> NDArray[np.floating] | float:
1122        """
1123        Equilibrium bulk modulus B0(T) from BM3 fit at each temperature.
1124    
1125        Parameters
1126        ----------
1127        T : array-like
1128            Temperatures, {T} [K]. Can be a scalar or array-like.
1129    
1130        Returns
1131        -------
1132        float or numpy.ndarray
1133            - If T is scalar: scalar float B0(T) [GPa]
1134            - If T is array-like: array of B0(T) [GPa] with shape (N,)
1135        """
1136        eosfits = self.calc_equilibrium(T)
1137        if eosfits.ndim == 1:
1138            return float(eosfits[2])
1139        return eosfits[:, 2]
1140
1141    
1142    def get_pressure(self, V: float) -> float:
1143        """
1144        Third-order Birch–Murnaghan pressure, P(V), as a function of per-atom volume.
1145    
1146        Computes the pressure using the BM equation of state:
1147        $$
1148        P(V) \;=\; \frac{3}{2} B_0
1149          \left[
1150            \left(\frac{V_0}{V}\right)^{7/3} - \left(\frac{V_0}{V}\right)^{5/3}
1151          \right]
1152          \left\{
1153            1 + \frac{3}{4}\,(B'_0 - 4)\left[\left(\frac{V_0}{V}\right)^{2/3} - 1\right]
1154          \right\}.
1155        $$
1156    
1157        Parameters
1158        ----------
1159        V : float
1160            Volume per atom, $V$ [Å^3/atom].
1161    
1162        Returns
1163        -------
1164        float
1165            Pressure, $P(V)$ [GPa].
1166    
1167        Notes
1168        -----
1169        - Units:
1170          - ``V0`` in [Å^3/atom], ``B0`` in [GPa], ``Bp`` dimensionless.
1171          - The returned pressure is in [GPa], consistent with ``B0``.
1172        """
1173        V0 = self.V0
1174        B0 = self.B0
1175        Bp = self.Bp
1176    
1177        eta = V0 / V
1178        return (3.0 * B0 / 2.0) * (eta**(7.0/3.0) - eta**(5.0/3.0)) * (1.0 + 0.75 * (Bp - 4.0) * (eta**(2.0/3.0) - 1.0))
1179
1180
1181    def get_bulk_modulus(self, V: float) -> float:
1182        """
1183        Volume-dependent (isothermal) bulk modulus, B(V), in the BM3 framework.
1184    
1185        Uses the linearized relation based on the definition of the pressure derivative
1186        at zero pressure:
1187        $$
1188        B(V) \;\approx\; B_0 \;+\; B'_0\,P(V),
1189        $$
1190        where $P(V)$ is the Birch–Murnaghan (third-order) pressure and $B'_0 = (dB/dP)_{P=0}$.
1191    
1192        Parameters
1193        ----------
1194        V : float
1195            Volume per atom, $V$ [Å^3/atom].
1196    
1197        Returns
1198        -------
1199        float
1200            Bulk modulus, $B(V)$ [GPa].
1201    
1202        Notes
1203        -----
1204        - This is a linear approximation in pressure around $P=0$; for large compressions
1205          or expansions, a full BM3 expression for $B(V)$ may be preferred.
1206        - Units: ``B0`` in [GPa], ``Bp`` dimensionless, returned $B(V)$ in [GPa].
1207        """
1208        p = self.get_pressure(V)
1209        return self.B0 + self.Bp * p
def get_rho(V: float, material: Union[str, ase.atoms.Atoms]) -> float: 
18def get_rho(V: float, material: Union[str, Atoms]) -> float:
19    """
20    Compute mass density rho [g/cm^3] from per-atom volume V [Å^3/atom]
21    and either an ASE Atoms object or a chemical formula.
22
23    Parameters
24    ----------
25    V : float
26        Per-atom volume in Å^3/atom.
27    material : str or ase.Atoms
28        Chemical formula string (e.g., 'Fe2O3', 'Si', 'AlN') or an ASE Atoms object.
29
30    Returns
31    -------
32    float
33        Mass density in g/cm^3.
34
35    Notes
36    -----
37    - Conversions: 1 amu = 1.66054e-24 g and 1 Å^3 = 1e-24 cm^3,
38      so rho = 1.66054 * (average_mass_in_amu) / V.
39    """
40    if isinstance(material, str):
41        atoms = Atoms(material)
42    elif isinstance(material, Atoms):
43        atoms = material
44    else:
45        raise TypeError("material must be a chemical formula string or an ase.Atoms object")
46
47    M_amu = atoms.get_masses().mean()  # average mass per atom (amu)
48    return 1.66054 * M_amu / V

Compute mass density rho [g/cm^3] from per-atom volume V [Å^3/atom] and either an ASE Atoms object or a chemical formula.

Parameters
  • V (float): Per-atom volume in Å^3/atom.
  • material (str or ase.Atoms): Chemical formula string (e.g., 'Fe2O3', 'Si', 'AlN') or an ASE Atoms object.
Returns
  • float: Mass density in g/cm^3.
Notes
  • Conversions: 1 amu = 1.66054e-24 g and 1 Å^3 = 1e-24 cm^3, so rho = 1.66054 * (average_mass_in_amu) / V.
class Debye: 
 53class Debye:
 54    """
 55    Implementation of the harmonic Debye approximation.
 56
 57    This class stores and provides accessors for thermodynamic and elastic
 58    descriptors required to evaluate properties in the harmonic Debye model,
 59    such as the Helmholtz free energy, heat capacity, and entropy.
 60
 61    Parameters
 62    ----------
 63    at_per_fu : int, optional
 64        Number of atoms in the primitive cell, $N_\\mathrm{at}$, by default 1.
 65
 66    Attributes
 67    ----------
 68    data : dict
 69        Internal storage for material fields:
 70        - 'V': float or None
 71            Volume per atom, $V$ [Å^3/atom].
 72        - 'rho': float or None
 73            Mass density, $\\rho$ [g/cm^3].
 74        - 'Nat': int
 75            Number of atoms in the primitive cell, $N_\\mathrm{at}$ [1].
 76        - 'B': float or None
 77            Bulk modulus, $B$ [GPa].
 78        - 'nu': float or None
 79            Poisson's ratio, $\\nu$ [1].
 80        - 'v': float or None
 81            Mean (Debye) sound velocity, $v$ [m/s].
 82        - 'ThetaD': float or None
 83            Debye temperature, $\\Theta_D$ [K].
 84        - 'E0': float or None
 85            Equilibrium total energy per atom, $E_0$ [eV/atom].
 86
 87    Notes
 88    -----
 89    In the harmonic Debye approximation, vibrational contributions depend
 90    primarily on the Debye temperature $\\Theta_D$ (or equivalently the
 91    Debye frequency), which can be related to the elastic constants (here
 92    represented by $B$ and $\\nu$) and the density $\\rho$ through the mean
 93    sound velocity $v$. The model is commonly used to approximate the
 94    vibrational free energy
 95    $$
 96    F_\\mathrm{vib}(T) = k_B T
 97    \\left[
 98      3 N_\\mathrm{at} \\ln\\!\\left(1 - e^{-\\Theta_D/T}\\right)
 99      - 9 N_\\mathrm{at} \\frac{T}{\\Theta_D} \\int_0^{\\Theta_D/T} \\frac{x^3}{e^x - 1}\\,dx
100    \\right],
101    $$
102    and the total Helmholtz free energy can be expressed as
103    $F(T) = E_0 + F_\\mathrm{vib}(T)$ at fixed volume in the harmonic limit.
104
105    Unit conventions are per-atom (where applicable) to avoid ambiguity when
106    comparing different structures or formula units.
107    """
108
109    def __init__(self, at_per_fu: int = 1):
110        self.data = {
111            'V': None,        # float: Å^3/atom
112            'rho': None,      # float: g/cm^3
113            'Nat': at_per_fu, # int: atoms in primitive cell
114            'B': None,        # float: GPa
115            'nu': None,       # float: dimensionless
116            'v': None,        # float: m/s
117            'ThetaD': None,   # float: K
118            'E0': None        # float: eV/atom
119        }
120
121    # -------------------------
122    # Volume V [Å^3/atom]
123    # -------------------------
124    @property
125    def V(self) -> Optional[float]:
126        """
127        Volume per atom, $V$.
128
129        Returns
130        -------
131        float or None
132            Volume per atom $V$ [Å^3/atom], or ``None`` if unset.
133
134        Notes
135        -----
136        - Use per-atom volume to avoid dependence on the choice of the
137          crystallographic cell or formula unit size.
138        """
139        return self.data['V']
140
141    @V.setter
142    def V(self, V: float) -> None:
143        """
144        Set the volume per atom.
145
146        Parameters
147        ----------
148        V : float
149            Volume per atom $V$ [Å^3/atom].
150
151        Notes
152        -----
153        - No validation is performed here; consider enforcing $V > 0$ upstream.
154        """
155        self.data['V'] = V
156
157    # -------------------------
158    # Density rho [g/cm^3]
159    # -------------------------
160    @property
161    def rho(self) -> Optional[float]:
162        """
163        Mass density, $\\rho$.
164
165        Returns
166        -------
167        float or None
168            Mass density $\\rho$ [g/cm^3], or ``None`` if unset.
169        """
170        return self.data['rho']
171
172    @rho.setter
173    def rho(self, rho: float) -> None:
174        """
175        Set the mass density.
176
177        Parameters
178        ----------
179        rho : float
180            Mass density $\\rho$ [g/cm^3].
181
182        Notes
183        -----
184        - In Debye-model workflows, $\\rho$ and elastic properties determine
185          the mean sound velocity $v$ (and thus $\\Theta_D$).
186        """
187        self.data['rho'] = rho
188
189    # -------------------------
190    # Number of atoms Nat [1]
191    # -------------------------
192    @property
193    def Nat(self) -> int:
194        """
195        Number of atoms in the primitive cell, $N_\\mathrm{at}$.
196
197        Returns
198        -------
199        int
200            $N_\\mathrm{at}$ [1].
201        """
202        return self.data['Nat']
203
204    @Nat.setter
205    def Nat(self, Nat: int) -> None:
206        """
207        Set the number of atoms in the primitive cell.
208
209        Parameters
210        ----------
211        Nat : int
212            $N_\\mathrm{at}$ [1].
213
214        Notes
215        -----
216        - Should be a positive integer.
217        """
218        self.data['Nat'] = Nat
219
220    # -------------------------
221    # Bulk modulus B [GPa]
222    # -------------------------
223    @property
224    def B(self) -> Optional[float]:
225        """
226        Bulk modulus, $B$.
227
228        Returns
229        -------
230        float or None
231            Bulk modulus $B$ [GPa], or ``None`` if unset.
232        """
233        return self.data['B']
234
235    @B.setter
236    def B(self, B: float) -> None:
237        """
238        Set the bulk modulus.
239
240        Parameters
241        ----------
242        B : float
243            Bulk modulus $B$ [GPa].
244        """
245        self.data['B'] = B
246
247    # -------------------------
248    # Poisson's ratio nu [1]
249    # -------------------------
250    @property
251    def nu(self) -> Optional[float]:
252        """
253        Poisson's ratio, $\\nu$.
254
255        Returns
256        -------
257        float or None
258            Poisson's ratio $\\nu$ [dimensionless], or ``None`` if unset.
259
260        Notes
261        -----
262        - For mechanically stable, isotropic media, $-1 < \\nu < 0.5$.
263        """
264        return self.data['nu']
265
266    @nu.setter
267    def nu(self, nu: float) -> None:
268        """
269        Set Poisson's ratio.
270
271        Parameters
272        ----------
273        nu : float
274            Poisson's ratio $\\nu$ [dimensionless].
275        """
276        self.data['nu'] = nu
277
278    # -------------------------
279    # Mean sound velocity v [m/s]
280    # -------------------------
281    @property
282    def v(self) -> Optional[float]:
283        """
284        Mean (Debye) sound velocity, $v$.
285
286        Returns
287        -------
288        float or None
289            Mean sound velocity $v$ [m/s], or ``None`` if unset.
290
291        Notes
292        -----
293        - Often derived from elastic constants and density
294          (e.g., via averages of longitudinal and transverse speeds).
295        """
296        return self.data['v']
297
298    @v.setter
299    def v(self, v: float) -> None:
300        """
301        Set the mean (Debye) sound velocity.
302
303        Parameters
304        ----------
305        v : float
306            Mean sound velocity $v$ [m/s].
307        """
308        self.data['v'] = v
309
310    # -------------------------
311    # Debye temperature ThetaD [K]
312    # -------------------------
313    @property
314    def ThetaD(self) -> Optional[float]:
315        """
316        Debye temperature, $\\Theta_D$.
317
318        Returns
319        -------
320        float or None
321            Debye temperature $\\Theta_D$ [K], or ``None`` if unset.
322
323        Notes
324        -----
325        - $\\Theta_D$ can be computed from $v$ and number density.
326        - In the harmonic Debye model, $\\Theta_D$ sets the vibrational
327          energy scale that governs $F_\\mathrm{vib}(T)$ and $C_V(T)$.
328        """
329        return self.data['ThetaD']
330
331    @ThetaD.setter
332    def ThetaD(self, ThetaD: float) -> None:
333        """
334        Set the Debye temperature.
335
336        Parameters
337        ----------
338        ThetaD : float
339            Debye temperature $\\Theta_D$ [K].
340        """
341        self.data['ThetaD'] = ThetaD
342
343    # -------------------------
344    # Equilibrium energy E0 [eV/atom]
345    # -------------------------
346    @property
347    def E0(self) -> Optional[float]:
348        """
349        Equilibrium total energy per atom, $E_0$.
350
351        Returns
352        -------
353        float or None
354            $E_0$ [eV/atom], or ``None`` if unset.
355
356        Notes
357        -----
358        - In the harmonic limit at fixed volume, the Helmholtz free energy is
359          $F(T) = E_0 + F_\\mathrm{vib}(T)$.
360        """
361        return self.data['E0']
362
363    @E0.setter
364    def E0(self, E0: float) -> None:
365        """
366        Set the equilibrium total energy per atom.
367
368        Parameters
369        ----------
370        E0 : float
371            $E_0$ [eV/atom].
372        """
373        self.data['E0'] = E0
374
375    
376    def calculate_mean_sound_velocity(self) -> None:
377        """
378        Compute and store the mean (Debye) sound velocity, v.
379    
380        The mean sound velocity is estimated from Poisson's ratio and the
381        bulk modulus and density via:
382        $$
383        v = f(\nu)\,\sqrt{B/\rho},
384        $$
385        where
386        $$
387        f(\nu) =
388        \left[
389            \frac{3}{
390                2\left(\tfrac{2}{3}\tfrac{1+\nu}{1-2\nu}\right)^{3/2}
391              + \left(\tfrac{1}{3}\tfrac{1+\nu}{1-\nu}\right)^{3/2}
392            }
393        \right]^{1/3}.
394        $$
395    
396        Parameters
397        ----------
398        None
399    
400        Returns
401        -------
402        None
403            The computed mean sound velocity is stored in ``self.data['v']`` [m/s].
404    
405        Notes
406        -----
407        - Units:
408          - Input: ``B`` in [GPa], ``rho`` in [g/cm^3], ``nu`` dimensionless.
409          - Output: ``v`` in [m/s].
410          - The factor ``1e3`` accounts for the conversion
411            $\\sqrt{\\text{GPa}/(\\text{g cm}^{-3})} \\to \\text{m s}^{-1}$,
412            since $1\\,\\text{GPa}=10^9\\,\\text{Pa}$ and
413            $1\\,\\text{g cm}^{-3}=10^3\\,\\text{kg m}^{-3}$.
414        - This method updates ``self.data['v']`` in place.
415        - Requires that ``nu``, ``B``, and ``rho`` are set beforehand.
416        """
417        v = self.data['nu']
418        fnu = (3/(2*(2/3*(1+v)/(1-2*v))**(3/2) + (1/3*(1+v)/(1-v))**(3/2)))**(1/3)
419        self.data['v'] = 1e3 * fnu * (self.data['B']/self.data['rho'])**0.5
420    
421    
422    def calculate_Debye_temperature(self) -> None:
423        """
424        Compute and store the Debye temperature, Θ_D.
425    
426        The Debye temperature is evaluated as
427        $$
428        \\Theta_D \\,=\\, \\frac{h}{k_B}\\, v \\,\\left(\\frac{3}{4\\pi\\,\\Omega}\\right)^{1/3},
429        $$
430        where $\\Omega = N_\\mathrm{at}\\,V$ is the primitive-cell volume and
431        $v$ is the mean sound velocity.
432    
433        Parameters
434        ----------
435        None
436    
437        Returns
438        -------
439        None
440            The computed Debye temperature is stored in ``self.data['ThetaD']`` [K].
441    
442        Notes
443        -----
444        - Units:
445          - Input: ``V`` in [Å^3/atom], ``Nat`` dimensionless (atoms/cell),
446            so $\\Omega$ is in [Å^3/cell]; ``v`` in [m/s].
447          - Constants: ``h`` in [J·s], ``kB`` in [J·K^-1].
448          - Output: ``ThetaD`` in [K].
449          - The factor ``1e10`` converts [Å^-1] to [m^-1].
450        - This method calls ``calculate_mean_sound_velocity()`` to ensure ``v`` is up to date.
451        - Requires that ``V``, ``Nat``, and the inputs needed for ``v`` are set.
452        """
453        kB = 1.38064852e-23  # J K^-1  (m2 kg s-2 K-1)
454        h = 6.62607015e-34   # J s     (m2 kg s^-1)
455        from math import pi
456    
457        self.calculate_mean_sound_velocity()
458        Omega = self.data['Nat'] * self.data['V']  # Å^3 per cell
459        self.data['ThetaD'] = 1e10 * h/kB * (3/(4*pi*Omega))**(1/3) * self.data['v']
460    
461    
462    def DebyeIntegral(self, TD: float | list[float], order: int = 3) -> float | list[float]:
463        """
464        Evaluate the Debye integral of order n.
465    
466        This computes the dimensionless Debye function
467        $$
468        D_n(x) \\,=\\, \\frac{n}{x^n} \\int_0^x \\frac{t^n}{e^t - 1}\\,dt,
469        $$
470        where typically $n=3$ and $x=\\Theta_D/T$.
471    
472        Parameters
473        ----------
474        TD : float or array-like
475            Upper bound $x = \\Theta_D/T$ (dimensionless).
476        order : int, optional
477            Debye integral order $n$ (default is 3).
478    
479        Returns
480        -------
481        float or list of float
482            Value(s) of $D_n(x)$ at the provided bound(s). Returns a float if
483            ``TD`` is scalar, otherwise a list of floats.
484    
485        Notes
486        -----
487        - Numerical evaluation is performed via ``scipy.integrate.quad``.
488        - This helper is used in the vibrational free-energy expression of
489          the harmonic Debye model.
490        """
491    
492        def __f(x):
493            return x**order/(exp(x)-1)
494    
495        if hasattr(TD, '__len__'):
496            return [order/(TDD**order) * integrate.quad(__f, 0, TDD)[0] for TDD in TD]
497        else:
498            return order/(TD**order) * integrate.quad(__f, 0, TD)[0]
499    
500    
501    def Fvib_harm(self, T: float) -> float:
502        """
503        Vibrational Helmholtz free energy in the harmonic Debye model, per atom.
504    
505        The vibrational contribution is computed as
506        $$
507        F_\\mathrm{vib}(T)
508          \\,=\\, \\frac{9}{8} k_B \\Theta_D
509          \\, - \\, k_B T
510          \\left[
511            D_3\\!\\left(\\frac{\\Theta_D}{T}\\right)
512            \\, -\\, 3\\ln\\!\\left(1 - e^{-\\Theta_D/T}\\right)
513          \\right],
514        $$
515        where $D_3(x)$ is the Debye function of order 3.
516    
517        Parameters
518        ----------
519        T : float
520            Temperature $T$ [K].
521    
522        Returns
523        -------
524        float
525            Vibrational free energy per atom, $F_\\mathrm{vib}(T)$ [eV/atom].
526    
527        Notes
528        -----
529        - Uses ``kB = 8.617333262145e-5`` eV/K.
530        - Calls ``calculate_Debye_temperature()`` internally; ensure the material
531          parameters needed for $\\Theta_D$ are set.
532        - The result is per atom, consistent with the class unit conventions.
533        """
534        kB = 8.617333262145e-5  # eV K^-1
535        from numpy import log, exp
536    
537        self.calculate_Debye_temperature()
538        ThetaD = self.data['ThetaD']
539        debInt = self.DebyeIntegral(ThetaD / T)
540    
541        return 9/8 * kB * ThetaD - kB * T * (debInt - 3 * log(1 - exp(-ThetaD / T)))
542    
543    
544    def F_harm(self, T: float) -> float:
545        """
546        Total Helmholtz free energy in the harmonic Debye model, per atom.
547    
548        The total free energy at fixed volume in the harmonic limit is
549        $$
550        F(T) \\,=\\, E_0 \\, + \\, F_\\mathrm{vib}(T).
551        $$
552    
553        Parameters
554        ----------
555        T : float
556            Temperature $T$ [K].
557    
558        Returns
559        -------
560        float
561            Total Helmholtz free energy per atom, $F(T)$ [eV/atom].
562    
563        Notes
564        -----
565        - Requires that ``E0`` is set (equilibrium energy per atom), and that the
566          parameters needed to compute $\\Theta_D$ are available.
567        - Internally calls ``Fvib_harm(T)``.
568        """
569        return self.data['E0'] + self.Fvib_harm(T)

Implementation of the harmonic Debye approximation.

This class stores and provides accessors for thermodynamic and elastic descriptors required to evaluate properties in the harmonic Debye model, such as the Helmholtz free energy, heat capacity, and entropy.

Parameters
  • at_per_fu (int, optional): Number of atoms in the primitive cell, $N_\mathrm{at}$, by default 1.
Attributes
  • data (dict): Internal storage for material fields:
    • 'V': float or None Volume per atom, $V$ [Å^3/atom].
    • 'rho': float or None Mass density, $\rho$ [g/cm^3].
    • 'Nat': int Number of atoms in the primitive cell, $N_\mathrm{at}$ [1].
    • 'B': float or None Bulk modulus, $B$ [GPa].
    • 'nu': float or None Poisson's ratio, $\nu$ [1].
    • 'v': float or None Mean (Debye) sound velocity, $v$ [m/s].
    • 'ThetaD': float or None Debye temperature, $\Theta_D$ [K].
    • 'E0': float or None Equilibrium total energy per atom, $E_0$ [eV/atom].
Notes

In the harmonic Debye approximation, vibrational contributions depend primarily on the Debye temperature $\Theta_D$ (or equivalently the Debye frequency), which can be related to the elastic constants (here represented by $B$ and $\nu$) and the density $\rho$ through the mean sound velocity $v$. The model is commonly used to approximate the vibrational free energy $$ F_\mathrm{vib}(T) = k_B T \left[ 3 N_\mathrm{at} \ln!\left(1 - e^{-\Theta_D/T}\right)

  • 9 N_\mathrm{at} \frac{T}{\Theta_D} \int_0^{\Theta_D/T} \frac{x^3}{e^x - 1}\,dx \right], $$ and the total Helmholtz free energy can be expressed as $F(T) = E_0 + F_\mathrm{vib}(T)$ at fixed volume in the harmonic limit.

Unit conventions are per-atom (where applicable) to avoid ambiguity when comparing different structures or formula units.

Debye(at_per_fu: int = 1) 
109    def __init__(self, at_per_fu: int = 1):
110        self.data = {
111            'V': None,        # float: Å^3/atom
112            'rho': None,      # float: g/cm^3
113            'Nat': at_per_fu, # int: atoms in primitive cell
114            'B': None,        # float: GPa
115            'nu': None,       # float: dimensionless
116            'v': None,        # float: m/s
117            'ThetaD': None,   # float: K
118            'E0': None        # float: eV/atom
119        }
data
V: Optional[float] 
124    @property
125    def V(self) -> Optional[float]:
126        """
127        Volume per atom, $V$.
128
129        Returns
130        -------
131        float or None
132            Volume per atom $V$ [Å^3/atom], or ``None`` if unset.
133
134        Notes
135        -----
136        - Use per-atom volume to avoid dependence on the choice of the
137          crystallographic cell or formula unit size.
138        """
139        return self.data['V']

Volume per atom, $V$.

Returns
  • float or None: Volume per atom $V$ [Å^3/atom], or None if unset.
Notes
  • Use per-atom volume to avoid dependence on the choice of the crystallographic cell or formula unit size.
rho: Optional[float] 
160    @property
161    def rho(self) -> Optional[float]:
162        """
163        Mass density, $\\rho$.
164
165        Returns
166        -------
167        float or None
168            Mass density $\\rho$ [g/cm^3], or ``None`` if unset.
169        """
170        return self.data['rho']

Mass density, $\rho$.

Returns
  • float or None: Mass density $\rho$ [g/cm^3], or None if unset.
Nat: int 
192    @property
193    def Nat(self) -> int:
194        """
195        Number of atoms in the primitive cell, $N_\\mathrm{at}$.
196
197        Returns
198        -------
199        int
200            $N_\\mathrm{at}$ [1].
201        """
202        return self.data['Nat']

Number of atoms in the primitive cell, $N_\mathrm{at}$.

Returns
  • int: $N_\mathrm{at}$ [1].
B: Optional[float] 
223    @property
224    def B(self) -> Optional[float]:
225        """
226        Bulk modulus, $B$.
227
228        Returns
229        -------
230        float or None
231            Bulk modulus $B$ [GPa], or ``None`` if unset.
232        """
233        return self.data['B']

Bulk modulus, $B$.

Returns
  • float or None: Bulk modulus $B$ [GPa], or None if unset.
nu: Optional[float] 
250    @property
251    def nu(self) -> Optional[float]:
252        """
253        Poisson's ratio, $\\nu$.
254
255        Returns
256        -------
257        float or None
258            Poisson's ratio $\\nu$ [dimensionless], or ``None`` if unset.
259
260        Notes
261        -----
262        - For mechanically stable, isotropic media, $-1 < \\nu < 0.5$.
263        """
264        return self.data['nu']

Poisson's ratio, $\nu$.

Returns
  • float or None: Poisson's ratio $\nu$ [dimensionless], or None if unset.
Notes
  • For mechanically stable, isotropic media, $-1 < \nu < 0.5$.
v: Optional[float] 
281    @property
282    def v(self) -> Optional[float]:
283        """
284        Mean (Debye) sound velocity, $v$.
285
286        Returns
287        -------
288        float or None
289            Mean sound velocity $v$ [m/s], or ``None`` if unset.
290
291        Notes
292        -----
293        - Often derived from elastic constants and density
294          (e.g., via averages of longitudinal and transverse speeds).
295        """
296        return self.data['v']

Mean (Debye) sound velocity, $v$.

Returns
  • float or None: Mean sound velocity $v$ [m/s], or None if unset.
Notes
  • Often derived from elastic constants and density (e.g., via averages of longitudinal and transverse speeds).
ThetaD: Optional[float] 
313    @property
314    def ThetaD(self) -> Optional[float]:
315        """
316        Debye temperature, $\\Theta_D$.
317
318        Returns
319        -------
320        float or None
321            Debye temperature $\\Theta_D$ [K], or ``None`` if unset.
322
323        Notes
324        -----
325        - $\\Theta_D$ can be computed from $v$ and number density.
326        - In the harmonic Debye model, $\\Theta_D$ sets the vibrational
327          energy scale that governs $F_\\mathrm{vib}(T)$ and $C_V(T)$.
328        """
329        return self.data['ThetaD']

Debye temperature, $\Theta_D$.

Returns
  • float or None: Debye temperature $\Theta_D$ [K], or None if unset.
Notes
  • $\Theta_D$ can be computed from $v$ and number density.
  • In the harmonic Debye model, $\Theta_D$ sets the vibrational energy scale that governs $F_\mathrm{vib}(T)$ and $C_V(T)$.
E0: Optional[float] 
346    @property
347    def E0(self) -> Optional[float]:
348        """
349        Equilibrium total energy per atom, $E_0$.
350
351        Returns
352        -------
353        float or None
354            $E_0$ [eV/atom], or ``None`` if unset.
355
356        Notes
357        -----
358        - In the harmonic limit at fixed volume, the Helmholtz free energy is
359          $F(T) = E_0 + F_\\mathrm{vib}(T)$.
360        """
361        return self.data['E0']

Equilibrium total energy per atom, $E_0$.

Returns
  • float or None: $E_0$ [eV/atom], or None if unset.
Notes
  • In the harmonic limit at fixed volume, the Helmholtz free energy is $F(T) = E_0 + F_\mathrm{vib}(T)$.
def calculate_mean_sound_velocity(self) -> None: 
376    def calculate_mean_sound_velocity(self) -> None:
377        """
378        Compute and store the mean (Debye) sound velocity, v.
379    
380        The mean sound velocity is estimated from Poisson's ratio and the
381        bulk modulus and density via:
382        $$
383        v = f(\nu)\,\sqrt{B/\rho},
384        $$
385        where
386        $$
387        f(\nu) =
388        \left[
389            \frac{3}{
390                2\left(\tfrac{2}{3}\tfrac{1+\nu}{1-2\nu}\right)^{3/2}
391              + \left(\tfrac{1}{3}\tfrac{1+\nu}{1-\nu}\right)^{3/2}
392            }
393        \right]^{1/3}.
394        $$
395    
396        Parameters
397        ----------
398        None
399    
400        Returns
401        -------
402        None
403            The computed mean sound velocity is stored in ``self.data['v']`` [m/s].
404    
405        Notes
406        -----
407        - Units:
408          - Input: ``B`` in [GPa], ``rho`` in [g/cm^3], ``nu`` dimensionless.
409          - Output: ``v`` in [m/s].
410          - The factor ``1e3`` accounts for the conversion
411            $\\sqrt{\\text{GPa}/(\\text{g cm}^{-3})} \\to \\text{m s}^{-1}$,
412            since $1\\,\\text{GPa}=10^9\\,\\text{Pa}$ and
413            $1\\,\\text{g cm}^{-3}=10^3\\,\\text{kg m}^{-3}$.
414        - This method updates ``self.data['v']`` in place.
415        - Requires that ``nu``, ``B``, and ``rho`` are set beforehand.
416        """
417        v = self.data['nu']
418        fnu = (3/(2*(2/3*(1+v)/(1-2*v))**(3/2) + (1/3*(1+v)/(1-v))**(3/2)))**(1/3)
419        self.data['v'] = 1e3 * fnu * (self.data['B']/self.data['rho'])**0.5

Compute and store the mean (Debye) sound velocity, v.

    The mean sound velocity is estimated from Poisson's ratio and the
    bulk modulus and density via:
    $$
    v = f(

u)\,\sqrt{B/ ho}, $$ where $$ f( u) = \left[ rac{3}{ 2\left( frac{2}{3} frac{1+ u}{1-2 u} ight)^{3/2} + \left( frac{1}{3} frac{1+ u}{1- u} ight)^{3/2} }

ight]^{1/3}. $$

    Parameters
    ----------
    None

    Returns
    -------
    None
        The computed mean sound velocity is stored in ``self.data['v']`` [m/s].

    Notes
    -----
    - Units:
      - Input: ``B`` in [GPa], ``rho`` in [g/cm^3], ``nu`` dimensionless.
      - Output: ``v`` in [m/s].
      - The factor ``1e3`` accounts for the conversion
        $\sqrt{\text{GPa}/(\text{g cm}^{-3})} \to \text{m s}^{-1}$,
        since $1\,\text{GPa}=10^9\,\text{Pa}$ and
        $1\,\text{g cm}^{-3}=10^3\,\text{kg m}^{-3}$.
    - This method updates ``self.data['v']`` in place.
    - Requires that ``nu``, ``B``, and ``rho`` are set beforehand.
def calculate_Debye_temperature(self) -> None: 
422    def calculate_Debye_temperature(self) -> None:
423        """
424        Compute and store the Debye temperature, Θ_D.
425    
426        The Debye temperature is evaluated as
427        $$
428        \\Theta_D \\,=\\, \\frac{h}{k_B}\\, v \\,\\left(\\frac{3}{4\\pi\\,\\Omega}\\right)^{1/3},
429        $$
430        where $\\Omega = N_\\mathrm{at}\\,V$ is the primitive-cell volume and
431        $v$ is the mean sound velocity.
432    
433        Parameters
434        ----------
435        None
436    
437        Returns
438        -------
439        None
440            The computed Debye temperature is stored in ``self.data['ThetaD']`` [K].
441    
442        Notes
443        -----
444        - Units:
445          - Input: ``V`` in [Å^3/atom], ``Nat`` dimensionless (atoms/cell),
446            so $\\Omega$ is in [Å^3/cell]; ``v`` in [m/s].
447          - Constants: ``h`` in [J·s], ``kB`` in [J·K^-1].
448          - Output: ``ThetaD`` in [K].
449          - The factor ``1e10`` converts [Å^-1] to [m^-1].
450        - This method calls ``calculate_mean_sound_velocity()`` to ensure ``v`` is up to date.
451        - Requires that ``V``, ``Nat``, and the inputs needed for ``v`` are set.
452        """
453        kB = 1.38064852e-23  # J K^-1  (m2 kg s-2 K-1)
454        h = 6.62607015e-34   # J s     (m2 kg s^-1)
455        from math import pi
456    
457        self.calculate_mean_sound_velocity()
458        Omega = self.data['Nat'] * self.data['V']  # Å^3 per cell
459        self.data['ThetaD'] = 1e10 * h/kB * (3/(4*pi*Omega))**(1/3) * self.data['v']

Compute and store the Debye temperature, Θ_D.

The Debye temperature is evaluated as $$ \Theta_D \,=\, \frac{h}{k_B}\, v \,\left(\frac{3}{4\pi\,\Omega}\right)^{1/3}, $$ where $\Omega = N_\mathrm{at}\,V$ is the primitive-cell volume and $v$ is the mean sound velocity.

Parameters
  • None
Returns
  • None: The computed Debye temperature is stored in self.data['ThetaD'] [K].
Notes
  • Units:
    • Input: V in [Å^3/atom], Nat dimensionless (atoms/cell), so $\Omega$ is in [Å^3/cell]; v in [m/s].
    • Constants: h in [J·s], kB in [J·K^-1].
    • Output: ThetaD in [K].
    • The factor 1e10 converts [Å^-1] to [m^-1].
  • This method calls calculate_mean_sound_velocity() to ensure v is up to date.
  • Requires that V, Nat, and the inputs needed for v are set.
def DebyeIntegral(self, TD: float | list[float], order: int = 3) -> float | list[float]: 
462    def DebyeIntegral(self, TD: float | list[float], order: int = 3) -> float | list[float]:
463        """
464        Evaluate the Debye integral of order n.
465    
466        This computes the dimensionless Debye function
467        $$
468        D_n(x) \\,=\\, \\frac{n}{x^n} \\int_0^x \\frac{t^n}{e^t - 1}\\,dt,
469        $$
470        where typically $n=3$ and $x=\\Theta_D/T$.
471    
472        Parameters
473        ----------
474        TD : float or array-like
475            Upper bound $x = \\Theta_D/T$ (dimensionless).
476        order : int, optional
477            Debye integral order $n$ (default is 3).
478    
479        Returns
480        -------
481        float or list of float
482            Value(s) of $D_n(x)$ at the provided bound(s). Returns a float if
483            ``TD`` is scalar, otherwise a list of floats.
484    
485        Notes
486        -----
487        - Numerical evaluation is performed via ``scipy.integrate.quad``.
488        - This helper is used in the vibrational free-energy expression of
489          the harmonic Debye model.
490        """
491    
492        def __f(x):
493            return x**order/(exp(x)-1)
494    
495        if hasattr(TD, '__len__'):
496            return [order/(TDD**order) * integrate.quad(__f, 0, TDD)[0] for TDD in TD]
497        else:
498            return order/(TD**order) * integrate.quad(__f, 0, TD)[0]

Evaluate the Debye integral of order n.

This computes the dimensionless Debye function $$ D_n(x) \,=\, \frac{n}{x^n} \int_0^x \frac{t^n}{e^t - 1}\,dt, $$ where typically $n=3$ and $x=\Theta_D/T$.

Parameters
  • TD (float or array-like): Upper bound $x = \Theta_D/T$ (dimensionless).
  • order (int, optional): Debye integral order $n$ (default is 3).
Returns
  • float or list of float: Value(s) of $D_n(x)$ at the provided bound(s). Returns a float if TD is scalar, otherwise a list of floats.
Notes
  • Numerical evaluation is performed via scipy.integrate.quad.
  • This helper is used in the vibrational free-energy expression of the harmonic Debye model.
def Fvib_harm(self, T: float) -> float: 
501    def Fvib_harm(self, T: float) -> float:
502        """
503        Vibrational Helmholtz free energy in the harmonic Debye model, per atom.
504    
505        The vibrational contribution is computed as
506        $$
507        F_\\mathrm{vib}(T)
508          \\,=\\, \\frac{9}{8} k_B \\Theta_D
509          \\, - \\, k_B T
510          \\left[
511            D_3\\!\\left(\\frac{\\Theta_D}{T}\\right)
512            \\, -\\, 3\\ln\\!\\left(1 - e^{-\\Theta_D/T}\\right)
513          \\right],
514        $$
515        where $D_3(x)$ is the Debye function of order 3.
516    
517        Parameters
518        ----------
519        T : float
520            Temperature $T$ [K].
521    
522        Returns
523        -------
524        float
525            Vibrational free energy per atom, $F_\\mathrm{vib}(T)$ [eV/atom].
526    
527        Notes
528        -----
529        - Uses ``kB = 8.617333262145e-5`` eV/K.
530        - Calls ``calculate_Debye_temperature()`` internally; ensure the material
531          parameters needed for $\\Theta_D$ are set.
532        - The result is per atom, consistent with the class unit conventions.
533        """
534        kB = 8.617333262145e-5  # eV K^-1
535        from numpy import log, exp
536    
537        self.calculate_Debye_temperature()
538        ThetaD = self.data['ThetaD']
539        debInt = self.DebyeIntegral(ThetaD / T)
540    
541        return 9/8 * kB * ThetaD - kB * T * (debInt - 3 * log(1 - exp(-ThetaD / T)))

Vibrational Helmholtz free energy in the harmonic Debye model, per atom.

The vibrational contribution is computed as $$ F_\mathrm{vib}(T) \,=\, \frac{9}{8} k_B \Theta_D \, - \, k_B T \left[ D_3!\left(\frac{\Theta_D}{T}\right) \, -\, 3\ln!\left(1 - e^{-\Theta_D/T}\right) \right], $$ where $D_3(x)$ is the Debye function of order 3.

Parameters
  • T (float): Temperature $T$ [K].
Returns
  • float: Vibrational free energy per atom, $F_\mathrm{vib}(T)$ [eV/atom].
Notes
  • Uses kB = 8.617333262145e-5 eV/K.
  • Calls calculate_Debye_temperature() internally; ensure the material parameters needed for $\Theta_D$ are set.
  • The result is per atom, consistent with the class unit conventions.
def F_harm(self, T: float) -> float: 
544    def F_harm(self, T: float) -> float:
545        """
546        Total Helmholtz free energy in the harmonic Debye model, per atom.
547    
548        The total free energy at fixed volume in the harmonic limit is
549        $$
550        F(T) \\,=\\, E_0 \\, + \\, F_\\mathrm{vib}(T).
551        $$
552    
553        Parameters
554        ----------
555        T : float
556            Temperature $T$ [K].
557    
558        Returns
559        -------
560        float
561            Total Helmholtz free energy per atom, $F(T)$ [eV/atom].
562    
563        Notes
564        -----
565        - Requires that ``E0`` is set (equilibrium energy per atom), and that the
566          parameters needed to compute $\\Theta_D$ are available.
567        - Internally calls ``Fvib_harm(T)``.
568        """
569        return self.data['E0'] + self.Fvib_harm(T)

Total Helmholtz free energy in the harmonic Debye model, per atom.

The total free energy at fixed volume in the harmonic limit is $$ F(T) \,=\, E_0 \, + \, F_\mathrm{vib}(T). $$

Parameters
  • T (float): Temperature $T$ [K].
Returns
  • float: Total Helmholtz free energy per atom, $F(T)$ [eV/atom].
Notes
  • Requires that E0 is set (equilibrium energy per atom), and that the parameters needed to compute $\Theta_D$ are available.
  • Internally calls Fvib_harm(T).
class DebyeQHA: 
 574class DebyeQHA:
 575    """
 576    Quasi-harmonic Debye model container for equation-of-state parameters.
 577
 578    This class stores per-atom reference quantities and material information
 579    needed for quasi-harmonic (volume-dependent) thermodynamic calculations,
 580    typically combining the Debye model for vibrational contributions with an
 581    equation of state (e.g., Birch–Murnaghan) for the static lattice energy.
 582
 583    Parameters
 584    ----------
 585    E0 : float
 586        Reference (equilibrium) energy per atom, $E_0$ [eV/atom].
 587    V0 : float
 588        Reference (equilibrium) volume per atom, $V_0$ [Å^3/atom].
 589    B0 : float
 590        Bulk modulus at $V_0$, $B_0$ [GPa].
 591    Bp : float
 592        Pressure derivative of the bulk modulus at $V_0$, $B'_0$ [dimensionless].
 593    material : str or ase.Atoms
 594        Chemical formula string (e.g., 'Al', 'Fe2O3') or an ASE ``Atoms`` object
 595        describing the material for which thermodynamic properties are computed.
 596    at_per_fu : int, optional
 597        Number of atoms in the primitive cell, $N_\\mathrm{at}$ [1], by default 1.
 598
 599    Attributes
 600    ----------
 601    data : dict
 602        Internal storage for model parameters:
 603        - 'E0': float
 604            $E_0$ [eV/atom].
 605        - 'V0': float
 606            $V_0$ [Å^3/atom].
 607        - 'B0': float
 608            $B_0$ [GPa].
 609        - 'Bp': float
 610            $B'_0$ [1].
 611        - 'material': str or ase.Atoms
 612            Material descriptor (formula or ``Atoms``).
 613        - 'Nat': int
 614            $N_\\mathrm{at}$ [1].
 615        - 'volumes': sequence of float or None
 616            Mesh of volumes per atom, $\\{V\\}$ [Å^3/atom], for QHA scans.
 617
 618    Notes
 619    -----
 620    - Units are per atom throughout to keep a consistent convention.
 621    - The pair $(B_0, B'_0)$ is commonly used with the third-order Birch–Murnaghan
 622      equation of state to describe $E(V)$ around $V_0$.
 623    """
 624
 625    def __init__(
 626        self,
 627        E0: float,
 628        V0: float,
 629        B0: float,
 630        Bp: float,
 631        material: Union[str, Atoms],
 632        at_per_fu: int = 1
 633    ) -> None:
 634        """
 635        Initialize the DebyeQHA container with equation-of-state parameters.
 636
 637        See class docstring for a detailed description of parameters and units.
 638        """
 639        self.data = {
 640            'E0': E0,           # float: eV/atom
 641            'V0': V0,           # float: Å^3/atom
 642            'B0': B0,           # float: GPa
 643            'Bp': Bp,           # float: 1 (dimensionless)
 644            'nu': None,         # float: 1 (dimensionless)
 645            'material': material,
 646            'Nat': at_per_fu,   # int: atoms in primitive cell
 647            'volumes': None     # Sequence[float]: Å^3/atom
 648        }
 649        self._deb = []          # array of Debye (HA) objects
 650
 651    # Nat [atoms in primitive cell]
 652    @property
 653    def Nat(self) -> int:
 654        """
 655        Number of atoms in the primitive cell, $N_\\mathrm{at}$.
 656    
 657        Returns
 658        -------
 659        int
 660            $N_\\mathrm{at}$ [1].
 661        """
 662        return self.data['Nat']
 663    
 664    @Nat.setter
 665    def Nat(self, Nat: int) -> None:
 666        """
 667        Set the number of atoms in the primitive cell.
 668    
 669        Parameters
 670        ----------
 671        Nat : int
 672            $N_\\mathrm{at}$ [1].
 673    
 674        Returns
 675        -------
 676        None
 677    
 678        Raises
 679        ------
 680        ValueError
 681            If ``Nat`` is not a positive integer.
 682        """
 683        if not isinstance(Nat, int) or Nat <= 0:
 684            raise ValueError("Nat must be a positive integer")
 685        self.data['Nat'] = Nat
 686
 687    # E0 [eV/atom]
 688    @property
 689    def E0(self) -> float:
 690        """
 691        Equilibrium energy per atom, $E_0$.
 692    
 693        Returns
 694        -------
 695        float
 696            $E_0$ [eV/atom].
 697        """
 698        return self.data['E0']
 699    
 700    @E0.setter
 701    def E0(self, E0: float) -> None:
 702        """
 703        Set the equilibrium energy per atom.
 704    
 705        Parameters
 706        ----------
 707        E0 : float
 708            $E_0$ [eV/atom].
 709    
 710        Returns
 711        -------
 712        None
 713        """
 714        self.data['E0'] = float(E0)
 715
 716
 717    # V0 [Å^3/atom]
 718    @property
 719    def V0(self) -> float:
 720        """
 721        Equilibrium volume per atom, $V_0$.
 722    
 723        Returns
 724        -------
 725        float
 726            $V_0$ [Å^3/atom].
 727        """
 728        return self.data['V0']
 729    
 730    @V0.setter
 731    def V0(self, V0: float) -> None:
 732        """
 733        Set the equilibrium volume per atom.
 734    
 735        Parameters
 736        ----------
 737        V0 : float
 738            $V_0$ [Å^3/atom].
 739    
 740        Returns
 741        -------
 742        None
 743        """
 744        self.data['V0'] = float(V0)
 745
 746
 747    # B0 [GPa]
 748    @property
 749    def B0(self) -> float:
 750        """
 751        Bulk modulus at $V_0$, $B_0$.
 752    
 753        Returns
 754        -------
 755        float
 756            $B_0$ [GPa].
 757        """
 758        return self.data['B0']
 759    
 760    @B0.setter
 761    def B0(self, B0: float) -> None:
 762        """
 763        Set the bulk modulus at $V_0$.
 764    
 765        Parameters
 766        ----------
 767        B0 : float
 768            $B_0$ [GPa].
 769    
 770        Returns
 771        -------
 772        None
 773        """
 774        self.data['B0'] = float(B0)
 775
 776
 777    # Bp [dimensionless]
 778    @property
 779    def Bp(self) -> float:
 780        """
 781        Pressure derivative of the bulk modulus at $V_0$, $B'_0$.
 782    
 783        Returns
 784        -------
 785        float
 786            $B'_0$ [dimensionless].
 787        """
 788        return self.data['Bp']
 789    
 790    @Bp.setter
 791    def Bp(self, Bp: float) -> None:
 792        """
 793        Set the pressure derivative of the bulk modulus at $V_0$.
 794    
 795        Parameters
 796        ----------
 797        Bp : float
 798            $B'_0$ [dimensionless].
 799    
 800        Returns
 801        -------
 802        None
 803        """
 804        self.data['Bp'] = float(Bp)
 805
 806
 807    # material [formula or ASE Atoms]
 808    @property
 809    def material(self) -> Union[str, Atoms]:
 810        """
 811        Material descriptor.
 812    
 813        Returns
 814        -------
 815        str or ase.Atoms
 816            Chemical formula string or an ASE ``Atoms`` object.
 817        """
 818        return self.data['material']
 819    
 820    @material.setter
 821    def material(self, material: Union[str, Atoms]) -> None:
 822        """
 823        Set the material descriptor.
 824    
 825        Parameters
 826        ----------
 827        material : str or ase.Atoms
 828            Chemical formula string (e.g., 'Al', 'Fe2O3') or an ASE ``Atoms`` object.
 829    
 830        Returns
 831        -------
 832        None
 833    
 834        Raises
 835        ------
 836        TypeError
 837            If ``material`` is not a string or ``ase.Atoms`` instance.
 838        """
 839        if not isinstance(material, (str, Atoms)):
 840            raise TypeError("material must be a chemical formula string or an ase.Atoms object")
 841        self.data['material'] = material
 842
 843    # Poisson's ratio
 844    @property
 845    def nu(self) -> Optional[float]:
 846        """
 847        Poisson's ratio, $\\nu$.
 848
 849        Returns
 850        -------
 851        float or None
 852            Poisson's ratio $\\nu$ [dimensionless], or ``None`` if unset.
 853
 854        Notes
 855        -----
 856        - For mechanically stable, isotropic media, $-1 < \\nu < 0.5$.
 857        """
 858        return self.data['nu']
 859
 860    @nu.setter
 861    def nu(self, nu: float) -> None:
 862        """
 863        Set Poisson's ratio.
 864
 865        Parameters
 866        ----------
 867        nu : float
 868            Poisson's ratio $\\nu$ [dimensionless].
 869        """
 870        self.data['nu'] = nu
 871
 872    # volumes [Å^3/atom] (array-like)
 873    @property
 874    def volumes(self) -> np.ndarray:
 875        """
 876        Volume mesh per atom, {V}.
 877    
 878        Returns
 879        -------
 880        numpy.ndarray or None
 881            1D array of volumes per atom [Å^3/atom], or ``None`` if unset.
 882    
 883        Notes
 884        -----
 885        - Used to sample the quasi-harmonic free energy across different volumes.
 886        - Stored internally in ``self.data['volumes']``.
 887        """
 888        return self.data['volumes']
 889    
 890    
 891    @volumes.setter
 892    def volumes(self, volumes: Any) -> None:
 893        """
 894        Set the volume mesh per atom.
 895    
 896        Parameters
 897        ----------
 898        volumes : array-like
 899            Any array-like of volumes per atom [Å^3/atom]. Examples include
 900            list, tuple, numpy.ndarray, pandas Series.
 901    
 902        Returns
 903        -------
 904        None
 905    
 906        Raises
 907        ------
 908        ValueError
 909            If the provided mesh cannot be converted to a 1D float array or is empty.
 910    
 911        Notes
 912        -----
 913        - The input is converted to a 1D numpy array of dtype float and stored in
 914          ``self.data['volumes']``.
 915        - Ensure values are per atom and, typically, positive.
 916        """
 917        arr = np.asarray(volumes, dtype=float)
 918        if arr.ndim != 1 or arr.size == 0:
 919            raise ValueError("volumes must be a non-empty 1D array-like of floats.")
 920        self.data['volumes'] = arr
 921        self._deb = []
 922        for V in self.data['volumes']:
 923            d = Debye()
 924            d.nu = self.nu
 925            d.V = V
 926            d.E0 = self.get_E0(V)
 927            d.B = self.get_bulk_modulus(V)
 928            d.rho = get_rho(V, self.material)
 929            self._deb.append(deepcopy(d))
 930
 931    
 932    def get_E0(self, V: float) -> float:
 933        """
 934        Energy at volume from stored Birch–Murnaghan parameters.
 935    
 936        Evaluates the energy per atom at a given per-atom volume using the
 937        Birch–Murnaghan equation of state with the parameters stored in the instance.
 938    
 939        Parameters
 940        ----------
 941        V : float
 942            Volume per atom, $V$ [Å^3/atom].
 943    
 944        Returns
 945        -------
 946        float
 947            Energy per atom, $E(V)$ [eV/atom].
 948    
 949        Notes
 950        -----
 951        - Uses the instance's ``E0``, ``V0``, ``B0`` (in GPa), and ``Bp`` values.
 952        - Delegates to ``BM_EoS(V, E0, V0, B0, Bp)``.
 953        """
 954        return self.BM_EoS(V, self.E0, self.V0, self.B0, self.Bp)
 955    
 956    
 957    def BM_EoS(self, V: float, E0: float, V0: float, B0: float, B0p: float) -> float:
 958        """
 959        Third-order Birch–Murnaghan equation of state, energy vs. volume per atom.
 960    
 961        Computes the per-atom energy as a function of per-atom volume using:
 962        $$
 963        E(V) \;=\; E_0 \;+\; \frac{9\,V_0\,B_0}{16}
 964        \\left\\{
 965          \\left[\\left(\\frac{V_0}{V}\\right)^{2/3}-1\\right]^3 B'_0
 966          +
 967          \\left[\\left(\\frac{V_0}{V}\\right)^{2/3}-1\\right]^2
 968          \\left[6 - 4 \\left(\\frac{V_0}{V}\\right)^{2/3}\\right]
 969        \\right\\}.
 970        $$
 971    
 972        Parameters
 973        ----------
 974        V : float
 975            Volume per atom, $V$ [Å^3/atom].
 976        E0 : float
 977            Reference energy per atom, $E_0$ [eV/atom].
 978        V0 : float
 979            Reference volume per atom, $V_0$ [Å^3/atom].
 980        B0 : float
 981            Bulk modulus at $V_0$, $B_0$ [GPa].
 982        B0p : float
 983            Pressure derivative of the bulk modulus at $V_0$, $B'_0$ [dimensionless].
 984    
 985        Returns
 986        -------
 987        float
 988            Energy per atom, $E(V)$ [eV/atom].
 989    
 990        Notes
 991        -----
 992        - If a consistent energy unit is required, convert $B_0$ to [eV/Å^3] via
 993          $1\\,\\text{eV}/\\text{Å}^3 = 160.21766208\\,\\text{GPa}$ before use.
 994        """
 995        return E0 + 9*V0*B0/16 * (((V0/V)**(2/3)-1)**3 * B0p + ((V0/V)**(2/3)-1)**2* (6-4*(V0/V)**(2/3)))
 996    
 997    
 998    def fit_BM_EoS(self, vols, enes):
 999        """
1000        Fit Birch–Murnaghan (third-order) EoS parameters to energy–volume data.
1001    
1002        Performs a nonlinear least-squares fit of the BM equation of state to the
1003        provided per-atom volume and energy data, returning fitted parameters.
1004    
1005        Parameters
1006        ----------
1007        vols : array-like
1008            Volumes per atom, $\\{V\\}$ [Å^3/atom].
1009        enes : array-like
1010            Energies per atom, $\\{E\\}$ [eV/atom].
1011    
1012        Returns
1013        -------
1014        numpy.ndarray
1015            Array of fitted parameters with shape (4,), ordered as:
1016            [E0, V0, B0, Bp], where
1017            - E0: float, E0 [eV/atom]
1018            - V0: float, V0 [Å^3/atom]
1019            - B0: float, B0 [GPa]
1020            - Bp: float, B0' [dimensionless]
1021    
1022        Notes
1023        -----
1024        - Initial guesses used for the fit are:
1025          $E_0 = \\min(E)$,
1026          $V_0 = \\tfrac{1}{2}(\\min(V)+\\max(V))$,
1027          $B_0 \\approx 1.5\\,\\text{eV}/\\text{Å}^3$,
1028          $B'_0 = 4$.
1029        - Converts the fitted bulk modulus from [eV/Å^3] to [GPa] using
1030          $1\\,\\text{eV}/\\text{Å}^3 \\approx 160.2\\,\\text{GPa}$.
1031        - Assumes a callable ``EoS(V, E0, V0, B, Bp)`` compatible with ``curve_fit``.
1032        """
1033        eosfit, _ = curve_fit(self.BM_EoS, vols, enes, [min(enes), 0.5*(min(vols)+max(vols)), 1.5, 4])
1034        eosfit[2] *= 160.2
1035        return eosfit
1036
1037        
1038    import numpy as np
1039    from numpy.typing import ArrayLike, NDArray
1040    
1041    def calc_equilibrium(self, T: ArrayLike) -> NDArray[np.floating]:
1042        """
1043        Fit BM3 EoS at each temperature and return the fitted parameters.
1044    
1045        For each temperature in T, computes the harmonic free energies across the
1046        stored volume mesh (via self._deb[i].F_harm(T_i)), fits the third-order
1047        Birch–Murnaghan EoS to E(V), and returns the fitted parameters.
1048    
1049        Parameters
1050        ----------
1051        T : array-like
1052            Temperatures, {T} [K]. Can be a scalar or array-like.
1053    
1054        Returns
1055        -------
1056        numpy.ndarray
1057            - If T is scalar: shape (4,) ordered as [E0, V0, B0, Bp]
1058            - If T is array-like: shape (N, 4) where N = len(T)
1059              ordered as [E0, V0, B0, Bp] for each temperature.
1060    
1061        Notes
1062        -----
1063        - Assumes self.volumes is the per-atom volume mesh [Å^3/atom].
1064        - Assumes self._deb is an iterable of Debye objects aligned with self.volumes,
1065          each providing F_harm(T) -> energy [eV/atom] at temperature T.
1066        """
1067        is_scalar = np.isscalar(T) or (np.ndim(T) == 0)
1068        T_arr = np.atleast_1d(T).astype(float)
1069    
1070        fits = []
1071        for t in T_arr:
1072            enes = [d.F_harm(t) for d in self._deb]
1073            eosfit = self.fit_BM_EoS(self.volumes, enes)  # returns [E0, V0, B0, Bp]
1074            fits.append(eosfit)
1075    
1076        fits = np.asarray(fits, dtype=float)
1077        return fits[0] if is_scalar else fits
1078    
1079    
1080    def calc_F(self, T: ArrayLike) -> NDArray[np.floating] | float:
1081        """
1082        Equilibrium energy E0(T) from BM3 fit at each temperature.
1083    
1084        Parameters
1085        ----------
1086        T : array-like
1087            Temperatures, {T} [K]. Can be a scalar or array-like.
1088    
1089        Returns
1090        -------
1091        float or numpy.ndarray
1092            - If T is scalar: scalar float E0(T) [eV/atom]
1093            - If T is array-like: array of E0(T) [eV/atom] with shape (N,)
1094        """
1095        eosfits = self.calc_equilibrium(T)
1096        if eosfits.ndim == 1:
1097            return float(eosfits[0])
1098        return eosfits[:, 0]
1099    
1100    
1101    def calc_V(self, T: ArrayLike) -> NDArray[np.floating] | float:
1102        """
1103        Equilibrium volume V0(T) from BM3 fit at each temperature.
1104    
1105        Parameters
1106        ----------
1107        T : array-like
1108            Temperatures, {T} [K]. Can be a scalar or array-like.
1109    
1110        Returns
1111        -------
1112        float or numpy.ndarray
1113            - If T is scalar: scalar float V0(T) [Å^3/atom]
1114            - If T is array-like: array of V0(T) [Å^3/atom] with shape (N,)
1115        """
1116        eosfits = self.calc_equilibrium(T)
1117        if eosfits.ndim == 1:
1118            return float(eosfits[1])
1119        return eosfits[:, 1]
1120    
1121    
1122    def calc_B(self, T: ArrayLike) -> NDArray[np.floating] | float:
1123        """
1124        Equilibrium bulk modulus B0(T) from BM3 fit at each temperature.
1125    
1126        Parameters
1127        ----------
1128        T : array-like
1129            Temperatures, {T} [K]. Can be a scalar or array-like.
1130    
1131        Returns
1132        -------
1133        float or numpy.ndarray
1134            - If T is scalar: scalar float B0(T) [GPa]
1135            - If T is array-like: array of B0(T) [GPa] with shape (N,)
1136        """
1137        eosfits = self.calc_equilibrium(T)
1138        if eosfits.ndim == 1:
1139            return float(eosfits[2])
1140        return eosfits[:, 2]
1141
1142    
1143    def get_pressure(self, V: float) -> float:
1144        """
1145        Third-order Birch–Murnaghan pressure, P(V), as a function of per-atom volume.
1146    
1147        Computes the pressure using the BM equation of state:
1148        $$
1149        P(V) \;=\; \frac{3}{2} B_0
1150          \left[
1151            \left(\frac{V_0}{V}\right)^{7/3} - \left(\frac{V_0}{V}\right)^{5/3}
1152          \right]
1153          \left\{
1154            1 + \frac{3}{4}\,(B'_0 - 4)\left[\left(\frac{V_0}{V}\right)^{2/3} - 1\right]
1155          \right\}.
1156        $$
1157    
1158        Parameters
1159        ----------
1160        V : float
1161            Volume per atom, $V$ [Å^3/atom].
1162    
1163        Returns
1164        -------
1165        float
1166            Pressure, $P(V)$ [GPa].
1167    
1168        Notes
1169        -----
1170        - Units:
1171          - ``V0`` in [Å^3/atom], ``B0`` in [GPa], ``Bp`` dimensionless.
1172          - The returned pressure is in [GPa], consistent with ``B0``.
1173        """
1174        V0 = self.V0
1175        B0 = self.B0
1176        Bp = self.Bp
1177    
1178        eta = V0 / V
1179        return (3.0 * B0 / 2.0) * (eta**(7.0/3.0) - eta**(5.0/3.0)) * (1.0 + 0.75 * (Bp - 4.0) * (eta**(2.0/3.0) - 1.0))
1180
1181
1182    def get_bulk_modulus(self, V: float) -> float:
1183        """
1184        Volume-dependent (isothermal) bulk modulus, B(V), in the BM3 framework.
1185    
1186        Uses the linearized relation based on the definition of the pressure derivative
1187        at zero pressure:
1188        $$
1189        B(V) \;\approx\; B_0 \;+\; B'_0\,P(V),
1190        $$
1191        where $P(V)$ is the Birch–Murnaghan (third-order) pressure and $B'_0 = (dB/dP)_{P=0}$.
1192    
1193        Parameters
1194        ----------
1195        V : float
1196            Volume per atom, $V$ [Å^3/atom].
1197    
1198        Returns
1199        -------
1200        float
1201            Bulk modulus, $B(V)$ [GPa].
1202    
1203        Notes
1204        -----
1205        - This is a linear approximation in pressure around $P=0$; for large compressions
1206          or expansions, a full BM3 expression for $B(V)$ may be preferred.
1207        - Units: ``B0`` in [GPa], ``Bp`` dimensionless, returned $B(V)$ in [GPa].
1208        """
1209        p = self.get_pressure(V)
1210        return self.B0 + self.Bp * p

Quasi-harmonic Debye model container for equation-of-state parameters.

This class stores per-atom reference quantities and material information needed for quasi-harmonic (volume-dependent) thermodynamic calculations, typically combining the Debye model for vibrational contributions with an equation of state (e.g., Birch–Murnaghan) for the static lattice energy.

Parameters
  • E0 (float): Reference (equilibrium) energy per atom, $E_0$ [eV/atom].
  • V0 (float): Reference (equilibrium) volume per atom, $V_0$ [Å^3/atom].
  • B0 (float): Bulk modulus at $V_0$, $B_0$ [GPa].
  • Bp (float): Pressure derivative of the bulk modulus at $V_0$, $B'_0$ [dimensionless].
  • material (str or ase.Atoms): Chemical formula string (e.g., 'Al', 'Fe2O3') or an ASE Atoms object describing the material for which thermodynamic properties are computed.
  • at_per_fu (int, optional): Number of atoms in the primitive cell, $N_\mathrm{at}$ [1], by default 1.
Attributes
  • data (dict): Internal storage for model parameters:
    • 'E0': float $E_0$ [eV/atom].
    • 'V0': float $V_0$ [Å^3/atom].
    • 'B0': float $B_0$ [GPa].
    • 'Bp': float $B'_0$ [1].
    • 'material': str or ase.Atoms Material descriptor (formula or Atoms).
    • 'Nat': int $N_\mathrm{at}$ [1].
    • 'volumes': sequence of float or None Mesh of volumes per atom, ${V}$ [Å^3/atom], for QHA scans.
Notes
  • Units are per atom throughout to keep a consistent convention.
  • The pair $(B_0, B'_0)$ is commonly used with the third-order Birch–Murnaghan equation of state to describe $E(V)$ around $V_0$.
DebyeQHA( E0: float, V0: float, B0: float, Bp: float, material: Union[str, ase.atoms.Atoms], at_per_fu: int = 1) 
625    def __init__(
626        self,
627        E0: float,
628        V0: float,
629        B0: float,
630        Bp: float,
631        material: Union[str, Atoms],
632        at_per_fu: int = 1
633    ) -> None:
634        """
635        Initialize the DebyeQHA container with equation-of-state parameters.
636
637        See class docstring for a detailed description of parameters and units.
638        """
639        self.data = {
640            'E0': E0,           # float: eV/atom
641            'V0': V0,           # float: Å^3/atom
642            'B0': B0,           # float: GPa
643            'Bp': Bp,           # float: 1 (dimensionless)
644            'nu': None,         # float: 1 (dimensionless)
645            'material': material,
646            'Nat': at_per_fu,   # int: atoms in primitive cell
647            'volumes': None     # Sequence[float]: Å^3/atom
648        }
649        self._deb = []          # array of Debye (HA) objects

Initialize the DebyeQHA container with equation-of-state parameters.

See class docstring for a detailed description of parameters and units.

data
Nat: int 
652    @property
653    def Nat(self) -> int:
654        """
655        Number of atoms in the primitive cell, $N_\\mathrm{at}$.
656    
657        Returns
658        -------
659        int
660            $N_\\mathrm{at}$ [1].
661        """
662        return self.data['Nat']

Number of atoms in the primitive cell, $N_\mathrm{at}$.

Returns
  • int: $N_\mathrm{at}$ [1].
E0: float 
688    @property
689    def E0(self) -> float:
690        """
691        Equilibrium energy per atom, $E_0$.
692    
693        Returns
694        -------
695        float
696            $E_0$ [eV/atom].
697        """
698        return self.data['E0']

Equilibrium energy per atom, $E_0$.

Returns
  • float: $E_0$ [eV/atom].
V0: float 
718    @property
719    def V0(self) -> float:
720        """
721        Equilibrium volume per atom, $V_0$.
722    
723        Returns
724        -------
725        float
726            $V_0$ [Å^3/atom].
727        """
728        return self.data['V0']

Equilibrium volume per atom, $V_0$.

Returns
  • float: $V_0$ [Å^3/atom].
B0: float 
748    @property
749    def B0(self) -> float:
750        """
751        Bulk modulus at $V_0$, $B_0$.
752    
753        Returns
754        -------
755        float
756            $B_0$ [GPa].
757        """
758        return self.data['B0']

Bulk modulus at $V_0$, $B_0$.

Returns
  • float: $B_0$ [GPa].
Bp: float 
778    @property
779    def Bp(self) -> float:
780        """
781        Pressure derivative of the bulk modulus at $V_0$, $B'_0$.
782    
783        Returns
784        -------
785        float
786            $B'_0$ [dimensionless].
787        """
788        return self.data['Bp']

Pressure derivative of the bulk modulus at $V_0$, $B'_0$.

Returns
  • float: $B'_0$ [dimensionless].
material: Union[str, ase.atoms.Atoms] 
808    @property
809    def material(self) -> Union[str, Atoms]:
810        """
811        Material descriptor.
812    
813        Returns
814        -------
815        str or ase.Atoms
816            Chemical formula string or an ASE ``Atoms`` object.
817        """
818        return self.data['material']

Material descriptor.

Returns
  • str or ase.Atoms: Chemical formula string or an ASE Atoms object.
nu: Optional[float] 
844    @property
845    def nu(self) -> Optional[float]:
846        """
847        Poisson's ratio, $\\nu$.
848
849        Returns
850        -------
851        float or None
852            Poisson's ratio $\\nu$ [dimensionless], or ``None`` if unset.
853
854        Notes
855        -----
856        - For mechanically stable, isotropic media, $-1 < \\nu < 0.5$.
857        """
858        return self.data['nu']

Poisson's ratio, $\nu$.

Returns
  • float or None: Poisson's ratio $\nu$ [dimensionless], or None if unset.
Notes
  • For mechanically stable, isotropic media, $-1 < \nu < 0.5$.
volumes: numpy.ndarray 
873    @property
874    def volumes(self) -> np.ndarray:
875        """
876        Volume mesh per atom, {V}.
877    
878        Returns
879        -------
880        numpy.ndarray or None
881            1D array of volumes per atom [Å^3/atom], or ``None`` if unset.
882    
883        Notes
884        -----
885        - Used to sample the quasi-harmonic free energy across different volumes.
886        - Stored internally in ``self.data['volumes']``.
887        """
888        return self.data['volumes']

Volume mesh per atom, {V}.

Returns
  • numpy.ndarray or None: 1D array of volumes per atom [Å^3/atom], or None if unset.
Notes
  • Used to sample the quasi-harmonic free energy across different volumes.
  • Stored internally in self.data['volumes'].
def get_E0(self, V: float) -> float: 
932    def get_E0(self, V: float) -> float:
933        """
934        Energy at volume from stored Birch–Murnaghan parameters.
935    
936        Evaluates the energy per atom at a given per-atom volume using the
937        Birch–Murnaghan equation of state with the parameters stored in the instance.
938    
939        Parameters
940        ----------
941        V : float
942            Volume per atom, $V$ [Å^3/atom].
943    
944        Returns
945        -------
946        float
947            Energy per atom, $E(V)$ [eV/atom].
948    
949        Notes
950        -----
951        - Uses the instance's ``E0``, ``V0``, ``B0`` (in GPa), and ``Bp`` values.
952        - Delegates to ``BM_EoS(V, E0, V0, B0, Bp)``.
953        """
954        return self.BM_EoS(V, self.E0, self.V0, self.B0, self.Bp)

Energy at volume from stored Birch–Murnaghan parameters.

Evaluates the energy per atom at a given per-atom volume using the Birch–Murnaghan equation of state with the parameters stored in the instance.

Parameters
  • V (float): Volume per atom, $V$ [Å^3/atom].
Returns
  • float: Energy per atom, $E(V)$ [eV/atom].
Notes
  • Uses the instance's E0, V0, B0 (in GPa), and Bp values.
  • Delegates to BM_EoS(V, E0, V0, B0, Bp).
def BM_EoS(self, V: float, E0: float, V0: float, B0: float, B0p: float) -> float: 
957    def BM_EoS(self, V: float, E0: float, V0: float, B0: float, B0p: float) -> float:
958        """
959        Third-order Birch–Murnaghan equation of state, energy vs. volume per atom.
960    
961        Computes the per-atom energy as a function of per-atom volume using:
962        $$
963        E(V) \;=\; E_0 \;+\; \frac{9\,V_0\,B_0}{16}
964        \\left\\{
965          \\left[\\left(\\frac{V_0}{V}\\right)^{2/3}-1\\right]^3 B'_0
966          +
967          \\left[\\left(\\frac{V_0}{V}\\right)^{2/3}-1\\right]^2
968          \\left[6 - 4 \\left(\\frac{V_0}{V}\\right)^{2/3}\\right]
969        \\right\\}.
970        $$
971    
972        Parameters
973        ----------
974        V : float
975            Volume per atom, $V$ [Å^3/atom].
976        E0 : float
977            Reference energy per atom, $E_0$ [eV/atom].
978        V0 : float
979            Reference volume per atom, $V_0$ [Å^3/atom].
980        B0 : float
981            Bulk modulus at $V_0$, $B_0$ [GPa].
982        B0p : float
983            Pressure derivative of the bulk modulus at $V_0$, $B'_0$ [dimensionless].
984    
985        Returns
986        -------
987        float
988            Energy per atom, $E(V)$ [eV/atom].
989    
990        Notes
991        -----
992        - If a consistent energy unit is required, convert $B_0$ to [eV/Å^3] via
993          $1\\,\\text{eV}/\\text{Å}^3 = 160.21766208\\,\\text{GPa}$ before use.
994        """
995        return E0 + 9*V0*B0/16 * (((V0/V)**(2/3)-1)**3 * B0p + ((V0/V)**(2/3)-1)**2* (6-4*(V0/V)**(2/3)))

Third-order Birch–Murnaghan equation of state, energy vs. volume per atom.

Computes the per-atom energy as a function of per-atom volume using: $$ E(V) \;=\; E_0 \;+\; rac{9\,V_0\,B_0}{16} \left{ \left[\left(\frac{V_0}{V}\right)^{2/3}-1\right]^3 B'_0 + \left[\left(\frac{V_0}{V}\right)^{2/3}-1\right]^2 \left[6 - 4 \left(\frac{V_0}{V}\right)^{2/3}\right] \right}. $$

Parameters
  • V (float): Volume per atom, $V$ [Å^3/atom].
  • E0 (float): Reference energy per atom, $E_0$ [eV/atom].
  • V0 (float): Reference volume per atom, $V_0$ [Å^3/atom].
  • B0 (float): Bulk modulus at $V_0$, $B_0$ [GPa].
  • B0p (float): Pressure derivative of the bulk modulus at $V_0$, $B'_0$ [dimensionless].
Returns
  • float: Energy per atom, $E(V)$ [eV/atom].
Notes
  • If a consistent energy unit is required, convert $B_0$ to [eV/Å^3] via $1\,\text{eV}/\text{Å}^3 = 160.21766208\,\text{GPa}$ before use.
def fit_BM_EoS(self, vols, enes): 
 998    def fit_BM_EoS(self, vols, enes):
 999        """
1000        Fit Birch–Murnaghan (third-order) EoS parameters to energy–volume data.
1001    
1002        Performs a nonlinear least-squares fit of the BM equation of state to the
1003        provided per-atom volume and energy data, returning fitted parameters.
1004    
1005        Parameters
1006        ----------
1007        vols : array-like
1008            Volumes per atom, $\\{V\\}$ [Å^3/atom].
1009        enes : array-like
1010            Energies per atom, $\\{E\\}$ [eV/atom].
1011    
1012        Returns
1013        -------
1014        numpy.ndarray
1015            Array of fitted parameters with shape (4,), ordered as:
1016            [E0, V0, B0, Bp], where
1017            - E0: float, E0 [eV/atom]
1018            - V0: float, V0 [Å^3/atom]
1019            - B0: float, B0 [GPa]
1020            - Bp: float, B0' [dimensionless]
1021    
1022        Notes
1023        -----
1024        - Initial guesses used for the fit are:
1025          $E_0 = \\min(E)$,
1026          $V_0 = \\tfrac{1}{2}(\\min(V)+\\max(V))$,
1027          $B_0 \\approx 1.5\\,\\text{eV}/\\text{Å}^3$,
1028          $B'_0 = 4$.
1029        - Converts the fitted bulk modulus from [eV/Å^3] to [GPa] using
1030          $1\\,\\text{eV}/\\text{Å}^3 \\approx 160.2\\,\\text{GPa}$.
1031        - Assumes a callable ``EoS(V, E0, V0, B, Bp)`` compatible with ``curve_fit``.
1032        """
1033        eosfit, _ = curve_fit(self.BM_EoS, vols, enes, [min(enes), 0.5*(min(vols)+max(vols)), 1.5, 4])
1034        eosfit[2] *= 160.2
1035        return eosfit

Fit Birch–Murnaghan (third-order) EoS parameters to energy–volume data.

Performs a nonlinear least-squares fit of the BM equation of state to the provided per-atom volume and energy data, returning fitted parameters.

Parameters
  • vols (array-like): Volumes per atom, ${V}$ [Å^3/atom].
  • enes (array-like): Energies per atom, ${E}$ [eV/atom].
Returns
  • numpy.ndarray: Array of fitted parameters with shape (4,), ordered as: [E0, V0, B0, Bp], where
    • E0: float, E0 [eV/atom]
    • V0: float, V0 [Å^3/atom]
    • B0: float, B0 [GPa]
    • Bp: float, B0' [dimensionless]
Notes
  • Initial guesses used for the fit are: $E_0 = \min(E)$, $V_0 = \tfrac{1}{2}(\min(V)+\max(V))$, $B_0 \approx 1.5\,\text{eV}/\text{Å}^3$, $B'_0 = 4$.
  • Converts the fitted bulk modulus from [eV/Å^3] to [GPa] using $1\,\text{eV}/\text{Å}^3 \approx 160.2\,\text{GPa}$.
  • Assumes a callable EoS(V, E0, V0, B, Bp) compatible with curve_fit.
def calc_equilibrium( self, T: Union[numpy._typing._array_like._SupportsArray[numpy.dtype[Any]], numpy._typing._nested_sequence._NestedSequence[numpy._typing._array_like._SupportsArray[numpy.dtype[Any]]], bool, int, float, complex, str, bytes, numpy._typing._nested_sequence._NestedSequence[Union[bool, int, float, complex, str, bytes]]]) -> numpy.ndarray[typing.Any, numpy.dtype[numpy.floating]]: 
1041    def calc_equilibrium(self, T: ArrayLike) -> NDArray[np.floating]:
1042        """
1043        Fit BM3 EoS at each temperature and return the fitted parameters.
1044    
1045        For each temperature in T, computes the harmonic free energies across the
1046        stored volume mesh (via self._deb[i].F_harm(T_i)), fits the third-order
1047        Birch–Murnaghan EoS to E(V), and returns the fitted parameters.
1048    
1049        Parameters
1050        ----------
1051        T : array-like
1052            Temperatures, {T} [K]. Can be a scalar or array-like.
1053    
1054        Returns
1055        -------
1056        numpy.ndarray
1057            - If T is scalar: shape (4,) ordered as [E0, V0, B0, Bp]
1058            - If T is array-like: shape (N, 4) where N = len(T)
1059              ordered as [E0, V0, B0, Bp] for each temperature.
1060    
1061        Notes
1062        -----
1063        - Assumes self.volumes is the per-atom volume mesh [Å^3/atom].
1064        - Assumes self._deb is an iterable of Debye objects aligned with self.volumes,
1065          each providing F_harm(T) -> energy [eV/atom] at temperature T.
1066        """
1067        is_scalar = np.isscalar(T) or (np.ndim(T) == 0)
1068        T_arr = np.atleast_1d(T).astype(float)
1069    
1070        fits = []
1071        for t in T_arr:
1072            enes = [d.F_harm(t) for d in self._deb]
1073            eosfit = self.fit_BM_EoS(self.volumes, enes)  # returns [E0, V0, B0, Bp]
1074            fits.append(eosfit)
1075    
1076        fits = np.asarray(fits, dtype=float)
1077        return fits[0] if is_scalar else fits

Fit BM3 EoS at each temperature and return the fitted parameters.

For each temperature in T, computes the harmonic free energies across the stored volume mesh (via self._deb[i].F_harm(T_i)), fits the third-order Birch–Murnaghan EoS to E(V), and returns the fitted parameters.

Parameters
  • T (array-like): Temperatures, {T} [K]. Can be a scalar or array-like.
Returns
  • numpy.ndarray: - If T is scalar: shape (4,) ordered as [E0, V0, B0, Bp]
    • If T is array-like: shape (N, 4) where N = len(T) ordered as [E0, V0, B0, Bp] for each temperature.
Notes
  • Assumes self.volumes is the per-atom volume mesh [Å^3/atom].
  • Assumes self._deb is an iterable of Debye objects aligned with self.volumes, each providing F_harm(T) -> energy [eV/atom] at temperature T.
def calc_F( self, T: Union[numpy._typing._array_like._SupportsArray[numpy.dtype[Any]], numpy._typing._nested_sequence._NestedSequence[numpy._typing._array_like._SupportsArray[numpy.dtype[Any]]], bool, int, float, complex, str, bytes, numpy._typing._nested_sequence._NestedSequence[Union[bool, int, float, complex, str, bytes]]]) -> numpy.ndarray[typing.Any, numpy.dtype[numpy.floating]] | float: 
1080    def calc_F(self, T: ArrayLike) -> NDArray[np.floating] | float:
1081        """
1082        Equilibrium energy E0(T) from BM3 fit at each temperature.
1083    
1084        Parameters
1085        ----------
1086        T : array-like
1087            Temperatures, {T} [K]. Can be a scalar or array-like.
1088    
1089        Returns
1090        -------
1091        float or numpy.ndarray
1092            - If T is scalar: scalar float E0(T) [eV/atom]
1093            - If T is array-like: array of E0(T) [eV/atom] with shape (N,)
1094        """
1095        eosfits = self.calc_equilibrium(T)
1096        if eosfits.ndim == 1:
1097            return float(eosfits[0])
1098        return eosfits[:, 0]

Equilibrium energy E0(T) from BM3 fit at each temperature.

Parameters
  • T (array-like): Temperatures, {T} [K]. Can be a scalar or array-like.
Returns
  • float or numpy.ndarray: - If T is scalar: scalar float E0(T) [eV/atom]
    • If T is array-like: array of E0(T) [eV/atom] with shape (N,)
def calc_V( self, T: Union[numpy._typing._array_like._SupportsArray[numpy.dtype[Any]], numpy._typing._nested_sequence._NestedSequence[numpy._typing._array_like._SupportsArray[numpy.dtype[Any]]], bool, int, float, complex, str, bytes, numpy._typing._nested_sequence._NestedSequence[Union[bool, int, float, complex, str, bytes]]]) -> numpy.ndarray[typing.Any, numpy.dtype[numpy.floating]] | float: 
1101    def calc_V(self, T: ArrayLike) -> NDArray[np.floating] | float:
1102        """
1103        Equilibrium volume V0(T) from BM3 fit at each temperature.
1104    
1105        Parameters
1106        ----------
1107        T : array-like
1108            Temperatures, {T} [K]. Can be a scalar or array-like.
1109    
1110        Returns
1111        -------
1112        float or numpy.ndarray
1113            - If T is scalar: scalar float V0(T) [Å^3/atom]
1114            - If T is array-like: array of V0(T) [Å^3/atom] with shape (N,)
1115        """
1116        eosfits = self.calc_equilibrium(T)
1117        if eosfits.ndim == 1:
1118            return float(eosfits[1])
1119        return eosfits[:, 1]

Equilibrium volume V0(T) from BM3 fit at each temperature.

Parameters
  • T (array-like): Temperatures, {T} [K]. Can be a scalar or array-like.
Returns
  • float or numpy.ndarray: - If T is scalar: scalar float V0(T) [Å^3/atom]
    • If T is array-like: array of V0(T) [Å^3/atom] with shape (N,)
def calc_B( self, T: Union[numpy._typing._array_like._SupportsArray[numpy.dtype[Any]], numpy._typing._nested_sequence._NestedSequence[numpy._typing._array_like._SupportsArray[numpy.dtype[Any]]], bool, int, float, complex, str, bytes, numpy._typing._nested_sequence._NestedSequence[Union[bool, int, float, complex, str, bytes]]]) -> numpy.ndarray[typing.Any, numpy.dtype[numpy.floating]] | float: 
1122    def calc_B(self, T: ArrayLike) -> NDArray[np.floating] | float:
1123        """
1124        Equilibrium bulk modulus B0(T) from BM3 fit at each temperature.
1125    
1126        Parameters
1127        ----------
1128        T : array-like
1129            Temperatures, {T} [K]. Can be a scalar or array-like.
1130    
1131        Returns
1132        -------
1133        float or numpy.ndarray
1134            - If T is scalar: scalar float B0(T) [GPa]
1135            - If T is array-like: array of B0(T) [GPa] with shape (N,)
1136        """
1137        eosfits = self.calc_equilibrium(T)
1138        if eosfits.ndim == 1:
1139            return float(eosfits[2])
1140        return eosfits[:, 2]

Equilibrium bulk modulus B0(T) from BM3 fit at each temperature.

Parameters
  • T (array-like): Temperatures, {T} [K]. Can be a scalar or array-like.
Returns
  • float or numpy.ndarray: - If T is scalar: scalar float B0(T) [GPa]
    • If T is array-like: array of B0(T) [GPa] with shape (N,)
def get_pressure(self, V: float) -> float: 
1143    def get_pressure(self, V: float) -> float:
1144        """
1145        Third-order Birch–Murnaghan pressure, P(V), as a function of per-atom volume.
1146    
1147        Computes the pressure using the BM equation of state:
1148        $$
1149        P(V) \;=\; \frac{3}{2} B_0
1150          \left[
1151            \left(\frac{V_0}{V}\right)^{7/3} - \left(\frac{V_0}{V}\right)^{5/3}
1152          \right]
1153          \left\{
1154            1 + \frac{3}{4}\,(B'_0 - 4)\left[\left(\frac{V_0}{V}\right)^{2/3} - 1\right]
1155          \right\}.
1156        $$
1157    
1158        Parameters
1159        ----------
1160        V : float
1161            Volume per atom, $V$ [Å^3/atom].
1162    
1163        Returns
1164        -------
1165        float
1166            Pressure, $P(V)$ [GPa].
1167    
1168        Notes
1169        -----
1170        - Units:
1171          - ``V0`` in [Å^3/atom], ``B0`` in [GPa], ``Bp`` dimensionless.
1172          - The returned pressure is in [GPa], consistent with ``B0``.
1173        """
1174        V0 = self.V0
1175        B0 = self.B0
1176        Bp = self.Bp
1177    
1178        eta = V0 / V
1179        return (3.0 * B0 / 2.0) * (eta**(7.0/3.0) - eta**(5.0/3.0)) * (1.0 + 0.75 * (Bp - 4.0) * (eta**(2.0/3.0) - 1.0))

Third-order Birch–Murnaghan pressure, P(V), as a function of per-atom volume.

Computes the pressure using the BM equation of state: $$ P(V) \;=\; rac{3}{2} B_0 \left[ \left( rac{V_0}{V} ight)^{7/3} - \left( rac{V_0}{V} ight)^{5/3}

ight] \left{ 1 + rac{3}{4}\,(B'_0 - 4)\left[\left( rac{V_0}{V} ight)^{2/3} - 1 ight]

ight}. $$

Parameters
  • V (float): Volume per atom, $V$ [Å^3/atom].
Returns
  • float: Pressure, $P(V)$ [GPa].
Notes
  • Units:
    • V0 in [Å^3/atom], B0 in [GPa], Bp dimensionless.
    • The returned pressure is in [GPa], consistent with B0.
def get_bulk_modulus(self, V: float) -> float: 
1182    def get_bulk_modulus(self, V: float) -> float:
1183        """
1184        Volume-dependent (isothermal) bulk modulus, B(V), in the BM3 framework.
1185    
1186        Uses the linearized relation based on the definition of the pressure derivative
1187        at zero pressure:
1188        $$
1189        B(V) \;\approx\; B_0 \;+\; B'_0\,P(V),
1190        $$
1191        where $P(V)$ is the Birch–Murnaghan (third-order) pressure and $B'_0 = (dB/dP)_{P=0}$.
1192    
1193        Parameters
1194        ----------
1195        V : float
1196            Volume per atom, $V$ [Å^3/atom].
1197    
1198        Returns
1199        -------
1200        float
1201            Bulk modulus, $B(V)$ [GPa].
1202    
1203        Notes
1204        -----
1205        - This is a linear approximation in pressure around $P=0$; for large compressions
1206          or expansions, a full BM3 expression for $B(V)$ may be preferred.
1207        - Units: ``B0`` in [GPa], ``Bp`` dimensionless, returned $B(V)$ in [GPa].
1208        """
1209        p = self.get_pressure(V)
1210        return self.B0 + self.Bp * p

Volume-dependent (isothermal) bulk modulus, B(V), in the BM3 framework.

Uses the linearized relation based on the definition of the pressure derivative at zero pressure: $$ B(V) \;pprox\; B_0 \;+\; B'_0\,P(V), $$ where $P(V)$ is the Birch–Murnaghan (third-order) pressure and $B'_0 = (dB/dP)_{P=0}$.

Parameters
  • V (float): Volume per atom, $V$ [Å^3/atom].
Returns
  • float: Bulk modulus, $B(V)$ [GPa].
Notes
  • This is a linear approximation in pressure around $P=0$; for large compressions or expansions, a full BM3 expression for $B(V)$ may be preferred.
  • Units: B0 in [GPa], Bp dimensionless, returned $B(V)$ in [GPa].
ArrayLike = typing.Union[numpy._typing._array_like._SupportsArray[numpy.dtype[typing.Any]], numpy._typing._nested_sequence._NestedSequence[numpy._typing._array_like._SupportsArray[numpy.dtype[typing.Any]]], bool, int, float, complex, str, bytes, numpy._typing._nested_sequence._NestedSequence[typing.Union[bool, int, float, complex, str, bytes]]]
NDArray = numpy.ndarray[typing.Any, numpy.dtype[+_ScalarType_co]]