diff --git a/burnman/classes/anisotropicsolution.py b/burnman/classes/anisotropicsolution.py
index 480cedcdc..99cdbd499 100644
--- a/burnman/classes/anisotropicsolution.py
+++ b/burnman/classes/anisotropicsolution.py
@@ -178,8 +178,10 @@ def set_composition(self, molar_fractions):
         )
         # Calculate all other required properties
         Solution.set_composition(self, molar_fractions)
-        if self._scalar_solution.pressure is not None:
+        try:  # if pressure and temperature are already set, compute base properties
             self.compute_base_properties()
+        except AttributeError:
+            pass
 
     @cached_property
     def ones(self):
@@ -349,9 +351,12 @@ def set_state(self, pressure, temperature, relaxed=True):
         """
         self.unrelaxed.set_state(pressure, temperature)
 
-        if hasattr(self.unrelaxed, "molar_fractions") and relaxed:
-            self._relax_at_PTX()
-            AnisotropicSolution.set_state(self, pressure, temperature)
+        if relaxed:
+            try:  # relax the structure if composition has been set
+                self._relax_at_PTX()
+                AnisotropicSolution.set_state(self, pressure, temperature)
+            except AttributeError:
+                pass
 
     def set_composition(self, molar_fractions, q_initial=None, relaxed=True):
         """
@@ -383,11 +388,14 @@ def set_composition(self, molar_fractions, q_initial=None, relaxed=True):
 
         self.unrelaxed.set_composition(n)
 
-        if self.unrelaxed.pressure is not None and relaxed:
-            AnisotropicSolution.set_composition(self, n)
-            self._relax_at_PTX()
-            self.unrelaxed._scalar_solution.set_composition(self.molar_fractions)
-            self.unrelaxed.set_composition(self.molar_fractions)
+        if relaxed:
+            try:  # relax the structure if state has been set
+                AnisotropicSolution.set_composition(self, n)
+                self._relax_at_PTX()
+                self.unrelaxed._scalar_solution.set_composition(self.molar_fractions)
+                self.unrelaxed.set_composition(self.molar_fractions)
+            except AttributeError:
+                pass
 
     def _relax_at_PTX(self):
         """
diff --git a/burnman/classes/composite.py b/burnman/classes/composite.py
index fc83ed86e..3024fbad3 100644
--- a/burnman/classes/composite.py
+++ b/burnman/classes/composite.py
@@ -214,7 +214,7 @@ def phase_volumes(self):
         """
         volumes = np.array(
             [
-                phase.molar_volume * molar_fraction
+                phase.volume * molar_fraction * self.number_of_moles
                 for (phase, molar_fraction) in zip(self.phases, self.molar_fractions)
             ]
         )
@@ -227,7 +227,7 @@ def phase_masses(self):
         """
         masses = np.array(
             [
-                phase.molar_mass * molar_fraction
+                phase.mass * molar_fraction * self.number_of_moles
                 for (phase, molar_fraction) in zip(self.phases, self.molar_fractions)
             ]
         )
@@ -236,9 +236,10 @@ def phase_masses(self):
     @material_property
     def formula(self):
         """
-        Returns molar chemical formula of the composite
+        Returns extensive chemical formula of the composite
         """
-        return sum_formulae([ph.formula for ph in self.phases], self.molar_fractions)
+        molar_amounts = self.molar_fractions * self.number_of_moles
+        return sum_formulae([ph.formula for ph in self.phases], molar_amounts)
 
     @material_property
     def molar_internal_energy(self):
@@ -282,14 +283,14 @@ def molar_volume(self):
         Returns molar volume of the composite [m^3/mol]
         Aliased with self.V
         """
-        return np.sum(self.phase_volumes)
+        return np.sum(self.phase_volumes) / self.number_of_moles
 
     @material_property
     def molar_mass(self):
         """
         Returns molar mass of the composite [kg/mol]
         """
-        return sum(self.phase_masses)
+        return sum(self.phase_masses) / self.number_of_moles
 
     @material_property
     def density(self):
@@ -438,7 +439,7 @@ def effective_isentropic_bulk_modulus(self):
         Returns the effective isotropic, isentropic bulk modulus of the composite [Pa]
         Aliased with self.K_eff
         """
-        V_frac = self.phase_volumes
+        V_frac = self.phase_volumes / sum(self.phase_volumes)
         K_ph = np.array([phase.isentropic_bulk_modulus_reuss for phase in self.phases])
         G_ph = np.array([phase.shear_modulus for phase in self.phases])
 
@@ -450,7 +451,7 @@ def effective_shear_modulus(self):
         Returns the effective isotropic shear modulus of the mineral [Pa]
         Aliased with self.G_eff
         """
-        V_frac = self.phase_volumes
+        V_frac = self.phase_volumes / sum(self.phase_volumes)
         K_ph = np.array([phase.isentropic_bulk_modulus_reuss for phase in self.phases])
         G_ph = np.array([phase.shear_modulus for phase in self.phases])
 
diff --git a/burnman/classes/elasticsolution.py b/burnman/classes/elasticsolution.py
index ebf9fbfd2..553b14de0 100644
--- a/burnman/classes/elasticsolution.py
+++ b/burnman/classes/elasticsolution.py
@@ -111,8 +111,10 @@ def set_composition(self, molar_fractions):
         self.reset()
         self.molar_fractions = np.array(molar_fractions)
 
-        if self.temperature is not None:
+        try:  # set the volume if pressure and temperature are already set
             _ = self.molar_volume
+        except AttributeError:
+            pass
 
     def set_method(self, method):
         for i in range(self.n_endmembers):
diff --git a/burnman/classes/material.py b/burnman/classes/material.py
index a9b158740..aa19cb9ff 100644
--- a/burnman/classes/material.py
+++ b/burnman/classes/material.py
@@ -55,8 +55,7 @@ class Material(object):
     """
 
     def __init__(self):
-        self._pressure = None
-        self._temperature = None
+        self._number_of_moles = 1.0
         if not hasattr(self, "name"):
             # if a derived class decides to set .name before calling this
             # constructor (I am looking at you, SLB_2011.py!), do not
@@ -130,13 +129,7 @@ def set_state(self, pressure, temperature):
         :param temperature: The desired temperature in [K].
         :type temperature: float
         """
-        if not hasattr(self, "_pressure"):
-            raise Exception(
-                "Material.set_state() could not find class member _pressure. "
-                "Did you forget to call Material.__init__(self) in __init___?"
-            )
         self.reset()
-
         self._pressure = pressure
         self._temperature = temperature
 
@@ -182,7 +175,7 @@ def _delta_volume(pressure, volume, temperature):
             # and brentq will try a larger pressure.
             try:
                 self.set_state(pressure, temperature)
-                return volume - self.molar_volume
+                return volume - self.volume
             except Exception:
                 return -np.inf
 
@@ -201,6 +194,28 @@ def _delta_volume(pressure, volume, temperature):
         pressure = brentq(_delta_volume, sol[0], sol[1], args=args)
         self.set_state(pressure, temperature)
 
+    @property
+    def number_of_moles(self):
+        """
+        Returns the number of moles in this material.
+
+        :returns: Number of moles.
+        :rtype: float
+        """
+        return self._number_of_moles
+
+    @number_of_moles.setter
+    def number_of_moles(self, value):
+        """
+        Sets the number of moles in this material.
+
+        :param value: Number of moles.
+        :type value: float
+        """
+        if value < 0:
+            raise ValueError("The number of moles cannot be negative")
+        self._number_of_moles = value
+
     def reset(self):
         """
         Resets all cached material properties.
@@ -250,8 +265,13 @@ def evaluate(self, vars_list, pressures, temperatures, molar_fractions=None):
             have different shapes.
         :rtype: list or :class:`numpy.array`, n-dimensional
         """
-        old_pressure = self.pressure
-        old_temperature = self.temperature
+        try:
+            old_pressure = self.pressure
+            old_temperature = self.temperature
+        except AttributeError:
+            old_pressure = None
+            old_temperature = None
+
         try:
             old_molar_fractions = self.molar_fractions
         except AttributeError:
@@ -286,9 +306,11 @@ def evaluate(self, vars_list, pressures, temperatures, molar_fractions=None):
             for j in range(len(vars_list)):
                 output[j][i] = getattr(self, vars_list[j])
         if old_pressure is None or old_temperature is None:
-            # do not set_state if old values were None. Just reset to None
-            # manually
-            self._pressure = self._temperature = None
+            # delete the pressure and temperature attributes
+            # if they were not set before
+            del self._pressure
+            del self._temperature
+
             self.reset()
         else:
             self.set_state(old_pressure, old_temperature)
@@ -328,8 +350,13 @@ def evaluate_with_volumes(
             have different shapes.
         :rtype: list or :class:`numpy.array`, n-dimensional
         """
-        old_pressure = self.pressure
-        old_temperature = self.temperature
+        try:
+            old_pressure = self.pressure
+            old_temperature = self.temperature
+        except AttributeError:
+            old_pressure = None
+            old_temperature = None
+
         try:
             old_molar_fractions = self.molar_fractions
         except AttributeError:
@@ -363,10 +390,10 @@ def evaluate_with_volumes(
             for j in range(len(vars_list)):
                 output[j][i] = getattr(self, vars_list[j])
         if old_pressure is None or old_temperature is None:
-            # do not set_state if old values were None. Just reset to None
-            # manually
-            self._pressure = self._temperature = None
-            self.reset()
+            # delete the pressure and temperature attributes
+            # if they were not set before
+            del self._pressure
+            del self._temperature
         else:
             self.set_state(old_pressure, old_temperature)
 
@@ -392,7 +419,12 @@ def pressure(self):
         :returns: Pressure in [Pa].
         :rtype: float
         """
-        return self._pressure
+        try:
+            return self._pressure
+        except AttributeError:
+            raise AttributeError(
+                "Material.pressure: pressure is not set. Did you forget to call set_state()?"
+            )
 
     @property
     def temperature(self):
@@ -405,7 +437,12 @@ def temperature(self):
         :returns: Temperature in [K].
         :rtype: float
         """
-        return self._temperature
+        try:
+            return self._temperature
+        except AttributeError:
+            raise AttributeError(
+                "Material.temperature: temperature is not set. Did you forget to call set_state()?"
+            )
 
     @material_property
     def molar_internal_energy(self):
@@ -462,6 +499,16 @@ def molar_mass(self):
         """
         raise NotImplementedError("need to implement molar_mass() in derived class!")
 
+    @material_property
+    def mass(self):
+        """
+        Returns the mass of this material.
+
+        :returns: Mass in [kg].
+        :rtype: float
+        """
+        return self.molar_mass * self.number_of_moles
+
     @material_property
     def molar_volume(self):
         """
@@ -739,7 +786,47 @@ def isentropic_thermal_gradient(self):
             "need to implement isentropic_thermal_gradient() in derived class!"
         )
 
-    #
+    # Extensive properties
+    @property
+    def internal_energy(self):
+        """Extensive internal energy of the material (J)."""
+        return self.molar_internal_energy * self.number_of_moles
+
+    @property
+    def helmholtz(self):
+        """Extensive Helmholtz free energy of the material (J)."""
+        return self.molar_helmholtz * self.number_of_moles
+
+    @property
+    def gibbs(self):
+        """Extensive Gibbs free energy of the material (J)."""
+        return self.molar_gibbs * self.number_of_moles
+
+    @property
+    def enthalpy(self):
+        """Extensive enthalpy of the material (J)."""
+        return self.molar_enthalpy * self.number_of_moles
+
+    @property
+    def volume(self):
+        """Extensive volume of the material (m^3)."""
+        return self.molar_volume * self.number_of_moles
+
+    @property
+    def entropy(self):
+        """Extensive entropy of the material (J/K)."""
+        return self.molar_entropy * self.number_of_moles
+
+    @property
+    def heat_capacity_v(self):
+        """Extensive heat capacity at constant volume of the material (J/K)."""
+        return self.molar_heat_capacity_v * self.number_of_moles
+
+    @property
+    def heat_capacity_p(self):
+        """Extensive heat capacity at constant pressure of the material (J/K)."""
+        return self.molar_heat_capacity_p * self.number_of_moles
+
     # Aliased properties
     @property
     def P(self):
@@ -751,25 +838,10 @@ def T(self):
         """Alias for :func:`~burnman.Material.temperature`"""
         return self.temperature
 
-    @property
-    def energy(self):
-        """Alias for :func:`~burnman.Material.molar_internal_energy`"""
-        return self.molar_internal_energy
-
-    @property
-    def helmholtz(self):
-        """Alias for :func:`~burnman.Material.molar_helmholtz`"""
-        return self.molar_helmholtz
-
-    @property
-    def gibbs(self):
-        """Alias for :func:`~burnman.Material.molar_gibbs`"""
-        return self.molar_gibbs
-
     @property
     def V(self):
-        """Alias for :func:`~burnman.Material.molar_volume`"""
-        return self.molar_volume
+        """Alias for :func:`~burnman.Material.volume`"""
+        return self.volume
 
     @property
     def rho(self):
@@ -778,13 +850,13 @@ def rho(self):
 
     @property
     def S(self):
-        """Alias for :func:`~burnman.Material.molar_entropy`"""
-        return self.molar_entropy
+        """Alias for :func:`~burnman.Material.entropy`"""
+        return self.entropy
 
     @property
     def H(self):
-        """Alias for :func:`~burnman.Material.molar_enthalpy`"""
-        return self.molar_enthalpy
+        """Alias for :func:`~burnman.Material.enthalpy`"""
+        return self.enthalpy
 
     @property
     def K_T(self):
@@ -848,10 +920,10 @@ def alpha(self):
 
     @property
     def C_v(self):
-        """Alias for :func:`~burnman.Material.molar_heat_capacity_v`"""
-        return self.molar_heat_capacity_v
+        """Alias for :func:`~burnman.Material.heat_capacity_v`"""
+        return self.heat_capacity_v
 
     @property
     def C_p(self):
-        """Alias for :func:`~burnman.Material.molar_heat_capacity_p`"""
-        return self.molar_heat_capacity_p
+        """Alias for :func:`~burnman.Material.heat_capacity_p`"""
+        return self.heat_capacity_p
diff --git a/burnman/classes/relaxedcomposite.py b/burnman/classes/relaxedcomposite.py
index e4aaae4f0..289b312fb 100644
--- a/burnman/classes/relaxedcomposite.py
+++ b/burnman/classes/relaxedcomposite.py
@@ -79,6 +79,8 @@ def __init__(self, composite, relaxation_vectors):
             name=composite.name,
         )
 
+        self.number_of_moles = composite.number_of_moles
+
     def set_state(self, pressure, temperature, bulk_composition=None):
         """
         Sets the state of the composite. Also relaxes the
@@ -97,27 +99,19 @@ def set_state(self, pressure, temperature, bulk_composition=None):
         :type bulk_composition: dict
         """
         if bulk_composition is None:
-            try:
-                bulk_composition = Counter(
-                    {
-                        el: count * self.unrelaxed.n_moles
-                        for el, count in self.unrelaxed.formula.items()
-                    }
-                )
-            except AttributeError:
-                bulk_composition = self.unrelaxed.formula
+            bulk_composition = self.unrelaxed.formula
 
         self.unrelaxed.set_state(pressure, temperature)
 
         sol, _ = equilibrate(
             bulk_composition, self.unrelaxed, [["P", pressure], ["T", temperature]]
         )
-        assert sol.success
+        assert sol.success, "Equilibration failed in RelaxedComposite.set_state()"
 
         # Set the relaxed composite to the equilibrated state
         Composite.set_state(self, pressure, temperature)
         Composite.set_fractions(self, self.unrelaxed.molar_fractions)
-        self.n_moles = self.unrelaxed.n_moles
+        self.number_of_moles = self.unrelaxed.number_of_moles
 
     @material_property
     def _d2Gdqdq_fixed_PT(self):
@@ -126,7 +120,7 @@ def _d2Gdqdq_fixed_PT(self):
         calculated at constant pressure and temperature.
         """
         phases = self.unrelaxed.phases
-        molar_amounts = self.unrelaxed.molar_fractions * self.unrelaxed.n_moles
+        molar_amounts = self.unrelaxed.molar_fractions * self.unrelaxed.number_of_moles
         hessians = [
             (
                 phase.gibbs_hessian / molar_amounts[i]
@@ -176,9 +170,9 @@ def _d2Gdzdz_Q(self):
         """
         beta_TR = self.unrelaxed.isothermal_compressibility_reuss
         alpha = self.unrelaxed.thermal_expansivity
-        c_p = self.unrelaxed.molar_heat_capacity_p * self.n_moles
-        V = self.molar_volume * self.n_moles
-        T = self.temperature
+        c_p = self.unrelaxed.heat_capacity_p
+        V = self.unrelaxed.volume
+        T = self.unrelaxed.temperature
         return np.array([[-V * beta_TR, V * alpha], [V * alpha, -c_p / T]])
 
     @material_property
@@ -199,17 +193,17 @@ def _d2Gdzdz(self):
     @material_property
     @copy_documentation(Composite.isothermal_compressibility_reuss)
     def isothermal_compressibility_reuss(self):
-        return -self._d2Gdzdz[0, 0] / self.V / self.n_moles
+        return -self._d2Gdzdz[0, 0] / self.volume
 
     @material_property
     @copy_documentation(Composite.thermal_expansivity)
     def thermal_expansivity(self):
-        return self._d2Gdzdz[0, 1] / self.V / self.n_moles
+        return self._d2Gdzdz[0, 1] / self.volume
 
     @material_property
     @copy_documentation(Composite.molar_heat_capacity_p)
     def molar_heat_capacity_p(self):
-        return -self._d2Gdzdz[1, 1] * self.T / self.n_moles
+        return -self._d2Gdzdz[1, 1] * self.T / self.number_of_moles
 
     @material_property
     @copy_documentation(Composite.isothermal_bulk_modulus_reuss)
diff --git a/burnman/classes/solution.py b/burnman/classes/solution.py
index 4dea25886..14cd025cf 100644
--- a/burnman/classes/solution.py
+++ b/burnman/classes/solution.py
@@ -865,10 +865,13 @@ def set_state(self, pressure, temperature, relaxed=True):
         :type relaxed: bool, optional
         """
         self.unrelaxed.set_state(pressure, temperature)
+        Solution.set_state(self, pressure, temperature)
 
-        if hasattr(self.unrelaxed, "molar_fractions") and relaxed:
-            self._relax_at_PTX()
-            Solution.set_state(self, pressure, temperature)
+        if relaxed:
+            try:  # relax if composition has already been set
+                self._relax_at_PTX()
+            except AttributeError:
+                pass
 
     def set_composition(self, molar_fractions, q_initial=None, relaxed=True):
         """
@@ -900,10 +903,12 @@ def set_composition(self, molar_fractions, q_initial=None, relaxed=True):
 
         self.unrelaxed.set_composition(n)
 
-        if self.unrelaxed.pressure is not None and relaxed:
-            Solution.set_composition(self, n)
-            self._relax_at_PTX()
-            self.unrelaxed.set_composition(self.molar_fractions)
+        if relaxed:
+            try:  # relax if state has already been set
+                Solution.set_composition(self, n)
+                self._relax_at_PTX()
+            except AttributeError:
+                pass
 
     def _relax_at_PTX(self):
         """
@@ -926,6 +931,7 @@ def G_func(dq):
         n = n0 + np.einsum("ij, j", self.dndq, sol.x)
         self.unrelaxed.set_composition(n)
         Solution.set_composition(self, n)
+        self.unrelaxed.set_composition(self.molar_fractions)
 
     def set_state_with_volume(
         self, volume, temperature, pressure_guesses=[0.0e9, 10.0e9]
diff --git a/burnman/tools/eos.py b/burnman/tools/eos.py
index 994de1a41..73a475899 100644
--- a/burnman/tools/eos.py
+++ b/burnman/tools/eos.py
@@ -59,35 +59,29 @@ def equilibration_function(mineral):
     dT = 1.0
     dP = 1000.0
 
-    def n_moles(m):
-        try:
-            return m.n_moles
-        except AttributeError:
-            return 1.0
-
     m.set_state(P, T)
     equilibration_function(m)
-    G0 = m.gibbs * n_moles(m)
-    S0 = m.S * n_moles(m)
-    V0 = m.V * n_moles(m)
+    G0 = m.gibbs
+    S0 = m.S
+    V0 = m.V
 
     expr = ["G = F + PV", "G = H - TS", "G = E - TS + PV"]
     eq = [
         [m.gibbs, (m.helmholtz + P * m.V)],
         [m.gibbs, (m.H - T * m.S)],
-        [m.gibbs, (m.molar_internal_energy - T * m.S + P * m.V)],
+        [m.gibbs, (m.internal_energy - T * m.S + P * m.V)],
     ]
 
     m.set_state(P, T + dT)
     equilibration_function(m)
-    G1 = m.gibbs * n_moles(m)
-    S1 = m.S * n_moles(m)
-    V1 = m.V * n_moles(m)
+    G1 = m.gibbs
+    S1 = m.S
+    V1 = m.V
 
     m.set_state(P + dP, T)
     equilibration_function(m)
-    G2 = m.gibbs * n_moles(m)
-    V2 = m.V * n_moles(m)
+    G2 = m.gibbs
+    V2 = m.V
 
     # T derivatives
     m.set_state(P, T + 0.5 * dT)
@@ -95,9 +89,12 @@ def n_moles(m):
     expr.extend(["S = -dG/dT", "alpha = 1/V dV/dT", "C_p = T dS/dT"])
     eq.extend(
         [
-            [m.S * n_moles(m), -(G1 - G0) / dT],
+            [m.S, -(G1 - G0) / dT],
             [m.alpha, (V1 - V0) / dT / (V1 + V0) * 2.0],
-            [m.molar_heat_capacity_p * n_moles(m), (T + 0.5 * dT) * (S1 - S0) / dT],
+            [
+                m.heat_capacity_p,
+                (T + 0.5 * dT) * (S1 - S0) / dT,
+            ],
         ]
     )
 
@@ -107,7 +104,7 @@ def n_moles(m):
     expr.extend(["V = dG/dP", "K_T = -V dP/dV"])
     eq.extend(
         [
-            [m.V * n_moles(m), (G2 - G0) / dP],
+            [m.V, (G2 - G0) / dP],
             [m.isothermal_bulk_modulus_reuss, -0.5 * (V2 + V0) * dP / (V2 - V0)],
         ]
     )
@@ -118,22 +115,17 @@ def n_moles(m):
     eq.extend(
         [
             [
-                m.molar_heat_capacity_v,
-                m.molar_heat_capacity_p
+                m.heat_capacity_v,
+                m.heat_capacity_p
                 - m.alpha * m.alpha * m.isothermal_bulk_modulus_reuss * m.V * T,
             ],
             [
                 m.isentropic_bulk_modulus_reuss,
-                m.isothermal_bulk_modulus_reuss
-                * m.molar_heat_capacity_p
-                / m.molar_heat_capacity_v,
+                m.isothermal_bulk_modulus_reuss * m.heat_capacity_p / m.heat_capacity_v,
             ],
             [
                 m.gr,
-                m.alpha
-                * m.isothermal_bulk_modulus_reuss
-                * m.V
-                / m.molar_heat_capacity_v,
+                m.alpha * m.isothermal_bulk_modulus_reuss * m.V / m.heat_capacity_v,
             ],
         ]
     )
@@ -248,7 +240,7 @@ def equilibration_function(mineral):
     eq = [
         [m.gibbs, (m.helmholtz + P * m.V)],
         [m.gibbs, (m.H - T * m.S)],
-        [m.gibbs, (m.molar_internal_energy - T * m.S + P * m.V)],
+        [m.gibbs, (m.internal_energy - T * m.S + P * m.V)],
     ]
 
     m.set_state(P, T + dT)
@@ -270,7 +262,7 @@ def equilibration_function(mineral):
         [
             [m.S, -(G1 - G0) / dT],
             [m.alpha, (V1 - V0) / dT / m.V],
-            [m.molar_heat_capacity_p, (T + 0.5 * dT) * (S1 - S0) / dT],
+            [m.heat_capacity_p, (T + 0.5 * dT) * (S1 - S0) / dT],
         ]
     )
 
@@ -289,15 +281,13 @@ def equilibration_function(mineral):
     eq.extend(
         [
             [
-                m.molar_heat_capacity_v,
-                m.molar_heat_capacity_p
+                m.heat_capacity_v,
+                m.heat_capacity_p
                 - m.alpha * m.alpha * m.isothermal_bulk_modulus_reuss * m.V * T,
             ],
             [
                 m.isentropic_bulk_modulus_reuss,
-                m.isothermal_bulk_modulus_reuss
-                * m.molar_heat_capacity_p
-                / m.molar_heat_capacity_v,
+                m.isothermal_bulk_modulus_reuss * m.heat_capacity_p / m.heat_capacity_v,
             ],
         ]
     )
diff --git a/burnman/tools/equilibration.py b/burnman/tools/equilibration.py
index bdacd2d5b..072b30973 100644
--- a/burnman/tools/equilibration.py
+++ b/burnman/tools/equilibration.py
@@ -101,7 +101,7 @@ def get_parameters(assemblage, n_free_compositional_vectors=0):
     :rtype: numpy.array
     """
     params = np.zeros(assemblage.n_endmembers + 2 + n_free_compositional_vectors)
-    n_moles_phase = assemblage.n_moles * np.array(assemblage.molar_fractions)
+    n_moles_phase = assemblage.number_of_moles * np.array(assemblage.molar_fractions)
 
     try:
         params[:2] = [
@@ -133,7 +133,7 @@ def get_endmember_amounts(assemblage):
     :returns: The current amounts of all the endmembers.
     :rtype: numpy.array
     """
-    phase_amounts = assemblage.n_moles * assemblage.molar_fractions
+    phase_amounts = assemblage.number_of_moles * assemblage.molar_fractions
     amounts = np.empty(assemblage.n_endmembers)
     j = 0
     for i, ph in enumerate(assemblage.phases):
@@ -177,8 +177,8 @@ def set_compositions_and_state_from_parameters(assemblage, parameters):
 
     assert np.all(phase_amounts > -1.0e-8)
     phase_amounts = np.abs(phase_amounts)
-    assemblage.n_moles = sum(phase_amounts)
-    assemblage.set_fractions(phase_amounts / assemblage.n_moles)
+    assemblage.number_of_moles = sum(phase_amounts)
+    assemblage.set_fractions(phase_amounts / assemblage.number_of_moles)
     return None
 
 
@@ -268,9 +268,9 @@ def F(
         elif type_c == "T":
             eqns[i] = (assemblage.temperature - eq_c) / T_scaling
         elif type_c == "S":
-            eqns[i] = (assemblage.molar_entropy * assemblage.n_moles - eq_c) / S_scaling
+            eqns[i] = (assemblage.entropy - eq_c) / S_scaling
         elif type_c == "V":
-            eqns[i] = (assemblage.molar_volume * assemblage.n_moles - eq_c) / V_scaling
+            eqns[i] = (assemblage.volume - eq_c) / V_scaling
         elif type_c == "X":
             eqns[i] = (np.dot(eq_c[0], x) - eq_c[1]) / X_scaling  # i.e. Ax = b
         else:
@@ -340,21 +340,15 @@ def jacobian(x, assemblage, equality_constraints, reduced_free_composition_vecto
         elif type_c == "S":  # dS/dx
             # dS/dP = -aV, dS/dT = Cp/T
             jacobian[ic, 0:2] = [
-                -assemblage.n_moles
-                * assemblage.alpha
-                * assemblage.molar_volume
-                / S_scaling,
-                assemblage.n_moles
-                * assemblage.molar_heat_capacity_p
-                / assemblage.temperature
-                / S_scaling,
+                -assemblage.alpha * assemblage.volume / S_scaling,
+                assemblage.heat_capacity_p / assemblage.temperature / S_scaling,
             ]
             j = 2
             for k, n in enumerate(assemblage.endmembers_per_phase):
                 jacobian[ic, j] = assemblage.phases[k].molar_entropy / S_scaling
                 if n > 1:  # for solutions with >1 endmember
                     jacobian[ic, j + 1 : j + n] = (
-                        assemblage.n_moles
+                        assemblage.number_of_moles
                         * assemblage.molar_fractions[k]
                         * (
                             assemblage.phases[k].partial_entropies[1:]
@@ -366,21 +360,17 @@ def jacobian(x, assemblage, equality_constraints, reduced_free_composition_vecto
         elif type_c == "V":  # dV/dx
             # dV/dP = -V/K_T, dV/dT = aV
             jacobian[ic, 0:2] = [
-                -assemblage.n_moles
-                * assemblage.molar_volume
+                -assemblage.volume
                 / assemblage.isothermal_bulk_modulus_reuss
                 / V_scaling,
-                assemblage.n_moles
-                * assemblage.molar_volume
-                * assemblage.alpha
-                / V_scaling,
+                assemblage.volume * assemblage.alpha / V_scaling,
             ]
             j = 2
             for k, n in enumerate(assemblage.endmembers_per_phase):
                 jacobian[ic, j] = assemblage.phases[k].molar_volume / V_scaling
                 if n > 1:  # for solutions with >1 stable endmember
                     jacobian[ic, j + 1 : j + n] = (
-                        assemblage.n_moles
+                        assemblage.number_of_moles
                         * assemblage.molar_fractions[k]
                         * (
                             assemblage.phases[k].partial_volumes[1:]
@@ -422,7 +412,7 @@ def jacobian(x, assemblage, equality_constraints, reduced_free_composition_vecto
 
     # Finally, let's build the compositional Hessian d2G/dfidfj = dmui/dfj
     # where fj is the fraction of endmember j in a phase
-    phase_amounts = np.array(assemblage.molar_fractions) * assemblage.n_moles
+    phase_amounts = np.array(assemblage.molar_fractions) * assemblage.number_of_moles
     comp_hessian = np.zeros((assemblage.n_endmembers, assemblage.n_endmembers))
     dfi_dxj = np.zeros((assemblage.n_endmembers, assemblage.n_endmembers))
     dpi_dxj = np.zeros((assemblage.n_endmembers, assemblage.n_endmembers))
@@ -1000,7 +990,7 @@ def equilibrate(
         assemblage.set_fractions([f for i in range(n_phases)])
 
     n_atoms = sum(composition.values())
-    assemblage.n_moles = n_atoms / sum(assemblage.formula.values())
+    assemblage.number_of_moles = n_atoms / sum(assemblage.formula.values())
 
     n_equality_constraints = len(equality_constraints)
     n_free_compositional_vectors = len(free_compositional_vectors)
@@ -1036,8 +1026,10 @@ def equilibrate(
     # Set up solves
     nc = [len(eq_constraint_list) for eq_constraint_list in eq_constraint_lists]
 
-    # Find the initial state (could be none here)
-    initial_state = [assemblage.pressure, assemblage.temperature]
+    try:
+        initial_state = [assemblage.pressure, assemblage.temperature]
+    except AttributeError:
+        initial_state = [None, None]
 
     # Reset initial state if equality constraints
     # are related to pressure or temperature
diff --git a/examples/example_equilibrate.py b/examples/example_equilibrate.py
index d1697b3e6..c6055f4e6 100644
--- a/examples/example_equilibrate.py
+++ b/examples/example_equilibrate.py
@@ -374,7 +374,7 @@
     equality_constraints = [("P", P0), ("T", T0)]
     sol, prm = equilibrate(composition, assemblage, equality_constraints)
 
-    S = np.array([assemblage.molar_entropy * assemblage.n_moles])
+    S = np.array([assemblage.entropy])
 
     # Now we calculate the assemblages and states at the post-perovskite-in and
     # bridgmanite-out boundaries at the desired entropy.
diff --git a/examples/example_optimal_thermobarometry.py b/examples/example_optimal_thermobarometry.py
index 1dbb49768..a0e22e699 100644
--- a/examples/example_optimal_thermobarometry.py
+++ b/examples/example_optimal_thermobarometry.py
@@ -63,7 +63,7 @@
     temperatures = np.linspace(600, 800, 5)
     equality_constraints = [["T", temperatures], ["phase_fraction", (sill, 0.0)]]
     sols, prm = equilibrate(
-        assemblage.formula, assemblage, equality_constraints=equality_constraints
+        ky.formula, assemblage, equality_constraints=equality_constraints
     )
     Ps = np.array([sol.assemblage.pressure for sol in sols])
     Ts = np.array([sol.assemblage.temperature for sol in sols])
diff --git a/misc/benchmarks/slb_2024_benchmarks.py b/misc/benchmarks/slb_2024_benchmarks.py
index 789545afe..37948bfc7 100644
--- a/misc/benchmarks/slb_2024_benchmarks.py
+++ b/misc/benchmarks/slb_2024_benchmarks.py
@@ -50,6 +50,7 @@ def check_fper_entropy():
     print("Endmember entropies (should all be zero):")
     print(np.abs(fper.solution_model.endmember_configurational_entropies))
 
+    fper.set_state(1.0e5, 300.0)
     fper.set_composition([0.0, 0.5, 0.0, 0.0, 0.5])
     Sxs = -2 * gas_constant * np.log(0.5)
     print(
diff --git a/tests/test_composite.py b/tests/test_composite.py
index 675000de6..b727f6249 100644
--- a/tests/test_composite.py
+++ b/tests/test_composite.py
@@ -502,7 +502,7 @@ def test_single_phase_composite_aliases(self):
         properties = [
             "P",
             "T",
-            "energy",
+            "internal_energy",
             "helmholtz",
             "gibbs",
             "V",
@@ -525,6 +525,8 @@ def test_single_phase_composite_aliases(self):
 
     def test_relaxed_composite(self):
         assemblage = setup_assemblage()
+        old_formula = assemblage.formula.copy()
+        old_number_of_moles = assemblage.number_of_moles
         relaxed_assemblage = RelaxedComposite(assemblage, assemblage.reaction_basis)
         relaxed_assemblage2 = RelaxedComposite(
             assemblage, assemblage.reaction_basis[:2, :]
@@ -536,6 +538,15 @@ def test_relaxed_composite(self):
         relaxed_assemblage.set_state(P, T)
         relaxed_assemblage2.set_state(P, T)
 
+        # Check that the formula is unchanged by isochemical reaction
+        for key in old_formula:
+            self.assertAlmostEqual(old_formula[key], relaxed_assemblage.formula[key])
+            self.assertAlmostEqual(old_formula[key], relaxed_assemblage2.formula[key])
+
+        # Check that the number of moles has changed by isochemical reaction
+        self.assertFalse(old_number_of_moles == relaxed_assemblage.number_of_moles)
+        self.assertFalse(old_number_of_moles == relaxed_assemblage2.number_of_moles)
+
         self.assertTrue(
             check_eos_consistency(
                 assemblage, P=P, T=T, including_shear_properties=False
@@ -549,7 +560,8 @@ def test_relaxed_composite(self):
             )
         )
 
-        # The properties for this assemblage are calculated by relaxing only a subset of reactions,
+        # The properties for this assemblage are calculated by relaxing
+        # only a subset of reactions,
         # so its properties should not be EOS consistent
         self.assertFalse(
             check_eos_consistency(
@@ -557,6 +569,182 @@ def test_relaxed_composite(self):
             )
         )
 
+    def test_formula_extensivity_composite(self):
+        min1 = minerals.SLB_2011.mg_perovskite()
+        min2 = minerals.SLB_2011.periclase()
+        molar_fractions = [0.5, 0.5]
+        c1 = burnman.Composite([min1, min2], molar_fractions, fraction_type="molar")
+        c2 = burnman.Composite([min1, min2], molar_fractions, fraction_type="molar")
+        c2.number_of_moles = 2.0
+
+        self.assertFloatEqual(c1.number_of_moles, c2.number_of_moles / 2.0)
+        self.assertFloatEqual(c1.formula["Mg"], c2.formula["Mg"] / 2.0)
+        self.assertFloatEqual(c1.formula["Si"], c2.formula["Si"] / 2.0)
+        self.assertFloatEqual(c1.formula["O"], c2.formula["O"] / 2.0)
+        self.assertFloatEqual(c1.molar_mass, c2.molar_mass)
+        self.assertFloatEqual(c1.mass, c2.mass / 2.0)
+
+    def test_property_extensivity_composite(self):
+        min1 = minerals.SLB_2011.mg_perovskite()
+        min2 = minerals.SLB_2011.periclase()
+        molar_fractions = [0.5, 0.5]
+        c1 = burnman.Composite([min1, min2], molar_fractions, fraction_type="molar")
+        c2 = burnman.Composite([min1, min2], molar_fractions, fraction_type="molar")
+        c2.number_of_moles = 2.0
+
+        c1.set_state(40.0e9, 2000.0)
+        c2.set_state(40.0e9, 2000.0)
+
+        properties_extensive = [
+            "internal_energy",
+            "helmholtz",
+            "gibbs",
+            "enthalpy",
+            "entropy",
+            "volume",
+            "heat_capacity_p",
+            "heat_capacity_v",
+            "mass",
+            "V",
+            "S",
+            "H",
+            "C_v",
+            "C_p",
+        ]
+
+        for prop in properties_extensive:
+            val1 = getattr(c1, prop)
+            val2 = getattr(c2, prop)
+            self.assertFloatEqual(val1, val2 / 2.0)
+
+        properties_intensive = [
+            "molar_internal_energy",
+            "molar_helmholtz",
+            "molar_gibbs",
+            "molar_enthalpy",
+            "molar_entropy",
+            "molar_volume",
+            "molar_heat_capacity_p",
+            "molar_heat_capacity_v",
+            "isothermal_bulk_modulus_reuss",
+            "isentropic_bulk_modulus_reuss",
+            "isothermal_compressibility_reuss",
+            "isentropic_compressibility_reuss",
+            "grueneisen_parameter",
+            "thermal_expansivity",
+            "molar_mass",
+            "density",
+        ]
+
+        for prop in properties_intensive:
+            val1 = getattr(c1, prop)
+            val2 = getattr(c2, prop)
+            self.assertFloatEqual(val1, val2)
+
+    def test_composite_state_not_set(self):
+        min1 = minerals.SLB_2011.mg_perovskite()
+        min2 = minerals.SLB_2011.periclase()
+        molar_fractions = [0.5, 0.5]
+        c = burnman.Composite([min1, min2], molar_fractions, fraction_type="molar")
+
+        with self.assertRaises(AttributeError):
+            _ = c.density
+
+    def test_formula_extensivity_relaxed_composite(self):
+        fper = minerals.SLB_2024.ferropericlase()
+        bdg = minerals.SLB_2024.bridgmanite()
+        c1_unrelaxed = burnman.Composite([fper, bdg], [0.5, 0.5])
+        c2_unrelaxed = burnman.Composite([fper, bdg], [0.5, 0.5])
+        c2_unrelaxed.number_of_moles = 2.0
+
+        fper.set_composition([0.35, 0.3, 0.2, 0.1, 0.05])
+        bdg.set_composition([0.7, 0.1, 0.1, 0.05, 0.03, 0.01, 0.01])
+
+        formula = c1_unrelaxed.formula.copy()
+        formula2 = c2_unrelaxed.formula.copy()
+
+        self.assertFloatEqual(formula["Mg"] * 2.0, formula2["Mg"])
+
+        c1 = RelaxedComposite(c1_unrelaxed, c1_unrelaxed.reaction_basis)
+        c2 = RelaxedComposite(c2_unrelaxed, c2_unrelaxed.reaction_basis)
+
+        c1.set_state(40.0e9, 2000.0)
+        c2.set_state(40.0e9, 2000.0)
+
+        self.assertFloatEqual(c1.number_of_moles, c2.number_of_moles / 2.0)
+        self.assertFloatEqual(c1.formula["Mg"], c2.formula["Mg"] / 2.0)
+        self.assertFloatEqual(c1.formula["Si"], c2.formula["Si"] / 2.0)
+        self.assertFloatEqual(c1.formula["O"], c2.formula["O"] / 2.0)
+        self.assertFloatEqual(c1.molar_mass, c2.molar_mass)
+        self.assertFloatEqual(c1.mass, c2.mass / 2.0)
+
+    def test_property_extensivity_relaxed_composite(self):
+        fper = minerals.SLB_2024.ferropericlase()
+        bdg = minerals.SLB_2024.bridgmanite()
+        c1_unrelaxed = burnman.Composite([fper, bdg], [0.5, 0.5])
+        c2_unrelaxed = burnman.Composite([fper, bdg], [0.5, 0.5])
+        c2_unrelaxed.number_of_moles = 2.0
+
+        fper.set_composition([0.35, 0.3, 0.2, 0.1, 0.05])
+        bdg.set_composition([0.7, 0.1, 0.1, 0.05, 0.03, 0.01, 0.01])
+
+        formula = c1_unrelaxed.formula.copy()
+        formula2 = c2_unrelaxed.formula.copy()
+
+        self.assertFloatEqual(formula["Mg"] * 2.0, formula2["Mg"])
+
+        c1 = RelaxedComposite(c1_unrelaxed, c1_unrelaxed.reaction_basis)
+        c2 = RelaxedComposite(c2_unrelaxed, c2_unrelaxed.reaction_basis)
+
+        c1.set_state(40.0e9, 2000.0)
+        c2.set_state(40.0e9, 2000.0)
+
+        properties_extensive = [
+            "internal_energy",
+            "helmholtz",
+            "gibbs",
+            "enthalpy",
+            "entropy",
+            "volume",
+            "heat_capacity_p",
+            "heat_capacity_v",
+            "mass",
+            "H",
+            "S",
+            "V",
+            "C_p",
+            "C_v",
+        ]
+
+        for prop in properties_extensive:
+            val1 = getattr(c1, prop)
+            val2 = getattr(c2, prop)
+            self.assertFloatEqual(val1, val2 / 2.0)
+
+        properties_intensive = [
+            "molar_internal_energy",
+            "molar_helmholtz",
+            "molar_gibbs",
+            "molar_enthalpy",
+            "molar_entropy",
+            "molar_volume",
+            "molar_heat_capacity_p",
+            "molar_heat_capacity_v",
+            "isothermal_bulk_modulus_reuss",
+            "isentropic_bulk_modulus_reuss",
+            "isothermal_compressibility_reuss",
+            "isentropic_compressibility_reuss",
+            "grueneisen_parameter",
+            "thermal_expansivity",
+            "molar_mass",
+            "density",
+        ]
+
+        for prop in properties_intensive:
+            val1 = getattr(c1, prop)
+            val2 = getattr(c2, prop)
+            self.assertFloatEqual(val1, val2)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tests/test_equilibration.py b/tests/test_equilibration.py
index bbdc48ac1..49b3df0fb 100644
--- a/tests/test_equilibration.py
+++ b/tests/test_equilibration.py
@@ -99,6 +99,32 @@ def test_univariant_line_fo(self):
         ]
         self.assertArraysAlmostEqual(Ps, Ps_ref)
 
+    def test_conservation_of_elements(self):
+        assemblage = make_ol_wad_assemblage()
+        ol = assemblage.phases[0]
+        wad = assemblage.phases[1]
+        assemblage.set_state(10.0e9, 1200.0)
+        equality_constraints = [
+            ("P", 10.0e9),
+            (
+                "phase_composition",
+                (ol, (["Mg_A", "Fe_A"], [0.0, 1.0], [1.0, 1.0], 0.45)),
+            ),
+        ]
+        composition = {"Mg": 1.0, "Fe": 1.0, "Si": 1.0, "O": 4.0}
+
+        sol, prm = equilibrate(composition, assemblage, equality_constraints)
+        for key in composition:
+            self.assertAlmostEqual(assemblage.formula[key], composition[key], places=5)
+        old_number_of_moles = assemblage.number_of_moles
+
+        # now double the composition and see that number of moles doubles
+        composition = {key: 2.0 * composition[key] for key in composition}
+        sol, prm = equilibrate(composition, assemblage, equality_constraints)
+        self.assertAlmostEqual(
+            assemblage.number_of_moles, 2.0 * old_number_of_moles, places=5
+        )
+
     def test_ol_wad_eqm(self):
         assemblage = make_ol_wad_assemblage()
         ol = assemblage.phases[0]
@@ -198,7 +224,7 @@ def test_ol_wad_eqm_entropy(self):
 
         P = assemblage.pressure
         T = assemblage.temperature
-        S = assemblage.molar_entropy * assemblage.n_moles
+        S = assemblage.entropy
 
         assemblage = make_ol_wad_assemblage()
         equality_constraints = [("P", P), ("S", S)]
@@ -220,7 +246,7 @@ def test_ol_wad_eqm_volume(self):
 
         P = assemblage.pressure
         T = assemblage.temperature
-        V = assemblage.molar_volume * assemblage.n_moles
+        V = assemblage.molar_volume * assemblage.number_of_moles
 
         assemblage = make_ol_wad_assemblage()
         equality_constraints = [("P", P), ("V", V)]
diff --git a/tests/test_material.py b/tests/test_material.py
index 324ca8adf..8d2f48605 100644
--- a/tests/test_material.py
+++ b/tests/test_material.py
@@ -133,6 +133,15 @@ def test_set_state_with_volume(self):
         m.set_state(P1, T1)  # return to new state
         self.assertFloatEqual(V, m.V)
 
+    def test_unset_state_raises(self):
+        m = self.min_with_name()
+        with self.assertRaises(AttributeError):
+            _ = m.pressure
+        with self.assertRaises(AttributeError):
+            _ = m.temperature
+        with self.assertRaises(AttributeError):
+            _ = m.V
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tutorial/tutorial_05_equilibrium.ipynb b/tutorial/tutorial_05_equilibrium.ipynb
index 6f6ee1829..ed14a9e43 100644
--- a/tutorial/tutorial_05_equilibrium.ipynb
+++ b/tutorial/tutorial_05_equilibrium.ipynb
@@ -112,8 +112,9 @@
         "print(sol.assemblage)\n",
         "\n",
         "# The total entropy of the assemblage is the molar entropy\n",
-        "# multiplied by the number of moles in the assemblage\n",
-        "entropy = sol.assemblage.S*sol.assemblage.n_moles"
+        "# multiplied by the number of moles in the assemblage. This can be\n",
+        "# accessed either by sol.assemblage.entropy, or by the shorter alias:\n",
+        "entropy = sol.assemblage.S"
       ]
     },
     {