diff --git a/burnman/classes/material.py b/burnman/classes/material.py
index fd32de37..3e7c36ac 100644
--- a/burnman/classes/material.py
+++ b/burnman/classes/material.py
@@ -283,6 +283,11 @@ def evaluate(self, vars_list, pressures, temperatures, molar_fractions=None):
         """
         old_pressure = self.pressure
         old_temperature = self.temperature
+        try:
+            old_molar_fractions = self.molar_fractions
+        except AttributeError:
+            old_molar_fractions = None
+
         pressures = np.array(pressures)
         temperatures = np.array(temperatures)
 
@@ -290,6 +295,7 @@ def evaluate(self, vars_list, pressures, temperatures, molar_fractions=None):
 
         if molar_fractions is not None:
             molar_fractions = np.array(molar_fractions)
+            self.set_composition(molar_fractions[0])
             assert temperatures.shape == molar_fractions.shape[:-1]
 
         # First, check the output types of all the requested variables:
@@ -316,6 +322,11 @@ def evaluate(self, vars_list, pressures, temperatures, molar_fractions=None):
             self.reset()
         else:
             self.set_state(old_pressure, old_temperature)
+            if old_molar_fractions is not None:
+                try:
+                    self.set_composition(old_molar_fractions)
+                except AttributeError:
+                    pass
 
         try:
             output = np.array(output)
@@ -349,6 +360,11 @@ def evaluate_with_volumes(
         """
         old_pressure = self.pressure
         old_temperature = self.temperature
+        try:
+            old_molar_fractions = self.molar_fractions
+        except AttributeError:
+            old_molar_fractions = None
+
         volumes = np.array(volumes)
         temperatures = np.array(temperatures)
 
@@ -356,6 +372,7 @@ def evaluate_with_volumes(
 
         if molar_fractions is not None:
             molar_fractions = np.array(molar_fractions)
+            self.set_composition(molar_fractions[0])
             assert temperatures.shape == molar_fractions.shape[:-1]
 
         # First, check the output types of all the requested variables:
@@ -383,6 +400,12 @@ def evaluate_with_volumes(
         else:
             self.set_state(old_pressure, old_temperature)
 
+            if old_molar_fractions is not None:
+                try:
+                    self.set_composition(old_molar_fractions)
+                except AttributeError:
+                    pass
+
         try:
             output = np.array(output)
         except ValueError:  # if the lists are different shapes