diff --git a/modules/ct_equil/equilibrate_T.m b/modules/ct_equil/equilibrate_T.m
index aa65789c..dc3af873 100644
--- a/modules/ct_equil/equilibrate_T.m
+++ b/modules/ct_equil/equilibrate_T.m
@@ -48,7 +48,7 @@
     % Compute property matrix of the species at chemical equilibrium
     % NOTE: If the ind variable is removed from the inputs, the set_species 
     % routine will completely fill the properties matrix
-    M0 = set_species(self, self.S.LS, N(:, 1), TP, ind);
+    M0 = set_species(self, self.S.LS, N(:, 1), TP, [ind, self.S.ind_frozen]);
     % Compute properties of final mixture
     mix2 = set_properties(self, mix1, M0, pP, TP, STOP, STOP_ions);
 end
diff --git a/modules/ct_equil/equilibrium_gibbs.m b/modules/ct_equil/equilibrium_gibbs.m
index 19ed2ddc..f0155040 100644
--- a/modules/ct_equil/equilibrium_gibbs.m
+++ b/modules/ct_equil/equilibrium_gibbs.m
@@ -40,7 +40,7 @@
     R0TP = self.C.R0 * TP; % [J/(mol)]
     
     % Initialization
-    NatomE = mix1.NatomE;
+    NatomE = mix1.NatomE_react;
     max_NatomE = max(NatomE);
     NP = 0.1;
     SIZE = -log(self.TN.tolN);
diff --git a/modules/ct_equil/equilibrium_gibbs_eos.m b/modules/ct_equil/equilibrium_gibbs_eos.m
index edc259ed..b8aa2dca 100644
--- a/modules/ct_equil/equilibrium_gibbs_eos.m
+++ b/modules/ct_equil/equilibrium_gibbs_eos.m
@@ -41,7 +41,7 @@
     A0 = C.A0.value;
     R0TP = C.R0 * TP; % [J/mol]
     % Initialization
-    NatomE = mix1.NatomE_react';
+    NatomE = mix1.NatomE_react;
     max_NatomE = max(NatomE);
     NP = 0.1;
     SIZE = -log(TN.tolN);
diff --git a/modules/ct_equil/equilibrium_gibbs_large.m b/modules/ct_equil/equilibrium_gibbs_large.m
index 6bafebb5..63fb37a3 100644
--- a/modules/ct_equil/equilibrium_gibbs_large.m
+++ b/modules/ct_equil/equilibrium_gibbs_large.m
@@ -38,7 +38,7 @@
     R0TP = self.C.R0 * TP; % [J/mol]
 
     % Initialization
-    NatomE = mix1.NatomE_react';
+    NatomE = mix1.NatomE_react;
     max_NatomE = max(NatomE);
     NP = 0.1;
     SIZE = -log(self.TN.tolN);
diff --git a/modules/ct_equil/equilibrium_helmholtz.m b/modules/ct_equil/equilibrium_helmholtz.m
index 641c0acd..da548343 100644
--- a/modules/ct_equil/equilibrium_helmholtz.m
+++ b/modules/ct_equil/equilibrium_helmholtz.m
@@ -41,7 +41,7 @@
     R0TP = self.C.R0 * TP; % [J/(mol)]
 
     % Initialization
-    NatomE = mix1.NatomE;
+    NatomE = mix1.NatomE_react;
     max_NatomE = max(NatomE);
     NP = 0.1;
     SIZE = -log(self.TN.tolN);
diff --git a/modules/ct_equil/equilibrium_helmholtz_eos.m b/modules/ct_equil/equilibrium_helmholtz_eos.m
index 7220bd62..27e243fd 100644
--- a/modules/ct_equil/equilibrium_helmholtz_eos.m
+++ b/modules/ct_equil/equilibrium_helmholtz_eos.m
@@ -41,7 +41,7 @@
     A0 = C.A0.value;
     R0TP = C.R0 * TP; % [J/mol]
     % Initialization
-    NatomE = mix1.NatomE_react';
+    NatomE = mix1.NatomE_react;
     max_NatomE = max(NatomE);
     NP = 0.1;
     SIZE = -log(TN.tolN);
diff --git a/modules/ct_equil/equilibrium_helmholtz_large.m b/modules/ct_equil/equilibrium_helmholtz_large.m
index b3170621..af833430 100644
--- a/modules/ct_equil/equilibrium_helmholtz_large.m
+++ b/modules/ct_equil/equilibrium_helmholtz_large.m
@@ -38,7 +38,7 @@
     R0TP = self.C.R0 * TP; % [J/mol]
 
     % Initialization
-    NatomE = mix1.NatomE_react';
+    NatomE = mix1.NatomE_react;
     max_NatomE = max(NatomE);
     NP = 0.1;
     SIZE = -log(self.TN.tolN);
diff --git a/utils/compute_properties.m b/utils/compute_properties.m
index 023fdf00..c475fe02 100644
--- a/utils/compute_properties.m
+++ b/utils/compute_properties.m
@@ -17,8 +17,8 @@
     % Definitions
     R0 = self.C.R0; % [J/(K mol)] Universal gas constant
     % Initialization
-    mix.error_moles = [];
-    mix.error_moles_ions = [];
+    mix.error_moles = 0;
+    mix.error_moles_ions = 0;
     % Inputs
     mix.p = p; % [bar]
     mix.T = T; % [K]