asaf.mpd

Module for handling Macrostate Probability Distribution (MPD) data.

MPD

MPD(
    dataframe: DataFrame,
    temperature: float,
    beta_mu: Optional[float] = None,
    fugacity: Optional[float] = None,
    metadata: Optional[dict[str, Any]] = None,
    order: int = 50,
    tolerance: float = 10.0,
)

Class for storing and processing macrostate probability distribution.

Parameters:

• dataframe (DataFrame) –

  a pandas dataframe with state specific data

• temperature (float) –

  temperature (in K) at which the simulation was performed

• beta_mu (Optional[float], default: None ) –

  beta_mu (unitless) at which the simulation was performed. At least one of beta_mu or fugacity must be specified

• fugacity (Optional[float], default: None ) –

  fugacity (in Pa) at which the simulation was performed. At least one of beta_mu or fugacity must be specified

• metadata (Optional[dict[str, Any]], default: None ) –

  a dictionary with the simulation metadata

• order (int, default: 50 ) –

  how many points on each side use to find minimum in lnp

• tolerance (float, default: 10.0 ) –

  used when checking the probability at lnp tail

Source code in src/asaf/mpd.py

def __init__(
    self,
    dataframe: pd.DataFrame,
    temperature: float,
    beta_mu: Optional[float] = None,
    fugacity: Optional[float] = None,
    metadata: Optional[dict[str, Any]] = None,
    order: int = 50,
    tolerance: float = 10.0,
) -> None:
    """Initialize the MPD class.

    Parameters
    ----------
    dataframe
        a pandas dataframe with state specific data
    temperature
        temperature (in K) at which the simulation was performed
    beta_mu
        beta_mu (unitless) at which the simulation was performed. At least one of beta_mu
        or fugacity must be specified
    fugacity
        fugacity (in Pa) at which the simulation was performed. At least one of beta_mu
        or fugacity must be specified
    metadata
        a dictionary with the simulation metadata
    order
        how many points on each side use to find minimum in lnp
    tolerance
        used when checking the probability at lnp tail
    """
    self.temperature = temperature

    if (beta_mu is None) and (fugacity is None):
        raise ValueError("Must provide `beta_mu` or/and `fugacity`.")
    if beta_mu is None:
        self.fugacity = fugacity
        self._beta_mu = self.beta * self.mu
    else:
        self._beta_mu = beta_mu
        self.mu = beta_mu / self.beta

    lnp_headers = ["macrostate", "lnp"]
    prob_headers = ["macrostate", "P_up", "P_down"]

    have_lnp = set(lnp_headers).issubset(dataframe.columns)
    have_prob = set(prob_headers).issubset(dataframe.columns)

    if not (have_lnp or have_prob):
        all_required = set(lnp_headers) | set(prob_headers)
        missing = all_required - set(dataframe.columns)
        raise ValueError(f"Some of the columns names {missing} are missing.")

    self._dataframe = dataframe

    if not have_lnp:
        lnp_df = calculate_lnp(dataframe[prob_headers])
        merged = lnp_df.merge(
            dataframe, on="macrostate", how="left", suffixes=("", "_inp")
        )
        self._dataframe = merged

    self._metadata = metadata or {}
    self.order = order
    self.tolerance = tolerance

    self._system_size_prod = 1

    if "system_size" in self.metadata:
        self.system_size = self.metadata["system_size"]
    else:
        self.system_size = [1, 1, 1]

    self.check_tail(order, tolerance)

beta property writable

beta: float

Return the beta (in J^-1).

beta_mu property

beta_mu: float

Return the beta_mu (unitless).

fugacity property writable

fugacity: float

Return the fugacity (in Pa).

lnp property

lnp: DataFrame

Return a dataframe with the natural logarithm of the macrostate probability.

metadata property writable

metadata: Dict[str, Any]

Return the metadata dictionary.

mu property writable

mu: float

Return the chemical potential (in J A^-3).

order property writable

order: int

Return the order used to find minimum in lnp.

system_size property writable

system_size: List[int]

Return the system size as a list of integers.

temperature property writable

temperature: float

Return the temperature (in K).

tolerance property writable

tolerance: float

Return the tolerance used when checking the probability at lnp tail.

average_macrostate

average_macrostate(
    lnp: Optional[DataFrame] = None,
) -> float

Calculate the average macrostate from the MPD data.

Note that this function does not check for multiple phases. Use average_macrostate_at_fugacity to calculate the average macrostate at a given fugacity, which checks for multiple phases.

Source code in src/asaf/mpd.py

def average_macrostate(self, lnp: Optional[pd.DataFrame] = None) -> float:
    """Calculate the average macrostate from the MPD data.

    Note that this function does not check for multiple phases. Use `average_macrostate_at_fugacity`
    to calculate the average macrostate at a given fugacity, which checks for multiple phases.
    """
    if lnp is None:
        lnp = self.lnp
    return (np.exp(lnp["lnp"]) * lnp["macrostate"]).sum()

average_macrostate_at_fugacity

average_macrostate_at_fugacity(
    fug: float, order: Optional[int] = None
) -> List[float]

Calculate the average macrostate at a given fugacity.

Source code in src/asaf/mpd.py

def average_macrostate_at_fugacity(
    self, fug: float, order: Optional[int] = None
) -> List[float]:
    """Calculate the average macrostate at a given fugacity."""
    beta_0 = self._beta
    mu_0 = self._mu
    mu = fugacity_to_mu(fug, beta_0)
    delta_beta_mu = beta_0 * (mu - mu_0)
    lnp_rw = self.reweight(delta_beta_mu)
    if order is None:
        order = self.order
    mins = self.minimums(lnp=lnp_rw["lnp"], order=order)

    if len(mins) == 0:
        return [self.average_macrostate(lnp_rw) / self._system_size_prod]
    else:
        minn = mins[mins.lnp == mins.lnp.min()].index[0]
        lnp_a = lnp_rw[:minn].copy()
        lnp_b = lnp_rw[minn + 1 :].copy()

        p_a = np.exp(lnp_a["lnp"]).sum()
        p_b = np.exp(lnp_b["lnp"]).sum()

        lnp_a["lnp"] = normalize(lnp_a["lnp"])
        lnp_b["lnp"] = normalize(lnp_b["lnp"])

        if p_a > p_b:
            return [
                self.average_macrostate(lnp_a) / self._system_size_prod,
                self.average_macrostate(lnp_b) / self._system_size_prod,
            ]
        else:
            return [
                self.average_macrostate(lnp_b) / self._system_size_prod,
                self.average_macrostate(lnp_a) / self._system_size_prod,
            ]

calculate_isotherm

calculate_isotherm(
    fugacity: ArrayLike,
    saturation_fugacity: Optional[float] = None,
    saturation_pressure: Optional[float] = None,
    pressure: Optional[ArrayLike] = None,
    order: Optional[int] = None,
    return_dataframe: bool = True,
) -> Union[DataFrame | Isotherm]

Calculate the adsorption isotherm.

Parameters:

• fugacity (ArrayLike) –

  Array of fugacities.

• saturation_fugacity (Optional[float], default: None ) –

  Saturation fugacity to calculate relative fugacity (f/f0).

• saturation_pressure (Optional[float], default: None ) –

  Saturation pressure to calculate relative pressure (p/p0).

• pressure (Optional[ArrayLike], default: None ) –

  Array of pressures corresponding to the fugacities.

• order (Optional[int], default: None ) –

  How many points on each side use to find minimum in lnp.

• return_dataframe (bool, default: True ) –

  Whether to return the adsorption isotherm as a dataframe or Isotherm instance.

Returns:

• DataFrame or Isotherm –

  DataFrame containing the adsorption isotherm or Isotherm instance if return_dataframe is False.

• Args ( Union[DataFrame | Isotherm] ) –

  pressure:

Source code in src/asaf/mpd.py

def calculate_isotherm(
    self,
    fugacity: ArrayLike,
    saturation_fugacity: Optional[float] = None,
    saturation_pressure: Optional[float] = None,
    pressure: Optional[ArrayLike] = None,
    order: Optional[int] = None,
    return_dataframe: bool = True,
) -> Union[pd.DataFrame | Isotherm]:
    """Calculate the adsorption isotherm.

    Parameters
    ----------
    fugacity
        Array of fugacities.
    saturation_fugacity
        Saturation fugacity to calculate relative fugacity (f/f0).
    saturation_pressure
        Saturation pressure to calculate relative pressure (p/p0).
    pressure
        Array of pressures corresponding to the fugacities.
    order
        How many points on each side use to find minimum in lnp.
    return_dataframe
        Whether to return the adsorption isotherm as a dataframe or Isotherm instance.

    Returns
    -------
    pd.DataFrame or Isotherm
        DataFrame containing the adsorption isotherm or Isotherm instance if return_dataframe is False.

    Args:
        pressure:
    """
    from asaf import Isotherm

    stable_phase = []
    metastable_gas = []
    metastable_liq = []

    if order is None:
        order = self.order

    for fug in fugacity:
        uptake = self.average_macrostate_at_fugacity(fug, order=order)
        if len(uptake) > 1:
            if uptake[0] > uptake[1]:
                stable_phase.append([fug, uptake[0]])
                metastable_gas.append([fug, uptake[1]])
            else:
                stable_phase.append([fug, uptake[0]])
                metastable_liq.append([fug, uptake[1]])
        else:
            stable_phase.append([fug, uptake[0]])

    isotherm = pd.DataFrame(stable_phase, columns=["fugacity", "uptake"])

    if len(metastable_gas) > 0:
        iso_metastable_gas = pd.DataFrame(
            metastable_gas, columns=["fugacity", "metastable_gas"]
        )

        isotherm = pd.merge(
            isotherm, iso_metastable_gas, on="fugacity", how="outer"
        )

    if len(metastable_liq) > 0:
        iso_metastable_liq = pd.DataFrame(
            metastable_liq, columns=["fugacity", "metastable_liq"]
        )

        isotherm = pd.merge(
            isotherm, iso_metastable_liq, on="fugacity", how="outer"
        )

    if saturation_fugacity is not None:
        isotherm.insert(1, "f/f0", isotherm["fugacity"] / saturation_fugacity)

    if pressure is not None:
        isotherm.insert(1, "pressure", np.array(pressure))
        if saturation_pressure is not None:
            isotherm.insert(2, "p/p0", isotherm["pressure"] / saturation_pressure)
    elif saturation_pressure is not None:
        raise ValueError(
            "`pressure` must be provided to calculate relative pressure (p/p0)."
        )

    if return_dataframe:
        return isotherm
    else:
        return Isotherm(
            data=isotherm,
            saturation_fugacity=saturation_fugacity,
            saturation_pressure=saturation_pressure,
            metadata=self.metadata,
        )

check_tail

check_tail(
    order: int,
    tolerance: float,
    lnp: Optional[DataFrame] = None,
) -> None

Check the probability at the tail of the lnp distribution.

Source code in src/asaf/mpd.py

def check_tail(
    self, order: int, tolerance: float, lnp: Optional[pd.DataFrame] = None
) -> None:
    """Check the probability at the tail of the lnp distribution."""
    if getattr(self, "_suppress_check_tail", False):
        return
    if lnp is None:
        lnp = self.lnp
    mins = self.minimums(lnp=lnp["lnp"], order=order)
    if len(mins) == 0:
        difference = lnp["lnp"].max() - lnp["lnp"].iloc[-1]
    else:
        minn = mins[mins.lnp == mins.lnp.min()].index[0]
        lnp_b = lnp[minn + 1 :].copy()
        difference = lnp_b["lnp"].max() - lnp_b["lnp"].iloc[-1]

    if difference < tolerance:
        print(
            f"WARNING! lnPi at N_max has a relative value higher ({difference:.1f}) than tolerance ({tolerance:.1f})."
        )
        print(
            "The results may be erroneous. Provide data for higher macrostate values."
        )

dataframe

dataframe() -> DataFrame

Return dataframe.

Source code in src/asaf/mpd.py

def dataframe(self) -> pd.DataFrame:
    """Return dataframe."""
    return self._dataframe

extrapolate

extrapolate(
    temperature: float,
    energy: Optional[DataFrame | Series] = None,
    terms: int = 1,
) -> "MPD"

Extrapolates the MPD to a new temperature.

Parameters:

• temperature (float) –

  Temperature (in K) to which to extrapolate MPD.

• energy (Optional[DataFrame | Series], default: None ) –

  Energy fluctuation data. If None ASAF will look for data in prob_df. Unit must be J.

• terms (int, default: 1 ) –

  Number of Taylor series terms used for extrapolation. Note that energy must contain columns named term_1, term_2, ..., term_n where n is the number of terms.

Returns:

• MPD –

  Extrapolated MPD.

Source code in src/asaf/mpd.py

def extrapolate(
    self,
    temperature: float,
    energy: Optional[pd.DataFrame | pd.Series] = None,
    terms: int = 1,
) -> "MPD":
    """Extrapolates the MPD to a new temperature.

    Parameters
    ----------
    temperature
        Temperature (in K) to which to extrapolate MPD.
    energy
        Energy fluctuation data. If None ASAF will look for data in prob_df. Unit must be J.
    terms
        Number of Taylor series terms used for extrapolation. Note that `energy` must contain columns
        named `term_1`, `term_2`, ..., `term_n` where n is the number of terms.

    Returns
    -------
    MPD
        Extrapolated MPD.
    """
    from math import factorial

    if terms < 1:
        raise ValueError("Number of terms must be at least 1.")

    if energy is None:
        if "term_1" in self._dataframe.columns:
            energy = self._dataframe[["macrostate", "term_1"]].copy()
        else:
            raise ValueError("Energy related data is missing.")

    beta = temperature_to_beta(temperature)
    delta_beta = beta - self.beta
    lnp_extrapolated = self.lnp.copy()
    lnp_extrapolated["lnp"] += (
        self.mu * lnp_extrapolated["macrostate"] - energy["term_1"]
    ) * delta_beta
    lnp_extrapolated["lnp"] = normalize(lnp_extrapolated["lnp"])

    for i in range(2, terms + 1):
        lnp_extrapolated["lnp"] += (
            1 / factorial(i) * energy[f"term_{i}"] * np.power(delta_beta, i)
        )
        lnp_extrapolated["lnp"] = normalize(lnp_extrapolated["lnp"])

    return MPD(
        dataframe=lnp_extrapolated,
        temperature=temperature,
        fugacity=mu_to_fugacity(self.mu, beta),
        metadata=self.metadata,
    )

find_phase_equilibrium

find_phase_equilibrium(
    delta_beta_mu_guess: Optional[float] = None,
    tolerance: float = 1e-06,
    return_probabilities: bool = False,
    max_delta_beta_mu: float = 8.0,
) -> Union[Tuple[float, float, float], float]

Find the fugacity at which the two phases are in equilibrium.

Uses a bounded phase-balance search. The search follows the physical direction of the current MPD: if the low-density phase dominates it shifts toward higher fugacity, and if the high-density phase dominates it shifts toward lower fugacity.

Parameters:

• delta_beta_mu_guess (Optional[float], default: None ) –

  Initial hint for the shift in beta*mu that brings the distribution closer to equilibrium. If provided, its sign is tried first and its magnitude is used as the first bracketing step. If None, the direction is auto-detected from the current distribution shape.

• tolerance (float, default: 1e-06 ) –

  Tolerance for the phase-probability balance.

• return_probabilities (bool, default: False ) –

  Whether to return the probabilities of the two phases at equilibrium.

• max_delta_beta_mu (float, default: 8.0 ) –

  Maximum absolute shift in beta*mu explored during bracketing. Since fugacity scales as exp(delta_beta_mu), this bounds the fugacity search and prevents excursions to unphysical values.

Returns:

• float or Tuple[float, float, float] –

  The fugacity at which the two phases are in equilibrium. If return_probabilities is True, returns (fugacity, p_low, p_high).

Raises:

• RuntimeError –

  If no phase equilibrium is found (distribution remains unimodal inside the bounded fugacity search).

Source code in src/asaf/mpd.py

def find_phase_equilibrium(
    self,
    delta_beta_mu_guess: Optional[float] = None,
    tolerance: float = 1e-6,
    return_probabilities: bool = False,
    max_delta_beta_mu: float = 8.0,
) -> Union[Tuple[float, float, float], float]:
    """Find the fugacity at which the two phases are in equilibrium.

    Uses a bounded phase-balance search.  The search follows the physical
    direction of the current MPD: if the low-density phase dominates it
    shifts toward higher fugacity, and if the high-density phase dominates
    it shifts toward lower fugacity.

    Parameters
    ----------
    delta_beta_mu_guess
        Initial hint for the shift in beta*mu that brings the distribution
        closer to equilibrium.  If provided, its sign is tried first and
        its magnitude is used as the first bracketing step.  If ``None``,
        the direction is auto-detected from the current distribution shape.
    tolerance
        Tolerance for the phase-probability balance.
    return_probabilities
        Whether to return the probabilities of the two phases at equilibrium.
    max_delta_beta_mu
        Maximum absolute shift in beta*mu explored during bracketing.  Since
        fugacity scales as ``exp(delta_beta_mu)``, this bounds the fugacity
        search and prevents excursions to unphysical values.

    Returns
    -------
    float or Tuple[float, float, float]
        The fugacity at which the two phases are in equilibrium. If
        ``return_probabilities`` is True, returns
        ``(fugacity, p_low, p_high)``.

    Raises
    ------
    RuntimeError
        If no phase equilibrium is found (distribution remains unimodal
        inside the bounded fugacity search).
    """
    from scipy.special import logsumexp

    if not np.isfinite(max_delta_beta_mu) or max_delta_beta_mu <= 0:
        raise ValueError("`max_delta_beta_mu` must be finite and positive.")

    macrostate = self.lnp["macrostate"].to_numpy(dtype=float)
    lnp = self.lnp["lnp"].to_numpy(dtype=float)
    macrostate_span = macrostate[-1] - macrostate[0]

    # Gas-like phases in adsorption MPDs sit near the low-N boundary.  A
    # unimodal distribution outside that region is treated as the dense
    # phase and nudged to lower fugacity.
    single_phase_low_fraction = 0.25

    def phase_balance(
        delta_beta_mu: float,
    ) -> dict[str, bool | float | None]:
        lnp_rw = lnp + delta_beta_mu * macrostate
        lnp_rw = lnp_rw - logsumexp(lnp_rw)
        min_loc = argrelextrema(lnp_rw, np.less, order=self._order)[0]
        min_loc = min_loc[(10 < min_loc) & (min_loc < lnp_rw.shape[0] - 10)]

        if len(min_loc) > 0:
            min_idx = int(min_loc[np.argmin(lnp_rw[min_loc])])
            p_low = float(np.exp(logsumexp(lnp_rw[:min_idx])))
            p_high = float(np.exp(logsumexp(lnp_rw[min_idx + 1 :])))
            return {
                "delta": delta_beta_mu,
                "balance": p_high - p_low,
                "two_phase": True,
                "p_low": p_low,
                "p_high": p_high,
            }

        mode_idx = int(np.argmax(lnp_rw))
        mode_fraction = (
            (macrostate[mode_idx] - macrostate[0]) / macrostate_span
            if macrostate_span > 0
            else 0.0
        )
        balance = -1.0 if mode_fraction <= single_phase_low_fraction else 1.0
        return {
            "delta": delta_beta_mu,
            "balance": balance,
            "two_phase": False,
            "p_low": None,
            "p_high": None,
        }

    def best_two_phase(
        best: dict[str, bool | float | None] | None,
        candidate: dict[str, bool | float | None],
        prefer_low_side: bool = False,
    ) -> dict[str, bool | float | None] | None:
        if not candidate["two_phase"]:
            return best
        if prefer_low_side and float(candidate["balance"]) > 0:
            return best
        if best is None:
            return candidate
        if abs(float(candidate["balance"])) < abs(float(best["balance"])):
            return candidate
        return best

    def bracket_in_direction(
        direction: float,
        start: dict[str, bool | float | None],
    ) -> (
        tuple[
            dict[str, bool | float | None],
            dict[str, bool | float | None],
            dict[str, bool | float | None] | None,
            dict[str, bool | float | None] | None,
        ]
        | None
    ):
        previous = start
        best = best_two_phase(None, start)
        best_low_side = best_two_phase(None, start, prefer_low_side=True)
        step = abs(delta_beta_mu_guess) if delta_beta_mu_guess else 0.05
        step = min(max(step, 0.05), max_delta_beta_mu)

        while step <= max_delta_beta_mu:
            delta = direction * step
            current = phase_balance(delta)
            best = best_two_phase(best, current)
            best_low_side = best_two_phase(
                best_low_side, current, prefer_low_side=True
            )

            if float(previous["balance"]) * float(current["balance"]) <= 0:
                return previous, current, best, best_low_side

            if step == max_delta_beta_mu:
                break

            previous = current
            step = min(step * 1.6, max_delta_beta_mu)

        return None

    start = phase_balance(0.0)
    if (
        start["two_phase"]
        and abs(float(start["balance"])) <= tolerance
        and float(start["balance"]) <= 0
    ):
        best = start
    else:
        if delta_beta_mu_guess is None or delta_beta_mu_guess == 0:
            direction = 1.0 if float(start["balance"]) < 0 else -1.0
        else:
            direction = float(np.sign(delta_beta_mu_guess))
        bracket = bracket_in_direction(direction, start)
        if bracket is None:
            bracket = bracket_in_direction(-direction, start)
        if bracket is None:
            raise RuntimeError(
                "No phase equilibrium found within the bounded fugacity search."
            )

        left, right, best, best_low_side = bracket
        if float(left["delta"]) > float(right["delta"]):
            left, right = right, left

        for _ in range(80):
            if (
                best_low_side is not None
                and abs(float(best_low_side["balance"])) <= tolerance
            ):
                best = best_low_side
                break

            midpoint = 0.5 * (float(left["delta"]) + float(right["delta"]))
            middle = phase_balance(midpoint)
            best = best_two_phase(best, middle)
            best_low_side = best_two_phase(
                best_low_side, middle, prefer_low_side=True
            )

            if float(left["balance"]) * float(middle["balance"]) <= 0:
                right = middle
            else:
                left = middle

        if best is None:
            raise RuntimeError(
                "No phase equilibrium found: distribution remains unimodal."
            )
        if (
            best_low_side is not None
            and abs(float(best_low_side["balance"])) <= tolerance
        ):
            best = best_low_side

    delta_beta_mu_eq = float(best["delta"])

    # Final validation with check_tail enabled
    lnp_eq = self.reweight(delta_beta_mu_eq)
    mins = self.minimums(order=self._order, lnp=lnp_eq["lnp"])

    if len(mins) == 0:
        raise RuntimeError(
            "No phase equilibrium found: distribution remains unimodal."
        )

    min_idx = int(mins[mins.lnp == mins.lnp.min()].index[0])
    p_low = float(np.exp(logsumexp(lnp_eq["lnp"].iloc[:min_idx])))
    p_high = float(np.exp(logsumexp(lnp_eq["lnp"].iloc[min_idx + 1 :])))

    equilibrium_beta_mu = self.beta_mu + delta_beta_mu_eq
    equilibrium_fugacity = mu_to_fugacity(
        equilibrium_beta_mu / self._beta, self._beta
    )

    if return_probabilities:
        return equilibrium_fugacity, p_low, p_high
    return equilibrium_fugacity

free_energy_at_fugacity

free_energy_at_fugacity(fug: float) -> DataFrame

Calculate the free energy profile at a given fugacity.

Source code in src/asaf/mpd.py

def free_energy_at_fugacity(self, fug: float) -> pd.DataFrame:
    """Calculate the free energy profile at a given fugacity."""
    beta_0 = self._beta
    mu_0 = self._mu
    mu = fugacity_to_mu(fug, beta_0)
    delta_beta_mu = beta_0 * (mu - mu_0)
    lnp_rw = self.reweight(delta_beta_mu)
    free_en = (
        -0.001
        * _BOLTZMANN_CONSTANT
        * _AVOGADRO_CONSTANT
        * self._temperature
        * lnp_rw["lnp"]
    )
    free_en -= free_en.min()
    free_energy = pd.DataFrame(
        {"macrostate": lnp_rw["macrostate"].copy(), "free_energy_kJ/mol": free_en}
    )

    return free_energy

from_csv classmethod

from_csv(file_name: str, **kwargs: object) -> MPD

Read natural logarithm of macrostates probability or transition probabilities from a csv file.

Source code in src/asaf/mpd.py

@classmethod
def from_csv(cls, file_name: str, **kwargs: object) -> MPD:
    """Read natural logarithm of macrostates probability or transition probabilities from a csv file."""
    df = pd.read_csv(file_name, **kwargs)

    metadata_filename = file_name.removesuffix(".csv") + ".metadata.json"
    with open(metadata_filename) as f:
        metadata = json.load(f)

    temperature = metadata.get("temperature")

    if temperature is None:
        raise ValueError("Metadata must contain 'temperature'.")

    beta_mu = metadata.get("beta_mu")
    fugacity = metadata.get("fugacity")

    if (beta_mu is None) and (fugacity is None):
        raise ValueError("Metadata must contain 'beta_mu' or/and 'fugacity'.")

    return cls(
        dataframe=df,
        temperature=temperature,
        beta_mu=beta_mu,
        fugacity=fugacity,
        metadata=metadata,
    )

minimums

minimums(
    order: int, lnp: Optional[Series] = None
) -> DataFrame

Find the local minimums in the lnp data.

Source code in src/asaf/mpd.py

def minimums(self, order: int, lnp: Optional[pd.Series] = None) -> pd.DataFrame:
    """Find the local minimums in the lnp data."""
    if lnp is None:
        lnp = self._dataframe["lnp"]
    elif isinstance(lnp, pd.DataFrame):
        lnp = lnp["lnp"]
    min_loc = argrelextrema(lnp.values, np.less, order=order)[0]
    min_loc = min_loc[(10 < min_loc) & (min_loc < lnp.shape[0] - 10)]
    minima = self._dataframe.iloc[min_loc][["macrostate"]].copy()
    minima["lnp"] = lnp.iloc[min_loc].to_numpy()
    return minima

plot

plot(
    fig: Optional[Figure] = None,
    name: Optional[str] = None,
    show: bool = True,
) -> None

Plot the MPD data using plotly.

Source code in src/asaf/mpd.py

def plot(
    self,
    fig: Optional[go.Figure] = None,
    name: Optional[str] = None,
    show: bool = True,
) -> None:
    """Plot the MPD data using plotly."""
    font = {"family": "Helvetica Neue", "size": 14, "color": "black"}

    axes = {
        "showline": True,
        "linewidth": 1,
        "linecolor": "black",
        "gridcolor": "lightgrey",
        "mirror": True,
        "zeroline": False,
        "ticks": "inside",
    }

    if fig is None:
        fig = go.Figure()

    fig.add_trace(
        go.Scatter(
            x=self.lnp["macrostate"],
            y=self.lnp["lnp"],
            mode="lines",
            name=name,
        )
    )

    xaxis_title = "Macrostate"
    yaxis_title = "lnΠ"

    fig.update_layout(
        font=font,
        xaxis=axes,
        xaxis_title=xaxis_title,
        yaxis=axes,
        yaxis_title=yaxis_title,
        plot_bgcolor="white",
        width=700,
        height=500,
        margin=dict(l=30, r=30, t=30, b=30),
    )

    if show:
        fig.show()

reweight

reweight(delta_beta_mu: float) -> DataFrame

Reweight the MPD to a new mu / fugacity value using delta_beta_mu.

Source code in src/asaf/mpd.py

def reweight(self, delta_beta_mu: float) -> pd.DataFrame:
    """Reweight the MPD to a new mu / fugacity value using `delta_beta_mu`."""
    lnp_rw = self.lnp.copy()
    lnp_rw["lnp"] += delta_beta_mu * lnp_rw["macrostate"]
    lnp_rw["lnp"] = normalize(lnp_rw["lnp"])
    self.check_tail(lnp=lnp_rw, order=self.order, tolerance=self.tolerance)

    return lnp_rw

reweight_to_fug

reweight_to_fug(
    fugacity: float, inplace: bool = True
) -> None | DataFrame

Reweight the MPD to a new mu / fugacity value using desired fugacity.

Source code in src/asaf/mpd.py

def reweight_to_fug(
    self, fugacity: float, inplace: bool = True
) -> None | pd.DataFrame:
    """Reweight the MPD to a new mu / fugacity value using desired fugacity."""
    beta_0 = self.beta
    mu_0 = self.mu
    mu = fugacity_to_mu(fugacity, beta_0)
    delta_beta_mu = beta_0 * (mu - mu_0)
    lnp_rw = self.reweight(delta_beta_mu)
    if inplace:
        self._dataframe["lnp"] = lnp_rw["lnp"]
        self.fugacity = fugacity
        return None
    else:
        return lnp_rw

options: filters: ["!^_"]