Module pychnosz.core.unicurve

CHNOSZ unicurve() function - Calculate univariant curves for geothermometry/geobarometry.

This module implements functions to solve for temperatures or pressures of equilibration for a given logK value, producing univariant curves useful for aqueous geothermometry and geobarometry applications.

Author: Based on pyCHNOSZ univariant.r by Grayson Boyer Optimized: Uses scipy.optimize.brentq() for efficient root-finding

Functions

def unicurve(logK: int | float | List[int | float],
species: str | List[str] | int | List[int],
coeff: int | float | List[int | float],
state: str | List[str],
pressures: float | List[float] = 1,
temperatures: float | List[float] = 25,
IS: float = 0,
minT: float = 0.1,
maxT: float = 100,
minP: float = 1,
maxP: float = 500,
tol: float | None = None,
solve: str = 'T',
messages: bool = True,
show: bool = True,
plot_it: bool = True,
width: int = 600,
height: int = 400,
res: int = 200) ‑> UnivariantResult | List[UnivariantResult]

Expand source code

def unicurve(logK: Union[float, int, List[Union[float, int]]],
             species: Union[str, List[str], int, List[int]],
             coeff: Union[int, float, List[Union[int, float]]],
             state: Union[str, List[str]],
             pressures: Union[float, List[float]] = 1,
             temperatures: Union[float, List[float]] = 25,
             IS: float = 0,
             minT: float = 0.1,
             maxT: float = 100,
             minP: float = 1,
             maxP: float = 500,
             tol: Optional[float] = None,
             solve: str = "T",
             messages: bool = True,
             show: bool = True,
             plot_it: bool = True,
             width: int = 600,
             height: int = 400,
             res: int = 200) -> Union[UnivariantResult, List[UnivariantResult]]:
    """
    Solve for temperatures or pressures of equilibration for a given logK value(s).

    This function calculates univariant curves useful for aqueous geothermometry
    and geobarometry. Given a measured equilibrium constant (logK) for a reaction,
    it solves for the temperatures (at specified pressures) or pressures (at
    specified temperatures) where the reaction would produce that logK value.

    The solver uses scipy.optimize.brentq (Brent's method), which combines
    bisection, secant, and inverse quadratic interpolation for efficient and
    robust convergence. This is ~100x faster than the original binary search
    algorithm while maintaining identical numerical accuracy.

    Parameters
    ----------
    logK : float, int, or list of float or int
        Logarithm (base 10) of the equilibrium constant(s). When a list is
        provided, each logK value is processed separately and a list of results
        is returned.
    species : str, int, or list of str or int
        Name, formula, or database index of species involved in the reaction
    coeff : int, float, or list
        Reaction stoichiometric coefficients (negative for reactants, positive for products)
    state : str or list of str
        Physical state(s) of species: "aq", "cr", "gas", "liq"
    pressures : float or list of float, default 1
        Pressure(s) in bars (used when solving for temperature)
    temperatures : float or list of float, default 25
        Temperature(s) in °C (used when solving for pressure)
    IS : float, default 0
        Ionic strength for activity corrections (mol/kg)
    minT : float, default 0.1
        Minimum temperature (°C) to search (ignored when solving for pressure)
    maxT : float, default 100
        Maximum temperature (°C) to search (ignored when solving for pressure)
    minP : float, default 1
        Minimum pressure (bar) to search (ignored when solving for temperature)
    maxP : float, default 500
        Maximum pressure (bar) to search (ignored when solving for temperature)
    tol : float, optional
        Tolerance for convergence. Default: 1/(10^(n+2)) where n is number of
        decimal places in logK, with maximum default of 1e-5
    solve : str, default "T"
        What to solve for: "T" for temperature or "P" for pressure
    messages : bool, default True
        Print informational messages
    show : bool, default True
        Display result table
    plot_it : bool, default True
        Display interactive plotly plot showing logK vs T (or P) with target logK
        as horizontal line and intersection points marked
    width : int, default 600
        Plot width in pixels (used if plot_it=True)
    height : int, default 400
        Plot height in pixels (used if plot_it=True)
    res : int, default 200
        Number of points to calculate for plotting the logK curve
        (used if plot_it=True)

    Returns
    -------
    UnivariantResult or list of UnivariantResult
        When logK is a single value: returns a UnivariantResult object.
        When logK is a list: returns a list of UnivariantResult objects.
        Each result contains:
        - reaction: DataFrame with reaction stoichiometry
        - out: DataFrame with solved T or P values and thermodynamic properties
        - warnings: List of warning messages

    Examples
    --------
    >>> from pychnosz import unicurve, reset
    >>> reset()
    >>>
    >>> # Solve for temperature: quartz dissolution
    >>> # SiO2(quartz) = SiO2(aq)
    >>> result = unicurve(logK=-2.71, species=["quartz", "SiO2"],
    ...                   state=["cr", "aq"], coeff=[-1, 1],
    ...                   pressures=200, minT=1, maxT=350)
    >>> print(result.out[["P", "T", "logK"]])
    >>>
    >>> # Solve for pressure: water dissociation
    >>> result = unicurve(logK=-14, species=["H2O", "H+", "OH-"],
    ...                   state=["liq", "aq", "aq"], coeff=[-1, 1, 1],
    ...                   temperatures=[25, 50, 75], solve="P",
    ...                   minP=1, maxP=1000)
    >>> print(result.out[["T", "P", "logK"]])

    Notes
    -----
    This function uses scipy.optimize.brentq for root-finding, which provides:
    - Guaranteed convergence if root is bracketed
    - Typical convergence in 5-15 function evaluations
    - ~100x speedup compared to custom binary search (1600 → 15 evaluations)
    - Identical numerical results to original implementation

    The algorithm also implements "warm start" optimization: when solving for
    multiple pressures/temperatures, previous solutions are used to intelligently
    bracket subsequent searches, further improving performance.

    References
    ----------
    Based on univariant.r from pyCHNOSZ by Grayson Boyer
    Optimized using Brent, R. P. (1973). Algorithms for Minimization without Derivatives.
    """
    # Track whether input was a single value or list
    single_logK_input = not isinstance(logK, list)

    # Ensure logK is a list for processing
    if single_logK_input:
        logK_list = [logK]
    else:
        logK_list = logK

    # Ensure species, state, and coeff are lists
    if not isinstance(species, list):
        species = [species]
    if not isinstance(state, list):
        state = [state]
    if not isinstance(coeff, list):
        coeff = [coeff]

    # Process each logK value
    results = []

    for this_logK in logK_list:
        result = UnivariantResult()

        # Set default tolerance based on logK precision
        if tol is None:
            # Count decimal places in logK
            logK_str = str(float(this_logK))
            if '.' in logK_str:
                n_decimals = len(logK_str.split('.')[1].rstrip('0'))
            else:
                n_decimals = 0
            this_tol = 10 ** (-(n_decimals + 2))
            if this_tol > 1e-5:
                this_tol = 1e-5
        else:
            this_tol = tol

        # Get reaction information from first subcrt call
        try:
            initial_calc = subcrt(species, coeff=coeff, state=state, T=25, P=1,
                                 exceed_Ttr=True, messages=False, show=False)
            result.reaction = initial_calc.reaction
        except Exception as e:
            if messages:
                warnings.warn(f"Error getting reaction information: {str(e)}")
            result.reaction = None

        if solve.upper() == "T":
            # Solve for temperature at given pressure(s)
            if not isinstance(pressures, list):
                pressures = [pressures]

            results_list = []
            prev_T = None  # For warm start optimization

            for i, pressure in enumerate(pressures):
                if messages:
                    print(f"Solving for T at P = {pressure} bar (logK = {this_logK})...")

                # Warm start: use previous solution to narrow search range if available
                current_minT = minT
                current_maxT = maxT
                if prev_T is not None and minT < prev_T < maxT:
                    # Center search around previous solution with a safety margin
                    # logK typically changes by ~0.006 per °C, so ±50°C should be safe
                    margin = 50
                    current_minT = max(minT, prev_T - margin)
                    current_maxT = min(maxT, prev_T + margin)
                    if messages:
                        print(f"  Using warm start: searching {current_minT:.1f} to {current_maxT:.1f}°C")

                result_dict = _solve_T_for_pressure(this_logK, species, state, coeff, pressure,
                                           IS, current_minT, current_maxT, this_tol,
                                           initial_guess=prev_T, messages=messages)

                # If warm start failed, try full range
                if result_dict['T'] is None and prev_T is not None:
                    if messages:
                        print(f"  Warm start failed, searching full range...")
                    result_dict = _solve_T_for_pressure(this_logK, species, state, coeff, pressure,
                                               IS, minT, maxT, this_tol, messages=messages)

                results_list.append(result_dict)

                # Update for next warm start
                if result_dict['T'] is not None:
                    prev_T = result_dict['T']

            result.out = pd.DataFrame(results_list)

        elif solve.upper() == "P":
            # Solve for pressure at given temperature(s)
            if not isinstance(temperatures, list):
                temperatures = [temperatures]

            results_list = []
            prev_P = None  # For warm start optimization

            for i, temperature in enumerate(temperatures):
                if messages:
                    print(f"Solving for P at T = {temperature} °C (logK = {this_logK})...")

                # Warm start: use previous solution to narrow search range if available
                current_minP = minP
                current_maxP = maxP
                if prev_P is not None and minP < prev_P < maxP:
                    # Center search around previous solution with a safety margin
                    # Pressure effects vary, use a generous ±500 bar margin
                    margin = 500
                    current_minP = max(minP, prev_P - margin)
                    current_maxP = min(maxP, prev_P + margin)
                    if messages:
                        print(f"  Using warm start: searching {current_minP:.0f} to {current_maxP:.0f} bar")

                result_dict = _solve_P_for_temperature(this_logK, species, state, coeff, temperature,
                                              IS, current_minP, current_maxP, this_tol,
                                              initial_guess=prev_P, messages=messages)

                # If warm start failed, try full range
                if result_dict['P'] is None and prev_P is not None:
                    if messages:
                        print(f"  Warm start failed, searching full range...")
                    result_dict = _solve_P_for_temperature(this_logK, species, state, coeff, temperature,
                                                  IS, minP, maxP, this_tol, messages=messages)

                results_list.append(result_dict)

                # Update for next warm start
                if result_dict['P'] is not None:
                    prev_P = result_dict['P']

            result.out = pd.DataFrame(results_list)

        else:
            raise ValueError(f"solve must be 'T' or 'P', got '{solve}'")

        # Create interactive plot if requested
        if plot_it:
            if not PLOTLY_AVAILABLE:
                warnings.warn("plotly is not installed. Set plot_it=False to suppress this warning, "
                             "or install plotly with: pip install plotly")
            else:
                result.fig = _create_unicurve_plot(this_logK, species, state, coeff, result, solve,
                                                   minT, maxT, minP, maxP, IS, width, height, res, messages)

        # Display result if requested
        if show and result.out is not None:
            try:
                from IPython.display import display
                if result.reaction is not None:
                    print("\nReaction:")
                    display(result.reaction)
                print(f"\nResults (logK = {this_logK}):")
                display(result.out)
            except ImportError:
                # Not in Jupyter, just print
                if result.reaction is not None:
                    print("\nReaction:")
                    print(result.reaction)
                print(f"\nResults (logK = {this_logK}):")
                print(result.out)

        # Add this result to the list
        results.append(result)

    # Return single result or list based on input
    if single_logK_input:
        return results[0]
    else:
        return results

Solve for temperatures or pressures of equilibration for a given logK value(s).

This function calculates univariant curves useful for aqueous geothermometry and geobarometry. Given a measured equilibrium constant (logK) for a reaction, it solves for the temperatures (at specified pressures) or pressures (at specified temperatures) where the reaction would produce that logK value.

The solver uses scipy.optimize.brentq (Brent's method), which combines bisection, secant, and inverse quadratic interpolation for efficient and robust convergence. This is ~100x faster than the original binary search algorithm while maintaining identical numerical accuracy.

Parameters

logK : float, int, or list of float or int

Logarithm (base 10) of the equilibrium constant(s). When a list is provided, each logK value is processed separately and a list of results is returned.

species : str, int, or list of str or int

Name, formula, or database index of species involved in the reaction

coeff : int, float, or list

Reaction stoichiometric coefficients (negative for reactants, positive for products)

state : str or list of str

Physical state(s) of species: "aq", "cr", "gas", "liq"

pressures : float or list of float, default 1

Pressure(s) in bars (used when solving for temperature)

temperatures : float or list of float, default 25

Temperature(s) in °C (used when solving for pressure)

IS : float, default 0

Ionic strength for activity corrections (mol/kg)

minT : float, default 0.1

Minimum temperature (°C) to search (ignored when solving for pressure)

maxT : float, default 100

Maximum temperature (°C) to search (ignored when solving for pressure)

minP : float, default 1

Minimum pressure (bar) to search (ignored when solving for temperature)

maxP : float, default 500

Maximum pressure (bar) to search (ignored when solving for temperature)

tol : float, optional

Tolerance for convergence. Default: 1/(10^(n+2)) where n is number of decimal places in logK, with maximum default of 1e-5

solve : str, default "T"

What to solve for: "T" for temperature or "P" for pressure

messages : bool, default True

Print informational messages

show : bool, default True

Display result table

plot_it : bool, default True

Display interactive plotly plot showing logK vs T (or P) with target logK as horizontal line and intersection points marked

width : int, default 600

Plot width in pixels (used if plot_it=True)

height : int, default 400

Plot height in pixels (used if plot_it=True)

res : int, default 200

Number of points to calculate for plotting the logK curve (used if plot_it=True)

Returns

UnivariantResult or list of UnivariantResult

When logK is a single value: returns a UnivariantResult object. When logK is a list: returns a list of UnivariantResult objects. Each result contains: - reaction: DataFrame with reaction stoichiometry - out: DataFrame with solved T or P values and thermodynamic properties - warnings: List of warning messages

Examples

>>> from pychnosz import unicurve, reset
>>> reset()
>>>
>>> # Solve for temperature: quartz dissolution
>>> # SiO2(quartz) = SiO2(aq)
>>> result = unicurve(logK=-2.71, species=["quartz", "SiO2"],
...                   state=["cr", "aq"], coeff=[-1, 1],
...                   pressures=200, minT=1, maxT=350)
>>> print(result.out[["P", "T", "logK"]])
>>>
>>> # Solve for pressure: water dissociation
>>> result = unicurve(logK=-14, species=["H2O", "H+", "OH-"],
...                   state=["liq", "aq", "aq"], coeff=[-1, 1, 1],
...                   temperatures=[25, 50, 75], solve="P",
...                   minP=1, maxP=1000)
>>> print(result.out[["T", "P", "logK"]])

Notes

This function uses scipy.optimize.brentq for root-finding, which provides: - Guaranteed convergence if root is bracketed - Typical convergence in 5-15 function evaluations - ~100x speedup compared to custom binary search (1600 → 15 evaluations) - Identical numerical results to original implementation

The algorithm also implements "warm start" optimization: when solving for multiple pressures/temperatures, previous solutions are used to intelligently bracket subsequent searches, further improving performance.

References

Based on univariant.r from pyCHNOSZ by Grayson Boyer Optimized using Brent, R. P. (1973). Algorithms for Minimization without Derivatives.

def univariant_TP(logK: int | float | List[int | float],
species: str | List[str] | int | List[int],
coeff: int | float | List[int | float],
state: str | List[str],
Trange: List[float],
Prange: List[float],
IS: float = 0,
xlim: List[float] | None = None,
ylim: List[float] | None = None,
line_type: str = 'markers+lines',
tol: float | None = None,
title: str | None = None,
res: int = 10,
width: int = 500,
height: int = 400,
save_as: str | None = None,
save_format: str = 'png',
save_scale: float = 1,
show: bool = False,
messages: bool = False,
parallel: bool = True,
plot_it: bool = True) ‑> List[UnivariantResult]

Expand source code

def univariant_TP(logK: Union[float, int, List[Union[float, int]]],
                  species: Union[str, List[str], int, List[int]],
                  coeff: Union[int, float, List[Union[int, float]]],
                  state: Union[str, List[str]],
                  Trange: List[float],
                  Prange: List[float],
                  IS: float = 0,
                  xlim: Optional[List[float]] = None,
                  ylim: Optional[List[float]] = None,
                  line_type: str = "markers+lines",
                  tol: Optional[float] = None,
                  title: Optional[str] = None,
                  res: int = 10,
                  width: int = 500,
                  height: int = 400,
                  save_as: Optional[str] = None,
                  save_format: str = "png",
                  save_scale: float = 1,
                  show: bool = False,
                  messages: bool = False,
                  parallel: bool = True,
                  plot_it: bool = True) -> List[UnivariantResult]:
    """
    Solve for temperatures and pressures of equilibration for given logK value(s)
    and produce an interactive T-P diagram.

    This function calculates univariant curves in temperature-pressure (T-P) space
    for one or more logK values. For each pressure in a range, it solves for the
    temperature where the reaction achieves the target logK. The resulting curves
    show phase boundaries or equilibrium conditions in T-P space.

    Parameters
    ----------
    logK : float, int, or list
        Logarithm (base 10) of equilibrium constant(s). Multiple values produce
        multiple curves on the same plot.
    species : str, int, or list of str or int
        Name, formula, or database index of species involved in the reaction
    coeff : int, float, or list
        Reaction stoichiometric coefficients (negative for reactants, positive for products)
    state : str or list of str
        Physical state(s) of species: "aq", "cr", "gas", "liq"
    Trange : list of two floats
        [min, max] temperature range (°C) to search for solutions
    Prange : list of two floats
        [min, max] pressure range (bar) to calculate along
    IS : float, default 0
        Ionic strength for activity corrections (mol/kg)
    xlim : list of two floats, optional
        [min, max] range for x-axis (temperature) in plot
    ylim : list of two floats, optional
        [min, max] range for y-axis (pressure) in plot
    line_type : str, default "markers+lines"
        Plotly line type: "markers+lines", "markers", or "lines"
    tol : float, optional
        Convergence tolerance. Default: 1/(10^(n+2)) where n is decimal places in logK
    title : str, optional
        Plot title. Default: auto-generated from reaction
    res : int, default 10
        Number of pressure points to calculate along the curve
    width : int, default 500
        Plot width in pixels
    height : int, default 400
        Plot height in pixels
    save_as : str, optional
        Filename to save plot (without extension)
    save_format : str, default "png"
        Save format: "png", "jpg", "jpeg", "webp", "svg", "pdf", "html"
    save_scale : float, default 1
        Scale factor for saved plot
    show : bool, default False
        Display subcrt result tables
    messages : bool, default False
        Print informational messages
    parallel : bool, default True
        Use parallel processing across multiple logK values for faster computation.
        Utilizes multiple CPU cores when processing multiple logK curves.
    plot_it : bool, default True
        Display the plot

    Returns
    -------
    list of UnivariantResult
        List of UnivariantResult objects, one for each logK value.
        Each contains reaction information and T-P curve data.

    Examples
    --------
    >>> from pychnosz import univariant_TP, reset
    >>> reset()
    >>>
    >>> # Calcite-aragonite phase boundary
    >>> result = univariant_TP(
    ...     logK=0,
    ...     species=["calcite", "aragonite"],
    ...     state=["cr", "cr"],
    ...     coeff=[-1, 1],
    ...     Trange=[0, 700],
    ...     Prange=[2000, 16000]
    ... )
    >>>
    >>> # Multiple curves for K-feldspar stability
    >>> result = univariant_TP(
    ...     logK=[-8, -6, -4, -2],
    ...     species=["K-feldspar", "kaolinite", "H2O", "SiO2", "muscovite"],
    ...     state=["cr", "cr", "liq", "aq", "cr"],
    ...     coeff=[-1, -1, 1, 2, 1],
    ...     Trange=[0, 350],
    ...     Prange=[1, 5000],
    ...     res=20
    ... )

    Notes
    -----
    This function creates T-P diagrams by:
    1. Generating a range of pressures from Prange[0] to Prange[1]
    2. For each pressure, solving for T where logK matches the target
    3. Plotting the resulting T-P points as a curve

    For multiple logK values, each curve represents a different equilibrium
    condition. This is useful for:
    - Phase diagrams (e.g., mineral stability fields)
    - Isopleths (lines of constant logK)
    - Reaction boundaries

    Requires plotly for interactive plotting. If plotly is not installed,
    set plot_it=False to just return the data without plotting.

    References
    ----------
    Based on univariant_TP from pyCHNOSZ by Grayson Boyer
    """

    # Check if plotly is available
    if plot_it and not PLOTLY_AVAILABLE:
        warnings.warn("plotly is not installed. Set plot_it=False to suppress this warning, "
                     "or install plotly with: pip install plotly")
        plot_it = False

    # Ensure logK is a list
    if not isinstance(logK, list):
        logK = [logK]

    # Create plotly figure
    if plot_it:
        fig = go.Figure()

    output = []

    # Generate pressure array
    pressures = np.linspace(Prange[0], Prange[1], res)

    # Process each logK value (in parallel if enabled)
    if parallel and len(logK) > 1:
        # Parallel processing
        max_workers = min(len(logK), multiprocessing.cpu_count())

        # Prepare arguments for each logK value
        args_list = [
            (this_logK, species, state, coeff, pressures, Trange, IS, tol, show, messages)
            for this_logK in logK
        ]

        # Process in parallel
        with ProcessPoolExecutor(max_workers=max_workers) as executor:
            # Submit all tasks
            future_to_logK = {
                executor.submit(_process_single_logK, args): args[0]
                for args in args_list
            }

            # Collect results as they complete (maintains order via logK list)
            results_dict = {}
            for future in as_completed(future_to_logK):
                this_logK = future_to_logK[future]
                try:
                    out = future.result()
                    results_dict[this_logK] = out
                except Exception as e:
                    if messages:
                        print(f"Error processing logK={this_logK}: {str(e)}")
                    # Create empty result
                    results_dict[this_logK] = None

            # Reorder results to match input logK order
            for this_logK in logK:
                out = results_dict.get(this_logK)
                if out is not None:
                    output.append(out)

                    # Add to plot if we have valid data
                    if plot_it and not out.out['T'].isnull().all():
                        fig.add_trace(go.Scatter(
                            x=out.out['T'],
                            y=out.out['P'],
                            mode=line_type,
                            name=f"logK={this_logK}",
                            text=[f"logK={this_logK}" for _ in range(len(out.out['T']))],
                            hovertemplate='%{text}<br>T, °C=%{x:.2f}<br>P, bar=%{y:.2f}<extra></extra>',
                        ))
                    elif out.out['T'].isnull().all():
                        if messages:
                            print(f"Could not find any T or P values in this range that correspond to a logK value of {this_logK}")

    else:
        # Sequential processing (original code)
        for this_logK in logK:
            # Set tolerance if not specified
            if tol is None:
                logK_str = str(float(this_logK))
                if '.' in logK_str:
                    n_decimals = len(logK_str.split('.')[1].rstrip('0'))
                else:
                    n_decimals = 0
                this_tol = 10 ** (-(n_decimals + 2))
                if this_tol > 1e-5:
                    this_tol = 1e-5
            else:
                this_tol = tol

            # Solve for T at each pressure
            out = unicurve(
                solve="T",
                logK=this_logK,
                species=species,
                state=state,
                coeff=coeff,
                pressures=list(pressures),
                minT=Trange[0],
                maxT=Trange[1],
                IS=IS,
                tol=this_tol,
                show=show,
                messages=messages,
                plot_it=False  # Don't plot individual curves - univariant_TP makes its own plot
            )

            # Add to plot if we have valid data
            if plot_it and not out.out['T'].isnull().all():
                fig.add_trace(go.Scatter(
                    x=out.out['T'],
                    y=out.out['P'],
                    mode=line_type,
                    name=f"logK={this_logK}",
                    text=[f"logK={this_logK}" for _ in range(len(out.out['T']))],
                    hovertemplate='%{text}<br>T, °C=%{x:.2f}<br>P, bar=%{y:.2f}<extra></extra>',
                ))
            elif out.out['T'].isnull().all():
                if messages:
                    print(f"Could not find any T or P values in this range that correspond to a logK value of {this_logK}")

            output.append(out)

    # Generate plot title if not specified
    if plot_it:
        if title is None and len(output) > 0 and output[0].reaction is not None:
            react_grid = output[0].reaction

            # Build reaction string
            reactants = []
            products = []
            for i, row in react_grid.iterrows():
                coeff_val = row['coeff']
                name = row['name'] if row['name'] != 'water' else 'H2O'

                if coeff_val < 0:
                    coeff_str = str(int(-coeff_val)) if -coeff_val != 1 else ""
                    reactants.append(f"{coeff_str} {name}".strip())
                elif coeff_val > 0:
                    coeff_str = str(int(coeff_val)) if coeff_val != 1 else ""
                    products.append(f"{coeff_str} {name}".strip())

            title = " + ".join(reactants) + " = " + " + ".join(products)

        # Update layout
        fig.update_layout(
            template="simple_white",
            title=str(title) if title else "",
            xaxis_title="T, °C",
            yaxis_title="P, bar",
            width=width,
            height=height,
            hoverlabel=dict(bgcolor="white"),
        )

        # Set axis limits if specified
        if xlim is not None:
            fig.update_xaxes(range=xlim)
        if ylim is not None:
            fig.update_yaxes(range=ylim)

        # Configure plot controls
        config = {
            'displaylogo': False,
            'modeBarButtonsToRemove': ['resetScale2d', 'toggleSpikelines'],
            'toImageButtonOptions': {
                'format': save_format,
                'filename': save_as if save_as else 'univariant_TP',
                'height': height,
                'width': width,
                'scale': save_scale,
            },
        }

        # Save plot if requested
        if save_as is not None:
            full_filename = f"{save_as}.{save_format}"
            if save_format == 'html':
                fig.write_html(full_filename)
            else:
                fig.write_image(full_filename, format=save_format,
                              width=width, height=height, scale=save_scale)
            if messages:
                print(f"Plot saved to {full_filename}")

        # Display plot
        fig.show(config=config)

        # Store figure in all result objects
        for out in output:
            out.fig = fig

    return output

Solve for temperatures and pressures of equilibration for given logK value(s) and produce an interactive T-P diagram.

This function calculates univariant curves in temperature-pressure (T-P) space for one or more logK values. For each pressure in a range, it solves for the temperature where the reaction achieves the target logK. The resulting curves show phase boundaries or equilibrium conditions in T-P space.

Parameters

logK : float, int, or list

Logarithm (base 10) of equilibrium constant(s). Multiple values produce multiple curves on the same plot.

species : str, int, or list of str or int

Name, formula, or database index of species involved in the reaction

coeff : int, float, or list

Reaction stoichiometric coefficients (negative for reactants, positive for products)

state : str or list of str

Physical state(s) of species: "aq", "cr", "gas", "liq"

Trange : list of two floats

[min, max] temperature range (°C) to search for solutions

Prange : list of two floats

[min, max] pressure range (bar) to calculate along

IS : float, default 0

Ionic strength for activity corrections (mol/kg)

xlim : list of two floats, optional

[min, max] range for x-axis (temperature) in plot

ylim : list of two floats, optional

[min, max] range for y-axis (pressure) in plot

line_type : str, default "markers+lines"

Plotly line type: "markers+lines", "markers", or "lines"

tol : float, optional

Convergence tolerance. Default: 1/(10^(n+2)) where n is decimal places in logK

title : str, optional

Plot title. Default: auto-generated from reaction

res : int, default 10

Number of pressure points to calculate along the curve

width : int, default 500

Plot width in pixels

height : int, default 400

Plot height in pixels

save_as : str, optional

Filename to save plot (without extension)

save_format : str, default "png"

Save format: "png", "jpg", "jpeg", "webp", "svg", "pdf", "html"

save_scale : float, default 1

Scale factor for saved plot

show : bool, default False

Display subcrt result tables

messages : bool, default False

Print informational messages

parallel : bool, default True

Use parallel processing across multiple logK values for faster computation. Utilizes multiple CPU cores when processing multiple logK curves.

plot_it : bool, default True

Display the plot

Returns

list of UnivariantResult

List of UnivariantResult objects, one for each logK value. Each contains reaction information and T-P curve data.

Examples

>>> from pychnosz import univariant_TP, reset
>>> reset()
>>>
>>> # Calcite-aragonite phase boundary
>>> result = univariant_TP(
...     logK=0,
...     species=["calcite", "aragonite"],
...     state=["cr", "cr"],
...     coeff=[-1, 1],
...     Trange=[0, 700],
...     Prange=[2000, 16000]
... )
>>>
>>> # Multiple curves for K-feldspar stability
>>> result = univariant_TP(
...     logK=[-8, -6, -4, -2],
...     species=["K-feldspar", "kaolinite", "H2O", "SiO2", "muscovite"],
...     state=["cr", "cr", "liq", "aq", "cr"],
...     coeff=[-1, -1, 1, 2, 1],
...     Trange=[0, 350],
...     Prange=[1, 5000],
...     res=20
... )

Notes

This function creates T-P diagrams by: 1. Generating a range of pressures from Prange[0] to Prange[1] 2. For each pressure, solving for T where logK matches the target 3. Plotting the resulting T-P points as a curve

For multiple logK values, each curve represents a different equilibrium condition. This is useful for: - Phase diagrams (e.g., mineral stability fields) - Isopleths (lines of constant logK) - Reaction boundaries

Requires plotly for interactive plotting. If plotly is not installed, set plot_it=False to just return the data without plotting.

References

Based on univariant_TP from pyCHNOSZ by Grayson Boyer

Classes

class UnivariantResult

Expand source code

class UnivariantResult:
    """Result structure for univariant curve calculations."""

    def __init__(self):
        self.reaction = None    # Reaction summary DataFrame
        self.out = None         # Results DataFrame with T/P and thermodynamic properties
        self.warnings = []      # Warning messages
        self.fig = None         # Plotly figure object (if plot_it=True)

    def __repr__(self):
        if self.out is not None:
            return f"UnivariantResult with {len(self.out)} points"
        return "UnivariantResult (no calculations performed)"

    def __getitem__(self, key):
        """Allow dictionary-style access to attributes."""
        return getattr(self, key)

Result structure for univariant curve calculations.