Merge pull request #655 from bobmyhill/improve_solution_testing

Remove "free" from energy argument and method names, improve solution testing

Author: bobmyhill

Readme.md

@@ -1,6 +1,6 @@
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.14238360.svg)](https://doi.org/10.5281/zenodo.14238360)
 [![DOI](https://joss.theoj.org/papers/10.21105/joss.05389/status.svg)](https://doi.org/10.21105/joss.05389)
-[![coverage](https://img.shields.io/badge/coverage-%3E92%25-green")](https://github.com/geodynamics/burnman/actions/workflows/coverage.yml)
+[![coverage](https://img.shields.io/badge/coverage-%3E93%25-green")](https://github.com/geodynamics/burnman/actions/workflows/coverage.yml)
 
 # BurnMan - a Python toolkit for planetary geophysics, geochemistry and thermodynamics
 

burnman/classes/combinedmineral.py

@@ -34,7 +34,7 @@ def __init__(
         self,
         mineral_list,
         molar_amounts,
-        free_energy_adjustment=[],
+        energy_adjustment=[],
         name="User-created endmember",
     ):
         model = MechanicalSolution(endmembers=[[m, ""] for m in mineral_list])
@@ -54,9 +54,9 @@ def __init__(
             "n": sum(self.mixture.formula.values()),
         }
 
-        if free_energy_adjustment != []:
-            assert len(free_energy_adjustment) == 3
-            dE, dS, dV = free_energy_adjustment
+        if energy_adjustment != []:
+            assert len(energy_adjustment) == 3
+            dE, dS, dV = energy_adjustment
             self.property_modifiers = [
                 ["linear", {"delta_E": dE, "delta_S": dS, "delta_V": dV}]
             ]

burnman/classes/elasticsolution.py

@@ -150,7 +150,7 @@ def activities(self):
         ]
 
         gibbs_pure = [
-            self.solution_model.endmembers[i][0].method.gibbs_free_energy(
+            self.solution_model.endmembers[i][0].method.gibbs_energy(
                 self.pressure,
                 self.temperature,
                 volumes[i],

burnman/classes/elasticsolutionmodel.py

@@ -296,7 +296,7 @@ def entropy_hessian(self, volume, temperature, molar_fractions):
     def pressure_hessian(self, volume, temperature, molar_fractions):
         return np.zeros((len(molar_fractions), len(molar_fractions)))
 
-    def _configurational_entropy(self, molar_fractions):
+    def configurational_entropy(self, molar_fractions):
         site_noccupancies = np.einsum(
             "i, ij", molar_fractions, self.endmember_noccupancies
         )

burnman/classes/mineral.py

@@ -224,7 +224,7 @@ def _delta_volume(pressure, volume, temperature):
     @copy_documentation(Material.molar_gibbs)
     def molar_gibbs(self):
         return (
-            self.method.gibbs_free_energy(
+            self.method.gibbs_energy(
                 self.pressure,
                 self.temperature,
                 self._molar_volume_unmodified,

burnman/classes/solution.py

@@ -201,7 +201,7 @@ def excess_partial_gibbs(self):
         Returns excess partial molar gibbs free energy [J/mol].
         Property specific to solutions.
         """
-        return self.solution_model.excess_partial_gibbs_free_energies(
+        return self.solution_model.excess_partial_gibbs_energies(
             self.pressure, self.temperature, self.molar_fractions
         )
 
@@ -264,7 +264,7 @@ def excess_gibbs(self):
         Returns molar excess gibbs free energy [J/mol].
         Property specific to solutions.
         """
-        return self.solution_model.excess_gibbs_free_energy(
+        return self.solution_model.excess_gibbs_energy(
             self.pressure, self.temperature, self.molar_fractions
         )
 
@@ -612,16 +612,25 @@ def f(i, j):
     @cached_property
     def stoichiometric_array(self):
         """
-        An array where each element arr[i,j] corresponds
+        A numpy array where each element arr[i,j] corresponds
         to the number of atoms of element[j] in endmember[i].
         """
         return np.array(self.stoichiometric_matrix)
 
     @cached_property
     def reaction_basis(self):
         """
-        An array where each element arr[i,j] corresponds
-        to the number of moles of endmember[j] involved in reaction[i].
+        Returns a basis for the reaction space, represented as the left
+        nullspace of the stoichiometric matrix.
+
+        Each row of the returned array corresponds to a linearly independent
+        chemical reaction that conserves mass and charge.
+        The value at arr[i, j] gives the number of moles of endmember[j]
+        involved in reaction[i].
+
+        :returns: An array where each row defines a balanced chemical reaction
+        in terms of molar amounts of endmembers.
+        :rtype: numpy.ndarray
         """
         reaction_basis = np.array(
             [v[:] for v in self.stoichiometric_matrix.T.nullspace()], dtype=float
@@ -641,11 +650,25 @@ def n_reactions(self):
 
     @cached_property
     def compositional_basis(self):
-        """_summary_
+        """
+        :returns: A basis for the compositional degrees of freedom,
+        orthogonal to the reaction space.
 
-        :return: _description_
-        :rtype: _type_
+        Each row of the returned array defines an independent
+        compositional constraint in terms of moles of endmembers.
+        These constraints span the row space of the stoichiometric
+        matrix and are complementary to the reaction basis
+        (which spans the left null space).
+
+        For example, endmembers with site-species occupancies
+        [Fe][Mg] and [Mg][Fe] have reaction basis vector [1, -1].
+        A compositional basis vector orthogonal to this is: [1, 1],
+        which satisfies the requirement that the Fe-Mg ratio of the
+        system remains fixed and the total amount of the two sites
+        remains equal.
+        :rtype: numpy.ndarray
         """
+
         return complete_basis(self.reaction_basis)[self.n_reactions :]
 
     @cached_property
@@ -672,8 +695,29 @@ def dependent_element_indices(self):
     @cached_property
     def compositional_null_basis(self):
         """
-        An array N such that N.b = 0 for all bulk compositions that can
-        be produced with a linear sum of the endmembers in the solution.
+        Returns a basis for the compositional nullspace (kernel)
+        of the stoichiometric matrix.
+
+        This basis consists of vectors that represent linear constraints
+        any valid bulk composition must satisfy to be formed from the
+        endmembers. In other words, for any bulk composition vector `b`
+        that can be expressed as a combination of the endmembers,
+        multiplying it by these vectors results in zero (`N @ b = 0`).
+
+        These constraints typically reflect conserved properties
+        such as elemental mass balance or charge neutrality.
+
+        Example: Suppose the stoichiometric matrix corresponds to
+        elements C and O in species, and the nullspace basis vector
+        is `[1, -2]`. Then for any bulk composition vector
+        `b = [b_C, b_O]`, the constraint `1 * b_C - 2 * b_O = 0`
+        must hold. This means the ratio of carbon to oxygen must
+        satisfy that relationship to be achievable from the given
+        endmembers.
+
+        :return: A 2D array where each row is a basis vector of the
+        nullspace.
+        :rtype: numpy.ndarray
         """
         null = self.stoichiometric_matrix.nullspace()
         null_basis = np.array([v[:] for v in null])
@@ -687,15 +731,15 @@ def compositional_null_basis(self):
     @cached_property
     def endmember_formulae(self):
         """
-        A list of formulae for all the endmember in the solution.
+        A list of formulae for all the endmembers in the solution.
         """
         mbrs = self.solution_model.endmembers
         return [mbr[0].params["formula"] for mbr in mbrs]
 
     @cached_property
     def endmember_names(self):
         """
-        A list of names for all the endmember in the solution.
+        A list of names for all the endmembers in the solution.
         """
         return [mbr[0].name for mbr in self.solution_model.endmembers]
 

burnman/classes/solutionmodel.py

@@ -185,7 +185,7 @@ def __init__(self):
         """
         pass
 
-    def excess_gibbs_free_energy(self, pressure, temperature, molar_fractions):
+    def excess_gibbs_energy(self, pressure, temperature, molar_fractions):
         """
         Given a list of molar fractions of different phases,
         compute the excess Gibbs free energy of the solution.
@@ -205,9 +205,7 @@ def excess_gibbs_free_energy(self, pressure, temperature, molar_fractions):
         """
         return np.dot(
             np.array(molar_fractions),
-            self.excess_partial_gibbs_free_energies(
-                pressure, temperature, molar_fractions
-            ),
+            self.excess_partial_gibbs_energies(pressure, temperature, molar_fractions),
         )
 
     def excess_volume(self, pressure, temperature, molar_fractions):
@@ -271,13 +269,11 @@ def excess_enthalpy(self, pressure, temperature, molar_fractions):
         :returns: The excess enthalpy of the solution.
         :rtype: float
         """
-        return self.excess_gibbs_free_energy(
+        return self.excess_gibbs_energy(
             pressure, temperature, molar_fractions
         ) + temperature * self.excess_entropy(pressure, temperature, molar_fractions)
 
-    def excess_partial_gibbs_free_energies(
-        self, pressure, temperature, molar_fractions
-    ):
+    def excess_partial_gibbs_energies(self, pressure, temperature, molar_fractions):
         """
         Given a list of molar fractions of different phases,
         compute the excess Gibbs free energy for each endmember
@@ -375,7 +371,7 @@ def __init__(self, endmembers):
         self.n_endmembers = len(endmembers)
         self.site_formulae = [e[1] for e in endmembers]
 
-    def excess_gibbs_free_energy(self, pressure, temperature, molar_fractions):
+    def excess_gibbs_energy(self, pressure, temperature, molar_fractions):
         return 0.0
 
     def excess_volume(self, pressure, temperature, molar_fractions):
@@ -387,9 +383,7 @@ def excess_entropy(self, pressure, temperature, molar_fractions):
     def excess_enthalpy(self, pressure, temperature, molar_fractions):
         return 0.0
 
-    def excess_partial_gibbs_free_energies(
-        self, pressure, temperature, molar_fractions
-    ):
+    def excess_partial_gibbs_energies(self, pressure, temperature, molar_fractions):
         return np.zeros_like(molar_fractions)
 
     def excess_partial_volumes(self, pressure, temperature, molar_fractions):
@@ -440,9 +434,7 @@ def _calculate_endmember_configurational_entropies(self):
         )
         self.endmember_configurational_entropies = S_conf
 
-    def excess_partial_gibbs_free_energies(
-        self, pressure, temperature, molar_fractions
-    ):
+    def excess_partial_gibbs_energies(self, pressure, temperature, molar_fractions):
         return self._ideal_excess_partial_gibbs(temperature, molar_fractions)
 
     def excess_partial_entropies(self, pressure, temperature, molar_fractions):
@@ -462,7 +454,7 @@ def entropy_hessian(self, pressure, temperature, molar_fractions):
     def volume_hessian(self, pressure, temperature, molar_fractions):
         return np.zeros((len(molar_fractions), len(molar_fractions)))
 
-    def _configurational_entropy(self, molar_fractions):
+    def configurational_entropy(self, molar_fractions):
         site_noccupancies = np.einsum(
             "i, ij", molar_fractions, self.endmember_noccupancies
         )
@@ -656,9 +648,7 @@ def _non_ideal_excess_partial_gibbs(self, pressure, temperature, molar_fractions
         Vint = self._non_ideal_interactions(self.Wv, molar_fractions)
         return Eint - temperature * Sint + pressure * Vint
 
-    def excess_partial_gibbs_free_energies(
-        self, pressure, temperature, molar_fractions
-    ):
+    def excess_partial_gibbs_energies(self, pressure, temperature, molar_fractions):
         ideal_gibbs = IdealSolution._ideal_excess_partial_gibbs(
             self, temperature, molar_fractions
         )
@@ -888,9 +878,7 @@ def _non_ideal_excess_partial_gibbs(self, pressure, temperature, molar_fractions
         Eint, Sint, Vint = self._non_ideal_interactions(molar_fractions)
         return Eint - temperature * Sint + pressure * Vint
 
-    def excess_partial_gibbs_free_energies(
-        self, pressure, temperature, molar_fractions
-    ):
+    def excess_partial_gibbs_energies(self, pressure, temperature, molar_fractions):
         ideal_gibbs = IdealSolution._ideal_excess_partial_gibbs(
             self, temperature, molar_fractions
         )
@@ -1031,9 +1019,7 @@ def volume_hess(pressure, temperature, molar_amounts):
 
         self.volume_hessian = volume_hess
 
-    def excess_partial_gibbs_free_energies(
-        self, pressure, temperature, molar_fractions
-    ):
+    def excess_partial_gibbs_energies(self, pressure, temperature, molar_fractions):
         ideal_gibbs = IdealSolution._ideal_excess_partial_gibbs(
             self, temperature, molar_fractions
         )
@@ -1251,8 +1237,11 @@ def _make_prefactors_and_exponents(self, a):
 
     def _make_interaction_arrays(self, Ws, n):
         """
-        A hidden convenience function that splits each
-        excess term into _summary_
+        A hidden convenience function that splits the user-defined
+        excess interaction terms into an array of interaction values
+        and a list of Interaction dataclass objects that facilitate
+        fast computation of first and second compositional derivatives
+        of thermodynamic properties.
 
         :param Ws: List of interactions in form input by the user.
         :type Ws: list of lists
@@ -1488,9 +1477,7 @@ def _non_ideal_excess_partial_gibbs(self, pressure, temperature, molar_fractions
         gibbs += np.einsum("i, ij->j", mbr_gibbs, mbr_gradients)
         return gibbs
 
-    def excess_partial_gibbs_free_energies(
-        self, pressure, temperature, molar_fractions
-    ):
+    def excess_partial_gibbs_energies(self, pressure, temperature, molar_fractions):
         ideal_gibbs = IdealSolution._ideal_excess_partial_gibbs(
             self, temperature, molar_fractions
         )

burnman/eos/aa.py

@@ -344,7 +344,7 @@ def molar_heat_capacity_p(self, pressure, temperature, volume, params):
 
         return C_p
 
-    def gibbs_free_energy(self, pressure, temperature, volume, params):
+    def gibbs_energy(self, pressure, temperature, volume, params):
         """
         Returns the Gibbs free energy at the pressure and temperature of the mineral [J/mol]
         F + PV

burnman/eos/birch_murnaghan.py

@@ -301,7 +301,7 @@ def _molar_internal_energy(self, pressure, temperature, volume, params):
 
         return -intPdV + params["E_0"]
 
-    def gibbs_free_energy(self, pressure, temperature, volume, params):
+    def gibbs_energy(self, pressure, temperature, volume, params):
         """
         Returns the Gibbs free energy :math:`\\mathcal{G}`
         of the mineral. :math:`[J/mol]`

burnman/eos/brosh_calphad.py

@@ -397,7 +397,7 @@ def _C_T(self, temperature, params):
 
         return C_T
 
-    def gibbs_free_energy(self, pressure, temperature, volume, params):
+    def gibbs_energy(self, pressure, temperature, volume, params):
         """
         Returns the Gibbs free energy of the mineral. :math:`[J/mol]`
         """