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,
    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 pressure to calculate the pressure in relative scale (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,
    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 pressure to calculate the pressure in relative scale (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 return_dataframe:
        return isotherm
    else:
        return Isotherm(
            data=isotherm,
            saturation_fugacity=saturation_fugacity,
            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,
) -> Union[Tuple[float, float, float], float]

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

Uses minimize_scalar (Brent's method) on an objective that is zero only when the low-density and high-density phase probabilities are equal. The objective is designed to guide the optimizer smoothly through unimodal, degenerate, and bimodal regimes.

Parameters:

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

  Initial hint for the shift in betamu that brings the distribution closer to equilibrium. Passed to minimize_scalar as a bracket starting point — the optimizer can search beyond* this value. If None, the direction is auto-detected from the current distribution shape.

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

  Tolerance for the minimization algorithm.

• return_probabilities (bool, default: False ) –

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

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 at the optimizer's solution).

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,
) -> Union[Tuple[float, float, float], float]:
    """Find the fugacity at which the two phases are in equilibrium.

    Uses ``minimize_scalar`` (Brent's method) on an objective that is zero
    only when the low-density and high-density phase probabilities are
    equal.  The objective is designed to guide the optimizer smoothly
    through unimodal, degenerate, and bimodal regimes.

    Parameters
    ----------
    delta_beta_mu_guess
        Initial hint for the shift in beta*mu that brings the distribution
        closer to equilibrium.  Passed to ``minimize_scalar`` as a bracket
        starting point — the optimizer **can search beyond** this value.
        If ``None``, the direction is auto-detected from the current
        distribution shape.
    tolerance
        Tolerance for the minimization algorithm.
    return_probabilities
        Whether to return the probabilities of the two phases at equilibrium.

    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
        at the optimizer's solution).
    """
    from scipy.optimize import minimize_scalar
    from scipy.special import logsumexp

    def objective(delta_beta_mu: float) -> float:
        lnp_rw = self.reweight(delta_beta_mu)
        mins = self.minimums(order=self._order, lnp=lnp_rw["lnp"])

        if len(mins) == 0:
            return 2.0 + float(
                lnp_rw["lnp"].iloc[0] - lnp_rw["lnp"].iloc[-1]
            ) ** 2

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

        if abs(p_low - 1) < 1e-6:
            return 1.1 + float(np.exp(-delta_beta_mu))
        if abs(p_high - 1) < 1e-6:
            return 1.1 + float(np.exp(delta_beta_mu))

        return abs(p_low - p_high)

    if delta_beta_mu_guess is None:
        edge_bias = float(self.lnp["lnp"].iloc[0] - self.lnp["lnp"].iloc[-1])
        delta_beta_mu_guess = 0.5 if edge_bias > 0 else -0.5

    self._suppress_check_tail = True
    try:
        result = minimize_scalar(
            objective,
            bracket=(0, delta_beta_mu_guess),
            tol=tolerance,
        )
    finally:
        self._suppress_check_tail = False

    delta_beta_mu_eq = result.x

    # 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: ["!^_"]