diff --git a/src/phreeqcpp/NameDouble.cxx b/src/phreeqcpp/NameDouble.cxx
index d9ceb6eb..ce485232 100644
--- a/src/phreeqcpp/NameDouble.cxx
+++ b/src/phreeqcpp/NameDouble.cxx
@@ -226,7 +226,7 @@ cxxNameDouble::add_intensive(const cxxNameDouble & addee, LDBLE f1,
 // Sums two name doubles, this*f1 + f2*nd2
 //
 {
-	assert(f1 >= 0 && f2 >= 0);
+	//assert(f1 >= 0 && f2 >= 0);
 	for (cxxNameDouble::const_iterator it = addee.begin(); it != addee.end();
 		 it++)
 	{
@@ -248,7 +248,7 @@ cxxNameDouble::add_log_activities(const cxxNameDouble & addee, LDBLE f1,
 // Sums two name doubles, this*f1 + f2*nd2, assuming log values
 //
 {
-	assert(f1 >= 0 && f2 >= 0);
+	//assert(f1 >= 0 && f2 >= 0);
 	for (cxxNameDouble::const_iterator it = addee.begin(); it != addee.end();
 		 it++)
 	{
diff --git a/src/phreeqcpp/Solution.cxx b/src/phreeqcpp/Solution.cxx
index b2f4187a..50a88427 100644
--- a/src/phreeqcpp/Solution.cxx
+++ b/src/phreeqcpp/Solution.cxx
@@ -129,24 +129,33 @@ cxxSolution::cxxSolution(std::map < int, cxxSolution > &solutions,
 		else
 			this->potV = cxxsoln_ptr1->potV;
 	}
+	//
+	// Sort to enable positive mixes first
+	const std::map < int, LDBLE >& mixcomps = mix.Get_mixComps();
+	std::set<std::pair<LDBLE, int> >s;
+	for (std::map < int, LDBLE >::const_iterator it = mixcomps.begin();
+		it != mixcomps.end(); it++)
+	{
+		std::pair<LDBLE, int> p(it->second, it->first);
+		s.insert(p); 
+	}
 //
 //   Mix solutions
 //
-	const std::map < int, LDBLE >&mixcomps = mix.Get_mixComps();
-	std::map < int, LDBLE >::const_iterator it;
-	for (it = mixcomps.begin(); it != mixcomps.end(); it++)
+	std::set < std::pair< LDBLE, int > >::const_reverse_iterator rit = s.rbegin();
+	for (rit = s.rbegin(); rit != s.rend(); rit++)
 	{
-		sol = solutions.find(it->first);
+		sol = solutions.find(rit->second);
 		if (sol == solutions.end())
 		{
 			std::ostringstream msg;
-			msg << "Solution " << it->first << " not found in mix_cxxSolutions.";
+			msg << "Solution " << rit->second << " not found in mix_cxxSolutions.";
 			error_msg(msg.str(), CONTINUE);
 		}
 		else
 		{
 			cxxsoln_ptr1 = &(sol->second);
-			this->add(*cxxsoln_ptr1, it->second);
+			this->add(*cxxsoln_ptr1, rit->first);
 		}
 	}
 
@@ -1418,8 +1427,25 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive)
 		return;
 	LDBLE ext1 = this->mass_water;
 	LDBLE ext2 = addee.mass_water * extensive;
-	LDBLE f1 = ext1 / (ext1 + ext2);
-	LDBLE f2 = ext2 / (ext1 + ext2);
+	this->mass_water += addee.mass_water * extensive;
+	if (this->mass_water <= 0.0)
+	{
+		std::ostringstream msg;
+		msg << "Negative mass of water when mixing solutions.";
+		error_msg(msg.str(), STOP);
+	}
+	LDBLE fconc = 0.0, f1 = 0.0, f2 = 0.0;
+	if (extensive > 0.0)
+	{
+
+		f1 = ext1 / (ext1 + ext2);
+		f2 = ext2 / (ext1 + ext2);
+	}
+	else
+	{
+		f1 = 1.0;
+		f2 = 0.0;
+	}
 	this->tc = f1 * this->tc + f2 * addee.tc;
 	this->ph = f1 * this->ph + f2 * addee.ph;
 	this->pe = f1 * this->pe + f2 * addee.pe;
@@ -1433,7 +1459,6 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive)
 	this->viscos_0 = f1 * this->viscos_0 + f2 * addee.viscos_0;
 	this->patm = f1 * this->patm + f2 * addee.patm;
 	// this->potV = f1 * this->potV + f2 * addee.potV; // appt
-	this->mass_water += addee.mass_water * extensive;
 	this->soln_vol += addee.soln_vol * extensive;
 	this->total_alkalinity += addee.total_alkalinity * extensive;
 	this->totals.add_extensive(addee.totals, extensive);
diff --git a/src/phreeqcpp/global_structures.h b/src/phreeqcpp/global_structures.h
index 3e7cc701..bedaee8f 100644
--- a/src/phreeqcpp/global_structures.h
+++ b/src/phreeqcpp/global_structures.h
@@ -1586,7 +1586,7 @@ public:
 typedef enum
 { TYPE_B0, TYPE_B1, TYPE_B2, TYPE_C0, TYPE_THETA, TYPE_LAMBDA, TYPE_ZETA,
   TYPE_PSI, TYPE_ETHETA, TYPE_ALPHAS, TYPE_MU, TYPE_ETA, TYPE_Other,
-  TYPE_SIT_EPSILON, TYPE_SIT_EPSILON_MU, TYPE_APHI
+  TYPE_SIT_EPSILON, TYPE_SIT_EPSILON_MU, TYPE_APHI, TYPE_SIT_EPSILON2
 } pitz_param_type;
 class pitz_param
 {
@@ -1625,6 +1625,7 @@ public:
 		LDBLE eps;
 		LDBLE eps1;
 		LDBLE aphi;
+		LDBLE eps2;
 	} U;
 	LDBLE a[6];
 	LDBLE alpha;
diff --git a/src/phreeqcpp/pitzer.cpp b/src/phreeqcpp/pitzer.cpp
index 9282b576..1f32a085 100644
--- a/src/phreeqcpp/pitzer.cpp
+++ b/src/phreeqcpp/pitzer.cpp
@@ -2690,7 +2690,7 @@ pitzer_make_lists(void)
 		/*
 		TYPE_B0, TYPE_B1, TYPE_B2, TYPE_C0, TYPE_THETA, TYPE_LAMBDA, TYPE_ZETA,
 		TYPE_PSI, TYPE_ETHETA, TYPE_ALPHAS, TYPE_MU, TYPE_ETA, TYPE_Other,
-		TYPE_SIT_EPSILON, TYPE_SIT_EPSILON_MU
+		TYPE_SIT_EPSILON, TYPE_SIT_EPSILON_MU, TYPE_APHI, TYPE_SIT_EPSILON2
 		*/
 		int i0 = pitz_params[i]->ispec[0];
 		int i1 = pitz_params[i]->ispec[1];
diff --git a/src/phreeqcpp/sit.cpp b/src/phreeqcpp/sit.cpp
index 4f3feb53..2d061f01 100644
--- a/src/phreeqcpp/sit.cpp
+++ b/src/phreeqcpp/sit.cpp
@@ -188,9 +188,10 @@ read_sit(void)
   const char* next_char;
   const char *opt_list[] = {
     "epsilon",					/* 0 */
-    "epsilon1"					/* 1 */
+    "epsilon1",					/* 1 */
+	"epsilon2"					/* 2 */
   };
-  int count_opt_list = 2;
+  int count_opt_list = 3;
   /*
    *   Read lines
    */
@@ -237,6 +238,11 @@ read_sit(void)
       n = 2;
       opt_save = OPTION_DEFAULT;
       break;
+	case 2:				/* epsilon2 */
+		pzp_type = TYPE_SIT_EPSILON2;
+		n = 2;
+		opt_save = OPTION_DEFAULT;
+		break;
     }
     if (return_value == EOF || return_value == KEYWORD)
       break;
@@ -275,6 +281,9 @@ calc_sit_param(class pitz_param *pz_ptr, LDBLE TK, LDBLE TR)
 	case TYPE_SIT_EPSILON_MU:
 		pz_ptr->U.eps1 = param;
 		break;
+	case TYPE_SIT_EPSILON2:
+		pz_ptr->U.eps2 = param;
+		break;
 	case TYPE_Other:
 	default:
 		error_msg("Should not be TYPE_Other in function calc_sit_param",
@@ -399,6 +408,7 @@ sit(void)
 	 *  Sums for sit_LGAMMA, and OSMOT
 	 *  epsilons are tabulated for log10 gamma (not ln gamma)
 	 */
+	LDBLE logmu = log10(I);
 	for (size_t j = 0; j < param_list.size(); j++)
 	{
 		int i = param_list[j];
@@ -425,7 +435,7 @@ sit(void)
 		case TYPE_SIT_EPSILON_MU:
 			sit_LGAMMA[i0] += sit_M[i1] * I * param;
 			sit_LGAMMA[i1] += sit_M[i0] * I * param;
-			OSMOT += sit_M[i0] * sit_M[i1] * param;
+			//OSMOT += sit_M[i0] * sit_M[i1] * param;
 			if (z0 == 0.0 && z1 == 0.0)
 			{
 				OSMOT += sit_M[i0] * sit_M[i1] * param * I / 2.0;
@@ -435,6 +445,19 @@ sit(void)
 				OSMOT += sit_M[i0] * sit_M[i1] * param * I;
 			}
 			break;
+		case TYPE_SIT_EPSILON2:
+			sit_LGAMMA[i0] += sit_M[i1] * logmu * param;
+			sit_LGAMMA[i1] += sit_M[i0] * logmu * param;
+			//OSMOT += sit_M[i0] * sit_M[i1] * param;
+			if (z0 == 0.0 && z1 == 0.0)
+			{
+				OSMOT += sit_M[i0] * sit_M[i1] * param * logmu / 2.0;
+			}
+			else
+			{
+				OSMOT += sit_M[i0] * sit_M[i1] * param * logmu;
+			}
+			break;
 		default:
 		case TYPE_Other:
 			error_msg("TYPE_Other in pitz_param list.", STOP);
diff --git a/src/phreeqcpp/step.cpp b/src/phreeqcpp/step.cpp
index 131aec3a..e7eb271e 100644
--- a/src/phreeqcpp/step.cpp
+++ b/src/phreeqcpp/step.cpp
@@ -621,70 +621,72 @@ add_mix(cxxMix *mix_ptr)
  *   calls add_solution to accumulate all data in master->totals
  *   and other variables.
  */
-	LDBLE sum_fractions, intensive, extensive;
-	LDBLE sum_fractions_water=0, sum_positive_water=0, intensive_water=0, extensive_water=0;
-	cxxSolution *solution_ptr;
-	int count_positive;
-	LDBLE sum_positive;
+	//LDBLE sum_fractions, intensive, extensive;
+	//LDBLE sum_fractions_water=0, sum_positive_water=0, intensive_water=0, extensive_water=0;
+	//cxxSolution *solution_ptr;
+	//int count_positive;
+	//LDBLE sum_positive;
 
 	if (mix_ptr == NULL)
 		return (OK);
 	if (mix_ptr->Get_mixComps().size() == 0)
 		return (OK);
-	sum_fractions = 0.0;
-	sum_positive = 0.0;
-	count_positive = 0;
-	std::map<int, LDBLE>::const_iterator it;
-	for (it = mix_ptr->Get_mixComps().begin(); it != mix_ptr->Get_mixComps().end(); it++)
-	{
-		solution_ptr = Utilities::Rxn_find(Rxn_solution_map, it->first);
-		if (solution_ptr == NULL)
-		{
-			error_string = sformatf("Mix solution not found, %d.",
-				it->first);
-			error_msg(error_string, CONTINUE);
-			input_error++;
-			continue;
-		}
-		sum_fractions += it->second;
-		sum_fractions_water += it->second * solution_ptr->Get_mass_water();
-		if (it->second > 0)
-		{
-			sum_positive += it->second;
-			sum_positive_water += it->second * solution_ptr->Get_mass_water();
-			count_positive++;
-		}
-	}
-	for (it = mix_ptr->Get_mixComps().begin(); it != mix_ptr->Get_mixComps().end(); it++)
-	{
-		solution_ptr = Utilities::Rxn_find(Rxn_solution_map, it->first);
-		if (solution_ptr == NULL)
-		{
-			error_string = sformatf( "Mix solution not found, %d.",
-						it->first);
-			error_msg(error_string, CONTINUE);
-			input_error++;
-			continue;
-		}
-		extensive = it->second;
-		extensive_water = it->second * solution_ptr->Get_mass_water();
-		intensive = extensive / sum_fractions;
-		intensive_water = extensive_water / sum_fractions_water;
-		
-		if (count_positive < (int) mix_ptr->Get_mixComps().size())
-		{
-			if (it->second > 0)
-			{
-				intensive = extensive / sum_positive;
-				intensive_water = extensive_water / sum_positive_water;
-			}
-			else
-			{
-				intensive = 0;
-			}
-		}
-		add_solution(solution_ptr, extensive, intensive_water);
-	}
+	//sum_fractions = 0.0;
+	//sum_positive = 0.0;
+	//count_positive = 0;
+	//std::map<int, LDBLE>::const_iterator it;
+	//for (it = mix_ptr->Get_mixComps().begin(); it != mix_ptr->Get_mixComps().end(); it++)
+	//{
+	//	solution_ptr = Utilities::Rxn_find(Rxn_solution_map, it->first);
+	//	if (solution_ptr == NULL)
+	//	{
+	//		error_string = sformatf("Mix solution not found, %d.",
+	//			it->first);
+	//		error_msg(error_string, CONTINUE);
+	//		input_error++;
+	//		continue;
+	//	}
+	//	sum_fractions += it->second;
+	//	sum_fractions_water += it->second * solution_ptr->Get_mass_water();
+	//	if (it->second > 0)
+	//	{
+	//		sum_positive += it->second;
+	//		sum_positive_water += it->second * solution_ptr->Get_mass_water();
+	//		count_positive++;
+	//	}
+	//}
+	//for (it = mix_ptr->Get_mixComps().begin(); it != mix_ptr->Get_mixComps().end(); it++)
+	//{
+	//	solution_ptr = Utilities::Rxn_find(Rxn_solution_map, it->first);
+	//	if (solution_ptr == NULL)
+	//	{
+	//		error_string = sformatf( "Mix solution not found, %d.",
+	//					it->first);
+	//		error_msg(error_string, CONTINUE);
+	//		input_error++;
+	//		continue;
+	//	}
+	//	extensive = it->second;
+	//	extensive_water = it->second * solution_ptr->Get_mass_water();
+	//	intensive = extensive / sum_fractions;
+	//	intensive_water = extensive_water / sum_fractions_water;
+	//	
+	//	if (count_positive < (int) mix_ptr->Get_mixComps().size())
+	//	{
+	//		if (it->second > 0)
+	//		{
+	//			intensive = extensive / sum_positive;
+	//			intensive_water = extensive_water / sum_positive_water;
+	//		}
+	//		else
+	//		{
+	//			intensive = 0;
+	//		}
+	//	}
+	//	add_solution(solution_ptr, extensive, intensive_water);
+	//}
+	cxxSolution s(Rxn_solution_map, *mix_ptr, 1);
+	add_solution(&s, 1.0, 1.0);
 	return (OK);
 }
 #ifdef SKIP_ERROR