diff --git a/phreeqcpp/Solution.cxx b/phreeqcpp/Solution.cxx
index d29a788f..85055223 100644
--- a/phreeqcpp/Solution.cxx
+++ b/phreeqcpp/Solution.cxx
@@ -114,6 +114,17 @@ cxxSolution::cxxSolution(std::map < int, cxxSolution > &solutions,
 	this->n_user = this->n_user_end = l_n_user;
 	this->new_def = false;
 	this->ah2o = 0;
+	// potV is an external variable, imposed in a given solution, not mixed.
+	std::map < int, cxxSolution >::const_iterator sol = solutions.find(mix.Get_n_user());
+	const cxxSolution *cxxsoln_ptr1;
+	if (sol != solutions.end())
+	{
+		cxxsoln_ptr1 = &(sol->second);
+		if (cxxsoln_ptr1->new_def)
+			this->potV = 0.0;
+		else
+			this->potV = cxxsoln_ptr1->potV;
+	}
 //
 //   Mix solutions
 //
@@ -121,8 +132,7 @@ cxxSolution::cxxSolution(std::map < int, cxxSolution > &solutions,
 	std::map < int, LDBLE >::const_iterator it;
 	for (it = mixcomps.begin(); it != mixcomps.end(); it++)
 	{
-		std::map < int, cxxSolution >::const_iterator sol =
-			solutions.find(it->first);
+		sol = solutions.find(it->first);
 		if (sol == solutions.end())
 		{
 			std::ostringstream msg;
@@ -131,7 +141,7 @@ cxxSolution::cxxSolution(std::map < int, cxxSolution > &solutions,
 		}
 		else
 		{
-			const cxxSolution *cxxsoln_ptr1 = &(sol->second);
+			cxxsoln_ptr1 = &(sol->second);
 			this->add(*cxxsoln_ptr1, it->second);
 		}
 	}
@@ -1388,7 +1398,7 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive)
 	this->cb += addee.cb * extensive;
 	this->density = f1 * this->density + f2 * addee.density;
 	this->patm = f1 * this->patm + f2 * addee.patm;
-	this->potV = f1 * this->potV + f2 * addee.potV;
+	// 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;
diff --git a/phreeqcpp/class_main.cpp b/phreeqcpp/class_main.cpp
index 9e6a0d0c..ee43e385 100644
--- a/phreeqcpp/class_main.cpp
+++ b/phreeqcpp/class_main.cpp
@@ -147,8 +147,8 @@ main_method(int argc, char *argv[])
 	}
 	Phreeqc MyCopy = *this;
 	this->clean_up();
-	//this->init();
-	//this->initialize();
+	// this->init();
+	// this->initialize();
 /*
  *   Read input data for simulation
  */
@@ -167,7 +167,7 @@ main_method(int argc, char *argv[])
  *   Display successful status
  */
 	pr.headings = TRUE;
-	//errors = do_status();
+	// errors = do_status();
 	if (errors != 0)
 	{
 		return errors;
@@ -292,7 +292,7 @@ write_banner(void)
 
 	/* version */
 #ifdef NPP
-	len = sprintf(buffer, "* PHREEQC-%s *", "3.6.5");
+	len = sprintf(buffer, "* PHREEQC-%s *", "3.7.1");
 #else
 	len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@");
 #endif
@@ -316,7 +316,7 @@ write_banner(void)
 
 	/* date */
 #ifdef NPP
-	len = sprintf(buffer, "%s", "February 24, 2021");
+	len = sprintf(buffer, "%s", "July 5, 2021");
 #else
 	len = sprintf(buffer, "%s", "@VER_DATE@");
 #endif
@@ -329,6 +329,7 @@ write_banner(void)
 
 	return 0;
 }
+
 /* ---------------------------------------------------------------------- */
 int Phreeqc::
 process_file_names(int argc, char *argv[], std::istream **db_cookie,
@@ -499,7 +500,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
 		output_msg(sformatf("   Input file: %s\n", in_file.c_str()));
 		output_msg(sformatf("  Output file: %s\n", out_file.c_str()));
 #ifdef NPP
-		output_msg(sformatf("Using PHREEQC: version 3.6.5, compiled February 24, 2021\n"));
+		output_msg(sformatf("Using PHREEQC: version 3.7.1, compiled July 5, 2021\n"));
 #endif
 		output_msg(sformatf("Database file: %s\n\n", token.c_str()));
 #ifdef NPP
@@ -508,6 +509,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie,
 		/*
 		*   local cleanup
 		*/
+							  
 		line = (char *) free_check_null(line);
 		line_save = (char *) free_check_null(line_save);
 
diff --git a/phreeqcpp/gases.cpp b/phreeqcpp/gases.cpp
index e2c5495c..d64d2c97 100644
--- a/phreeqcpp/gases.cpp
+++ b/phreeqcpp/gases.cpp
@@ -101,7 +101,7 @@ build_fixed_volume_gas(void)
 		}
 
 		/* All elements in gas */
-		for (j = 0; j < count_elts; j++)
+		for (j = 0; j < (int) count_elts; j++)
 		{
 			unknown_ptr = NULL;
 			if (strcmp(elt_list[j].elt->name, "H") == 0)
@@ -149,7 +149,7 @@ build_fixed_volume_gas(void)
 			output_msg(sformatf( "\n\tJacobian summations %s.\n\n",
 					   phase_ptr->name));
 		}
-		for (j = 0; j < count_elts; j++)
+		for (j = 0; j < (int) count_elts; j++)
 		{
 			unknown_ptr = NULL;
 			if (strcmp(elt_list[j].elt->name, "H") == 0)
@@ -436,23 +436,31 @@ calc_PR(void)
 			{
 				if (!strcmp(phase_ptr1->name, "CO2(g)"))
 					a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217
-				else if (!strcmp(phase_ptr1->name, "H2S(g)"))
+				else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)"))
 					a_aa *= 0.81;
-				else if (!strcmp(phase_ptr1->name, "CH4(g)"))
+				else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)") || !strcmp(phase_ptr1->name, "Methane(g)"))
 					a_aa *= 0.51;
-				else if (!strcmp(phase_ptr1->name, "N2(g)"))
+				else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)"))
 					a_aa *= 0.51;
+				else if (!strcmp(phase_ptr1->name, "Ethane(g)"))
+					a_aa *= 0.51;
+				else if (!strcmp(phase_ptr1->name, "Propane(g)"))
+					a_aa *= 0.45;
 			}
 			if (!strcmp(phase_ptr1->name, "H2O(g)"))
 			{
 				if (!strcmp(phase_ptr->name, "CO2(g)"))
 					a_aa *= 0.81;
-				else if (!strcmp(phase_ptr->name, "H2S(g)"))
+				else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)"))
 					a_aa *= 0.81;
-				else if (!strcmp(phase_ptr->name, "CH4(g)"))
+				else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)") || !strcmp(phase_ptr->name, "Methane(g)"))
 					a_aa *= 0.51;
-				else if (!strcmp(phase_ptr->name, "N2(g)"))
+				else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)"))
 					a_aa *= 0.51;
+				else if (!strcmp(phase_ptr->name, "Ethane(g)"))
+					a_aa *= 0.51;
+				else if (!strcmp(phase_ptr->name, "Propane(g)"))
+					a_aa *= 0.45;
 			}
 			a_aa_sum += phase_ptr->fraction_x * phase_ptr1->fraction_x * a_aa;
 			a_aa_sum2 += phase_ptr1->fraction_x * a_aa;
@@ -598,10 +606,10 @@ calc_PR(void)
 		{
 			phi = B_r * (rz - 1) - log(rz - B) + A / (2.828427 * B) * (B_r - 2.0 * phase_ptr->pr_aa_sum2 / a_aa_sum) *
 				  log((rz + 2.41421356 * B) / (rz - 0.41421356 * B));
-			phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi));
+			//phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi));
 		}
 		else
-			phi = -3.0; // fugacity coefficient > 0.05
+			phi = -3.0; // fugacity coefficient = 0.05
 		phase_ptr->pr_phi = exp(phi);
 		phase_ptr->pr_si_f = phi / LOG_10;										// pr_si_f updated
 		// ****
diff --git a/phreeqcpp/global_structures.h b/phreeqcpp/global_structures.h
index e66b374c..b217d965 100644
--- a/phreeqcpp/global_structures.h
+++ b/phreeqcpp/global_structures.h
@@ -41,6 +41,7 @@
 #define DISP 2
 #define STAG 3
 #define NOMIX 4
+#define MIX_BS 5 // mix boundary solutions in electromigration
 
 #define CONVERGED 2
 #define MASS_BALANCE 3
diff --git a/phreeqcpp/inverse.cpp b/phreeqcpp/inverse.cpp
index 7ce99fbb..81a33781 100644
--- a/phreeqcpp/inverse.cpp
+++ b/phreeqcpp/inverse.cpp
@@ -272,7 +272,7 @@ setup_inverse(class inverse *inv_ptr)
 /*
  *   Define inv_zero and inv_zero array, delta
  */
-	for (i = 0; i < max; i++)
+	for (i = 0; i < (int) max; i++)
 		inv_zero[i] = 0.0;
 
 	memcpy((void *) &(delta[0]), (void *) &(inv_zero[0]),
@@ -3595,7 +3595,7 @@ check_isotopes(class inverse *inv_ptr)
 				/* use solution-defined uncertainties second */
 			}
 #ifdef NPP
-			else if (inv_ptr->i_u[i].count_uncertainties > 0
+			else if (inv_ptr->i_u[i].uncertainties.size() > 0
 				&& !isnan(inv_ptr->i_u[i].uncertainties[inv_ptr->i_u[i].uncertainties.size() - 1]))
 #else
 			else if (inv_ptr->i_u[i].uncertainties.size() > 0
diff --git a/phreeqcpp/kinetics.cpp b/phreeqcpp/kinetics.cpp
index dab9ab62..27ce4bdc 100644
--- a/phreeqcpp/kinetics.cpp
+++ b/phreeqcpp/kinetics.cpp
@@ -1595,7 +1595,7 @@ set_transport(int i, int use_mix, int use_kinetics, int nsaver)
 /*
  *   i			--user number for soln, reaction, etc.
  *   use_mix	  --integer flag
-					state == TRANSPORT: DISP, STAG, NOMIX
+					state == TRANSPORT: DISP, STAG, NOMIX, MIX_BS
 					state == REACTION: TRUE, FALSE
  *   use_kinetics --true or false flag to calculate kinetic reactions
  *   nsaver	   --user number to store solution
@@ -1615,7 +1615,7 @@ set_transport(int i, int use_mix, int use_kinetics, int nsaver)
 		use.Set_n_mix_user(i);
 		use.Set_n_mix_user_orig(i);
 	}
-	else if (use_mix == STAG && multi_Dflag != TRUE)
+	else if ((use_mix == STAG && multi_Dflag != TRUE) || use_mix == MIX_BS)
 	{
 		use.Set_mix_ptr(Utilities::Rxn_find(Rxn_mix_map, i));
 		if (use.Get_mix_ptr() != NULL)
diff --git a/phreeqcpp/model.cpp b/phreeqcpp/model.cpp
index cebc94e8..1d31226a 100644
--- a/phreeqcpp/model.cpp
+++ b/phreeqcpp/model.cpp
@@ -2613,13 +2613,12 @@ calc_gas_pressures(void)
 					V_m = (2. * gas_phase_ptr->Get_v_m() + V_m) / 3;
 				else
 					V_m = (1. * gas_phase_ptr->Get_v_m() + V_m) / 2;
-				if (iterations > 99 && numerical_fixed_volume == false)
+				if ((pitzer_model || iterations > 99) && numerical_fixed_volume == false)
 				{
 					//V_m *= 1; /* debug */
 					numerical_fixed_volume = true;
 					//switch_numerical = true;
-					warning_msg
-						("Numerical method failed, switching to numerical derivatives.");
+					if (!pitzer_model) warning_msg("Numerical method failed, switching to numerical derivatives.");
 					prep();
 					//switch_numerical = false;
 				}
@@ -2649,8 +2648,8 @@ calc_gas_pressures(void)
 				lp += rxn_ptr->s->la * rxn_ptr->coef;
 			}
 			phase_ptr->p_soln_x = exp(LOG_10 * (lp - phase_ptr->pr_si_f));
-			if (!strcmp(phase_ptr->name, "H2O(g)") && phase_ptr->p_soln_x > 90)
-					phase_ptr->p_soln_x = 90;
+			//if (!strcmp(phase_ptr->name, "H2O(g)") && phase_ptr->p_soln_x > 90)
+			//		phase_ptr->p_soln_x = 90;
 
 			if (gas_phase_ptr->Get_type() == cxxGasPhase::GP_PRESSURE)
 			{
@@ -5305,7 +5304,8 @@ ss_f(LDBLE xb, LDBLE l_a0, LDBLE l_a1, LDBLE l_kc, LDBLE l_kb, LDBLE xcaq,
 	f = xcaq * (xb / r + xc) + xbaq * (xb + r * xc) - 1;
 	return (f);
 }
-
+//#define ORIGINAL
+#ifdef ORIGINAL
 /* ---------------------------------------------------------------------- */
 int Phreeqc::
 numerical_jacobian(void)
@@ -5425,8 +5425,8 @@ numerical_jacobian(void)
 		case GAS_MOLES:
 			if (gas_in == FALSE)
 				continue;
-
-			d2 = d * x[i]->moles;
+			d2 = (x[i]->moles > 1 ? 1 : 20);
+			d2 *= d * x[i]->moles;
 			if (d2 < 1e-14)
 				d2 = 1e-14;
 			x[i]->moles += d2;
@@ -5546,7 +5546,265 @@ numerical_jacobian(void)
 	calculating_deriv = FALSE;
 	return OK;
 }
+#else
+/* ---------------------------------------------------------------------- */
+int Phreeqc::
+numerical_jacobian(void)
+/* ---------------------------------------------------------------------- */
+{
+	std::vector<double> base;
+	LDBLE d, d1, d2;
+	int i, j;
+	cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
+	std::vector<class phase*> phase_ptrs;
+	std::vector<class phase> base_phases;
+	double base_mass_water_bulk_x = 0, base_moles_h2o = 0;
+	cxxGasPhase base_gas_phase;
+	cxxSurface base_surface;
 
+	if (!
+		(numerical_deriv ||
+			(use.Get_surface_ptr() != NULL && use.Get_surface_ptr()->Get_type() == cxxSurface::CD_MUSIC) ||
+			(gas_phase_ptr != NULL && gas_phase_ptr->Get_type() == cxxGasPhase::GP_VOLUME &&
+				(gas_phase_ptr->Get_pr_in() || force_numerical_fixed_volume) && numerical_fixed_volume)
+			))
+		return(OK);
+
+	//jacobian_sums();
+	if (use.Get_surface_ptr() != NULL)
+	{
+		base_surface = *use.Get_surface_ptr();
+	}
+	if (use.Get_gas_phase_ptr() != NULL)
+	{
+		//cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
+		base_gas_phase = *gas_phase_ptr;
+		base_phases.resize(gas_phase_ptr->Get_gas_comps().size());
+		for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++)
+		{
+			const cxxGasComp* gas_comp_ptr = &(gas_phase_ptr->Get_gas_comps()[i]);
+			class phase* phase_ptr = phase_bsearch(gas_comp_ptr->Get_phase_name().c_str(), &j, FALSE);
+			phase_ptrs.push_back(phase_ptr);
+			base_phases[i] = *phase_ptr;
+		}
+	}
+	calculating_deriv = TRUE;
+	gammas(mu_x);
+	molalities(TRUE);
+	mb_sums();
+	//mb_gases();
+	//mb_ss();
+	residuals();
+	/*
+	 *   Clear array, note residuals are in array[i, count_unknowns+1]
+	 */
+
+	//for (i = 0; i < count_unknowns; i++)
+	//{
+	//	my_array[i] = 0.0;
+	//}
+	//for (i = 1; i < count_unknowns; i++)
+	//{
+	//	memcpy((void*)&(my_array[(size_t)i * (count_unknowns + 1)]),
+	//		(void*)&(my_array[0]), count_unknowns * sizeof(LDBLE));
+	//}
+
+	base.resize(count_unknowns);
+	base = residual;
+	d = 0.0001;
+	d1 = d * LOG_10;
+	d2 = 0;
+	for (i = 0; i < count_unknowns; i++)
+	{
+		switch (x[i]->type)
+		{
+		case MB:
+		case ALK:
+		case CB:
+		case SOLUTION_PHASE_BOUNDARY:
+		case EXCH:
+		case SURFACE:
+		case SURFACE_CB:
+		case SURFACE_CB1:
+		case SURFACE_CB2:
+			x[i]->master[0]->s->la += d;
+			d2 = d;// *LOG_10;
+			break;
+		case MH:
+			s_eminus->la += d;
+			d2 = d;// *LOG_10;
+			break;
+		case AH2O:
+			x[i]->master[0]->s->la += d;
+			d2 = d;// *LOG_10;
+			break;
+		case PITZER_GAMMA:
+			x[i]->s->lg += d;
+			d2 = d;
+			break;
+		case MH2O:
+			//mass_water_aq_x *= (1 + d);
+			//x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
+			//d2 = log(1.0 + d);
+			//break;
+			// DL_pitz
+			d1 = mass_water_aq_x * d;
+			mass_water_aq_x += d1;
+			if (use.Get_surface_in() && dl_type_x == cxxSurface::DONNAN_DL)
+				mass_water_bulk_x += d1;
+			x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
+			d2 = d1;
+			break;
+		case MU:
+			d2 = d * mu_x;
+			mu_x += d2;
+			//gammas(mu_x);
+			break;
+		case PP:
+			for (j = 0; j < count_unknowns; j++)
+			{
+				delta[j] = 0.0;
+			}
+			d2 = -1e-8;
+			delta[i] = d2;
+			reset();
+			d2 = delta[i];
+			break;
+		case SS_MOLES:
+			if (x[i]->ss_in == FALSE)
+				continue;
+			for (j = 0; j < count_unknowns; j++)
+			{
+				delta[j] = 0.0;
+			}
+			/*d2 = -1e-8; */
+			d2 = d * 10 * x[i]->moles;
+			//d2 = -.1 * x[i]->moles;
+			/*
+			   if (d2 > -1e-10) d2 = -1e-10;
+			   calculating_deriv = FALSE;
+			 */
+			delta[i] = d2;
+			/*fprintf (stderr, "delta before reset %e\n", delta[i]); */
+			reset();
+			d2 = delta[i];
+			/*fprintf (stderr, "delta after reset %e\n", delta[i]); */
+			break;
+		case GAS_MOLES:
+			if (gas_in == FALSE)
+				continue;
+			d2 = (x[i]->moles > 1 ? 1 : 30);
+			d2 *= d * x[i]->moles;
+			d2 = (d2 < ineq_tol ? ineq_tol : d2);
+			//if (d2 < 1e-14)
+			//	d2 = 1e-14;
+			x[i]->moles += d2;
+			break;
+		}
+		gammas(mu_x);
+		molalities(TRUE);
+		mb_sums();
+		//mb_gases();
+		//mb_ss();
+		residuals();
+
+		for (j = 0; j < count_unknowns; j++)
+		{
+			my_array[(size_t)j * (count_unknowns + 1) + (size_t)i] = -(residual[j] - base[j]) / d2;
+			if (x[i]->type == MH2O) // DL_pitz
+				my_array[(size_t)j * (count_unknowns + 1) + (size_t)i] *= mass_water_aq_x;
+		}
+		switch (x[i]->type)
+		{
+		case MB:
+		case ALK:
+		case CB:
+		case SOLUTION_PHASE_BOUNDARY:
+		case EXCH:
+		case SURFACE:
+		case SURFACE_CB:
+		case SURFACE_CB1:
+		case SURFACE_CB2:
+		case AH2O:
+			x[i]->master[0]->s->la -= d2;
+			break;
+		case MH:
+			s_eminus->la -= d2;
+			if (my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] == 0)
+			{
+				/*output_msg(sformatf( "Zero diagonal for MH\n")); */
+				my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] =
+					under(s_h2->lm) * 2;
+			}
+			break;
+		case PITZER_GAMMA:
+			x[i]->s->lg -= d2;
+			break;
+		case MH2O:
+			//mass_water_aq_x /= (1 + d);
+			//x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
+			//break;
+			//DL_pitz
+			mass_water_aq_x -= d2;
+			if (use.Get_surface_in() && dl_type_x == cxxSurface::DONNAN_DL)
+				mass_water_bulk_x -= d2;
+			x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
+			break;
+		case MU:
+			mu_x -= d2;
+			//gammas(mu_x);
+			break;
+		case PP:
+			delta[i] = -d2;
+			reset();
+			break;
+		case SS_MOLES:
+			delta[i] = -d2;
+			reset();
+			break;
+		case GAS_MOLES:
+			x[i]->moles -= d2;
+			break;
+		}
+		if (use.Get_surface_ptr() != NULL)
+		{
+			*use.Get_surface_ptr() = base_surface;
+		}
+		if (use.Get_gas_phase_ptr() != NULL)
+		{
+			*use.Get_gas_phase_ptr() = base_gas_phase;
+			for (size_t g = 0; g < base_phases.size(); g++)
+			{
+				*phase_ptrs[g] = base_phases[g];
+			}
+		}
+		//gammas(mu_x);
+		//molalities(TRUE);
+		//mb_sums();
+		////mb_gases();
+		////mb_ss();
+		//residuals();
+	}
+	gammas(mu_x);
+	molalities(TRUE);
+	mb_sums();
+	//mb_gases();
+	//mb_ss();
+	residuals();
+	//for (i = 0; i < count_unknowns; i++)
+	//{
+	//	//Debugging
+	//	if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 1e-2 &&
+	//		fabs(residual[i]) + fabs(base[i]) > 1e-8)
+	//	{
+	//		std::cerr << iterations << ": " << x[i]->description << "  " << residual[i] << "  " << base[i] << std::endl;
+	//	}
+	//}
+	base.clear();
+	calculating_deriv = FALSE;
+	return OK;
+}
+#endif
 /* ---------------------------------------------------------------------- */
 void Phreeqc::
 set_inert_moles(void)
diff --git a/phreeqcpp/pitzer.cpp b/phreeqcpp/pitzer.cpp
index 59882063..68a57a39 100644
--- a/phreeqcpp/pitzer.cpp
+++ b/phreeqcpp/pitzer.cpp
@@ -1948,19 +1948,40 @@ jacobian_pz(void)
 /* ---------------------------------------------------------------------- */
 { // calculate the derivatives numerically
 	std::vector<double> base;
+	std::vector<class phase*> phase_ptrs;
+	std::vector<class phase> base_phases;
+	cxxGasPhase base_gas_phase;
+	cxxSurface base_surface;
 	LDBLE d, d1, d2;
 	int i, j;
-
-	calculating_deriv = 1;
 Restart:
-	size_t pz_max_unknowns = max_unknowns;
-	//k_temp(tc_x, patm_x);
+	if (use.Get_surface_ptr() != NULL)
+	{
+		base_surface = *use.Get_surface_ptr();
+	}
+	if (use.Get_gas_phase_ptr() != NULL)
+	{
+		cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
+		base_gas_phase = *gas_phase_ptr;
+		base_phases.resize(gas_phase_ptr->Get_gas_comps().size());
+		for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++)
+		{
+			const cxxGasComp* gas_comp_ptr = &(gas_phase_ptr->Get_gas_comps()[i]);
+			class phase* phase_ptr = phase_bsearch(gas_comp_ptr->Get_phase_name().c_str(), &j, FALSE);
+			phase_ptrs.push_back(phase_ptr);
+			base_phases[i] = *phase_ptr;
+		}
+	}
+	calculating_deriv = 1;
+	molalities(TRUE);
 	if (full_pitzer == TRUE)
 	{
-		molalities(TRUE);
 		pitzer();
-		residuals();
 	}
+	mb_sums();
+	residuals();
+
+	size_t pz_max_unknowns = max_unknowns;
 	base.resize(count_unknowns);
 	for (i = 0; i < count_unknowns; i++)
 	{
@@ -2026,9 +2047,11 @@ Restart:
 		case GAS_MOLES:
 			if (gas_in == FALSE)
 				continue;
-			d2 = d * x[i]->moles;
-			if (d2 < 1e-14)
-				d2 = 1e-14;
+			d2 = (x[i]->moles > 1 ? 1 : 30);
+			d2 *= d * x[i]->moles;
+			d2 = (d2 < ineq_tol ? ineq_tol : d2);
+			//if (d2 < 1e-14)
+			//	d2 = 1e-14;
 			x[i]->moles += d2;
 			break;
 		case MU:
@@ -2124,12 +2147,38 @@ Restart:
 			reset();
 			break;
 		}
+		if (use.Get_surface_ptr() != NULL)
+		{
+			*use.Get_surface_ptr() = base_surface;
+		}
+		if (use.Get_gas_phase_ptr() != NULL)
+		{
+			*use.Get_gas_phase_ptr() = base_gas_phase;
+			for (size_t g = 0; g < base_phases.size(); g++)
+			{
+				*phase_ptrs[g] = base_phases[g];
+			}
+		}
+		//molalities(TRUE);
+		//if (full_pitzer == TRUE)
+		//	pitzer();
+		//mb_sums();
+		//residuals();
 	}
 	molalities(TRUE);
 	if (full_pitzer == TRUE)
 		pitzer();
 	mb_sums();
 	residuals();
+	//for (i = 0; i < count_unknowns; i++)
+	//{
+	//	//Debugging
+	//	if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 1e-2 &&
+	//		fabs(residual[i]) + fabs(base[i]) > 1e-6)
+	//	{
+	//		std::cerr << iterations << ": " << x[i]->description << "  " << residual[i] << "  " << base[i] << std::endl;
+	//	}
+	//}
 	base.clear();
 	calculating_deriv = 0;
 	return OK;
diff --git a/phreeqcpp/prep.cpp b/phreeqcpp/prep.cpp
index 5791517b..82b3125c 100644
--- a/phreeqcpp/prep.cpp
+++ b/phreeqcpp/prep.cpp
@@ -3841,15 +3841,11 @@ calc_PR(std::vector<class phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
 			{
 				if (!strcmp(phase_ptr1->name, "CO2(g)"))
 					a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217
-				else if (!strcmp(phase_ptr1->name, "H2S(g)"))
+				else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)"))
 					a_aa *= 0.81;
-				else if (!strcmp(phase_ptr1->name, "CH4(g)"))
+				else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)") || !strcmp(phase_ptr1->name, "Methane(g)"))
 					a_aa *= 0.51;
-				else if (!strcmp(phase_ptr1->name, "Mtg(g)"))
-					a_aa *= 0.51;
-				else if (!strcmp(phase_ptr1->name, "Methane(g)"))
-					a_aa *= 0.51;
-				else if (!strcmp(phase_ptr1->name, "N2(g)"))
+				else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)"))
 					a_aa *= 0.51;
 				else if (!strcmp(phase_ptr1->name, "Ethane(g)"))
 					a_aa *= 0.51;
@@ -3860,15 +3856,11 @@ calc_PR(std::vector<class phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
 			{
 				if (!strcmp(phase_ptr->name, "CO2(g)"))
 					a_aa *= 0.81;
-				else if (!strcmp(phase_ptr->name, "H2S(g)"))
+				else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)"))
 					a_aa *= 0.81;
-				else if (!strcmp(phase_ptr->name, "CH4(g)"))
+				else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)") || !strcmp(phase_ptr->name, "Methane(g)"))
 					a_aa *= 0.51;
-				else if (!strcmp(phase_ptr->name, "Mtg(g)"))
-					a_aa *= 0.51;
-				else if (!strcmp(phase_ptr->name, "Methane(g)"))
-					a_aa *= 0.51;
-				else if (!strcmp(phase_ptr->name, "N2(g)"))
+				else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)"))
 					a_aa *= 0.51;
 				else if (!strcmp(phase_ptr->name, "Ethane(g)"))
 					a_aa *= 0.51;
@@ -4002,17 +3994,17 @@ calc_PR(std::vector<class phase *> phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m)
 		{
 			phi = B_r * (rz - 1) - log(rz - B) + A / (2.828427 * B) * (B_r - 2.0 * phase_ptr->pr_aa_sum2 / a_aa_sum) *
 				log((rz + 2.41421356 * B) / (rz - 0.41421356 * B));
-			phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi));
+			//phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi));
 			//if (phi > 4.44)
 			//	phi = 4.44;
 		}
 		else
-			phi = -3.0; // fugacity coefficient > 0.05
-		if (/*!strcmp(phase_ptr->name, "H2O(g)") && */phi < -3)
-		{
-			// avoid such phi...
-			phi = -3;
-		}
+			phi = -3.0; // fugacity coefficient = 0.05
+		//if (/*!strcmp(phase_ptr->name, "H2O(g)") && */phi < -3)
+		//{
+		////	 avoid such phi...
+		//	phi = -3;
+		//}
 		phase_ptr->pr_phi = exp(phi);
 		phase_ptr->pr_si_f = phi / LOG_10;
 		// for initial equilibrations, adapt log_k of the gas phase...
diff --git a/phreeqcpp/print.cpp b/phreeqcpp/print.cpp
index cb01ea9f..d4ad7b79 100644
--- a/phreeqcpp/print.cpp
+++ b/phreeqcpp/print.cpp
@@ -708,8 +708,8 @@ print_gas_phase(void)
 				   (double) initial_moles,
 				   (double) moles,
 				   (double) delta_moles));
-		if (!strcmp(phase_ptr->name, "H2O(g)") && phase_ptr->p_soln_x == 90)
-			output_msg("       WARNING: The pressure of H2O(g) is above the program limit: use the polynomial for log_k.\n");
+		//if (!strcmp(phase_ptr->name, "H2O(g)") && phase_ptr->p_soln_x == 90)
+		//	output_msg("       WARNING: The pressure of H2O(g) is fixed to the program limit.\n");
 
 	}
 	output_msg("\n");
diff --git a/phreeqcpp/read.cpp b/phreeqcpp/read.cpp
index fb9e78e1..4ade4fe4 100644
--- a/phreeqcpp/read.cpp
+++ b/phreeqcpp/read.cpp
@@ -2619,7 +2619,7 @@ read_aq_species_vm_parms(const char* cptr, LDBLE * delta_v)
 	/*
 	*   Read supcrt parms and Ionic strength terms
 	*/
-	for (j = 0; j < 9; j++)
+	for (j = 0; j < 11; j++)
 	{
 		delta_v[j] = 0.0;
 	}
@@ -2627,7 +2627,7 @@ read_aq_species_vm_parms(const char* cptr, LDBLE * delta_v)
 /* Vmax, dmax...
 	delta_v[10] = 999.0;
 	delta_v[11] = 1.0; */
-	j = sscanf(cptr, SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT /*SCANFORMAT SCANFORMAT */,
+	j = sscanf(cptr, SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT ,
 		/* a1..a4 */
 		&(delta_v[0]), &(delta_v[1]), &(delta_v[2]), &(delta_v[3]),
 		/* wref */
@@ -2635,9 +2635,9 @@ read_aq_species_vm_parms(const char* cptr, LDBLE * delta_v)
 		/* b_Av */
 		&(delta_v[5]),
 		/* c1..c4 */
-		&(delta_v[6]), &(delta_v[7]), &(delta_v[8]), &(delta_v[9]));
-		/* vmax, dmax 
-		&(delta_v[10]), &(delta_v[11])); */
+		&(delta_v[6]), &(delta_v[7]), &(delta_v[8]), &(delta_v[9]),
+		 //vmP, exP
+		&(delta_v[10]), &(delta_v[11])); 
 	if (j < 1)
 	{
 		input_error++;
@@ -2665,8 +2665,8 @@ read_vm_only(const char* cptr, LDBLE * delta_v, DELTA_V_UNIT * units)
 	char token[MAX_LENGTH];
 	/*
 	*   Read analytical expression
-	*/
-	for (j = 0; j < 8; j++)
+	*/			
+	for (j = 0; j < 9; j++)
 	{
 		delta_v[j] = 0.0;
 	}
diff --git a/phreeqcpp/sit.cpp b/phreeqcpp/sit.cpp
index 27825c26..fb372c3b 100644
--- a/phreeqcpp/sit.cpp
+++ b/phreeqcpp/sit.cpp
@@ -798,7 +798,8 @@ sit_revise_guesses(void)
 	/*gammas(mu_x); */
 	return (OK);
 }
-
+//#define ORIGINAL
+#ifdef ORIGINAL
 /* ---------------------------------------------------------------------- */
 int Phreeqc::
 jacobian_sit(void)
@@ -870,7 +871,8 @@ Restart:
 		case GAS_MOLES:
 			if (gas_in == FALSE)
 				continue;
-			d2 = d * x[i]->moles;
+			d2 = (x[i]->moles > 1 ? 1 : 20);
+			d2 *= d * x[i]->moles;
 			if (d2 < 1e-14)
 				d2 = 1e-14;
 			x[i]->moles += d2;
@@ -951,7 +953,209 @@ Restart:
 	residuals();
 	return OK;
 }
-
+#else
+/* ---------------------------------------------------------------------- */
+int Phreeqc::
+jacobian_sit(void)
+/* ---------------------------------------------------------------------- */
+{
+	std::vector<double> base;
+	LDBLE d, d1, d2;
+	int i, j;
+	cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
+	std::vector<class phase*> phase_ptrs;
+	std::vector<class phase> base_phases;
+	double base_mass_water_bulk_x = 0, base_moles_h2o = 0;
+	cxxGasPhase base_gas_phase;
+	cxxSurface base_surface;
+Restart:
+	if (use.Get_surface_ptr() != NULL)
+	{
+		base_surface = *use.Get_surface_ptr();
+	}
+	if (use.Get_gas_phase_ptr() != NULL)
+	{
+		cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
+		base_gas_phase = *gas_phase_ptr;
+		base_phases.resize(gas_phase_ptr->Get_gas_comps().size());
+		for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++)
+		{
+			const cxxGasComp* gas_comp_ptr = &(gas_phase_ptr->Get_gas_comps()[i]);
+			class phase* phase_ptr = phase_bsearch(gas_comp_ptr->Get_phase_name().c_str(), &j, FALSE);
+			phase_ptrs.push_back(phase_ptr);
+			base_phases[i] = *phase_ptr;
+		}
+	}
+	calculating_deriv = 1;
+	size_t pz_max_unknowns = max_unknowns;
+	//k_temp(tc_x, patm_x);
+	molalities(TRUE);
+	if (full_pitzer == TRUE)
+	{
+		
+		sit();
+	}
+	mb_sums();
+	residuals();
+	base = residual; // std::vectors
+	d = 0.0001;
+	d1 = d * LOG_10;
+	d2 = 0;
+	for (i = 0; i < count_unknowns; i++)
+	{
+		switch (x[i]->type)
+		{
+		case MB:
+		case ALK:
+		case CB:
+		case SOLUTION_PHASE_BOUNDARY:
+		case EXCH:
+		case SURFACE:
+		case SURFACE_CB:
+		case SURFACE_CB1:
+		case SURFACE_CB2:
+			x[i]->master[0]->s->la += d;
+			d2 = d1;
+			break;
+		case AH2O:
+			x[i]->master[0]->s->la += d;
+			d2 = d1;
+			break;
+		case PITZER_GAMMA:
+			if (!full_pitzer)
+				continue;
+			x[i]->s->lg += d;
+			d2 = d;
+			break;
+		case MH2O:
+			mass_water_aq_x *= (1.0 + d);
+			x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
+			d2 = log(1.0 + d);
+			break;
+		case MH:
+			s_eminus->la += d;
+			d2 = d1;
+			break;
+			/*
+			if (pitzer_pe == TRUE)
+			{
+				s_eminus->la += d;
+				d2 = d1;
+				break;
+			}
+			else
+			{
+				continue;
+			}
+			*/
+		case GAS_MOLES:
+			if (gas_in == FALSE)
+				continue;
+			d2 = (x[i]->moles > 1 ? 1 : 20);
+			d2 *= d * x[i]->moles;
+			if (d2 < 1e-14)
+				d2 = 1e-14;
+			x[i]->moles += d2;
+			break;
+		case MU:
+			//continue;
+			d2 = d * mu_x;
+			mu_x += d2;
+			//k_temp(tc_x, patm_x);
+			gammas_sit();
+			break;
+		case PP:
+		case SS_MOLES:
+			continue;
+			break;
+		}
+		molalities(TRUE);
+		if (max_unknowns > pz_max_unknowns)
+		{
+			gammas_sit();
+			jacobian_sums();
+			goto Restart;
+		}
+		if (full_pitzer == TRUE)
+			sit();
+		mb_sums();
+		residuals();
+		for (j = 0; j < count_unknowns; j++)
+		{
+			my_array[(size_t)j * (count_unknowns + 1) + (size_t)i] =
+				-(residual[j] - base[j]) / d2;
+		}
+		switch (x[i]->type)
+		{
+		case MB:
+		case ALK:
+		case CB:
+		case SOLUTION_PHASE_BOUNDARY:
+		case EXCH:
+		case SURFACE:
+		case SURFACE_CB:
+		case SURFACE_CB1:
+		case SURFACE_CB2:
+		case AH2O:
+			x[i]->master[0]->s->la -= d;
+			break;
+		case MH:
+			s_eminus->la -= d;
+			if (my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] == 0)
+			{
+				my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] =
+					exp(s_h2->lm * LOG_10) * 2;
+			}
+			break;
+		case PITZER_GAMMA:
+			x[i]->s->lg -= d;
+			break;
+		case MH2O:
+			mass_water_aq_x /= (1 + d);
+			x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
+			break;
+		case MU:
+			mu_x -= d2;
+			//k_temp(tc_x, patm_x);
+			gammas_sit();
+			break;
+		case GAS_MOLES:
+			if (gas_in == FALSE)
+				continue;
+			x[i]->moles -= d2;
+			break;
+		}
+		if (use.Get_surface_ptr() != NULL)
+		{
+			*use.Get_surface_ptr() = base_surface;
+		}
+		if (use.Get_gas_phase_ptr() != NULL)
+		{
+			*use.Get_gas_phase_ptr() = base_gas_phase;
+			for (size_t g = 0; g < base_phases.size(); g++)
+			{
+				*phase_ptrs[g] = base_phases[g];
+			}
+		}
+	}
+	molalities(TRUE);
+	if (full_pitzer == TRUE)
+		sit();
+	mb_sums();
+	residuals();
+	//for (i = 0; i < count_unknowns; i++)
+	//{
+	//	//Debugging
+	//	if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 1e-2 &&
+	//		fabs(residual[i]) + fabs(base[i]) > 1e-8)
+	//	{
+	//		std::cerr << i << ": " << x[i]->description << "  " << residual[i] << "  " << base[i] << std::endl;
+	//	}
+	//}
+	calculating_deriv = 0;
+	return OK;
+}
+#endif
 /* ---------------------------------------------------------------------- */
 int Phreeqc::
 model_sit(void)
diff --git a/phreeqcpp/step.cpp b/phreeqcpp/step.cpp
index b7204ecd..a5503f02 100644
--- a/phreeqcpp/step.cpp
+++ b/phreeqcpp/step.cpp
@@ -57,11 +57,19 @@ step(LDBLE step_fraction)
  */
 	if (use.Get_mix_ptr() != NULL)
 	{
-		add_mix(use.Get_mix_ptr());
+			add_mix(use.Get_mix_ptr());
+			int n = use.Get_n_mix_user_orig();
+			if (n == 0 || n == count_cells + 1)
+			{
+				cxxSolution *solution_ptr = Utilities::Rxn_find(Rxn_solution_map, n);
+				if (solution_ptr != NULL && !solution_ptr->Get_new_def())
+					potV_x = solution_ptr->Get_potV();
+			}
 	}
 	else if (use.Get_solution_ptr() != NULL)
 	{
 		add_solution(use.Get_solution_ptr(), 1.0, 1.0);
+		potV_x = use.Get_solution_ptr()->Get_potV();
 		cell_no = use.Get_n_solution_user();
 	}
 	else
@@ -367,7 +375,7 @@ add_solution(cxxSolution *solution_ptr, LDBLE extensive, LDBLE intensive)
 	tc_x += solution_ptr->Get_tc() * intensive;
 	ph_x += solution_ptr->Get_ph() * intensive;
 	patm_x += solution_ptr->Get_patm() * intensive;
-	potV_x += solution_ptr->Get_potV() * intensive;
+	//potV_x += solution_ptr->Get_potV() * intensive;
 	solution_pe_x += solution_ptr->Get_pe() * intensive;
 	mu_x += solution_ptr->Get_mu() * intensive;
 	ah2o_x += solution_ptr->Get_ah2o() * intensive;
diff --git a/phreeqcpp/transport.cpp b/phreeqcpp/transport.cpp
index ca2317c0..00fedd8a 100644
--- a/phreeqcpp/transport.cpp
+++ b/phreeqcpp/transport.cpp
@@ -1,4 +1,4 @@
-#include "Utils.h"
+#include "Utils.h"
 #include "Phreeqc.h"
 #include "phqalloc.h"
 #include "Exchange.h"
@@ -605,7 +605,10 @@ transport(void)
 
 						if (i == 0 || i == count_cells + 1)
 						{
-							run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i
+							if (dV_dcell)
+								run_reactions(i, kin_time, MIX_BS, step_fraction); // nsaver = i
+							else
+								run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i
 						}
 						else
 						{
@@ -868,7 +871,12 @@ transport(void)
 					status(0, token);
 
 					if (i == 0 || i == count_cells + 1)
-						run_reactions(i, kin_time, NOMIX, step_fraction);
+					{
+						if (dV_dcell)
+							run_reactions(i, kin_time, MIX_BS, step_fraction); // nsaver = i
+						else
+							run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i
+					}
 					else
 					{
 						run_reactions(i, kin_time, DISP, step_fraction);
@@ -1070,6 +1078,12 @@ print_punch(int i, boolean active)
 	if (!active)
 		run_reactions(i, 0, NOMIX, 0);
 	cell_no = i;
+	if (dV_dcell)
+	{
+		use.Set_n_solution_user(i);
+		use.Get_solution_ptr()->Set_potV(cell_data[i].potV);
+		potV_x = cell_data[i].potV;
+	}
 	use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, i));
 	if (use.Get_kinetics_ptr() != NULL)
 	{
@@ -1413,7 +1427,7 @@ mix_stag(int i, LDBLE kin_time, int l_punch, LDBLE step_fraction)
 				for (std::map < int, LDBLE >::const_iterator it = use.Get_mix_ptr()->Get_mixComps().begin();
 					it != use.Get_mix_ptr()->Get_mixComps().end(); it++)
 				{
-					if (it->first > i && it->first < all_cells)
+					if (it->first > i && it->first < all_cells && it->first != count_cells + 1)
 					{
 						k = it->first;
 						ptr_imm = Utilities::Rxn_find(Rxn_solution_map, k);
@@ -6053,6 +6067,7 @@ calc_vm_Cl(void)
 		V_Cl = s_ptr->logk[vma1] + s_ptr->logk[vma2] / pb_s +
 			(s_ptr->logk[vma3] + s_ptr->logk[vma4] / pb_s) / TK_s -
 			s_ptr->logk[wref] * QBrn;
+
 		/* the ionic strength term * I^0.5... */
 		if (s_ptr->logk[b_Av] < 1e-5)
 			V_Cl += s_ptr->z * s_ptr->z * 0.5 * DH_Av * sqrt_mu;