feat(calculators): garnet fe_basis + wt fractionation

src/phasetools/calculators/garnet.py

@@ -2,51 +2,94 @@
 import sys
 from .pt_grid import MAGEMinPTGridCalculator
 from ..core.phase_properties import phase_frac, extract_end_member, get_oxide_apfu
-from ..utils.bulk_rock import atomic_mass_dict, convert_mol_percent_to_wt_percent
+from ..utils.bulk_rock import atomic_mass_dict, atomic_frac_to_wt_frac
 from phasetools import MAGEMin_C
 from juliacall import Main as jl, convert as jlconvert
 
 class MAGEMinGarnetCalculator(MAGEMinPTGridCalculator):
     """High-level wrappers for garnet-focused MAGEMin calculations."""
-    def __init__(self, db="ig", dataset=636, verbose=False):
+
+    def __init__(
+        self, db: str = "ig", dataset: int = 636, verbose: bool = False, fe_basis: str = "FeOt"
+    ) -> None:
+        """
+        Parameters
+        ----------
+        db : str, default="ig"
+            Thermodynamic database label.
+        dataset : int, default=636
+            Thermodynamic dataset version.
+        verbose : bool, default=False
+            Whether to print progress information.
+        fe_basis : str, default="FeOt"
+            Iron basis used for garnet X-site fractions:
+
+            * ``'FeOt'`` -- all iron treated as Fe2+ (total iron,
+              ``FeO + 2*Fe2O3`` APFU) placed on the divalent site.  This is
+              the standard community convention for garnet end-members and
+              X-site fractions (e.g. Williams & Grambling 1990; Krogh Ravna
+              2000) and is consistent with the pyralspite garnet model used
+              by MAGEMin and Holland-Powell-type databases.
+            * ``'Fe2+'`` -- only the stoichiometrically estimated ferrous
+              iron (``_extract_fe_split_from_apfu``) is placed on the
+              divalent site, excluding Fe3+.  Use only when garnet Fe3+ is
+              known to be significant (oxidised eclogites, skarns) or
+              measured directly (XANES, Mössbauer).
+            * ``'Fe2'`` or ``'fe2'`` -- alias for ``'Fe2+'``.
+
+            Case-insensitive.
+        """
         super().__init__(db, dataset, verbose)
+        basis = str(fe_basis).casefold()
+        if basis not in ("feot", "fe2+", "fe2"):
+            raise ValueError(f"fe_basis must be 'FeOt' or 'Fe2+', got {fe_basis!r}")
+        self.fe_basis = "fe2+" if basis.startswith("fe2") else "feot"
 
     def _extract_garnet_elements_from_oxides(self, out, sys_in):
         """
         Extract garnet Mg-Mn-Fe-Ca cation fractions for the divalent (X) site.
-        Strictly isolates Fe2+ to ensure X-site fractions sum to 1.0.
+
+        The Fe basis is controlled by ``self.fe_basis``:
+
+        * ``'feot'`` (default) -- all iron treated as Fe2+ (total iron,
+          ``FeO + 2*Fe2O3`` APFU) placed on the divalent site.  The standard
+          community convention for garnet X-site fractions.
+        * ``'fe2+'`` -- only the stoichiometrically estimated ferrous iron
+          is placed on the divalent site, excluding Fe3+.
         """
         if 'g' not in out.ph:
             return 0.0, 0.0, 0.0, 0.0
 
-        # Isolating divalent iron using the robust Fe-split method
-        split = self._extract_fe_split_from_apfu(out, 'g')
-        fe2_moles = split['fe2']
-
-        elements = get_oxide_apfu(out, 'g', ['MgO', 'MnO', 'CaO'])
+        if self.fe_basis == 'feot':
+            # Total iron on the divalent site (all Fe as Fe2+)
+            elements = get_oxide_apfu(out, 'g', ['MgO', 'MnO', 'CaO', 'FeO', 'Fe2O3'])
+            fe_moles = elements.get("FeO", 0.0) + 2.0 * elements.get("Fe2O3", 0.0)
+        else:
+            # Ferrous iron only, from the excess-oxygen Fe2+/Fe3+ split
+            elements = get_oxide_apfu(out, 'g', ['MgO', 'MnO', 'CaO'])
+            fe_moles = self._extract_fe_split_from_apfu(out, 'g')['Fe2']
+
         mg_moles = elements.get("MgO", 0.0)
         mn_moles = elements.get("MnO", 0.0)
         ca_moles = elements.get("CaO", 0.0)
 
         # Total atoms in the X-site (should be approx 3.0)
-        total_x_site_moles = mg_moles + mn_moles + fe2_moles + ca_moles
-
+        total_x_site_moles = mg_moles + mn_moles + fe_moles + ca_moles
         if total_x_site_moles <= 0:
             return 0.0, 0.0, 0.0, 0.0
 
         # Normalise to 1.0 (mole fractions of the divalent site)
         Mg = mg_moles / total_x_site_moles
         Mn = mn_moles / total_x_site_moles
-        Fe = fe2_moles / total_x_site_moles
+        Fe = fe_moles / total_x_site_moles
         Ca = ca_moles / total_x_site_moles
 
         if sys_in.casefold() == 'wt':
-            wt_percent_list = convert_mol_percent_to_wt_percent(
-                [Mg, Mn, Fe, Ca],
-                ["Mg", "Mn", "Fe", "Ca"],
+            wt_frac = atomic_frac_to_wt_frac(
+                {"Mg": Mg, "Mn": Mn, "Fe": Fe, "Ca": Ca},
                 atomic_mass_dict,
             )
-            Mg, Mn, Fe, Ca = [val / 100.0 for val in wt_percent_list]
+            Mg, Mn, Fe, Ca = wt_frac["Mg"], wt_frac["Mn"], wt_frac["Fe"], wt_frac["Ca"]
 
         return Mg, Mn, Fe, Ca
 
@@ -127,7 +170,61 @@ def _gt_single_point_from_jl(self, P, T, X_jl, Xoxides_jl, sys_in, rm_list=None)
         return gt_frac, gt_wt, gt_vol, Mg, Mn, Fe, Ca, out
 
     def gt_along_path(self, P, T, fractionate=False, normalise_start=True):
-        """Calculate garnet fractions and element chemistry along a P-T path."""
+        """Calculate garnet fractions and element chemistry along a P-T path.
+
+        Parameters
+        ----------
+        P : array-like
+            Pressure values along the path (kbar).
+        T : array-like
+            Temperature values along the path (°C).
+        fractionate : bool, default=False
+            If True, fractionate garnet from the bulk composition as it grows.
+        normalise_start : bool, default=True
+            Controls how the first P-T point is treated:
+
+            * ``True`` — The first P-T point is treated as a nucleation
+              barrier: garnet volume starts at zero and only new growth is
+              modelled.  The initial garnet fraction (if any) is used as a
+              baseline for measuring incremental growth but is **not** removed
+              from the reactive bulk.  Use when the path starts outside the
+              garnet stability field or at the nucleation threshold.
+
+            * ``False`` — The first P-T point has an initial garnet volume
+              (overstepped nucleation).  That fraction is removed from the
+              bulk at step 0, and subsequent growth is measured relative to
+              it.  Use when the path starts well inside the garnet stability
+              field.
+
+        Returns
+        -------
+        gt_mol_frac : numpy.ndarray
+            Garnet molar fraction at each P-T point.
+        gt_wt_frac : numpy.ndarray
+            Garnet weight fraction at each P-T point.
+        gt_vol_frac : numpy.ndarray
+            Garnet volume fraction at each P-T point.
+        Mgi : numpy.ndarray
+            Garnet X-site Mg fraction at each P-T point.
+        Mni : numpy.ndarray
+            Garnet X-site Mn fraction at each P-T point.
+        Fei : numpy.ndarray
+            Garnet X-site Fe fraction at each P-T point.  The Fe basis
+            (total Fe as Fe2+ by default, or ferrous-only) is set by the
+            ``fe_basis`` argument passed to :meth:`__init__`.
+        Cai : numpy.ndarray
+            Garnet X-site Ca fraction at each P-T point.
+        X_along_path : numpy.ndarray
+            Bulk composition after each step's fractionation.  Each row is
+            normalised to sum to 1.  Row ``i`` is the bulk **after** step
+            ``i``'s fractionation has been applied.
+
+        Notes
+        -----
+        Fractionation uses the **current-step** garnet composition (not the
+        growth-increment composition) — a first-order approximation valid for
+        small P-T steps.
+        """
         from .phase_search import PhaseFunctions
 
         X = self.X
@@ -144,6 +241,7 @@ def gt_along_path(self, P, T, fractionate=False, normalise_start=True):
         X_along_path = np.zeros(shape=(n_points, len(self._Xoxides_py)) )
 
         gt_frac_max_previous = 0.
+        gt_wt_max_previous = 0.
         phase_functions = PhaseFunctions(db=self.db, dataset=self.dataset, verbose=self.verbose) if fractionate else None
         if phase_functions:
             # Sync standardised state to the helper instance
@@ -158,19 +256,28 @@ def gt_along_path(self, P, T, fractionate=False, normalise_start=True):
             gt_wt_frac[i] = gt_wt
             gt_vol_frac[i] = gt_vol
 
+            # Select the fraction basis consistent with sys_in
+            if self.sys_in.casefold() == 'wt':
+                gt_frac_for_fractionation = gt_wt
+                gt_frac_max_prev_for_fractionation = gt_wt_max_previous
+            else:
+                gt_frac_for_fractionation = gt_frac
+                gt_frac_max_prev_for_fractionation = gt_frac_max_previous
+
             if phase_functions is not None:
                 if i == 0 and not normalise_start:
-                    if gt_frac > 0:
-                        X_py = phase_functions.fractionate_phase('g', out, self.sys_in, frac_amount=gt_frac)
+                    if gt_frac_for_fractionation > 0:
+                        X_py = phase_functions.fractionate_phase('g', out, self.sys_in, frac_amount=gt_frac_for_fractionation)
                         X = jlconvert(jl.Vector[jl.Float64], X_py)
                 elif i > 0:
-                    frac_amount = max(gt_frac - gt_frac_max_previous, 0.0)
+                    frac_amount = max(gt_frac_for_fractionation - gt_frac_max_prev_for_fractionation, 0.0)
                     if frac_amount > 0:
                         X_py = phase_functions.fractionate_phase('g', out, self.sys_in, frac_amount=frac_amount)
                         X = jlconvert(jl.Vector[jl.Float64], X_py)
 
-            X_along_path[i] = np.array(X)
+            X_along_path[i] = np.array(X) / np.sum(X)
             gt_frac_max_previous = max(gt_frac_max_previous, gt_frac)
+            gt_wt_max_previous = max(gt_wt_max_previous, gt_wt)
 
             Mgi[i] = Mg
             Mni[i] = Mn

src/phasetools/calculators/phase_search.py

@@ -41,24 +41,72 @@ def liquidus_func(T):
         return result.root
 
     def fractionate_phase(self, phase, out, sys_in, frac_amount=None):
-        """Perform batch fractionation of a phase from the bulk rock composition."""
+        """Perform batch fractionation of a phase from the bulk rock composition.
+
+        Removes a specified fraction of a phase from the current bulk composition
+        and renormalises the remainder. The output is always normalised to sum to 1
+        (molar or weight fractions, depending on ``sys_in``).
+
+        MAGEMin orders ``out.ph`` with solution phases first (``SS_vec``), then
+        pure phases (``PP_vec``).  This method handles both indexing ranges
+        transparently.
+
+        .. note::
+
+            When a phase appears multiple times in ``out.ph`` (e.g. two
+            coexisting pyroxenes on a solvus), only the **first** instance is
+            removed.  Use ``phase_frac()`` to obtain the summed fraction before
+            calling this method if the total solvus fraction is needed.
+
+        Parameters
+        ----------
+        phase : str
+            Phase name as it appears in ``out.ph`` (e.g. ``'g'``, ``'liq'``).
+        out : MAGEMinOutput
+            Raw output from a MAGEMin single-point minimisation.
+        sys_in : str
+            ``'mol'`` or ``'wt'`` — determines which bulk composition and
+            phase composition vectors are used.
+        frac_amount : float or None, optional
+            Fraction of the phase to remove.  If ``None``, the full phase
+            fraction from the minimisation result is used.  Values of 0 or
+            less are treated as a no-op and return the current bulk unchanged.
+
+        Returns
+        -------
+        numpy.ndarray
+            Renormalised bulk composition after fractionation.
+        """
         if sys_in.casefold() == "wt":
             current_X = out.bulk_wt
         else:
             current_X = out.bulk
-        
+
         if phase in out.ph:
             phase_ind = out.ph.index(phase)
-            if sys_in.casefold() == "wt":
-                ph_comp = np.array(out.SS_vec[phase_ind].Comp_wt)
-                ph_frac = out.ph_frac_wt[phase_ind]
+            # MAGEMin orders: solution phases (SS_vec) first, then pure phases (PP_vec)
+            if phase_ind < out.n_SS:
+                if sys_in.casefold() == "wt":
+                    ph_comp = np.array(out.SS_vec[phase_ind].Comp_wt)
+                    ph_frac = out.ph_frac_wt[phase_ind]
+                else:
+                    ph_comp = np.array(out.SS_vec[phase_ind].Comp)
+                    ph_frac = out.ph_frac[phase_ind]
             else:
-                ph_comp = np.array(out.SS_vec[phase_ind].Comp)
-                ph_frac = out.ph_frac[phase_ind]
+                pp_ind = phase_ind - out.n_SS
+                if sys_in.casefold() == "wt":
+                    ph_comp = np.array(out.PP_vec[pp_ind].Comp_wt)
+                    ph_frac = out.ph_frac_wt[phase_ind]
+                else:
+                    ph_comp = np.array(out.PP_vec[pp_ind].Comp)
+                    ph_frac = out.ph_frac[phase_ind]
 
             if frac_amount is None:
                 frac_amount = ph_frac
 
+            if frac_amount <= 0:
+                return np.array(current_X, dtype=float)
+
             if frac_amount >= 1.0:
                 warnings.warn(f"fractionate_phase: requested frac_amount={frac_amount} >= 1.0; skipping.")
                 return np.array(current_X, dtype=float)

src/phasetools/calculators/pt_estimation.py

@@ -26,8 +26,8 @@ def _get_phase_composition(self, out, phase, components, comp_type='element'):
         # Handle special components Fe2 and Fe3 first
         if 'Fe2' in components or 'Fe3' in components:
             split = self._extract_fe_split_from_apfu(out, phase)
-            element_map['Fe2'] = split['fe2']
-            element_map['Fe3'] = split['fe3']
+            element_map['Fe2'] = split['Fe2']
+            element_map['Fe3'] = split['Fe3']
 
         for ox, stoichiometry in self._stoich_map.items():
             for el, mult in stoichiometry.items():

tests/test_fe_basis.py

@@ -0,0 +1,89 @@
+"""Mock-based tests for the garnet ``fe_basis`` option.
+
+Verifies that ``MAGEMinGarnetCalculator._extract_garnet_elements_from_oxides``
+honours the ``fe_basis`` flag:
+
+* ``'FeOt'`` (default) -- total Fe (FeO + 2*Fe2O3 APFU) placed on the
+  divalent X-site.  This is the standard community convention for garnet
+  end-members / X-site fractions.
+* ``'Fe2+'`` -- only the stoichiometrically estimated ferrous iron is
+  placed on the divalent site, excluding Fe3+.
+
+No live Julia runtime is needed -- ``get_oxide_apfu`` and the Fe2+/Fe3+
+split are mocked.
+"""
+
+import unittest
+from unittest.mock import MagicMock, patch
+
+from phasetools.calculators.garnet import MAGEMinGarnetCalculator
+
+
+class TestGarnetFeBasis(unittest.TestCase):
+    """Fe-basis flag on the garnet X-site extraction."""
+
+    def _make_calc(self, fe_basis):
+        calc = MAGEMinGarnetCalculator.__new__(MAGEMinGarnetCalculator)
+        calc.fe_basis = fe_basis  # canonical lowercase form set by __init__
+        return calc
+
+    def test_feot_uses_total_iron(self):
+        """FeOt puts FeO + 2*Fe2O3 on the X-site and normalises to 1."""
+        calc = self._make_calc('feot')
+        out = MagicMock()
+        out.ph = ['g', 'q']
+        apfu = {'MgO': 0.6, 'MnO': 0.1, 'CaO': 0.8, 'FeO': 1.2, 'Fe2O3': 0.15}
+        with patch('phasetools.calculators.garnet.get_oxide_apfu', return_value=apfu):
+            Mg, Mn, Fe, Ca = calc._extract_garnet_elements_from_oxides(out, 'mol')
+        fe_total = 1.2 + 2.0 * 0.15
+        total = 0.6 + 0.1 + 0.8 + fe_total
+        self.assertAlmostEqual(Fe, fe_total / total, places=6)
+        self.assertAlmostEqual(Mg, 0.6 / total, places=6)
+        self.assertAlmostEqual(Mn, 0.1 / total, places=6)
+        self.assertAlmostEqual(Ca, 0.8 / total, places=6)
+        self.assertAlmostEqual(Mg + Mn + Fe + Ca, 1.0, places=6)
+
+    def test_fe2_uses_split(self):
+        """Fe2+ uses only the ferrous split, excluding Fe3+."""
+        calc = self._make_calc('fe2+')
+        out = MagicMock()
+        out.ph = ['g', 'q']
+        apfu = {'MgO': 0.6, 'MnO': 0.1, 'CaO': 0.8}
+        split = {'Fe2': 1.35, 'Fe3': 0.15}
+        with patch('phasetools.calculators.garnet.get_oxide_apfu', return_value=apfu), \
+             patch.object(calc, '_extract_fe_split_from_apfu', return_value=split):
+            Mg, Mn, Fe, Ca = calc._extract_garnet_elements_from_oxides(out, 'mol')
+        total = 0.6 + 0.1 + 0.8 + 1.35
+        self.assertAlmostEqual(Fe, 1.35 / total, places=6)
+        self.assertAlmostEqual(Mg + Mn + Fe + Ca, 1.0, places=6)
+        self.assertAlmostEqual(Mn, 0.1 / total, places=6)
+
+    def test_absent_garnet_returns_zeros(self):
+        """No garnet in the assemblage -> all-zero X-site fractions."""
+        calc = self._make_calc('feot')
+        out = MagicMock()
+        out.ph = ['q', 'dio']
+        Mg, Mn, Fe, Ca = calc._extract_garnet_elements_from_oxides(out, 'mol')
+        self.assertEqual((Mg, Mn, Fe, Ca), (0.0, 0.0, 0.0, 0.0))
+
+    def test_invalid_basis_raises(self):
+        """Unsupported fe_basis values must raise ValueError at construction."""
+        with patch('phasetools.calculators.garnet.MAGEMinPTGridCalculator.__init__',
+                   return_value=None):
+            with self.assertRaises(ValueError):
+                MAGEMinGarnetCalculator(db='ig', fe_basis='Fe3')
+
+    def test_default_is_feot(self):
+        """The default fe_basis is 'FeOt' (community convention)."""
+        with patch('phasetools.calculators.garnet.MAGEMinPTGridCalculator.__init__',
+                   return_value=None):
+            calc = MAGEMinGarnetCalculator(db='ig')
+            self.assertEqual(calc.fe_basis, 'feot')
+            calc = MAGEMinGarnetCalculator(db='ig', fe_basis='Fe2+')
+            self.assertEqual(calc.fe_basis, 'fe2+')
+            calc = MAGEMinGarnetCalculator(db='ig', fe_basis='feot')
+            self.assertEqual(calc.fe_basis, 'feot')
+
+
+if __name__ == '__main__':
+    unittest.main()