diff --git a/espei/error_functions/activity_error.py b/espei/error_functions/activity_error.py
index a2d5e239..690f3b0c 100644
--- a/espei/error_functions/activity_error.py
+++ b/espei/error_functions/activity_error.py
@@ -19,6 +19,7 @@
 from pycalphad.core.utils import filter_phases, unpack_species
 from scipy.stats import norm
 from pycalphad import Workspace
+from pycalphad.property_framework import JanssonDerivative
 
 from espei.core_utils import ravel_conditions
 from espei.error_functions.residual_base import ResidualFunction, residual_function_registry
@@ -29,7 +30,7 @@
 _log = logging.getLogger(__name__)
 
 
-def target_chempots_from_activity(component, target_activity, temperatures, wks_ref):
+def target_chempots_from_activity(component, parameters, target_activity, temperatures, wks_ref):
     """
     Return an array of experimental chemical potentials for the component
 
@@ -49,15 +50,21 @@ def target_chempots_from_activity(component, target_activity, temperatures, wks_
     numpy.ndarray
         Array of experimental chemical potentials
     """
-    # acr_i = exp((mu_i - mu_i^{ref})/(RT))
-    # so mu_i = R*T*ln(acr_i) + mu_i^{ref}
+    
     ref_chempot = wks_ref.get(v.MU(component))
     exp_chem_pots = v.R * temperatures * np.log(target_activity) + ref_chempot
-    return exp_chem_pots
+    
+    gradient_props = [JanssonDerivative(v.MU(component), key) for key in parameters]
+    gradients = wks_ref.get(*gradient_props)
+    if type(gradients) is list:
+        ref_grads = [float(element) for element in gradients]
+    else:
+        ref_grads = gradients
+    return exp_chem_pots, ref_grads
 
 
 # TODO: roll this function into ActivityResidual
-def calculate_activity_residuals(dbf, comps, phases, datasets, parameters=None, phase_models=None, callables=None, data_weight=1.0) -> Tuple[List[float], List[float]]:
+def calculate_activity_residuals(dbf, comps, phases, datasets, parameters=None, phase_models=None, callables=None, data_weight=1.0) -> Tuple[List[float], List[float], List[float]]:
     """
     Notes
     -----
@@ -75,12 +82,23 @@ def calculate_activity_residuals(dbf, comps, phases, datasets, parameters=None,
     if parameters is None:
         parameters = {}
 
+    params_keys = []
+
+    # This mutates the global pycalphad namespace
+    for key in parameters.keys():
+        if not hasattr(v, key):
+            setattr(v, key, v.IndependentPotential(key))
+        params_keys.append(getattr(v, key))
+        # Mutates argument to function
+        dbf.symbols.pop(key,None)
+        
     activity_datasets = datasets.search(
         (tinydb.where('output').test(lambda x: 'ACR' in x)) &
         (tinydb.where('components').test(lambda x: set(x).issubset(comps))))
 
     residuals = []
     weights = []
+    gradients = []
     for ds in activity_datasets:
         acr_component = ds['output'].split('_')[1]  # the component of interest
         # calculate the reference state equilibrium
@@ -90,7 +108,9 @@ def calculate_activity_residuals(dbf, comps, phases, datasets, parameters=None,
         data_comps = ds['components']
         data_phases = filter_phases(dbf, unpack_species(dbf, data_comps), candidate_phases=phases)
         ref_conditions = {_map_coord_to_variable(coord): val for coord, val in ref['conditions'].items()}
-        wks_ref = Workspace(database=dbf, components=data_comps, phases=ref['phases'], conditions=ref_conditions, parameters=parameters)
+        # removed parameter assignment from wks_ref
+        ref_conditions.update(parameters)
+        wks_ref = Workspace(database=dbf, components=data_comps, phases=ref['phases'], conditions=ref_conditions)
 
         # calculate current chemical potentials
         # get the conditions
@@ -102,6 +122,7 @@ def calculate_activity_residuals(dbf, comps, phases, datasets, parameters=None,
         # we will ravel each composition individually, since they all must have the same shape
         dataset_computed_chempots = []
         dataset_weights = []
+        dataset_gradients = []
         for comp_name, comp_x in conds_list:
             P, T, X = ravel_conditions(ds['values'], ds['conditions']['P'], ds['conditions']['T'], comp_x)
             conditions[v.P] = P
@@ -113,22 +134,34 @@ def calculate_activity_residuals(dbf, comps, phases, datasets, parameters=None,
         # assume now that the ravelled conditions all have the same size
         conditions_list = [{c: conditions[c][i] for c in conditions.keys()} for i in range(len(conditions[v.T]))]
         for conds in conditions_list:
-            wks_sample = Workspace(database=dbf, components=data_comps, phases=data_phases, conditions=conds, parameters=parameters)
+            conds.update(parameters)
+            wks_sample = Workspace(database=dbf, components=data_comps, phases=data_phases, conditions=conds)
             dataset_computed_chempots.append(wks_sample.get(v.MU(acr_component)))
             dataset_weights.append(std_dev / data_weight / ds.get("weight", 1.0))
+            gradient_props = [JanssonDerivative(v.MU(acr_component), key) for key in parameters]
+            grads = wks_sample.get(*gradient_props)
+            if type(grads) is list:
+                sample_grads = [float(element) for element in grads]
+            else:
+                sample_grads = grads
+            dataset_gradients.append(sample_grads)
 
         # calculate target chempots
         dataset_activities = np.array(ds['values']).flatten()
-        dataset_target_chempots = target_chempots_from_activity(acr_component, dataset_activities, conditions[v.T], wks_ref)
+        dataset_target_chempots, ref_grads = target_chempots_from_activity(acr_component, parameters, dataset_activities, conditions[v.T], wks_ref)
         dataset_residuals = (np.asarray(dataset_computed_chempots) - np.asarray(dataset_target_chempots, dtype=float)).tolist()
+        adjusted_gradient = []
+        for element in dataset_gradients:
+            adjusted_gradient.append((np.asarray(element) - np.asarray(ref_grads)).tolist())
         _log.debug('Data: %s, chemical potential difference: %s, reference: %s', dataset_activities, dataset_residuals, ds["reference"])
         residuals.extend(dataset_residuals)
         weights.extend(dataset_weights)
-    return residuals, weights
+        gradients.append(adjusted_gradient)
+    return residuals, weights, gradients
 
 
 # TODO: roll this function into ActivityResidual
-def calculate_activity_error(dbf, comps, phases, datasets, parameters=None, phase_models=None, callables=None, data_weight=1.0) -> float:
+def calculate_activity_error(dbf, comps, phases, datasets, parameters=None, phase_models=None, callables=None, data_weight=1.0) -> Tuple[float, List[float]]:
     """
     Return the sum of square error from activity data
 
@@ -160,14 +193,23 @@ def calculate_activity_error(dbf, comps, phases, datasets, parameters=None, phas
 
 
     """
-    residuals, weights = calculate_activity_residuals(dbf, comps, phases, datasets, parameters=parameters, phase_models=phase_models, callables=callables, data_weight=data_weight)
+    residuals, weights, gradients = calculate_activity_residuals(dbf, comps, phases, datasets, parameters=parameters, phase_models=phase_models, callables=callables, data_weight=data_weight)
     likelihood = np.sum(norm(0, scale=weights).logpdf(residuals))
+    if len(gradients) == 0:
+        likelihood_grads = []
+    else: 
+        gradients = np.concatenate(gradients)
+        derivative = -np.array(residuals)*np.array(gradients).T/np.array(weights)**2
+        if derivative.ndim == 1:
+            likelihood_grads = np.sum(derivative, axis=0)
+        else:
+            likelihood_grads = np.sum(derivative, axis=1)
     if np.isnan(likelihood):
         # TODO: revisit this case and evaluate whether it is resonable for NaN
         # to show up here. When this comment was written, the test
         # test_subsystem_activity_probability would trigger a NaN.
-        return -np.inf
-    return likelihood
+        return -np.inf, np.zeros(len(parameters))
+    return likelihood, likelihood_grads
 
 
 # TODO: the __init__ method should pre-compute Model and PhaseRecord objects
@@ -215,13 +257,14 @@ def __init__(
 
     def get_residuals(self, parameters: npt.ArrayLike) -> Tuple[List[float], List[float]]:
         parameters = {param_name: param for param_name, param in zip(self._symbols_to_fit, parameters.tolist())}
-        residuals, weights = calculate_activity_residuals(parameters=parameters, **self._activity_likelihood_kwargs)
+        residuals, weights, grads = calculate_activity_residuals(parameters=parameters, **self._activity_likelihood_kwargs)
         return residuals, weights
 
-    def get_likelihood(self, parameters: npt.NDArray) -> float:
+    def get_likelihood(self, parameters: npt.NDArray) -> Tuple[float, List[float]]:
         parameters = {param_name: param for param_name, param in zip(self._symbols_to_fit, parameters.tolist())}
-        likelihood = calculate_activity_error(parameters=parameters, **self._activity_likelihood_kwargs)
-        return likelihood
+        likelihood, gradients = calculate_activity_error(parameters=parameters, **self._activity_likelihood_kwargs)
+        return likelihood, gradients
 
 
 residual_function_registry.register(ActivityResidual)
+
diff --git a/espei/error_functions/equilibrium_thermochemical_error.py b/espei/error_functions/equilibrium_thermochemical_error.py
index 01d31e71..da9872f5 100644
--- a/espei/error_functions/equilibrium_thermochemical_error.py
+++ b/espei/error_functions/equilibrium_thermochemical_error.py
@@ -12,13 +12,13 @@
 from tinydb import where
 from scipy.stats import norm
 from pycalphad import Database, Model, ReferenceState, variables as v
-from pycalphad.core.utils import instantiate_models, filter_phases, extract_parameters, unpack_species, unpack_condition
+from pycalphad.core.utils import instantiate_models, filter_phases, unpack_species, unpack_condition
 from pycalphad.codegen.phase_record_factory import PhaseRecordFactory
 from pycalphad import Workspace, as_property
+from pycalphad.property_framework import JanssonDerivative
 
 from espei.error_functions.residual_base import ResidualFunction, residual_function_registry
 from espei.phase_models import PhaseModelSpecification
-from espei.shadow_functions import update_phase_record_parameters
 from espei.typing import SymbolName
 from espei.utils import PickleableTinyDB, database_symbols_to_fit
 
@@ -73,13 +73,21 @@ def build_eqpropdata(data: tinydb.database.Document,
         'SM':   0.2,  # J/K-mol
         'CPM':  0.2,  # J/K-mol
     }
+    
+    params_keys = []
 
-    params_keys, _ = extract_parameters(parameters)
+    # This mutates the global pycalphad namespace
+    for key in parameters.keys():
+        if not hasattr(v, key):
+            setattr(v, key, v.IndependentPotential(key))
+        params_keys.append(getattr(v, key))
+        # Mutates argument to function
+        dbf.symbols.pop(key,None)
 
     data_comps = list(set(data['components']).union({'VA'}))
     species = sorted(unpack_species(dbf, data_comps), key=str)
     data_phases = filter_phases(dbf, species, candidate_phases=data['phases'])
-    models = instantiate_models(dbf, species, data_phases, model=model, parameters=parameters)
+    models = instantiate_models(dbf, species, data_phases, model=model)
     output = data['output']
     property_output = output.split('_')[0]  # property without _FORM, _MIX, etc.
     samples = np.array(data['values']).flatten()
@@ -99,14 +107,7 @@ def build_eqpropdata(data: tinydb.database.Document,
     pot_conds = OrderedDict([(getattr(v, key), unpack_condition(data['conditions'][key])) for key in sorted(data['conditions'].keys()) if not key.startswith('X_')])
     comp_conds = OrderedDict([(v.X(key[2:]), unpack_condition(data['conditions'][key])) for key in sorted(data['conditions'].keys()) if key.startswith('X_')])
 
-    phase_record_factory = PhaseRecordFactory(dbf, species, {**pot_conds, **comp_conds}, models, parameters=parameters)
-
-    # Now we need to unravel the composition conditions
-    # (from Dict[v.X, Sequence[float]] to Sequence[Dict[v.X, float]]), since the
-    # composition conditions are only broadcast against the potentials, not
-    # each other. Each individual composition needs to be computed
-    # independently, since broadcasting over composition cannot be turned off
-    # in pycalphad.
+    phase_record_factory = PhaseRecordFactory(dbf, species, {**pot_conds, **comp_conds, **parameters}, models)
     rav_comp_conds = [OrderedDict(zip(comp_conds.keys(), pt_comps)) for pt_comps in zip(*comp_conds.values())]
 
     # Build weights, should be the same size as the values
@@ -173,7 +174,7 @@ def get_equilibrium_thermochemical_data(dbf: Database, comps: Sequence[str],
 def calc_prop_differences(eqpropdata: EqPropData,
                           parameters: np.ndarray,
                           approximate_equilibrium: Optional[bool] = False,
-                          ) -> Tuple[np.ndarray, np.ndarray]:
+                          ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
     """
     Calculate differences between the expected and calculated values for a property
 
@@ -196,42 +197,45 @@ def calc_prop_differences(eqpropdata: EqPropData,
         * weights for this dataset
 
     """
+        
     dbf = eqpropdata.dbf
     species = eqpropdata.species
     phases = eqpropdata.phases
     pot_conds = eqpropdata.potential_conds
     models = eqpropdata.models
     phase_record_factory = eqpropdata.phase_record_factory
-    update_phase_record_parameters(phase_record_factory, parameters)
     params_dict = OrderedDict(zip(map(str, eqpropdata.params_keys), parameters))
     output = as_property(eqpropdata.output)
     weights = np.array(eqpropdata.weight, dtype=np.float64)
     samples = np.array(eqpropdata.samples, dtype=np.float64)
-    wks = Workspace(database=dbf, components=species, phases=phases, models=models, phase_record_factory=phase_record_factory, parameters=params_dict)
+    wks = Workspace(database=dbf, components=species, phases=phases, models=models, phase_record_factory=phase_record_factory)
 
     calculated_data = []
+    gradient_data = []
     for comp_conds in eqpropdata.composition_conds:
-        cond_dict = OrderedDict(**pot_conds, **comp_conds)
+        cond_dict = OrderedDict(**pot_conds, **comp_conds, **params_dict)
         wks.conditions = cond_dict
-        wks.parameters = params_dict  # these reset models and phase_record_factory through depends_on -> lose Model.shift_reference_state, etc.
         wks.models = models
-        wks.phase_record_factory = phase_record_factory
         vals = wks.get(output)
         calculated_data.extend(np.atleast_1d(vals).tolist())
+        gradient_props = [JanssonDerivative(output, key) for key in params_dict]
+        gradients = wks.get(*gradient_props)
+        gradient_data.append(gradients)
 
     calculated_data = np.array(calculated_data, dtype=np.float64)
+    gradient_data = np.array(gradient_data, dtype=np.float64)
 
     assert calculated_data.shape == samples.shape, f"Calculated data shape {calculated_data.shape} does not match samples shape {samples.shape}"
     assert calculated_data.shape == weights.shape, f"Calculated data shape {calculated_data.shape} does not match weights shape {weights.shape}"
     differences = calculated_data - samples
     _log.debug('Output: %s differences: %s, weights: %s, reference: %s', output, differences, weights, eqpropdata.reference)
-    return differences, weights
+    return differences, weights, gradient_data
 
 
 def calculate_equilibrium_thermochemical_probability(eq_thermochemical_data: Sequence[EqPropData],
                                                      parameters: np.ndarray,
                                                      approximate_equilibrium: Optional[bool] = False,
-                                                     ) -> float:
+                                                     ) -> Tuple[float, List[float]]:
     """
     Calculate the total equilibrium thermochemical probability for all EqPropData
 
@@ -252,23 +256,31 @@ def calculate_equilibrium_thermochemical_probability(eq_thermochemical_data: Seq
 
     """
     if len(eq_thermochemical_data) == 0:
-        return 0.0
+        return 0.0, np.zeros(len(parameters))
 
     differences = []
     weights = []
+    gradients = []
     for eqpropdata in eq_thermochemical_data:
-        diffs, wts = calc_prop_differences(eqpropdata, parameters, approximate_equilibrium)
+        diffs, wts, grads = calc_prop_differences(eqpropdata, parameters, approximate_equilibrium)
         if np.any(np.isinf(diffs) | np.isnan(diffs)):
             # NaN or infinity are assumed calculation failures. If we are
             # calculating log-probability, just bail out and return -infinity.
-            return -np.inf
+            return -np.inf, np.zeros(len(parameters))
         differences.append(diffs)
         weights.append(wts)
+        gradients.append(grads)
 
     differences = np.concatenate(differences, axis=0)
     weights = np.concatenate(weights, axis=0)
+    gradients = np.concatenate(gradients, axis=0)
     probs = norm(loc=0.0, scale=weights).logpdf(differences)
-    return np.sum(probs)
+    grad_probs = -differences*gradients.T/weights**2
+    if grad_probs.ndim == 1:
+        grad_probs_sum = np.sum(grad_probs, axis=0)
+    else:
+        grad_probs_sum = np.sum(grad_probs, axis=1)
+    return np.sum(probs), grad_probs_sum
 
 
 class EquilibriumPropertyResidual(ResidualFunction):
@@ -304,14 +316,14 @@ def get_residuals(self, parameters: npt.ArrayLike) -> Tuple[List[float], List[fl
         residuals = []
         weights = []
         for data in self.property_data:
-            dataset_residuals, dataset_weights = calc_prop_differences(data, parameters)
+            dataset_residuals, dataset_weights, dataset_grads = calc_prop_differences(data, parameters)
             residuals.extend(dataset_residuals.tolist())
             weights.extend(dataset_weights.tolist())
         return residuals, weights
 
-    def get_likelihood(self, parameters) -> float:
-        likelihood = calculate_equilibrium_thermochemical_probability(self.property_data, parameters)
-        return likelihood
+    def get_likelihood(self, parameters) -> Tuple[float, List[float]]:
+        likelihood, gradients = calculate_equilibrium_thermochemical_probability(self.property_data, parameters)
+        return likelihood, gradients
 
 
-residual_function_registry.register(EquilibriumPropertyResidual)
\ No newline at end of file
+residual_function_registry.register(EquilibriumPropertyResidual)
diff --git a/espei/error_functions/non_equilibrium_thermochemical_error.py b/espei/error_functions/non_equilibrium_thermochemical_error.py
index 60dd806c..02d8f21d 100644
--- a/espei/error_functions/non_equilibrium_thermochemical_error.py
+++ b/espei/error_functions/non_equilibrium_thermochemical_error.py
@@ -16,7 +16,7 @@
 from pycalphad import Database, Model, ReferenceState, variables as v
 from pycalphad.core.utils import unpack_species, get_pure_elements, filter_phases
 from pycalphad import Workspace
-from pycalphad.property_framework import IsolatedPhase
+from pycalphad.property_framework import IsolatedPhase, ModelComputedProperty
 from pycalphad.property_framework.metaproperties import find_first_compset
 from pycalphad.core.solver import Solver, SolverResult
 
@@ -249,7 +249,7 @@ def get_sample_condition_dicts(calculate_dict: Dict[Any, Any], configuration_tup
     return sample_condition_dicts
 
 
-def get_thermochemical_data(dbf, comps, phases, datasets, model=None, weight_dict=None, symbols_to_fit=None):
+def get_thermochemical_data(dbf, comps, phases, datasets, parameters, model=None, weight_dict=None, symbols_to_fit=None):
     """
 
     Parameters
@@ -274,6 +274,17 @@ def get_thermochemical_data(dbf, comps, phases, datasets, model=None, weight_dic
     list
         List of data dictionaries to iterate over
     """
+    
+    params_keys = []
+
+    # This mutates the global pycalphad namespace
+    for key in parameters.keys():
+        if not hasattr(v, key):
+            setattr(v, key, v.IndependentPotential(key))
+        params_keys.append(getattr(v, key))
+        # Mutates argument to function
+        dbf.symbols.pop(key,None)
+        
     # phase by phase, then property by property, then by model exclusions
     if weight_dict is None:
         weight_dict = {}
@@ -330,7 +341,7 @@ def get_thermochemical_data(dbf, comps, phases, datasets, model=None, weight_dic
                 curr_data = filter_temperatures(curr_data)
                 calculate_dict = get_prop_samples(curr_data, constituents)
                 model_cls = model.get(phase_name, Model)
-                mod = model_cls(dbf, comps, phase_name, parameters=symbols_to_fit)
+                mod = model_cls(dbf, comps, phase_name)
                 if prop.endswith('_FORM'):
                     output = ''.join(prop.split('_')[:-1])+"R"
                     mod.shift_reference_state(ref_states, dbf, contrib_mods={e: symengine.S.Zero for e in exclusion})
@@ -368,7 +379,6 @@ def compute_fixed_configuration_property_differences(dbf, calc_data: FixedConfig
     phase_name = calc_data['phase_name']
     models = calc_data['model']  # Dict[PhaseName: Model]
     output = calc_data['output']
-    phase_record_factory = calc_data['phase_record_factory']
     sample_values = calc_data['calculate_dict']['values']
     str_statevar_dict = calc_data['str_statevar_dict']
 
@@ -382,8 +392,9 @@ def compute_fixed_configuration_property_differences(dbf, calc_data: FixedConfig
         counter = counter + 1
 
     differences = []
+    gradients = []
     for index in range(len(sample_values)):
-        cond_dict = {}
+        cond_dict = {**parameters}
         for sv_key, sv_val in str_statevar_dict.items():
             cond_dict.update({sv_key: sv_val[index]})
 
@@ -407,7 +418,7 @@ def compute_fixed_configuration_property_differences(dbf, calc_data: FixedConfig
         # Need to be careful here. Making a workspace erases the custom models that have some contributions excluded (which are passed in). Not sure exactly why.
         # The models themselves are preserved, but the ones inside the workspace's phase_record_factory get clobbered.
         # We workaround this by replacing the phase_record_factory models with ours, but this is definitely a hack we'd like to avoid.
-        wks = Workspace(database=dbf, components=species, phases=[phase_name], conditions={**cond_dict}, models=models, phase_record_factory=phase_record_factory, parameters=parameters, solver=NoSolveSolver())
+        wks = Workspace(database=dbf, components=species, phases=[phase_name], conditions={**cond_dict}, models=models, solver=NoSolveSolver())
         # We then get a composition set and we use a special "NoSolveSolver" to
         # ensure that we don't change from the data-specified DOF.
         compset = find_first_compset(phase_name, wks)
@@ -417,9 +428,18 @@ def compute_fixed_configuration_property_differences(dbf, calc_data: FixedConfig
         iso_phase = IsolatedPhase(compset, wks=wks)
         iso_phase.solver = NoSolveSolver()
         results = wks.get(iso_phase(output))
+        
+        prop = ModelComputedProperty(output)
+        # chemical potentials are an input to _compute_property_gradient, but they aren't used, so here they are just
+        # passed as zeros
+        grads = prop._compute_property_gradient([compset], wks.conditions,np.zeros(len(active_pure_elements)))
+        grads = grads[0][len(str_statevar_dict):len(str_statevar_dict)+len(parameters)]
+        
+                
         sample_differences = results - sample_values[index]
         differences.append(sample_differences)
-    return differences
+        gradients.append(grads)
+    return differences, gradients
 
 
 def calculate_non_equilibrium_thermochemical_probability(thermochemical_data: List[FixedConfigurationCalculationData], dbf, parameters=None):
@@ -443,16 +463,24 @@ def calculate_non_equilibrium_thermochemical_probability(thermochemical_data: Li
         parameters = {}
 
     prob_error = 0.0
+    overall_grad = np.zeros(len(parameters))
     for data in thermochemical_data:
         phase_name = data['phase_name']
         sample_values = data['calculate_dict']['values']
-        differences = compute_fixed_configuration_property_differences(dbf, data, parameters)
+        differences, gradients = compute_fixed_configuration_property_differences(dbf, data, parameters)
         differences = np.array(differences)
+        gradients = np.array(gradients)
         probabilities = norm.logpdf(differences, loc=0, scale=data['weights'])
         prob_sum = np.sum(probabilities)
         _log.debug("%s(%s) - probability sum: %0.2f, data: %s, differences: %s, probabilities: %s, references: %s", data['prop'], phase_name, prob_sum, sample_values, differences, probabilities, data['calculate_dict']['references'])
         prob_error += prob_sum
-    return prob_error
+        derivative = -differences*gradients.T/data['weights']**2
+        if derivative.ndim == 1:
+            grad_prob = np.sum(derivative, axis = 0)
+        else:
+            grad_prob = np.sum(derivative, axis=1)
+        overall_grad += grad_prob
+    return prob_error, overall_grad
 
 
 class FixedConfigurationPropertyResidual(ResidualFunction):
@@ -480,7 +508,8 @@ def __init__(
         phases = sorted(filter_phases(database, unpack_species(database, comps), database.phases.keys()))
         if symbols_to_fit is None:
             symbols_to_fit = database_symbols_to_fit(database)
-        self.thermochemical_data = get_thermochemical_data(database, comps, phases, datasets, model_dict, weight_dict=self.weight, symbols_to_fit=symbols_to_fit)
+        parameters = dict(zip(symbols_to_fit, [0]*len(symbols_to_fit)))
+        self.thermochemical_data = get_thermochemical_data(database, comps, phases, datasets, parameters, model_dict, weight_dict=self.weight, symbols_to_fit=symbols_to_fit)
         self._symbols_to_fit = symbols_to_fit
         self.dbf = database
 
@@ -488,7 +517,7 @@ def get_residuals(self, parameters: npt.ArrayLike) -> Tuple[List[float], List[fl
         residuals = []
         weights = []
         for data in self.thermochemical_data:
-            dataset_residuals = compute_fixed_configuration_property_differences(self.dbf, data, dict(zip(self._symbols_to_fit, parameters)))
+            dataset_residuals, grads = compute_fixed_configuration_property_differences(self.dbf, data, dict(zip(self._symbols_to_fit, parameters)))
             residuals.extend(dataset_residuals)
             dataset_weights = np.asarray(data["weights"], dtype=float).flatten().tolist()
             if len(dataset_weights) != len(dataset_residuals):
@@ -498,10 +527,10 @@ def get_residuals(self, parameters: npt.ArrayLike) -> Tuple[List[float], List[fl
             weights.extend(dataset_weights)
         return residuals, weights
 
-    def get_likelihood(self, parameters) -> float:
+    def get_likelihood(self, parameters) -> Tuple[float, List[float]]:
         parameters = {param_name: param for param_name, param in zip(self._symbols_to_fit, parameters.tolist())}
-        likelihood = calculate_non_equilibrium_thermochemical_probability(self.thermochemical_data, self.dbf, parameters)
-        return likelihood
+        likelihood, gradient = calculate_non_equilibrium_thermochemical_probability(self.thermochemical_data, self.dbf, parameters)
+        return likelihood, gradient
 
 
-residual_function_registry.register(FixedConfigurationPropertyResidual)
\ No newline at end of file
+residual_function_registry.register(FixedConfigurationPropertyResidual)
diff --git a/espei/error_functions/zpf_error.py b/espei/error_functions/zpf_error.py
index f6d961e2..121b4a8a 100644
--- a/espei/error_functions/zpf_error.py
+++ b/espei/error_functions/zpf_error.py
@@ -22,14 +22,14 @@
 import numpy as np
 from numpy.typing import ArrayLike
 from pycalphad import Database, Model, Workspace, variables as v
-from pycalphad.property_framework import IsolatedPhase
+from pycalphad.property_framework import IsolatedPhase, JanssonDerivative
 from pycalphad.core.utils import filter_phases, unpack_species
 from pycalphad.codegen.phase_record_factory import PhaseRecordFactory
 from scipy.stats import norm
 import tinydb
 
 from espei.phase_models import PhaseModelSpecification
-from espei.shadow_functions import calculate_, update_phase_record_parameters
+from espei.shadow_functions import calculate_
 from espei.typing import SymbolName
 from espei.utils import PickleableTinyDB, database_symbols_to_fit
 from .residual_base import ResidualFunction, residual_function_registry
@@ -141,9 +141,16 @@ def get_zpf_data(dbf: Database, comps: Sequence[str], phases: Sequence[str], dat
     desired_data = datasets.search((tinydb.where('output') == 'ZPF') &
                                    (tinydb.where('components').test(lambda x: set(x).issubset(comps))) &
                                    (tinydb.where('phases').test(lambda x: len(set(phases).intersection(x)) > 0)))
-    wks = Workspace(dbf, comps, phases, parameters=parameters)
+    wks = Workspace(dbf, comps, phases)
     if model is None:
         model = {}
+        
+    params_keys = []
+    for key in parameters.keys():
+        if not hasattr(v, key):
+            setattr(v, key, v.IndependentPotential(key))
+        params_keys.append(getattr(v, key))
+        dbf.symbols.pop(key,None)
 
     zpf_data = []  # 1:1 correspondence with each dataset
     for data in desired_data:
@@ -203,7 +210,7 @@ def get_zpf_data(dbf: Database, comps: Sequence[str], phases: Sequence[str], dat
     return zpf_data
 
 
-def estimate_hyperplane(phase_region: PhaseRegion, dbf: Database, parameters: np.ndarray, approximate_equilibrium: bool = False) -> np.ndarray:
+def estimate_hyperplane(phase_region: PhaseRegion, dbf: Database, parameter_dict, parameters: np.ndarray, approximate_equilibrium: bool = False) -> Tuple[np.ndarray, np.ndarray]:
     """
     Calculate the chemical potentials for the target hyperplane, one vertex at a time
 
@@ -217,25 +224,33 @@ def estimate_hyperplane(phase_region: PhaseRegion, dbf: Database, parameters: np
 
     """
     target_hyperplane_chempots = []
+    target_hyperplane_chempots_grads = []
     species = phase_region.species
     phases = phase_region.phases
     models = phase_region.models
-    param_keys = list(models.values())[0]._parameters_arg
-    parameters_dict = dict(zip(sorted(map(str, param_keys)), parameters))
     for vertex in phase_region.hyperplane_vertices:
-        update_phase_record_parameters(vertex.phase_record_factory, parameters)
         cond_dict = {**vertex.comp_conds, **phase_region.potential_conds}
+        params_keys = list(parameter_dict.keys())
+        for index in range(len(parameters)):
+            cond_dict.update({params_keys[index]: parameters[index]})
         if vertex.has_missing_comp_cond:
             # This composition is unknown -- it doesn't contribute to hyperplane estimation
             pass
         else:
             # Extract chemical potential hyperplane from multi-phase calculation
             # Note that we consider all phases in the system, not just ones in this tie region
-            wks = Workspace(database=dbf, components=species, phases=phases, models=models, phase_record_factory=vertex.phase_record_factory, conditions=cond_dict, parameters=parameters_dict)
+            wks = Workspace(database=dbf, components=species, phases=phases, models=models, phase_record_factory=vertex.phase_record_factory, conditions=cond_dict)
             # TODO: active_pure_elements should be replaced with wks.components when wks.components no longer includes phase constituent Species
             active_pure_elements = [list(x.constituents.keys()) for x in species]
             active_pure_elements = sorted(set(el.upper() for constituents in active_pure_elements for el in constituents) - {"VA"})
-            MU_values = [wks.get(v.MU(comp)) for comp in active_pure_elements]
+            MU_values = [wks.get(v.MU(comp)) for comp in active_pure_elements] 
+            gradient_params = [JanssonDerivative(v.MU(spc), key) for spc in active_pure_elements for key in parameter_dict]
+            gradients = wks.get(*gradient_params)
+            if type(gradients) is list:
+                gradients_magnitude = [float(element) for element in gradients]
+            else:
+                gradients_magnitude = gradients
+            gradients_magnitude = np.reshape(gradients_magnitude,(len(MU_values),len(parameter_dict)))
             num_phases = np.sum(wks.eq.Phase.squeeze() != '')
             Y_values = wks.eq.Y.squeeze()
             no_internal_dof = np.all((np.isclose(Y_values, 1.)) | np.isnan(Y_values))
@@ -243,31 +258,67 @@ def estimate_hyperplane(phase_region: PhaseRegion, dbf: Database, parameters: np
                 target_hyperplane_chempots.append(np.full_like(MU_values, np.nan))
             else:
                 target_hyperplane_chempots.append(MU_values)
+                target_hyperplane_chempots_grads.append(gradients_magnitude)
     target_hyperplane_mean_chempots = np.nanmean(target_hyperplane_chempots, axis=0, dtype=np.float64)
-    return target_hyperplane_mean_chempots
+    target_hyperplane_chempots_grads = np.nanmean(target_hyperplane_chempots_grads, axis=0, dtype=np.float64) 
+    return target_hyperplane_mean_chempots, target_hyperplane_chempots_grads
 
 
-def driving_force_to_hyperplane(target_hyperplane_chempots: np.ndarray,
+def driving_force_to_hyperplane(target_hyperplane_chempots: np.ndarray, target_hyperplane_chempots_grads: np.ndarray,
                                 phase_region: PhaseRegion, dbf: Database, parameter_dict, vertex: RegionVertex,
                                 parameters: np.ndarray, approximate_equilibrium: bool = False) -> Tuple[float,List[float]]:
     """Calculate the integrated driving force between the current hyperplane and target hyperplane.
     """
     species = phase_region.species
     models = phase_region.models
-    param_keys = list(models.values())[0]._parameters_arg
-    parameters_dict = dict(zip(sorted(map(str, param_keys)), parameters))
     current_phase = vertex.phase_name
     cond_dict = {**phase_region.potential_conds, **vertex.comp_conds}
+    params_keys = list(parameter_dict.keys())
+    for index in range(len(parameters)):
+        cond_dict.update({params_keys[index]: parameters[index]})
     str_statevar_dict = OrderedDict([(str(key),cond_dict[key]) for key in sorted(phase_region.potential_conds.keys(), key=str)])
     phase_record_factory = vertex.phase_record_factory
-    update_phase_record_parameters(phase_record_factory, parameters)
     if vertex.has_missing_comp_cond:
+        for index in range(len(parameters)):
+            str_statevar_dict.update({params_keys[index]: parameters[index]})
         # We don't have the phase composition here, so we estimate the driving force.
         # Can happen if one of the composition conditions is unknown or if the phase is
         # stoichiometric and the user did not specify a valid phase composition.
-        single_eqdata = calculate_(species, [current_phase], str_statevar_dict, models, phase_record_factory, pdens=50)
+        single_eqdata = calculate_(species, [current_phase], str_statevar_dict, models, phase_record_factory, pdens=50) ## SOMETHING WEIRD HAPPENS WHEN PDENS IS TOO HIGH!
         df = np.multiply(target_hyperplane_chempots, single_eqdata.X).sum(axis=-1) - single_eqdata.GM
-        driving_force = float(df.max())
+        
+        if np.squeeze(single_eqdata.X).ndim == 1:
+            vertex_comp_estimate = np.squeeze(single_eqdata.X)
+        elif (np.isnan(df).all()):
+            driving_force = np.nan
+            driving_force_gradient = np.full(len(parameters), np.nan)
+            return driving_force, driving_force_gradient
+        else:
+            vertex_comp_estimate = np.squeeze(single_eqdata.X)[np.nanargmax(df),:]
+            
+        counter = 0
+        for comp in species:
+            if v.Species(comp) != v.Species('VA'):
+                if v.X(comp) in cond_dict.keys():
+                    if vertex_comp_estimate[counter] < 5e-6:
+                        vertex_comp_estimate[counter] = 5e-6
+                    elif vertex_comp_estimate[counter] > 1 - 5e-6:
+                        vertex_comp_estimate[counter] = 1 - 5e-6
+                    cond_dict.update({v.X(comp): vertex_comp_estimate[counter]})
+                counter = counter + 1
+                
+        # local_conds = dict(zip(single_eqdata.components, single_eqdata.X))
+            
+        wks = Workspace(database=dbf, components=species, phases=current_phase,  phase_record_factory=phase_record_factory, conditions=cond_dict)
+        constrained_energy = wks.get(IsolatedPhase(current_phase,wks=wks)('GM'))
+        driving_force = np.dot(np.squeeze(target_hyperplane_chempots), vertex_comp_estimate) - constrained_energy
+        ip = IsolatedPhase(current_phase, wks=wks)
+        constrained_energy_gradient = []
+        for key in parameter_dict:
+            constrained_energy_gradient.append(wks.get(ip('GM.'+key)))
+        
+        driving_force_gradient = np.squeeze(np.matmul(vertex_comp_estimate,target_hyperplane_chempots_grads) - constrained_energy_gradient)
+        
     elif vertex.is_disordered:
         # Construct disordered sublattice configuration from composition dict
         # Compute energy
@@ -295,17 +346,21 @@ def driving_force_to_hyperplane(target_hyperplane_chempots: np.ndarray,
         driving_force = np.multiply(target_hyperplane_chempots, single_eqdata.X).sum(axis=-1) - single_eqdata.GM
         driving_force = float(np.squeeze(driving_force))
     else:
-        wks = Workspace(database=dbf, components=species, phases=current_phase, models=models, phase_record_factory=phase_record_factory, conditions=cond_dict, parameters=parameters_dict)
+        wks = Workspace(database=dbf, components=species, phases=current_phase, models=models, phase_record_factory=phase_record_factory, conditions=cond_dict)
         constrained_energy = wks.get(IsolatedPhase(current_phase,wks=wks)('GM'))
         driving_force = np.dot(np.squeeze(target_hyperplane_chempots), vertex.composition) - constrained_energy
-    return driving_force
+        constrained_energy_gradient = []
+        for key in parameter_dict:
+            constrained_energy_gradient.append(wks.get(IsolatedPhase(current_phase,wks=wks)('GM.'+key)))
+        driving_force_gradient = np.squeeze(np.matmul(np.reshape(vertex.composition,(1,-1)),target_hyperplane_chempots_grads) - constrained_energy_gradient)
+    return driving_force, driving_force_gradient
 
 
 def calculate_zpf_driving_forces(zpf_data: Sequence[Dict[str, Any]],
                                  parameters: ArrayLike = None,
                                  approximate_equilibrium: bool = False,
                                  short_circuit: bool = False
-                                 ) -> Tuple[List[List[float]], List[List[float]]]:
+                                 ) -> Tuple[List[List[float]], List[List[float]], List[List[float]]]:
     """
     Calculate error due to phase equilibria data
 
@@ -339,41 +394,47 @@ def calculate_zpf_driving_forces(zpf_data: Sequence[Dict[str, Any]],
         parameters = np.array([])
     driving_forces = []
     weights = []
+    gradients = []
     for data in zpf_data:
         data_driving_forces = []
         data_weights = []
+        data_gradients = []
         weight = data['weight']
         dataset_ref = data['dataset_reference']
         # for the set of phases and corresponding tie-line verticies in equilibrium
         for phase_region in data['phase_regions']:
             # 1. Calculate the average multiphase hyperplane
             eq_str = phase_region.eq_str()
-            target_hyperplane = estimate_hyperplane(phase_region, data['dbf'], parameters, approximate_equilibrium=approximate_equilibrium)
+            target_hyperplane, hyperplane_grads = estimate_hyperplane(phase_region, data['dbf'], data['parameter_dict'], parameters, approximate_equilibrium=approximate_equilibrium)
             if np.any(np.isnan(target_hyperplane)):
                 _log.debug('NaN target hyperplane. Equilibria: (%s), driving force: 0.0, reference: %s.', eq_str, dataset_ref)
                 data_driving_forces.extend([0]*len(phase_region.vertices))
                 data_weights.extend([weight]*len(phase_region.vertices))
+                if len(parameters) == 1:
+                    data_gradients.extend(np.zeros(len(phase_region.vertices)))
+                else:
+                    data_gradients.extend([[0]*len(parameters)]*len(phase_region.vertices))
                 continue
             # 2. Calculate the driving force to that hyperplane for each vertex
             for vertex in phase_region.vertices:
-                driving_force = driving_force_to_hyperplane(target_hyperplane, phase_region, data['dbf'], data['parameter_dict'], vertex, parameters,
-                                                            approximate_equilibrium=approximate_equilibrium,
-                                                            )
+                driving_force, df_grad = driving_force_to_hyperplane(target_hyperplane, hyperplane_grads, phase_region, data['dbf'], data['parameter_dict'], vertex, parameters, approximate_equilibrium=approximate_equilibrium)
                 if np.isinf(driving_force) and short_circuit:
                     _log.debug('Equilibria: (%s), current phase: %s, hyperplane: %s, driving force: %s, reference: %s. Short circuiting.', eq_str, vertex.phase_name, target_hyperplane, driving_force, dataset_ref)
                     return [[np.inf]], [[np.inf]]
                 data_driving_forces.append(driving_force)
                 data_weights.append(weight)
+                data_gradients.append(df_grad)
                 _log.debug('Equilibria: (%s), current phase: %s, hyperplane: %s, driving force: %s, reference: %s', eq_str, vertex.phase_name, target_hyperplane, driving_force, dataset_ref)
         driving_forces.append(data_driving_forces)
         weights.append(data_weights)
-    return driving_forces, weights
+        gradients.append(data_gradients)
+    return driving_forces, weights, gradients
 
 
 def calculate_zpf_error(zpf_data: Sequence[Dict[str, Any]],
                         parameters: np.ndarray = None,
                         data_weight: int = 1.0,
-                        approximate_equilibrium: bool = False) -> float:
+                        approximate_equilibrium: bool = False) -> Tuple[float, List[float]]:
     """
     Calculate the likelihood due to phase equilibria data.
 
@@ -386,16 +447,27 @@ def calculate_zpf_error(zpf_data: Sequence[Dict[str, Any]],
 
     """
     if len(zpf_data) == 0:
-        return 0.0
-    driving_forces, weights = calculate_zpf_driving_forces(zpf_data, parameters, approximate_equilibrium, short_circuit=True)
+        return 0.0, []
+    driving_forces, weights, gradients = calculate_zpf_driving_forces(zpf_data, parameters, approximate_equilibrium, short_circuit=True)
     # Driving forces and weights are 2D ragged arrays with the shape (len(zpf_data), len(zpf_data['values']))
     driving_forces = np.concatenate(driving_forces).T
     weights = np.concatenate(weights)
+    gradients = np.concatenate(gradients)
     if np.any(np.logical_or(np.isinf(driving_forces), np.isnan(driving_forces))):
-        return -np.inf
+        if len(parameters) == 1:
+            return -np.inf, -np.inf
+        else:
+            return -np.inf, np.ones(len(parameters))*(-np.inf)
+        
     log_probabilites = norm.logpdf(driving_forces, loc=0, scale=1000/data_weight/weights)
+    grad_log_probs = -driving_forces*gradients.T/(1000/data_weight/weights)**2
+
+    if grad_log_probs.ndim == 1:
+        grad_log_probs_sum = np.sum(grad_log_probs, axis =0)
+    else:
+        grad_log_probs_sum = np.sum(grad_log_probs, axis =1)
     _log.debug('Data weight: %s, driving forces: %s, weights: %s, probabilities: %s', data_weight, driving_forces, weights, log_probabilites)
-    return np.sum(log_probabilites)
+    return np.sum(log_probabilites), grad_log_probs_sum
 
 
 class ZPFResidual(ResidualFunction):
@@ -426,15 +498,15 @@ def __init__(
         self.zpf_data = get_zpf_data(database, comps, phases, datasets, parameters, model_dict)
 
     def get_residuals(self, parameters: ArrayLike) -> Tuple[List[float], List[float]]:
-        driving_forces, weights = calculate_zpf_driving_forces(self.zpf_data, parameters, short_circuit=True)
+        driving_forces, weights, grads = calculate_zpf_driving_forces(self.zpf_data, parameters, short_circuit=True)
         # Driving forces and weights are 2D ragged arrays with the shape (len(zpf_data), len(zpf_data['values']))
         residuals = np.concatenate(driving_forces).tolist()
         weights = np.concatenate(weights).tolist()
         return residuals, weights
 
-    def get_likelihood(self, parameters) -> float:
-        likelihood = calculate_zpf_error(self.zpf_data, parameters, data_weight=self.weight)
-        return likelihood
+    def get_likelihood(self, parameters) -> Tuple[float, List[float]]:
+        likelihood, gradients = calculate_zpf_error(self.zpf_data, parameters, data_weight=self.weight)
+        return likelihood, gradients
 
 
-residual_function_registry.register(ZPFResidual)
\ No newline at end of file
+residual_function_registry.register(ZPFResidual)
diff --git a/espei/optimizers/opt_mcmc.py b/espei/optimizers/opt_mcmc.py
index 5d0cc865..9b0d2229 100644
--- a/espei/optimizers/opt_mcmc.py
+++ b/espei/optimizers/opt_mcmc.py
@@ -290,7 +290,7 @@ def predict(params, **ctx):
         likelihoods = {}
         for residual_obj in ctx.get("residual_objs", []):
             residual_starttime = time.time()
-            likelihood = residual_obj.get_likelihood(params)
+            likelihood, grads = residual_obj.get_likelihood(params)
             residual_time = time.time() - residual_starttime
             likelihoods[type(residual_obj).__name__] = (likelihood, residual_time)
             lnlike += likelihood
diff --git a/tests/test_error_functions.py b/tests/test_error_functions.py
index 7af9e5c4..71abbefa 100644
--- a/tests/test_error_functions.py
+++ b/tests/test_error_functions.py
@@ -33,7 +33,7 @@ def test_activity_error(datasets_db):
     datasets_db.insert(CU_MG_EXP_ACTIVITY)
 
     dbf = Database(CU_MG_TDB)
-    error = calculate_activity_error(dbf, ['CU','MG','VA'], list(dbf.phases.keys()), datasets_db, {}, {}, {})
+    error, grads = calculate_activity_error(dbf, ['CU','MG','VA'], list(dbf.phases.keys()), datasets_db, parameters=None, phase_models=None, callables=None, data_weight=1.0)
     assert np.isclose(error, -257.41020886970756, rtol=1e-6)
 
 
@@ -47,7 +47,7 @@ def test_activity_residual_function(datasets_db):
     residuals, weights = residual_func.get_residuals(np.asarray([]))
     assert len(residuals) == len(weights)
     assert np.allclose(residuals, [6522.652187085958, -1890.1414208991046, -4793.211215856485, -3018.311675280318, -1062.6724585088668, -2224.814500229084, -2256.9820026771777, -1735.8692674535414, -805.219891012428, 0.0])
-    likelihood = residual_func.get_likelihood(np.asarray([]))
+    likelihood, grads = residual_func.get_likelihood(np.asarray([]))
     assert np.isclose(likelihood, -257.41020886970756, rtol=1e-6)
 
 
@@ -61,14 +61,14 @@ def test_subsystem_activity_probability(datasets_db):
     phases = list(dbf_tern.phases.keys())
 
     # Truth
-    bin_prob = calculate_activity_error(dbf_bin, ['CR','NI','VA'], phases, datasets_db, {}, {}, {})
+    bin_prob, bin_prob_grads = calculate_activity_error(dbf_bin, ['CR','NI','VA'], phases, datasets_db, {}, {}, {})
 
     # Getting binary subsystem data explictly (from binary input)
-    prob = calculate_activity_error(dbf_tern, ['CR','NI','VA'], phases, datasets_db, {}, {}, {})
+    prob, prob_grads = calculate_activity_error(dbf_tern, ['CR','NI','VA'], phases, datasets_db, {}, {}, {})
     assert np.isclose(prob, bin_prob)
 
     # Getting binary subsystem from ternary input
-    prob = calculate_activity_error(dbf_tern, ['CR', 'FE', 'NI', 'VA'], phases, datasets_db, {}, {}, {})
+    prob, prob_grads = calculate_activity_error(dbf_tern, ['CR', 'FE', 'NI', 'VA'], phases, datasets_db, {}, {}, {})
     assert np.isclose(prob, bin_prob)
 
 
@@ -78,11 +78,11 @@ def test_get_thermochemical_data_filters_invalid_sublattice_configurations(datas
     dbf = Database(CU_MG_TDB)
     comps = ["CU", "MG", "VA"]
     phases = ["CUMG2"]
-    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db, symbols_to_fit=[])
+    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db, parameters={}, symbols_to_fit=[])
     print('thermochemical data:', thermochemical_data)
     assert thermochemical_data[0]["calculate_dict"]["values"].shape == (2,)
 
-    error = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
+    error, grads = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
     assert np.isclose(error, -14.28729)
 
 
@@ -96,7 +96,7 @@ def test_fixed_configuration_residual_function(datasets_db):
     residuals, weights = residual_func.get_residuals(np.asarray([]))
     assert len(residuals) == len(weights)
     assert np.allclose(residuals, [-10.0, -100.0])
-    likelihood = residual_func.get_likelihood(np.asarray([]))
+    likelihood, grads = residual_func.get_likelihood(np.asarray([]))
     assert np.isclose(likelihood, -14.28729, rtol=1e-6)
 
 
@@ -119,11 +119,11 @@ def test_get_thermochemical_data_filters_configurations_when_all_configurations_
     dbf = Database(CU_MG_TDB)
     comps = ["CU", "MG", "VA"]
     phases = ["CUMG2"]
-    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db)
+    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db, parameters={})
     print('thermochemical data:', thermochemical_data)
     assert thermochemical_data[0]["calculate_dict"]["values"].shape == (0,)
 
-    error = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
+    error, grads = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
     assert np.isclose(error, 0)
 
 
@@ -134,9 +134,8 @@ def test_non_equilibrium_thermochemical_error_with_multiple_X_points(datasets_db
     dbf = Database(CU_MG_TDB)
     phases = list(dbf.phases.keys())
     comps = ['CU', 'MG', 'VA']
-    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db, symbols_to_fit=[])
-    error = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
-
+    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db, parameters={}, symbols_to_fit=[])
+    error, grads = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
     assert np.isclose(error, -4061.119001241541, rtol=1e-6)
 
 
@@ -147,8 +146,8 @@ def test_non_equilibrium_thermochemical_error_with_multiple_T_points(datasets_db
     dbf = Database(CU_MG_TDB)
     phases = list(dbf.phases.keys())
     comps = ['CU', 'MG', 'VA']
-    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db, symbols_to_fit=[])
-    error = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
+    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db, parameters={}, symbols_to_fit=[])
+    error, grads = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
     assert np.isclose(error,-14.287293263253728, rtol=1e-6)
 
 
@@ -159,11 +158,12 @@ def test_non_equilibrium_thermochemical_error_with_multiple_T_X_points(datasets_
     dbf = Database(CU_MG_TDB)
     phases = list(dbf.phases.keys())
     comps = ['CU', 'MG', 'VA']
-    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db)
     symbols_to_fit = database_symbols_to_fit(dbf)
-    initial_guess = np.array([unpack_piecewise(dbf.symbols[s]) for s in symbols_to_fit])
-    error = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf, parameters=dict(zip(symbols_to_fit, initial_guess)))
-
+    params = dict(zip(symbols_to_fit, [0]*len(symbols_to_fit)))
+    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db, parameters=params)
+    #initial_guess = np.array([unpack_piecewise(dbf.symbols[s]) for s in symbols_to_fit])
+    #error, grads = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf, parameters=dict(zip(symbols_to_fit, initial_guess)))
+    error, grads = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf, parameters=params)
     assert np.isclose(float(error), -3282497.2380024833, rtol=1e-6)
 
 def test_non_equilibrium_thermochemical_error_for_mixing_entropy_error_is_excess_only(datasets_db):
@@ -213,8 +213,8 @@ def test_non_equilibrium_thermochemical_error_for_mixing_entropy_error_is_excess
     # the dataset is excess only
     zero_error_prob = scipy.stats.norm(loc=0, scale=0.2).logpdf(0.0)  # SM weight = 0.2
     # Explicitly pass parameters={} to not try fitting anything
-    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db, symbols_to_fit=[])
-    error = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
+    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db, parameters = {}, symbols_to_fit=[])
+    error, grads = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
     assert np.isclose(error, zero_error_prob, atol=1e-6)
 
 
@@ -264,8 +264,8 @@ def test_non_equilibrium_thermochemical_error_for_of_enthalpy_mixing(datasets_db
     # the dataset is excess only
     zero_error_prob = scipy.stats.norm(loc=0, scale=500.0).logpdf(0.0)  # HM weight = 500
     # Explicitly pass parameters={} to not try fitting anything
-    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db, symbols_to_fit=[])
-    error = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
+    thermochemical_data = get_thermochemical_data(dbf, comps, phases, datasets_db, parameters = {},symbols_to_fit=[])
+    error, grads = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
     assert np.isclose(error, zero_error_prob, atol=1e-6)
 
 
@@ -279,17 +279,17 @@ def test_subsystem_non_equilibrium_thermochemcial_probability(datasets_db):
     phases = list(dbf_tern.phases.keys())
 
     # Truth
-    thermochemical_data = get_thermochemical_data(dbf_bin, ['CR', 'NI', 'VA'], phases, datasets_db)
-    bin_prob = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf_bin)
+    thermochemical_data = get_thermochemical_data(dbf_bin, ['CR', 'NI', 'VA'], phases, datasets_db,parameters={})
+    bin_prob, grads = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf_bin)
 
     # Getting binary subsystem data explictly (from binary input)
-    thermochemical_data = get_thermochemical_data(dbf_tern, ['CR', 'NI', 'VA'], phases, datasets_db)
-    prob = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf_tern)
+    thermochemical_data = get_thermochemical_data(dbf_tern, ['CR', 'NI', 'VA'], phases, datasets_db, parameters = {})
+    prob, grads = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf_tern)
     assert np.isclose(prob, bin_prob)
 
     # Getting binary subsystem from ternary input
-    thermochemical_data = get_thermochemical_data(dbf_tern, ['CR', 'FE', 'NI', 'VA'], phases, datasets_db)
-    prob = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf_tern)
+    thermochemical_data = get_thermochemical_data(dbf_tern, ['CR', 'FE', 'NI', 'VA'], phases, datasets_db, parameters = {})
+    prob, grads = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf_tern)
     assert np.isclose(prob, bin_prob)
 
 
@@ -305,7 +305,7 @@ def test_zpf_error_zero(datasets_db):
     zero_error_prob = 2 * scipy.stats.norm(loc=0, scale=1000.0).logpdf(0.0)
 
     zpf_data = get_zpf_data(dbf, comps, phases, datasets_db, {})
-    error = calculate_zpf_error(zpf_data, np.array([]))
+    error, grads = calculate_zpf_error(zpf_data, np.array([]))
     assert np.isclose(error, zero_error_prob, rtol=1e-6)
 
 
@@ -319,7 +319,7 @@ def test_zpf_residual_function(datasets_db):
     residuals, weights = residual_func.get_residuals(np.asarray([]))
     assert len(residuals) == len(weights)
     assert np.allclose(residuals, [0.0, 0.0], atol=1e-3)  # looser tolerance due to numerical instabilities
-    likelihood = residual_func.get_likelihood(np.asarray([]))
+    likelihood, grads = residual_func.get_likelihood(np.asarray([]))
     # ZPF weight = 1 kJ and there are two points in the tieline
     zero_error_prob = np.sum(scipy.stats.norm(loc=0, scale=1000.0).logpdf([0.0, 0.0]))
     assert np.isclose(likelihood, zero_error_prob, rtol=1e-6)
@@ -336,16 +336,16 @@ def test_subsystem_zpf_probability(datasets_db):
 
     # Truth
     zpf_data = get_zpf_data(dbf_bin, ['CR', 'NI', 'VA'], phases, datasets_db, {})
-    bin_prob = calculate_zpf_error(zpf_data, np.array([]))
+    bin_prob, grads = calculate_zpf_error(zpf_data, np.array([]))
 
     # Getting binary subsystem data explictly (from binary input)
     zpf_data = get_zpf_data(dbf_tern, ['CR', 'NI', 'VA'], phases, datasets_db, {})
-    prob = calculate_zpf_error(zpf_data, np.array([]))
+    prob, grads = calculate_zpf_error(zpf_data, np.array([]))
     assert np.isclose(prob, bin_prob)
 
     # Getting binary subsystem from ternary input
     zpf_data = get_zpf_data(dbf_tern, ['CR', 'FE', 'NI', 'VA'], phases, datasets_db, {})
-    prob = calculate_zpf_error(zpf_data, np.array([]))
+    prob, grads = calculate_zpf_error(zpf_data, np.array([]))
     assert np.isclose(prob, bin_prob)
 
 
@@ -367,7 +367,7 @@ def test_zpf_error_species(datasets_db):
     zero_error_probability = 2 * scipy.stats.norm(loc=0, scale=1000.0).logpdf(0.0)
 
     zpf_data = get_zpf_data(dbf, comps, phases, datasets_db, {})
-    likelihood = calculate_zpf_error(zpf_data)
+    likelihood, grads = calculate_zpf_error(zpf_data)
     assert np.isclose(likelihood, zero_error_probability)
 
 
@@ -384,7 +384,7 @@ def test_zpf_error_equilibrium_failure(datasets_db):
     zpf_data = get_zpf_data(dbf, comps, phases, datasets_db, {})
 
     with mock.patch('espei.error_functions.zpf_error.estimate_hyperplane', return_value=np.array([np.nan, np.nan])):
-        likelihood = calculate_zpf_error(zpf_data)
+        likelihood, grads = calculate_zpf_error(zpf_data)
         assert np.isclose(likelihood, zero_error_probability, rtol=1e-6)
 
 
@@ -400,7 +400,7 @@ def test_zpf_error_works_for_stoichiometric_cmpd_tielines(datasets_db):
     zero_error_probability = 2 * scipy.stats.norm(loc=0, scale=1000.0).logpdf(0.0)
 
     zpf_data = get_zpf_data(dbf, comps, phases, datasets_db, {})
-    likelihood = calculate_zpf_error(zpf_data)
+    likelihood, grads = calculate_zpf_error(zpf_data)
     assert np.isclose(likelihood, zero_error_probability, rtol=1e-6)
 
 
@@ -412,8 +412,8 @@ def test_non_equilibrium_thermochemcial_species(datasets_db):
     dbf = Database(LI_SN_TDB)
     phases = ['LIQUID']
 
-    thermochemical_data = get_thermochemical_data(dbf, ['LI', 'SN'], phases, datasets_db)
-    prob = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
+    thermochemical_data = get_thermochemical_data(dbf, ['LI', 'SN'], phases, datasets_db, parameters={})
+    prob, grads = calculate_non_equilibrium_thermochemical_probability(thermochemical_data, dbf)
     # Near zero error and non-zero error
     assert np.isclose(prob, (-7.13354663 + -22.43585011))
 
@@ -429,7 +429,7 @@ def test_equilibrium_thermochemcial_error_species(datasets_db):
     eqdata = get_equilibrium_thermochemical_data(dbf, ['LI', 'SN'], phases, datasets_db)
     # Thermo-Calc
     truth_values = np.array([0.0, -28133.588, -40049.995, 0.0])
-    residuals, weights = calc_prop_differences(eqdata[0], np.array([]))
+    residuals, weights, grads = calc_prop_differences(eqdata[0], np.array([]))
     assert np.all(np.isclose(residuals, truth_values, atol=3e-5))
 
 
@@ -443,7 +443,7 @@ def test_equilibrium_thermochemical_error_unsupported_property(datasets_db):
     phases = list(dbf.phases.keys())
 
     eqdata = get_equilibrium_thermochemical_data(dbf, ['CR', 'NI'], phases, datasets_db)
-    residuals, weights = calc_prop_differences(eqdata[0], np.array([]))
+    residuals, weights, grads = calc_prop_differences(eqdata[0], np.array([]))
     assert np.all(np.isclose(residuals, EXPECTED_VALUES, atol=1e-3))
 
 
@@ -457,7 +457,7 @@ def test_equilibrium_property_residual_function(datasets_db):
     assert len(residuals) == len(weights)
     assert np.allclose(residuals, [374.6625, 0.0, 0.0])
     # Regression test "truth" values - got values by running
-    likelihood = residual_func.get_likelihood(np.asarray([]))
+    likelihood, grads = residual_func.get_likelihood(np.asarray([]))
     assert np.isclose(likelihood, -70188.75126872442, rtol=1e-6)
 
 
@@ -470,17 +470,17 @@ def test_equilibrium_thermochemical_error_computes_correct_probability(datasets_
     # Test that errors update in response to changing parameters
     # no parameters
     eqdata = get_equilibrium_thermochemical_data(dbf, ['CU', 'MG'], phases, datasets_db)
-    errors, weights = calc_prop_differences(eqdata[0], np.array([]))
+    errors, weights, grads = calc_prop_differences(eqdata[0], np.array([]))
     expected_vals = [-31626.6*0.5*0.5]
     assert np.all(np.isclose(errors, expected_vals))
-
+    
     # VV0017 (LIQUID, L0)
     eqdata = get_equilibrium_thermochemical_data(dbf, ['CU', 'MG'], phases, datasets_db, parameters={'VV0017': -31626.6})
     # unchanged, should be the same as before
-    errors, weights = calc_prop_differences(eqdata[0], np.array([-31626.6]))
+    errors, weights, grads = calc_prop_differences(eqdata[0], np.array([-31626.6]))
     assert np.all(np.isclose(errors, [-31626.6*0.5*0.5]))
     # change to -40000
-    errors, weights = calc_prop_differences(eqdata[0], np.array([-40000], np.float64))
+    errors, weights, grads = calc_prop_differences(eqdata[0], np.array([-40000], np.float64))
     assert np.all(np.isclose(errors, [-40000*0.5*0.5]))
 
 
@@ -495,17 +495,17 @@ def test_driving_force_miscibility_gap(datasets_db):
 
     # Ideal solution case
     params = np.array([0.0])
-    prob = calculate_zpf_error(zpf_data, parameters=params)
+    prob, grads = calculate_zpf_error(zpf_data, parameters=params)
     assert np.isclose(prob, zero_error_prob)
 
     # Negative interaction case
     params = np.array([-10000.0])
-    prob = calculate_zpf_error(zpf_data, parameters=params)
+    prob, grads = calculate_zpf_error(zpf_data, parameters=params)
     assert np.isclose(prob, zero_error_prob)
 
     # Miscibility gap case
     params = np.array([10000.0])
-    prob = calculate_zpf_error(zpf_data, parameters=params)
+    prob, grads = calculate_zpf_error(zpf_data, parameters=params)
     # Remember these are log probabilities, so more negative means smaller probability and larger error
     assert prob < zero_error_prob
 
@@ -536,7 +536,7 @@ def test_setting_up_context_with_custom_models(datasets_db):
     with pytest.raises(ModelTestException):
         ctx = setup_context(dbf, datasets_db, phase_models=phase_models)
 
-
+#this one requires work on the activity error function
 def test_zpf_context_is_pickleable(datasets_db):
     """Test that the context for ZPF data is pickleable"""
     datasets_db.insert(CU_MG_DATASET_ZPF_ZERO_ERROR)
@@ -619,7 +619,7 @@ def test_zpf_error_for_prescribed_hyperplane_composition(datasets_db):
     datasets_db.insert(A_B_DATASET_ALPHA)
     dbf = Database(A_B_REGULAR_SOLUTION_TDB)  # Ideal solution case by default
     zpf_data = get_zpf_data(dbf, ["A", "B"], ["ALPHA"], datasets_db, {})
-    driving_forces, weights = calculate_zpf_driving_forces(zpf_data)
+    driving_forces, weights, grads = calculate_zpf_driving_forces(zpf_data)
     flat_driving_forces = np.asarray(driving_forces).flatten()
     assert len(flat_driving_forces) == 1
     assert np.isclose(flat_driving_forces[0], 0.0)
@@ -630,7 +630,7 @@ def test_zpf_error_hyperplane_with_null_phases(datasets_db):
     datasets_db.insert(CU_MG_DATASET_ZPF_HYPERPLANE_TWOPHASE)
     dbf = Database(CU_MG_TDB)  # Ideal solution case by default
     zpf_data = get_zpf_data(dbf, ["CU", "MG"], list(dbf.phases.keys()), datasets_db, {})
-    driving_forces, weights = calculate_zpf_driving_forces(zpf_data)
+    driving_forces, weights, grads = calculate_zpf_driving_forces(zpf_data)
     flat_driving_forces = np.asarray(driving_forces).flatten()
     assert len(flat_driving_forces) == 2  # One for each vertex, HCP_A3 and CUMG2
     assert np.allclose(flat_driving_forces, [-18.05883506, -780.50836135])