Merge pull request #571 from bobmyhill/reformat_solution

reformatted solution

burnman/classes/solution.py

@@ -1,6 +1,6 @@
 # This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
 # the Earth and Planetary Sciences
-# Copyright (C) 2012 - 2022 by the BurnMan team, released under the GNU
+# Copyright (C) 2012 - 2024 by the BurnMan team, released under the GNU
 # GPL v2 or later.
 
 
@@ -46,8 +46,8 @@ class Solution(Mineral):
     :param solution_model: The SolutionModel object defining the properties
         of the solution.
     :type solution_model: :class:`burnman.SolutionModel`
-    :param molar_fractions: The molar fractions of each endmember in the solution.
-        Can be reset using the set_composition() method.
+    :param molar_fractions: The molar fractions of each endmember
+        in the solution. Can be reset using the set_composition() method.
     :type molar_fractions: numpy.array
     """
 
@@ -67,10 +67,8 @@ def __init__(self, name=None, solution_model=None, molar_fractions=None):
             self.solution_model = solution_model
 
         # Equation of state
-        for i in range(self.n_endmembers):
-            self.solution_model.endmembers[i][0].set_method(
-                self.solution_model.endmembers[i][0].params["equation_of_state"]
-            )
+        for mbr in self.solution_model.endmembers:
+            mbr[0].set_method(mbr[0].params["equation_of_state"])
 
         # Molar fractions
         if molar_fractions is not None:
@@ -85,16 +83,21 @@ def set_composition(self, molar_fractions):
         Set the composition for this solution.
         Resets cached properties.
 
-        :param molar_fractions: Molar abundance for each endmember, needs to sum to one.
+        :param molar_fractions: Molar abundance for each endmember,
+            needs to sum to one.
         :type molar_fractions: list of float
         """
         assert len(self.solution_model.endmembers) == len(molar_fractions)
 
-        if type(self.solution_model) != MechanicalSolution:
-            assert sum(molar_fractions) > 0.9999
-            assert sum(molar_fractions) < 1.0001
+        if type(self.solution_model) is not MechanicalSolution:
+            if np.abs(sum(molar_fractions) - 1.0) > 1.0e-4:
+                raise ValueError(
+                    "Molar fractions do not sum to one for "
+                    "the current instance of "
+                    f"<{self.name}>: {molar_fractions}"
+                )
 
-        if type(self.solution_model) == PolynomialSolution:
+        if type(self.solution_model) is PolynomialSolution:
             self.solution_model.set_composition(molar_fractions)
 
         self.reset()
@@ -107,12 +110,12 @@ def set_method(self, method):
         self.reset()
 
     def set_state(self, pressure, temperature):
-        if type(self.solution_model) == PolynomialSolution:
+        if type(self.solution_model) is PolynomialSolution:
             self.solution_model.set_state(pressure, temperature)
 
         Mineral.set_state(self, pressure, temperature)
-        for i in range(self.n_endmembers):
-            self.solution_model.endmembers[i][0].set_state(pressure, temperature)
+        for mbr in self.solution_model.endmembers:
+            mbr[0].set_state(pressure, temperature)
 
     @material_property
     def formula(self):
@@ -128,23 +131,22 @@ def site_occupancies(self):
         :returns: The fractional occupancies of species on each site.
         :rtype: list of OrderedDicts
         """
-        occs = np.einsum(
-            "ij, i", self.solution_model.endmember_occupancies, self.molar_fractions
-        )
+        f = self.molar_fractions
+        occs = np.einsum("ij, i", self.solution_model.endmember_occupancies, f)
         site_occs = []
         k = 0
-        for i in range(self.solution_model.n_sites):
+        for site in self.solution_model.sites:
             site_occs.append(OrderedDict())
-            for j in range(len(self.solution_model.sites[i])):
-                site_occs[-1][self.solution_model.sites[i][j]] = occs[k]
+            for species in site:
+                site_occs[-1][species] = occs[k]
                 k += 1
 
         return site_occs
 
     def site_formula(self, precision=2):
         """
-        Returns the molar chemical formula of the solution with site occupancies.
-            For example, [Mg0.4Fe0.6]2SiO4.
+        Returns the molar chemical formula of the solution with
+            site occupancies. For example, [Mg0.4Fe0.6]2SiO4.
 
         :param precision: Precision with which to print the site occupancies
         :type precision: int
@@ -225,12 +227,7 @@ def partial_gibbs(self):
         Property specific to solutions.
         """
         return (
-            np.array(
-                [
-                    self.solution_model.endmembers[i][0].gibbs
-                    for i in range(self.n_endmembers)
-                ]
-            )
+            np.array([mbr[0].gibbs for mbr in self.solution_model.endmembers])
             + self.excess_partial_gibbs
         )
 
@@ -241,12 +238,7 @@ def partial_volumes(self):
         Property specific to solutions.
         """
         return (
-            np.array(
-                [
-                    self.solution_model.endmembers[i][0].molar_volume
-                    for i in range(self.n_endmembers)
-                ]
-            )
+            np.array([mbr[0].molar_volume for mbr in self.solution_model.endmembers])
             + self.excess_partial_volumes
         )
 
@@ -257,12 +249,7 @@ def partial_entropies(self):
         Property specific to solutions.
         """
         return (
-            np.array(
-                [
-                    self.solution_model.endmembers[i][0].molar_entropy
-                    for i in range(self.n_endmembers)
-                ]
-            )
+            np.array([mbr[0].molar_entropy for mbr in self.solution_model.endmembers])
             + self.excess_partial_entropies
         )
 
@@ -312,10 +299,11 @@ def molar_gibbs(self):
         Returns molar Gibbs free energy of the solution [J/mol].
         Aliased with self.gibbs.
         """
+        f = self.molar_fractions
         return (
             sum(
                 [
-                    self.solution_model.endmembers[i][0].gibbs * self.molar_fractions[i]
+                    self.solution_model.endmembers[i][0].gibbs * f[i]
                     for i in range(self.n_endmembers)
                 ]
             )
@@ -394,10 +382,11 @@ def molar_entropy(self):
         Returns molar entropy of the solution [J/K/mol].
         Aliased with self.S.
         """
+        f = self.molar_fractions
         return (
             sum(
                 [
-                    self.solution_model.endmembers[i][0].S * self.molar_fractions[i]
+                    self.solution_model.endmembers[i][0].S * f[i]
                     for i in range(self.n_endmembers)
                 ]
             )
@@ -420,10 +409,11 @@ def molar_enthalpy(self):
         Returns molar enthalpy of the solution [J/mol].
         Aliased with self.H.
         """
+        f = self.molar_fractions
         return (
             sum(
                 [
-                    self.solution_model.endmembers[i][0].H * self.molar_fractions[i]
+                    self.solution_model.endmembers[i][0].H * f[i]
                     for i in range(self.n_endmembers)
                 ]
             )
@@ -671,7 +661,8 @@ 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.
         """
-        null_basis = np.array([v[:] for v in self.stoichiometric_matrix.nullspace()])
+        null = self.stoichiometric_matrix.nullspace()
+        null_basis = np.array([v[:] for v in null])
 
         M = null_basis[:, self.dependent_element_indices]
         assert (M.shape[0] == M.shape[1]) and (M == np.eye(M.shape[0])).all()
@@ -683,7 +674,8 @@ def endmember_formulae(self):
         """
         A list of formulae for all the endmember in the solution.
         """
-        return [mbr[0].params["formula"] for mbr in self.solution_model.endmembers]
+        mbrs = self.solution_model.endmembers
+        return [mbr[0].params["formula"] for mbr in mbrs]
 
     @cached_property
     def endmember_names(self):